texts:
I have set of referenced rasters (.tif), which need to be tiled - but each of them need to be rotated by a certain degree before that using python (gdal mainly as other routines use it as well). In QGIS, for example, I can achieve the right view by changing map rotation value at the bottom right corner of the screen - it does exactly what i want, but does not affect tiling process in any way and can't be done outside of QGIS interface, which is important too.

While using North Pole Azimuthal projection (EPSG:102016) I tried to define custom CRS similar to it, but had no sucess as +alpha and +gamma parameters work only with +omerc projections.

I would be very grateful for pointing me in the right direction. Should i use custom wkt-defined CRS-es or, perhaps, in gdal tiling instrument there is some built-in way for rotating an image..?

P.S. Using azimuthal projection for tiling is kind of weird but it has its purpose.
Thanks in advance for your help!

UPDATE:
Seems like rotation can be performed by changing affine parameters via gdal.SetGeoTransform. However with formulae provided by @MortenSickel i got following result (1 is raster before changing affines, 2 - after). Btw the order of parameters is different in gdal compared to esri world file and that can lead to errors, but i double-checked that.


texts:
QGIS error message when in Georeferencer:
.jpg is not a supported raster data source.
libjpeg: Bogus Huffman table definition
How do I fix this?
Tried different .jpeg and still doesn't work.
Previous fixes have included coding stuff which I don't understand, so haven't used.
Anyone have a plain English explaintion/fix?
The current released version of QGIS is 3.14.15.
On Windows 7.


texts: Defining the boundary of a catchment area

texts:
I digitized contour lines from topographic maps and created a DEM. With the DEM, I estimated the boundary of a catchment area. Estimated positional accuracy of topographic maps is RMSE = 65 m, with a systematic bias towards East. Using an affine transformation (R-package vec2dtransf), I shifted the data towards West. The RMSE of the affine transformation is 30 m.

Visual check of the boundary using satellite images suggests that its actual position is more accurate than the former. However, I understand that there still is an uncertainty about the correct position and shape of the boundary.

I will use the estimated boundary to limit the area were to make predictions of soil properties. But I want to be more confident that predictions cover the real area of the catchment. Thus, I am planning to increase the area adding a buffer equal to the RMSE of the affine transformation (30 m).

I like to know if this approach is correct and if I am missing something.


texts: Creating Time-series (Temporal) Map of NDVI?

texts:
Is there a best practice for a cartographer to create a Time-series or Temporal Map ? I have 83 of the NDVI in an imagine format and should I converted them to polygons and populated them and display them in Monthly series as a field ?

This is the one I will be doing a similar to this but it will be mainly the NDVI and how do I put the legend ?

NDVI Map using Time-series


texts: Creating transparent marker/symbol in ArcGIS Desktop?

texts:
How can I create a transparent marker/symbol in ArcGIS Desktop 10?

I am trying to do this by setting up the transparency of the layer, but the markers seem unaffected by this change.

Any ideas?


texts:
I have the following query which gives me two columns store and customer_count...one row per store.

SELECT a.id AS store,count(b.*) AS customer_count FROM stores a, customers b
WHERE ST_DWithin(a.geom::geography,b.geom::geography,5000) GROUP BY a.id

This works perfectly.

What I need to add to this is another variable that represents the percentage of customers for each store. That is, customer_count/total number of rows in customer table. I'm just not sure how to structure this.


texts: Identifying streets of same name where split by osm

texts:
I have a large map of every street in the country.
This is a vector layer in QGIS
many streets are split by OSM , IE have a different OSM ID yet the street has the same name.
Often these streets are split at intersecting points.
Ideally I'd like to join them using Qgis and then export my shapefile
My software previously looked at points in a route and determined if the street the user was on, or turning onto, continued or ended at the intersection: hence allowing me to determine road rules for the intersection.
For example our main street is one road, but is split into about 10 sections.
Is there an "easy" way to join streets or determine if the street is the entire street or just a section of the street ?


texts:
I installed OSGEO4w (selecting Qgis Desktop 2.12, Qgis Server and Apache Web Server) on Windows 7 because I must publish some Qgis project on my pc as server.
I have a problem to set my pc as Server WMS.
When I try to execute "Add layer from WMS Server" and click on Search Server I do not find my localhost, if I try to create "new server" setting URL as http://localhost  and click on Connect button I receive following message:

  "The server you are trying to connect does not appear to be a WMS
  server . Check the URL address ."

and in "show detail" for error I obtain:

Instead of the string of capabilities, it received the following response :

<title>OSGeo4W - OSGeo for Windows</title>

OSGeo4W - OSGeo for Windows

 Installed Web Applications and Features

There are currently
   2 osgeo4w applications   installed.

Apache 2.2.8:

Manual

PHP 5.2.5:

phpinfo

What is the problem?


texts: Save GPX from drawn feature in OpenLayers

texts:
I am trying to modify the Vector Formats example to only display the GPX in the textarea rather than all the options in the original example but I can't get it working.

My code below currently displays [object Object] is the text area. The next stage would be to save it as a GPX file but I can't even get this part to work. Help appreciated.

Draw route to update the output below.

<textarea id="output"></textarea>


texts: Restricting ArcGIS network analysis to finding origins/destinations with common ID?

texts:
I have two point feature classes, patients and physicians. The patient feature class has a record of which physician each patient visits (a one to one relationship). I would like to run an OD Cost Matrix (or whatever works) and restrict analysis to finding routes where the patient's physician ID (origins) and the physician ID (destination) match.

The only solution that I can think of (which is far from elegant) is to find all destinations for each origin and then extract those where the two IDs match.

Am I missing something obvious?

P.s. The network layer does not have a time field, so VRP is unavailable.


texts:
I am using ArcGIS Desktop.

When using the MODIS MCD12Q1 IGBP classifications, is there any way to run accuracy assessment or Quality Control is the only way?


texts: Radiometric calibration of Sentinel-2 products

texts:
How can I calibrate Sentinel-2 products, to extract the reflectance?

I've downloaded few Sentinel-2 products from https://scihub.copernicus.eu/s2/#/home, and opened them with SNAP but its radiometric calibration is dedicated only to SAR products.

On ENVI, the Radiometric calibration tool is seemingly only capable of extracting the radiance. What is the difference between the radiance and the reflectance? Is the estimation of the radiance counted as a radiometric calibration of a product?


texts: Get attributes from one line feature to another

texts:
I have to use two different networks (e.g., Network A and Network B) to do analysis and they don't align with each other perfectly. I used the same points to split the lines in the two networks, so each line in the Network A can find another corresponding line in the Network B, as shown in the figure. How can I get some attributes from the red line between the black triangle points to the green line?
I am using ArcGIS Advanced License Desktop 10.3.1.


texts: Missing data and negative predictions when kriging daily rainfall?

texts:
I want to perform the interpolation of daily rainfall data for 7 years. I am now working on python using arcpy module. I face several problems while performing interpolation

Some day there is no rainfall in all the station so that I am getting error for this issue.
After doing the interpolation I am getting some absurd values with negative values which is far from the real observation.

Is there a way to identify the parameters used for kriging interpolation?

I found that I need to change the parameters of the variogram for each day rainfall. Is this true?


texts: Are there Sentinel crop labeled datasets available?

texts: I am starting a project where I am going to use Sentinel-2 imagery and machine learning for detecting sugarcane and coffee. In order to train the algorithm I would need labeled sentinel datasets. I've found EuroSat which is a labeled dataset for different land cover classes. I have also found the Brazilian Coffee Scenes Dataset. Is there a similar dataset for sugarcane, coffee and other crops with Sentinel?


texts: ArcMap grids with different coordinate systems problem

texts:
I'm doing an exercice in cartography. We should do the canvas to a topographic map (1:25000), we start with a given point, we make the other 15, into an Excel file, (map is 7'30" * 7'30") the goal is to make two grids :

First : Graticule grid, graticule and labels (interval is 2'30"), geographical coordinate system

Second : Measured grid, tick marks and labels, projected coordinate system

The problem is when I want to do the first grid, it doesn't show as intended due to the data frame's coordinate system being the projected one, and if I change to the geographical one the shape of the map changes to a square instead of a rectangle. The second grid I don't have a problem with.


texts: Correct shapefile for Montenegro municipalities?

texts:
I downloaded a shapefile for Montenegro from https://gadm.org/download_country_v3.html and created a table in PostGIS using shp2pgsql. The table has column geom of type geometry which is a multipolygon. I'm drawing these multipolygons as polygons using Leaflet.draw (converting multipolygon to single polygon is the reason why some municipalities in the south with islands have strange boundaries).

As you can see data from shapefile isn't correct (especially in the north). I'm using OpenStreetMap as basemap provider, so the map images are correct. Also to prove GADM data isn't correct, here is the image of the most eastern municipality on their website https://gadm.org/maps/MNE/rozaje.html.

It's clearly the same shape as the one on the map. I don't mind shapes not being too detailed or slightly incorrect, but I would like at least their positions to be accurate. Each of the municipalities in the north is moved like 10 km to the south.

Please share if you know a better source for this type of data.


texts: Using QField to connect to external Bluetooth GPS

texts: I am using QField for field data acquisition. Can QField connect to an external Bluetooth GPS?


texts: Unable to export time lapse video in GEE because of bit size of satellite imagery

texts:
I am trying to create a time lapse video in Google Earth Engine code editor of the Arikaree glacier in CO with the following code:
  .filter(ee.Filter.eq('WRS_PATH', 34))
  .filter(ee.Filter.eq('WRS_ROW', 32))
  .filter(ee.Filter.lt('CLOUD_COVER',25))
  .filterBounds(ee.Geometry.Rectangle(-105.647002, 40.042805, -105.625015, 40.060087))
  .select(['B3', 'B2', 'B1']);

var L7coll = ee.ImageCollection('LANDSAT/LE07/C01/T1_SR')
  .filter(ee.Filter.eq('WRS_PATH', 34))
  .filter(ee.Filter.eq('WRS_ROW', 32))
  .filter(ee.Filter.lt('CLOUD_COVER',25))
  .filterBounds(ee.Geometry.Rectangle(-105.647002, 40.042805, -105.625015, 40.060087))
  .select(['B3', 'B2', 'B1'])
  .map(function(image){
  var filled1a = image.focal_mean(2, 'square', 'pixels', 1);
  return filled1a.blend(image);
});

var L8coll = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR')
  .filter(ee.Filter.eq('WRS_PATH', 34))
  .filter(ee.Filter.eq('WRS_ROW', 32))
  .filterBounds(ee.Geometry.Rectangle(-105.647002, 40.042805, -105.625015, 40.060087))
  .filter(ee.Filter.lt('CLOUD_COVER',25))
  .map(function(image){
  return image.rename(['B0', 'B1', 'B2', 'B3', 'B4', 'B5', 'B6', 'B7', 'B8', 'B9', 'B10', 'B11']);
  })
  .select(['B3', 'B2', 'B1']);

var collection_merge = ee.ImageCollection(L5coll.merge(L7coll.merge(L8coll)));
print (collection_merge);

var years = ee.List.sequence(1984, 2020);
print (years);

// Export the image, specifying scale and region.
Export.video.toDrive({
  collection: collection_merge,
  description: 'Arikaree3',
  framesPerSecond: 5,
  scale : (0,255 )});

When I try to create the video, I get the error:
"Error: Data type for band [B3] is not supported by video export: Type<Short>. Video export requires 8-bit bands; consider scaling to 0-255 and casting to 'uint8'."

I'm unsure what to do. Any suggestions?
From my research, it seems that Landsat 8 may have a 16-bit band which I think maybe causing the problem, but when I try to run this without the Landsat 8 images I still get the same error.


texts: Can' update a WMS-source with changed params (updating cql-filter) in OL3 / Geoserver

texts:
I'm updating a WMS-Layer in a OL3-Map.

I tried to reload the wms-layer with updated clq_filter-params.

Can anyone correct my code or give me a working example?

// initialize filtervariable to show all states
var filterstring = "PERSONS>1";

// WMS Layer with cql-filter var filterstring
    var layers = [
            new ol.layer.Image({extent: [-13884991, 2870341, -7455066, 6338219],
            source: new ol.source.ImageWMS({url: 'http://demo.opengeo.org/geoserver/wms',
            params: {'LAYERS':'topp:states','cql_filter': filterstring},
            serverType: 'geoserver'})}
        )
    ];

// OL Karte mit WMS
var map = new ol.Map({
    layers: layers,
    target: 'map',
    view: new ol.View({
    center: [-10997148, 4569099],
    zoom: 4
    })
});

// get Filtervalue from form and replace clq-filter in ol.layer.tile params
document.getElementById("submitBtn").addEventListener("click", func1);
    function func1(){
        var a=document.getElementById("zahl").value;
        var filterstring = "PERSONS>" + a;
    layer[ol.layer.Image].getSource().updateParams({params});
};


texts: Creating simple maps to view on Mobile devices?

texts:
My organization wants to create a few maps that will be interactive and viewable via phones (android, iphones,etc) and ipads.

We already have a solution using Geocortex (Silverlight) for browser which cannot display on most devices.
I've been looking into HTML5, Javascript and ArcGIS Online.

What I need is something that displays nicely for mobile use and that is easy to set up.
Looking at Javascript is slightly intimidating to me and ArcGIS Online maps don't display well in the phone's browser.

I do not come from a Comp Sci or development background so if there's already an existing interface that makes it easily to display for mobile I'd prefer that rather than messing with code.
The map would contain point features with hopefully a pop-up containing a picture and a short description of the property.

I'm pretty computer savvy but I'd like this to be as painless as possible.

Are there any other interfaces that exist that make this easier?

The mobile map will be only for viewing purposing, no editing or capture.

I'm open to looking at other platforms outside of ArcGIS as well.

As an added bonus (but not mandatory), I've been looking into Story Maps and plan on using one in the future, specifically the Short List. However, I don't see this working out for a mobile device as its very picture heavy and would crowd the display. If anyone has any other alternative ideas I'm open to them!

I've been looking at Leaflet and Mapbox as a possible solution for mobile use - there is plenty of documentation - but still a lot of code!

For now we've gained access to HTML5 viewer for geocortex and will be using it for our mobile sites. If there are any other alternatives I haven't already listed above - that are about on-par with ArcGIS Online for ease of use for non-programmers, I would be interested to hear about them.


texts: Why does QGIS hang when I try to create a vector grid?

texts: I am working in OpenLayers plugin. I need to project the satellite data with respect to vector data shapefile I have. Andre Joost suggested to create vector grid with vector data reference system. But when I create the vector grid the system hangs, what is the problem?


texts: In MapBasic, how do I display a dialog box as the result of a menu item being clicked on?

texts:
I have created a menu item in MapBasic.  I need to call a procedure to open a dialog box as the result of one of the menu items being selected.  I believe I am having trouble declaring the sub procedure that will be called.  I have not included the sub procedure, because it is not right.  Here is my code:

Declare Sub Main
Sub Main
    Create Menu "Test" ID "12" As
        "Preferences",
        "Leads",
        "LeadsRemain",
        "Accts",
        "AcctsRemain"
    Alter Menu Bar Add "Test"
End Sub


texts: Converting multilinestring and multipolygon to multipoint in PostGIS?

texts:
In PostGIS I have 2 layers:'lines' (multilinestring, public, gid, geom) and 'polygons' (multipolygon, public, gid, geom). I try to find how to convert them to multipoint (or point). I suspect that ST_CollectionExtract (https://postgis.net/docs/ST_CollectionExtract.html) can do that, but I didn't find any sql example; I tried sth like:

SELECT ST_CollectionExtract (public.polygons as geom,1)

but it doesn't work.

I'm not good at SQL.  I usually try to find some sql example that uses all those things like gid, geom, table_name, column_name etc.


texts:
How are you doing ? I am Mehul Nahar working on a project which is similar to http://www.tableausoftware.com/learning/examples/gpa-analysis. I have used Polymaps platform, where I have used cloudmade API for the base map.

My problem -

I am unable to implement the Geojson file, in fact I am still unable to plot Lat&long. I was wondering if you guys could help me on this. If you can give me a right direction to work on and also if you have some sample code which you could share with will be real helpful.

Look forward to hearing from you.


texts: Imported from Excel XY coordinates into ArcMap 10.4.1 display at zero latitude

texts:
When importing spatial data from excel, the imported point layer is displayed onto 34 latitude and 0 longitude. It should be 47-48 latitude and 37-39 longitude.
What could be causing it?
I tried to change the spatial reference, it does not seem to help.


texts:
I am working with QGIS 3.4

Is there a way to easily copy the styles from one group of layers to another group of layers with the same names? The outcome should be that layers that share the same name, but which are in different groups, end up with the same style.


texts: How to change the Python Version in QGIS?

texts:
I have installed ImageMagick via Homebrew and apparently it's within Python 3.8. But my QGIS is using Python 3.7.3, therefore, it's telling me that I don't have ImageMagick installed (and therefore won't do what I want to do, which is to export an animated GIF using the TimeManager plugin). I'm running OS 10.15.4, QGIS 3.12.2. everything patched up.

So how do I tell QGIS to use Python 3.8 and not 3.7?

PS. I am not a programmer. I would appreciate if someone can educate me what it means to have "ImageMagick installed in Path". I'm assuming it's the path of the version of Python that QGIS is using. so how do I do that? I don't really understand the instruction on ImageMagick web


texts: What is PyQGIS equivalent code of expression "attribute(@atlas_feature, 'FIELD_NAME')"?

texts:
I have a pretty long expression filter so it would be more neat to put it all in a function expression instead. So far, I just found that equivalent PyQGIS within function body tends to be ironically even longer, though.

I'm trying to find equivalent expression in PyQGIS of

attribute(@atlas_feature, 'FIELD_NAME')


texts: QGIS - online Bing basemap shifts when map is exported as a PDF

texts:
Using QGIS 2.18 “Las Palmas”

Created a map to print using the “New Print Composer” extension.

When working on this all polygons are perfectly in place.

Once the map is saved as a PDF/Image the polygons shift.

Extra info: Bing map (Base Map from Open Layers Plugin) and all polygons are set to ESPG:29902 TM65/Irish Grid


texts: Differences between different georeferencing transformation algorithms in QGIS

texts:
Expanding on Transformation types in geo-referencing of QGIS?, I'm trying to understand the exact transformations used by QGIS, and how they differ. The georeferencer documentation is of course a good starting point, but I need more details :

Currently, the following Transformation types are available:
The Linear algorithm is used to create a world file and is different
from the other algorithms, as it does not actually transform the
raster. This algorithm likely won’t be sufficient if you are dealing
with scanned material.
The Helmert transformation performs simple scaling and rotation
transformations.
The Polynomial algorithms 1-3 are among the most widely used
algorithms introduced to match source and destination ground control
points. The most widely used polynomial algorithm is the second-order
polynomial transformation, which allows some curvature. First-order
polynomial transformation (affine) preserves collinearity and allows
scaling, translation and rotation only.
The Thin Plate Spline (TPS) algorithm is a more modern georeferencing
method, which is able to introduce local deformations in the data.
This algorithm is useful when very low quality originals are being
georeferenced.
The Projective transformation is a linear rotation and translation of
coordinates.

Specifically, what are the differences between linear, Helmert and projective?
The "linear" algorithm creates a world file, that includes 6 parameters allowing to define translation (lines 5 and 6), scaling (1 and 4) and rotation (2 and 3). Although, I remember reading somewhere that QGIS does not honour rotation in world files (?), but I cannot trace the source of this info - empirically, it certainly feels as if the "linear" fit in the georeferencer does not rotate the image.
The "Helmert" algorithm "performs simple scaling and rotation" (and translation, presumably, as you can use it to "move" an image to any coordinate)
The "Projective" transformation "is a linear rotation and translation of the coordinates".
I must confess, I cannot understand the differences between the above three.
In addition, a first-order polynomial "allows scaling, translation and rotation only".
In all cases, it seems to me that the algorithm uses fitting of the 6 parameters (2 each for translation, rotation and scale) to minimize the differences between GCP and target positions. Apart from different fitting strategies (plain least-square, vs. more robust fitting, outlier-tolerant or so), I do not see the differences, and in any case with 6 GCP the matrix being square, there should be an unique solution common to all.
On the other hand, I understand indeed how a second-order polynomial third-order polynomial and spline differ, of course. But this really gives only 4 strategies (5, if indeed linear does not do rotation): first order without rotation, first order, second order, third order and spline, for the 7 options available.


texts:
I am looking to convert some Geosoft .gdbs for use in QGIS. The .gdb files are in Geosoft .gdb format and not in Esri File Geodatabase format. I do not have Geosoft installed (only QGIS) so am looking for a conversion workaround. I get errors when trying to import into QGIS as a vector/directory.
data here: http://www.geosciencebc.com/major-projects/quest/


texts: Arcpy ListDatasets function to specific directories instead of workspace

texts:
Are there any alternatives to the arcpy list functions (ListDatasets, etc.) where the directory can be specified rather than using the env.workspace?

For example, what if I would like to get a list of datasets in several different directories?

Or if I have my workspace in_memory?

The os.walk might partially work, but won't recognize non-file-based data types.  Is another function or snippets that can accomplish this?


texts: Converting from Julian Date to Calendar Date?

texts: Does anyone know of a way to convert from a Julian date to a regular calendar date?


texts:
I would like to export an image from ArcMap to a text document. My problem is that, as you can see, my image

does not show all the colors and the areas, because the geographic area is too big and the polygons are too small. As a consequence, large part of the geographic area I am interest in looks just brown, because of the polygons borders overlapping.

My question is: is there a way to improve the image export in order to make all the colors visible, taking into account that this image must fit in a normal A4 paper page?


texts: I'm using a multipolygon PostGIS layer in QGIS 3 and created a rule-based labeling. I'd like to know if when a label is moved by move label and diagram will it be possible to update the label by changing the entry in the attribute table?


texts: GeoServer - Show Maps in a Application Java

texts:
How can I display a GeoServer's Map inside a Java's Desktop Application?

I would like to display it inside a Java App instead of a browser. Is this possible?


texts: Configurable Map Viewer (CMV) - Content Pane with CMV

texts:
I’m upgrading from the old viewer to the new CMV version.

Is there a better way for handling a content pane in the sidebar under the new format? In the previous version I did it thru the controller, but this seems to be against the best practices of the new CMV. Any suggestions?

This is what I was doing: (Last portion of the posted code from Controller.js creates the content pane pictured below circled in red and adds it to the 5th position.)

Controller.js

    initWidgets: function(evt) {
        var widgets = [];

        for (var key in config.widgets) {
            if (config.widgets.hasOwnProperty(key)) {
                var widget = lang.clone(config.widgets[key]);
                if (widget.include) {
                    widget.position = ('undefined' !== typeof(widget.position)) ? widget.position : 10000;
                    widgets.push(widget);
                }
            }
        }

        widgets.sort(function(a, b) {
            return a.position - b.position;
        });

        array.forEach(widgets, function(widget, i) {
            this.widgetLoader(widget, i);
        }, this);

        var wcCP = new ContentPane({
            content: "

Divison of Oil & Gas
2045 Morse Rd
Columbus, OH 43229
614-265-6922
oilandgas@dnr.state.oh.us
ODNR Oil & Gas Website

Download GIS Oil & Gas Well Data:
Download Well Data by County"
        }).placeAt(this.sidebar, 5);
    },


texts:
I'm using ArcObjects .NET to load WMS 1.3.0 service layers from our ArcGIS server into ArcMap 10.2.2. The process of connecting to the server and adding layers works fine, however none of the layers will draw - I simply have a blank canvas. I have tried two methods to add WMS layers to the map:

Using 'Add Data' button > GIS Server > Add WMS Server
With ArcObjects .NET (see code below)

    'setup connection properties
    Dim pMap As IMap = pMxDoc.ActiveView.FocusMap
    Dim pWMSLayer As IWMSGroupLayer = New WMSMapLayer
    Dim pConnectionName As IWMSConnectionName = New WMSConnectionNameClass
    Dim pPropSet As IPropertySet = New PropertySetClass
    pPropSet.SetProperty("URL", txtURL.Text)
    pConnectionName.ConnectionProperties = pPropSet
    Dim pDataLayer As IDataLayer = CType(pWMSLayer, IDataLayer)

    'connect to server
    Try
        pDataLayer.Connect(CType(pConnectionName, IName))
    Catch ex As Exception
        MsgBox("Error connecting to server: " & ex.Message)
    End Try

    'get top level group layer from server
    Dim pServiceDesc As IWMSServiceDescription = pWMSLayer.WMSServiceDescription
    Dim pLayerDesc As IWMSLayerDescription = pServiceDesc.LayerDescription(0)
    _WMSGroupLayer = pWMSLayer.CreateWMSGroupLayers(pLayerDesc)

    'add layers to map
    For i As Integer = 0 To _WMSGroupLayer.Count - 1
        Dim pLayer As ILayer = _WMSGroupLayer.Layer(i)
        pMap.AddLayer(pLayer)
    Next

Both methods add the layers to the map and they are visible in the table of contents with no errors or warnings. I can view layer source (seems correct), extent (also seems correct), some symbology properties, and I can also zoom to the extent of the layer (which changes map zoom), however nothing at all is visible on the map canvas. Scale range isn't the issue as I've tried zooming in far beyond the specified range.

I've come across this page relating to WMS layers not drawing in ArcMap 9, but it isn't offering me any working solutions for 10.2:

WMS layers won't draw in 9.3

I have tried changing coordinate system of data frame (and layers), no to effect. They are all in WGS1984 which is supported by the WMS service.

At this point I don't know what else I can try besides contacting ESRI for assistance. Does anyone have any ideas or similar issues with WMS layers? Thanks.

EDIT: also looks like this person had the same problem but no solution was posted.


texts: Bathymetry map (DEM) and the reservoir Volume from bathymetry points

texts:
Could anyone provide me with a workflow from POINTS --> Bathymetry map (DEM) and  Volume calculation in ArcGIS?

I would like to create a bathymetry map (DEM) and calculate the volume for a large riverine reservoir (Lake Diefenbaker). I have the bathymetry depth points.

The first issue is how can I restrict the interpolation area to the shoreline areas?

Because points are along the river valley which it divides into a fork at the end (T shape reservoir, attached photo), I can’t simply create a TIN surface without creating a boundary around the points.


texts: To Use or Not to Use Relationship Classes?

texts:
I am building a government database from a massive amount of data. Now I am following the standard approach of the ESRI Local Government Data Model, but modified to fits the needs and wants of a specific government.

I've noticed that in the Local Government Data Model ESRI uses a lot of relationship classes. I've used relationship classes in the past, but tend to steer clear of them. This is because I usually don't see a reason why I would want attributes about a shape in a separate table from the feature class.

For almost all of the data I am assembling I still cannot see a reason for this. For Parcel Fabric I can see some reasons where it would be helpful, but still I believe it would be best to still keep the attributes directly in the feature class.

Now to be clear I am going to use relationship classes to reference feature classes to feature classes, this question is more specified to why I would store data in a non-spatial table that could just as easily be store in a feature class.

Anyways thanks for any input folks!


texts:
There has been an error when uploading a shapefile

I want to start with some imported data. Upload a file.

Please, select a valid format file

# when upload a zip file:

  Oops! There has been an error
  Unknown
  Code 99999
      Begin:
      TRANSITION: preprocessing => uploading, 2014-01-07 02:59:15 -0800
      getaddrinfo: Name or service not known
      Table error, #


texts: Which OpenLayers highlights should I present to a non-specialist audience?

texts:
There is an event in my University that interests in Open Source Software. I am supposed to give a short presentation about anything related so I thought about giving a simple overview about OpenLayers Library.

Would you please give me some highlighted points that really matter non-specialist audience?


texts: Pythonic way to set extent and save maps in a cycle

texts:
I got number of QgsRectangles and i want to save images bounded by them.

Code written in python:

    for box in self.boxes:
        xs = [box[0][0],box[1][0],box[2][0],box[3][0]]
        ys = [box[0][1],box[1][1],box[2][1],box[3][1]]
        min_x = min(xs)
        min_y = min(ys)
        max_x = max(xs)
        max_y = max(ys)
        box = QgsRectangle(min_x,min_y,max_x,max_y)
        self.iface.mapCanvas().setExtent(box)
        self.iface.mapCanvas().refresh()
        time.sleep(1)
        self.utils.iface.mapCanvas().saveAsImage('%s.tiff' %(str(min_x)))

The problem is script walks over all my boxes, but extents are set only for the last one. I thought qgis will set extents, show location of my box and save an image associated with it.


texts: Dynamic defining of transformation type for ISpatialReferenceFactory.CreateGeoTransformation call

texts:
Suppose given a point (IPoint) in some coordinate system. Required to transform it to another coordinate system. If both of coordinate systems are predefined, we can use IGeometry2.ProjectEx with transformation specified by corresponded constant from esriSRGeoTransformationType, esriSRGeoTransformation2Type, or esriSRGeoTransformation3Type enumerations.

The code looks like (source)

package arcgissamples.geometry;

import com.esri.arcgis.geometry.*;
import com.esri.arcgis.system.*;

public class ProjectGeometry {

  public static void main(String[] args) {
    try {
      EngineInitializer.initializeEngine();
      AoInitialize aoInit = new AoInitialize();

      // Create a point with Geographic coordinates...
      Point point = new Point();
      point.putCoords(-100.0, 40.0);

      System.out.println("");
      System.out.println("Original coordinates: " + point.getX() + ","
                + point.getY());

      // Create the SpatialReferenceEnvironment...
      SpatialReferenceEnvironment spatialReferenceEnvironment = new SpatialReferenceEnvironment();

      // Apply the initial spatial reference...
      ISpatialReference geographicCoordinateSystem = spatialReferenceEnvironment
                .createGeographicCoordinateSystem(esriSRGeoCSType.esriSRGeoCS_NAD1927);
      point.setSpatialReferenceByRef(geographicCoordinateSystem);

      // Create the output projected coordinate system...
      ISpatialReference projectedCoordinateSystem = spatialReferenceEnvironment
                .createProjectedCoordinateSystem(esriSRProjCSType.esriSRProjCS_NAD1983UTM_13N);

      // Create the GeoTransformation...
      //***************************************************
      // !!!! How to define arg for createGeoTransformation??
      IGeoTransformation iGeoTransformation = (IGeoTransformation) spatialReferenceEnvironment
                .createGeoTransformation(esriSRGeoTransformationType.esriSRGeoTransformation_NAD1927_To_WGS1984_5);

      // Project the point...
      point.projectEx(projectedCoordinateSystem, esriTransformDirection.esriTransformForward,
                iGeoTransformation, false, 0.0, 0.0);
      System.out.println("Projected coordinates: " + point.getX() + " , "
                + point.getY());

      System.out.println("Done!");

      aoInit.shutdown();

    } catch (Exception ex) {
      ex.printStackTrace();
    }
  }
}

But problem is that both original coordinate system and output coordinate system are defining at the time of program execution. So I need a function to define constant for ISpatialReferenceFactory.CreateGeoTransformation method like

int defineTransformationType(ISpatialReference srcCoordSystem, ISpatialReference destCoordSys)

Is it possible? And how to implement this function?


texts: How do I convert a USGS Geotiff from NAD83 into ECEF?

texts:
I have downloaded elevation data from USGS, such as this file covering Atlanta.
The file is a Geotiff in the NAD83 coordinate system. I would like to do some geometry with this that requires ECEF. I have tried a variety of gdalwarp commands, including:
gdalwarp -t_srs EPSG:4978 USGS_1_n34w085.tif ecef_n34w085.tif

which gives this message:
Creating output file that is 4515P x 3206L.
Processing USGS_1_n34w085.tif [1/1] : 0Using internal nodata values (e.g. -999999) for image USGS_1_n34w085.tif.
Copying nodata values from source USGS_1_n34w085.tif to destination ecef.tif.
...10...20...30...40...50...60...70...80...90...100 - done.

Indeed, the output is a Geotiff where every point is just -999999. So I tried this instead to see if it's a problem outputting ECEF in Geotiff:
gdalwarp -t_srs EPSG:4978 USGS_1_n34w085.tif ecef_n34w085.xyz

But the result is the same. Here are the first few points in the XYZ file:
461186.393299947202 -5264153.32883292157 -999999
461208.248676233809 -5264153.32883292157 -999999
461230.104052520415 -5264153.32883292157 -999999
461251.959428806964 -5264153.32883292157 -999999
461273.814805093571 -5264153.32883292157 -999999

What am I doing wrong? Is there something inherently incompatible with the input data, or with gdalwarp trying to output ECEF? I don't really care whether I have to use gdalwarp, but if I manually iterate all 13 million points through the GDAL library, it takes a long time.


texts: Unable to load an OSM .bz2 file into a PostGIS database using osm2pgsql

texts:
Following is the query :
C:>osm2pgsql -c --slim -d osm_1 -U postgres -H localhost -S C:\default.style C:\mumbai_india.osm.pbf

Following is the error on which i am stuck:

Using projection SRS 900913 (Spherical Mercator)
Setting up table: planet_osm_point
NOTICE:  table "planet_osm_point" does not exist, skipping
NOTICE:  table "planet_osm_point_tmp" does not exist, skipping
PREPARE get_way (int4) AS SELECT AsText(way) FROM planet_osm_point WHERE osm_id
= $1;
 failed: ERROR:  function astext(geometry) does not exist
LINE 1: PREPARE get_way (int4) AS SELECT AsText(way) FROM planet_osm...
                                         ^
HINT:  No function matches the given name and argument types. You might need to add explicit type casts.

Error occurred, cleaning up

Following is the version of osm2pgsql:

C:\>osm2pgsql
osm2pgsql SVN version 0.69-21289M

I am using Windows 7 64bit PC with PostgreSql database installed with pgAdmin4..
I did pass various  osm.pbf files in the above query but still getting the same issue.
I went through various post and solutions and also created Extension postgis and Extension hstore in the database named osm_1 in pgAdmin4 postgresql but still the same issue..

I re-installed osm2pgsql with the latest version found on Google but still got the same error..


texts: QGIS Symbology: how to make a symbol a shape according to a set of rules and a color according to another set of rules

texts:
I'm using QGIS 3.14 with GRASS 7.8.3. I have a data set with point coordinates on a plan (artifacts on an archaeological site). I have two main variables that I want to show. So, I'd like to have the type of artifacts as a letter (that I managed with the set of rules) and a color depending of the raw materials used for the artifact.
For example:

Biface, Chert
Biface, Quartz
Arrowhead, Jasper
Flake, Chert

I'd like to have all the bifaces with a B, and all the artifacts in chert showing up in grey. And that for all my categories (I have a dozen for each variable) and my data base includes hundreds of entries. I'd like a solution that is not too tedious if it possible.


texts: Is there an automatic way to create a diverging colour scale centered around zero?

texts:
I'm trying to plot a map of the residuals of a linear model. My residuals centre around zero and can obviously go in either direction.

Is there an easy way to apply a graduated colour scale to these values in QGIS, one that centres on zero and then shades in red above zero and blue below (or whatever colour ramp you choose) without having to manually set the upper and lower limits of each category?


texts: ERROR 999999 in ModelBuilder because of Undefined Spatial Reference?

texts:
I have a folder with around 300 Shapefiles I want to use as masks in the tool Extract by Mask on the same Raster image. I'm using the following model :

But I get the following error :

ERROR 999999: Error executing function.
The table name is invalid.
Failed to  execute (Extract by Mask)

I'm pretty sure it's because the new raster image loses the original image's spatial reference :

The original image is in GCS WGS 1984.
How can I solve this error?


texts:
I am using a basemap (http://basemaps.cartocdn.com/light_nolabels/{z}/{x}/{y}.png) in QGIS to draw a map for a scientific paper.
I am trying to understand who is the "owner" of this map and who should I credit to.
Also I would like to know what license this basemap has (CC by 4.0,others?)
If you open the basemap on your browser at http://basemaps.cartocdn.com/#11/50.2342/-66.4254
you can see at the right bottom corner

"Leaflet| (c) OpenStreetmap contibutors, (c) rastertiles/voyager"

when clicking on "(c) OpenStreetmap" I am sent to OpenStreetmap website, clearly.
when clicking on "(c) rastertiles/voyager" I am sent to CARTO.com.
Is it CARTO.com or OpenStreetMap ?
PS: I wrote to CARTO.com via their website form but the only thing I got was a marketing mailing list inscription...


texts: Installing QGIS with Python 2 and Python 3 on OSX

texts:
I just installed qgis (and all the dependencies, etc.) on my Mac OSX 10.10 via homebrew.  Currently, I use both Python 2 and Python 3 on the machine.  Whenever I change the Python path [to open qgis], I can no longer use Python 3 in any of my py3 virtual environments until I unset the python path.

Has anyone found a good solution to having both Python2 + Python3 on their machine with qgis?


texts:
I'm writing a script in python that loops the LayerToKML conversion tool on a set of layers; the problem is that the first execution works fine, while the rest (hundreds) produce .kmz files like this...

...all the geometries are invisible and placed at the Equator or North Pole, while the attributes are correct.

I run the script with an IDE (PyScripter, 32 Bit, python 2.7.2): if I do it with the ArcGIS (10.1) python window, it works instead!

PS: athis is the code function that I use to convert the shapefiles:

def creatingKmz(fcKmz,fcName):

    arcpy.MakeFeatureLayer_management(fcKmz, fcName)
    arcpy.ApplySymbologyFromLayer_management(fcName,inFolder + os.sep + fc + ".lyr")
    arcpy.LayerToKML_conversion(fcName,outKMZDirectory + os.sep + fcName + ".kmz","1","false",arcpy.Describe(fcName).extent,"","","CLAMPED_TO_GROUND")
    arcpy.Delete_management(fcKmz)


texts: Summing values that are geographically close to each other using Carto?

texts:
I have a dataset of Boko Harams rampage in Nigeria uploaded to CartoDB.

I have columns like coordinates, dates, number of victims.

Would it be possible to write a query that would sum the victims within a given radius?

If attack x happens within 10 km from attack y, add the number of victims of attack x to the number of victims of attack y.


texts: Loop through time values in a raster brick in R

texts:
I have raster bricks of temperature and rainfall consisting of one layer for each month of the year.  The raster brick looks like this:
> tmx
class      : RasterBrick
dimensions : 3875, 8017, 31065875, 240  (nrow, ncol, ncell, nlayers)
resolution : 0.04490319, 0.04490319  (x, y)
extent     : -180.0393, 179.9495, -90.04088, 83.95898  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
source     : /scratch/bntjoa002/chelsa_cmip5_ts/tasmax/CHELSAcmip5ts_tasmax_ACCESS1-3_rcp8                                      5_2070-2089.nc
names      : X2070.01.15, X2070.02.15, X2070.03.15, X2070.04.15, X2070.05.15, X2070.06.15,                                       X2070.07.15, X2070.08.15, X2070.09.15, X2070.10.15, X2070.11.15, X2070.12.15, X2071.01.15                                      , X2071.02.15, X2071.03.15, ...
Date       : 2070-01-15, 2089-12-15 (min, max)
varname    : air_temperature

I have been using the following code to get the average monthly values for a specified period (e.g. resulting in 12 averaged layers, one for each month, over a X year period):
new_tmx <- subset(tmx, which(getZ(tmx) >= as.Date("2070-01-15") & getZ(tmx) <= as.Date("2089-12-15")))
indices <- format(as.Date(names(new_tmx), format = "X%Y.%m.%d"), format = "%m")
indices <- as.numeric(indices)
Monthlytmx<- stackApply(new_tmx, indices, fun = mean)

I do this separately for the maximum temperature, minimum temperature, and precipitation bricks.
What I need now is to do the same process, but instead to output one raster layer for each year consisting of the averaged values for one-year periods (i.e. for 2070, 2071, 2072...).
I've only done very simple loops in R, and I'm not sure how to deal with the time/date function in this example?


texts: Creating a distance map from OpenStreetMap data

texts:
I would like to create a distance map using OpenStreetMap objects as input objects. To explain what I mean with distance map, let's look at the wikipedia definition:

  A distance transform, also known as distance map or distance field, is a derived representation of a digital image. [...]
  The map labels each pixel of the image with the distance to the nearest obstacle pixel.

Imagine, for example, one wants to create a distance map for Sweden, indicating the distance for each pixel to the nearest alcohol takeway (OSM key: shop, OSM value: alcohol). This would require the following steps:

Gather all OSM data for the area inside the country borders of Sweden
Define pixel spacing of desired raster map (e.g. 1m or 1 km etc.)
Calculate distance for every pixel of that map to nearest alcohol shop OSM map feature using a to-be-defined distance metric (e.g. Euclidian distance)

While I have a clear conceptual idea of how to do it, I don't know where to start implementation-wise.

Is there a good Python API that allows to do such a thing without too much overhead?

Or is it easier to achieve my goal in QGIS?


texts: Investigating slow performance when opening tables in ArcGIS Desktop?

texts:
My ArcGIS 10.5 performance is very slow just in opening some thinks like: Options, tables, labels, arctoolbox, draw menu ect. It open the mentioned things in about 3 minutes.
I have tried different solutions like: l deleted Normal.mxt, remove and reinstalled the ArcGIS, restore my PC and tried more solutions that I found in Internet but no one was useful to solve the problem.
My Computer  is Dell 5559 cori7 2.5GHz 4cpus, 8GB ram and 4GB Graphic dedicated memory
This is bothering me. What is the problem?


texts: Performance of joining tables on a geometry column with SQL Server

texts:
I had a shapefile that I needed to integrate into my SQL Server DB. To do so, I used the very useful Shape2Sql.exe with which I chose 'Set SRID 4326' in the settings.

I now have a table that contains a geometry column while another table is showing to me the latitude and longitude of different addresses along with other useful information. This can be summed up as follows:

CREATE TABLE dbo.myTable1 (
       FULL_ADDRESS varchar(120) NULL,
       LONGITUDE varchar(50) NULL,
       LATITUDE varchar(50) NULL
                           )

CREATE TABLE dbo.myTable2 (
       POLYGON_NAME nvarchar(255) NULL,
       geom geometry NULL
                           )

What I would like is to join the two tables based on the sets of coordinates that match a certain polygon. I took some inspiration from here and I basically did this:

SELECT DISTINCT
      mt1.FULL_ADDRESS
    , mt1.LONGITUDE
    , mt1.LATITUDE
    , mt2.POLYGON_NAME
INTO #TEMP_TABLE
FROM myTable1 mt1
LEFT JOIN myTable2 mt2
  ON geom.STIntersects(geometry::STGeomFromWKB((geography::Point(mt2.LATITUDE, mt2.LONGITUDE, 4326)).STAsBinary(), 4326)) = 1

(the DISTINCT is useful because, for some reason that I can't explain yet, I get several matches for a single set)

The results thereof are exactly what I expect but, for 10,000 sets of coordinates, it almost took 3 hours. Problem is, I'm supposed to find the polygons of 25,000,000 sets (which means 300 days of computation...). I therefore need to find a trick to improve that query performance by several orders of magnitude.

BONUS QUESTION:

When I don't store the results of my query inside a temp table (i.e. when I want to see the results), I get the following error message:

  A .NET Framework error occurred during execution of user-defined routine or aggregate "No extended information available.": .

But everything is fine when using a SELECT... INTO #TEMP_TABLE instead of a plain and simple SELECT. I couldn't find any piece of information about this message. Any idea to solve it?


texts: Representing multiple polylines in unique vector layer?

texts:
I have imported a layer which contains a lot of lines. They overlap each other in a cross kind of pattern. When zoomed out it can look like this:

When zooming in on one of the broader lines it looks like below, where the thicker black line indicates my wanted output?

Is there a way of identifying what lines that overlap each other and merge them to one line that stretches across where they are? It is also OK with a straight line if the curved is not doable:

The background is that the lines comes from a tool called Kartapullautin and are representing cliffs from LIDAR data. I want them to be merged so I can easily add it to a map and not have thousands of small cliffs on top of each other, but rather have a few simplified one. Both from readability purpose, but also because the large amount of cliffs are growing my map project in file size so it is difficult to work with. It is worth to mention, I am QGIS newbie!


texts: Getting Area of Polygon Intersections Stratified by Categorical Variable in PostGIS

texts:
I've got two layers: counties and forest type area. The county layer is what you'd expect: a bunch of simple polygons representing counties in Texas. The forest type area layer looks like Tetris shapes, and consists of "geom" and "type" attributes. See the attached screenshot for an example of what it looks like.

What I'm trying to do is get the total area of each "type" from the forest type layer that falls into each county. I've already intersected the forest type and county layers so that forest type is split where it intersects county lines.

The results should look something like this:

|---------------------|------------------|----------------|
|County Name          |Type 1 Area       |Type 2 Area     |
|---------------------|------------------|----------------|
|          Tarrant    |         565.27   |     2737.28    |
|---------------------|------------------|----------------|

I've taken a stab at this before, before I had to put this project on hold for about a year. Now I'm revisiting it and cannot figure it out. Here's my previous attempt:

/*
Take forest fragmentation table and calculate the sum area of each patch type as it intersects with counties. Save the results to a new PostGIS table
*/

CREATE TABLE countyFragStats AS SELECT * FROM texascounties;

UPDATE countyFragStats
        SELECT texascounties.geom FROM texascounties,  coalesce(sum(area_1992_w3_decid.area),0) AS areaSum
        FROM area_1992_w3_decid
        LEFT OUTER JOIN texascounties
        ON ST_Intersects(area_1992_w3_decid.geom, texascounties.geom)
        GROUP area_w3_1992_decid.id
        LIMIT 1000;

It's a non-starter, giving me

LINE 2:  SELECT texascounties.geom FROM texascounties,  coalesce(sum...

Any ideas?


texts: How to obfuscate data when using Google Maps?

texts:
I've been working on a google map (v3 of the api) that is plotting 120 markers or so (or will do shortly) - see http://www.mediwales.com/mapping. If you view the source the map data is for all the see, is it possible to hide this?

I'm not worried about the code that is generating the map, just the data. The data is grabbed from a Wordpress cms.

Here's the code that's generating it all:


texts: How to merge polygons into a single item where an attribute is same?

texts:
I have a shapefile of an area. It consists of a list of fields in which sub-locality is also a field. Under the sub-locality field all the areas are common, say for example: perungudi, perungudi, perungudi.
When I export this shapefile in Post Gis in the sub-locality field each attribute is represented as a separate shapefile. So I want all the attributes in the sublocality field to have common boundary. How can I do this?


texts: OpenLayers Full Screen control not working

texts:
I have a problem with OpenLayers controls. I have two controls which are perfectly working,

var mousePositionControl = new ol.control.MousePosition({
        coordinateFormat: ol.coordinate.createStringXY(6),
        projection: 'EPSG:4326',
        // comment the following two lines to have the mouse position
        // be placed within the map.
        className: 'custom-mouse-position',
        target: document.getElementById('mouse-position'),
        undefinedHTML: 'Haritada degilsiniz.'
    });
    var scaleLine = new ol.control.ScaleLine();

...

controls: ol.control.defaults().extend([mousePositionControl, scaleLine])

But when I try add fullscreen control as a third one, like this:

var fullScreen = new ol.control.FullScreen();

and edit controls like this:

controls: ol.control.defaults().extend([mousePositionControl, scaleLine, fullScreen])

it's not working and also I lose map.on('click', function(event)){...} functions.

By the way I try to add only FullScreen control, it works. Have a limitation there of adding controls on OpenLayers 3?


texts: How to determine when Mapbox GL JS FlyTo has "arrived"

texts:
I'd like to show an overlay when the Mapbox flyto camera movement has completed its action, and is in place at the correct position and zoom level.  Is there a way to know when that action is complete?

So, fly to a position and then show the overlay, basically.


texts: How can I project this data correctly in node-mapnik?

texts:
I'm trying to use node-mapnik to rasterize vector data into a transparent image and overlay it onto a basemap.

I'm generating OpenStreetMap basemaps with tile-stitch.  That's working well.

I'm then using node-mapnik to read vector data from a geojson file and render to an image.  I'm getting the image okay, but the projection is off.

The coordinates are all EPSG:4326, decimal notation, derived from a slippy map.
My mapnik xml stylesheet looks like this

<Map background-color='#00000000' srs='+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs'>
    <Style name='test style'>
        <Rule><PolygonSymbolizer fill='#bb5acc66' /><LineSymbolizer stroke='#3ea3df66' stroke-width='1' /></Rule>
    </Style>
    <Layer name='test' srs='+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs'>
        <StyleName>test style</StyleName>
        <Datasource>
            <Parameter name='type'>geojson</Parameter>
            <Parameter name='file'>test.geojson</Parameter>
        </Datasource>
    </Layer>
</Map>

My geojson looks like this:

{"type": "FeatureCollection",
"crs": { "type": "name", "properties": { "name": "urn:ogc:def:crs:OGC:1.3:CRS84" } },
"features": [
{ "type": "Feature", "properties": { }, "geometry": { "type": "Polygon", "coordinates": [ [ [ -1.182891726493835, 51.630744982473551 ], [ -1.183235049247742, 51.630671725935755 ], [ -1.183497905731201, 51.630445295887 ], [ -1.183481812477112, 51.630258823233 ], [ -1.183358430862427, 51.630108978580999 ], [ -1.183213591575622, 51.630029061230999 ], [ -1.182989378008187, 51.630094380903088 ], [ -1.18303120136261, 51.630212204949991 ], [ -1.183058023452759, 51.630352059655998 ], [ -1.182972192764282, 51.630468604915002 ], [ -1.18278980255127, 51.630571830465001 ], [ -1.182720065116882, 51.630675055780998 ], [ -1.182703971862793, 51.630718343745997 ], [ -1.182891726493835, 51.630744982473551 ] ] ] } },
{ "type": "Feature", "properties": { }, "geometry": { "type": "Polygon", "coordinates": [ [ [ -1.181958317756653, 51.630367044063 ], [ -1.182041466236115, 51.630565170758999 ], [ -1.182301640510559, 51.630656741628989 ], [ -1.182315051555634, 51.630578490169995 ], [ -1.18239551782608, 51.630508563219998 ], [ -1.182478666305542, 51.630385358331999 ], [ -1.182626187801361, 51.630297116786998 ], [ -1.182596683502197, 51.63016059179499 ], [ -1.182666420936584, 51.629989102503011 ], [ -1.182325780391693, 51.629940818994001 ], [ -1.182111203670502, 51.629947478791003 ], [ -1.1820387840271, 51.630093994089002 ], [ -1.181958317756653, 51.630367044063 ] ] ] } },
{ "type": "Feature", "properties": { }, "geometry": { "type": "Polygon", "coordinates": [ [ [ -1.182666420936584, 51.629989102503011 ], [ -1.18283862812109, 51.630068380938702 ], [ -1.182989378008187, 51.630094380903088 ], [ -1.18303120136261, 51.630212204949991 ], [ -1.183058023452759, 51.630352059655998 ], [ -1.182972192764282, 51.630468604915002 ], [ -1.18278980255127, 51.630571830465001 ], [ -1.182720065116882, 51.630675055780998 ], [ -1.182703971862793, 51.630718343745997 ], [ -1.182395517689429, 51.630700862234136 ], [ -1.182301640510559, 51.630656741628989 ], [ -1.182315051555634, 51.630578490169995 ], [ -1.18239551782608, 51.630508563219998 ], [ -1.182478666305542, 51.630385358331999 ], [ -1.182626187801361, 51.630297116786998 ], [ -1.182596683502197, 51.63016059179499 ], [ -1.182666420936584, 51.629989102503011 ] ] ] } }]}

The data on the slippy map looks like this:

but the mapnik image looks like this:

I'm giving tile-stitch and mapnik the same bounds, and the same image dimensions.


texts: Controlling processing order in FME workspace?

texts:
I have created a workbench with 3 bookmarks:

a Creator/ProjectWise connector  which downloads an Excel file from ProjectWise
a processing workflow
a second Creator/ProjectWise connector which initiates the upload of the spreadsheet

The idea being that I can write information about the data transformation step (2) into the spreadsheet and then upload it. The problem is that I can't control the sequence of events very well which means that the final step (step(3) - upload) occurs before my reader/transform/writer process (2) has had a chance to update the Excel file. I've tried setting the Creator order (no impact) and the tried a decelerator on Creator 2 but with the latter the whole workbench gets stalled.

Is there a way to control the workbench processing order without using multiple workspaces?


texts: Fast minimum distance between two sets of polygons

texts:
I'm trying to calculate the minimum distance between a set a polygons, and a subset thereof. I've used a quick-and-dirty GeoPandas apply:

df['subset_distance'] = df.geometry.apply(lambda g: df_subset.distance(g).min())

which works, but it's pretty slow, even on my test dataset of around 700 polygons. I need to find the minimum distances for around 2 million polygons and a subset of 4k polygons, so I suspect I'm going to need a different strategy. Could anyone suggest a better approach? I'm comfortable with PostGIS, and any combination of Python spatial libraries.


texts: Looping Dataframe in Python doesn't Return Expected Results

texts:
I'm creating a python addin script to be used on some mxd's I have created as templates. Each map contains 2 dataframes. I need to loop through each dataframe and do some processing on each. I set my loop as follows:

for df in arcpy.mapping.ListDataFrames(mxd, "DFMAIN"):
                    for lyr in arcpy.mapping.ListLayers(mxd, "", df):
                           # do some processing here
for df in arcpy.mapping.ListDataFrames(mxd, "DFINSET"):
                    for lyr in arcpy.mapping.ListLayers(mxd, "", df):
                           # do some processing here

In this scenario, if the second dataframe is the 'active' dataframe, it gets chosen first and then some weird stuff happens. I get layers from my first dataframe copied into the second which baffles me and I am not even using any functionality that would copy layers (i.e. arcpy.mapping.AddLayer or anything similar) I just update some definition queries.

If I change the template so that the first dataframe is the active one, then everything works fine.

I tried indexing the dataframe loops:

for df in arcpy.mapping.ListDataFrames(mxd, "DFMAIN")[0]:

but I get the same weirdness if the second dataframe is active.

I guess the easy answer is, make sure the template has the first dataframe as the active one and save it that way, but I have users who may change this unknowingly. I want the geoprocessing to take care of the active dataframe issue. I can't seen anything in the help menu to set the dataframe as the active one.


texts: Loading raster in postgis via raster2pgsql

texts:
I need to import data file to my database.
I created a web interface in order to upload file
but I have already an error

  Warning: pg_query(): Query failed: ERROR: rt_raster_from_wkb: wkb size
  (0) < min size (61) LINE 1: INSERT INTO
  eau_param_hydro(rast,filename) VALUES('','testte... ^ in
  C:\wamp\www\test\test.php on line 32

My code is :

<form method="post" action="" enctype="multipart/form-data">
     Ajouter raster :

     <input type="file" name="fichier"  />
     <input type="submit" name="submit" value="Envoyer" /> </form>
       <?php $name = $_FILES['fichier']['name']; $db = pg_connect('host=localhost port=5432 dbname=gisdb user=postgres >password=******');           $output =shell_exec("raster2pgsql -a -s 32629 -I -C -F -M '".$_FILES["fichier"]["tmp_name"]."' eau_param_hydro | psql -U postgres -d gisdb -h localhost -p 5432 ");         $req = exec(" SELECT * FROM '".$output."' ");
    $result = pg_query($db,"INSERT INTO eau_param_hydro(rast,filename)  >VALUES('".$req."','".$name."');");
         if (!$result) {
         $errormessage = pg_last_error();
         echo "Error with query: " . $errormessage;
         exit();
        }  pg_close();  ?>

in PostGIS I have a table :

  CREATE TABLE eau_param_hydro (   rid serial NOT NULL,   rast raster,
  filename text,   CONSTRAINT eau_param_hydro_pkey PRIMARY KEY (rid ),
  CONSTRAINT enforce_num_bands_rast CHECK (st_numbands(rast) = 1),
  CONSTRAINT enforce_srid_rast CHECK (st_srid(rast) = 4326) ) WITH (
  OIDS=FALSE ); ALTER TABLE eau_param_hydro   OWNER TO postgres;

  -- Index: eau_param_hydro_st_convexhull_idx

  -- DROP INDEX eau_param_hydro_st_convexhull_idx;

  CREATE INDEX eau_param_hydro_st_convexhull_idx   ON eau_param_hydro
  USING gist   (st_convexhull(rast) );


texts: Iterating intersect between polylines and polygons in QGIS using Python script?

texts:
I have a shapefile made up of 127 separate polygons that represent different districts.  Overlaid on this layer is a polyline layer of 2800 streets, some which are completely bound by the polygons and some which crisscross two or more polygons.

My task is two find the total length of streets within each polygon (district).  I've already done this with one district using Vector> Geoprocessing Tools > Intersect and then performing summary statsitics.  However, in order to avoid doing this manually 127 times, I would like to write a script that will iterate through each polygon - performing the intersect with the polyline layer and leaving me with 127 separate shapefiles.

I not very familiar with Python - I only know some very basic commands.  I've attached images to help you visualize and conceptualize the problem.


texts: Why unable to edit ArcSDE from multiple computers on LAN using custom ArcGIS Engine application?

texts:
I cannot edit an ArcSDE database from multiple computers on a LAN network using a custom ArcGIS Engine application, I have tried both versioned and non-versioned feature classes, but they cannot be edited from more than one computer.

Inside the custom ArcGIS Engine application, for every single operation (Create, Remove), I open and close an edit session.

Any advice or experiences are welcome as why I may be unable to edit an ArcSDE database from multiple computers on a LAN network using a custom ArcGIS Engine application?

I tried to edit the same feature classe from 2 computers in the same time , but when i closed the edit session in the second computer, it gimme this error Your edits could not be saved because the target version has been updated since you started editing. These edits have been merged into your edit session. Review the changes and save again.

i'm editing feature classes like in this example :

    IWorkspaceEdit workspaceEdit = (IWorkspaceEdit)DB_Access.get_obj().ws;
    workspaceEdit.StartEditing(true);
    workspaceEdit.StartEditOperation();

        IFeatureClass fc = DB_Access.get_obj().ws.OpenFeatureClass("points");

        IFeature f = fc.CreateFeature();

        f.set_Value(fc.Fields.FindField("name"), "name_of_point");

        f.Shape = geometry;
        f.Store();

        workspaceEdit.StopEditOperation();

        workspaceEdit.StopEditing(true);


texts: Issues merging divided carriageways using Skeletron

texts:
I would like to merge divided roads from openstreetmap, and a tool that keeps getting mentioned is Skeletron, for example here.

I tried the tool on a small geographic area, and the results are fairly poor. I used zoom 18 as input, and received this:

Based on the documentation, and the linked paper, I was surprised to get a result like that. I feel like I must be doing something wrong. Is there a way to get a more high quality result?


texts: Edit Geometry with digitizing tools in QGIS of a postgis layer

texts: I have a multiline PostGIS layer added to my QGIS 2.10.1 Pisa. I'm able to add features using the digitizing tools, but I'm not able to use the delete or node tool. It remains gray even if I'm in editing mode. Every thing works fine on a temporary layer. Does anyone has a clue why I can't edit?


texts:
I have 1 shp line, I want to determine the intersection positions of the roads at 4 the intersection.

I used "intersection point at the road" but it gave an unexpected result (picture 1)

I want to results as picture 2

Is there any method to make it?


texts: Getting parameters of Python Toolbox tool to appear in ArcGIS Pro?

texts:
I can't get the parameters of my tool to show up in ArcGIS Pro 2.1 Because of my work I can't
put my exact code here, but I'll put the skeleton of what I'm doing and hopefully someone
can see what I'm doing wrong.

As a disclaimer, I know I can right click and add the parameters directly into my script tool,
through pro directly, but I want to do it strictly through python so that I can add my own
error messages. However, my python skills are not the strongest and I don't see where I am
going wrong.

import arcpy

class Toolbox(object):
    def __init__(self):
        self.label = "My label"
        self.alias = "The alias"
        self.description = "The description of the tool"
        self.tools = [MyTools]

# Functions

def firstFunction(table):
    """function stuff here"""
def secondFunction(param01):
    """function stuff here"""

# Parameter Info

class MyTools(object):

    def __init__(self):
        self.label = "test"
        self.description = "The description of the tool"
        self.canRunInBackground = True

    def getParameterInfo(self):
        table = arcpy.Parameter(
            displayName = "TABLE",
            name = "table",
            datatype = ["Table", "Feature Layer"],
            parameterType = "Required",
            direction = "input")
        param01 = arcpy.Parameter(
            displayName = "param01",
            name = "param01",
            datatype = "Long",
            parameterType = "Required",
            direction = "input")
        parameters = [table, param01]
        return parameters

    def execute(self, parameters):
        table = parameters[0].valueAsText
        param01 = parameters[1].valueAsText

        MyTools()
        firstFunction(table)
        secondFunction(param01)


texts:
I need to filter an OpenLayers layer using a CQL filter to an specific Bounding Box. Is it possible?

Something like this isn't working:

igh > 5000 AND bbox(igh,20786691,22051867,-305577260,-306370356)


texts: ogr2ogr large gdb file to SQLite warning: "organizePolygons() received a polygon with more than 100 parts"

texts:
I'm trying to write a .gdb file (could also be any other format, geojson, gpkg, etc.) to a Spatialite DB but it takes awhile and I am getting a warning.
Command:
ogr2ogr --config OGR_SQLITE_CACHE 1024 -f SQLite -dsco SPATIALITE=YES -nlt MULTIPOLYGON slite.db tmp_dls/data_d4e815db1/Universe.gdb -nln paqrcel_lines lotlines -gt unlimited

Warning message:
Warning 1: organizePolygons() received a polygon with more than 100 parts. The processing may be really slow.  You can skip the processing by setting METHOD=SKIP, or only make it analyze counter-clock wise parts by setting METHOD=ONLY_CCW if you can assume that the outline of holes is counter-clock wise defined

I know SQLite has low default settings for transactions and memory which is why I added the -gt unlimited and --config OGR_SQLITE_CACHE 1024 options but they dont have any affect on the speed of the command. With those 2 options removed it takes the same amount of time.
Is there a command in ogr2ogr I can use to make this run faster, or maybe I should switch to geopandas and use some shapely function to solve this?
Thanks any help is appreciated!


texts:
I have QGIS 2.18.3 installed. I'm trying to change angle units from degrees to Gon/gradians (90 degrees = 100 Gon).

I navigated to Settings--> Options --> Map Tools --> Preferred angle units, and changed it from "degrees" to "Gon/gradians".

It seems that Gons/Gradians are applied to measurement tools, but not for the Advanced Digitizing tools. Is there any way to get Gons/Gradians also to work in Advanced Digitizing tools as input angle units?

Bellow are the results (after hitting "OK" and restarting QGIS)


texts: Is the resolution of MYD13A2/MOD13A2 ( MODIS 16 day 1km Global Vegetation Indices) erroneous in Google Earth Engine?

texts:
The MODIS vegetation indices collection shows a resolution of 111319.4 metres in Google Earth Engine when, in fact, it should be 1000 metres.

Below is my code:

var Veg = ee.ImageCollection("MODIS/006/MOD13A2").filterDate ('2009-04-15',
'2009-04-30').max().select('NDVI');
print('Scale in meters:', Veg.projection().nominalScale());

I checked the resolution because I wanted to change resolution of my 90m Digital Elevation Model to 1000m using MYD13A2 data, but it gave the following error:

Layer 1: Tile error: Too many input pixels per output pixel. Need 1440001, but only 10000 allowed.


texts: Using SQL parameter when saving to GeoPackage

texts:
In QGIS 3.0, when I am saving a layer as a GeoPackage when using the Points to Path tool.

The 'Paths' parameter is pre-filled with something like this,

ogr:dbname='C:/Paths/Paths.gpkg' table="1119" (geom) sql=

Where do I find information on the parameters allowed in this field?  It would be great to use the SQL option to only choose the geometry I want here.  I haven't found a solution in the QGIS docs (the newest ones).  Where else could I look?


texts: How to fillup the holes at the time of object detection in GEE?

texts:
I permormed focal_max filter on my binary image on GEE but it not give me an output like python operator scipy.ndimage.morphology.binary_fill_holes.
var training = filt.sample({
  region: area,
  scale: 30,
  numPixels: 10000
});
var clusterer = ee.Clusterer.wekaKMeans(2).train(training);
var result = filt.cluster(clusterer);
print(result)
Map.addLayer(result.eq(1), {}, 'clusters');
var res = result.eq(1)
var kernel = ee.Kernel.square({radius: 3});
var opened = res.focal_max({kernel: kernel, iterations: 1})
Map.addLayer(opened)

How to perform this in GEE?


texts:
I have a line feature class that I would like to add to an AGOL map. To get the appropriate symbology, I need to have certain portions of the line to have two colors running parallel to each other.
In ArcMap, this is easy since you would just offset one of the lines. However there is no offset feature in AGOL. Also when I create the symbology in ArcMap and bring it into AGOL, I am unable to import symbology.
Is there something I should be doing? Or am I just not able to do this AGOL?


texts: Using Model Builder to find the overlap of any 2 routes

texts:
As you can see from the following picture, I want to use ModelBuilder to find the overlap segment between any two routes, and then I need to merge the overlap segment with the original route.shp with a feedback loop. In the end, I will have a shp file, including all the routes and all the overlap segment between any two routes. The model shown in the following picture works well without iteration since I just manually assign the Route 1 and Route 2 as the input for the INTERSECT.

Suppose I have 100 routes. I need to do the INTERSECT of Route 1 and Route 2, Route 1 and Route 3, ..., Route 1 and Route 100, and Route 2 and Route 3, Route 2 and Route 4, ..., Route 2 and Route 100, ..., Route 98 and Route 99, Route 98 and Route 99, and Route 99 and Route 100.

ModelBuilder doesn't allow us to do the nested iteration so I can't use two for loop in one model. I wonder if there is a way to do it in separate models to realize the thing I want.

I am not familiar with Python code.


texts: Creating random points along polyline in QGIS?

texts:
I am trying to create random points along a polyline shapefile in QGIS. Basically, like the 'random points' tool in 'research tools' but for polylines instead of polygons.

I have tried converting the line file to a polygon shapefile, but it fills in some areas with polygons, whereas other areas remain long line-type polygons.

I am pretty new to QGIS, and am not really familiar with Python code.


texts: Export classified land cover types as table - Google Earth Engine

texts:
My question relates to extracting the area of different land cover types in a region. So far, I classified a Landsat7 image using MODIS data as classifier. When I choose a point on the Map, the inspector shows me which landcover type is at this point. However, I do not know how to calculate the area of this landcover type for the region and export it.

- So, how do I calculate and export the area of all landcover types in my classified image?

- Also, how do I export the "upsampled" image properly? I only get black images with my current export function

    // Use the MCD12 land-cover as training data.
    var modisLandcover = ee.Image('MODIS/051/MCD12Q1/2011_01_01')
      .select('Land_Cover_Type_1');

    // A pallete to use for visualizing landcover images.
    var landcoverPalette = [
        'aec3d4', // water
        '152106', '225129', '369b47', '30eb5b', '387242', // forest
        '6a2325', 'c3aa69', 'b76031', 'd9903d', '91af40',  // shrub, grass, savanah
        '111149', // wetlands
        '8dc33b', // croplands
        'cc0013', // urban
        '6ca80d', // crop mosaic
        'd7cdcc', // snow and ice
        'f7e084', // barren
        '6f6f6f'  // tundra
    ];

    // A set of visualization parameters using the landcover palette.
    var landcoverVisualization = {palette: landcoverPalette, min: 0, max: 17, format: 'png'};
    // Center over our region of interest.
    Map.centerObject(geometry, 9);
    // Draw the MODIS landcover image.
    Map.addLayer(modisLandcover, landcoverVisualization, 'MODIS landcover');

    // Load and filter Landsat data.
    var l7 = ee.ImageCollection('LANDSAT/LE07/C01/T1')
        .filterBounds(geometry)
        .filterDate('2011-01-01', '2012-01-01')
        .sort('CLOUD_COVER', false);            //sort by cloud cover ('false' indivcates descending )

    // Draw the Landsat composite, visualizing true color bands.
    var landsatComposite = ee.Algorithms.Landsat.simpleComposite({
      collection: l7,
      asFloat: true
    });
    Map.addLayer(landsatComposite, {min: 0, max: 0.3, bands: ['B3','B2','B1']}, 'Landsat composite');

    // Make a training dataset by sampling the stacked images.
    var training = modisLandcover.addBands(landsatComposite).sample({
      region: geometry,
      scale: 30,
      numPixels: 5000,
      seed: 0
    });

    // Train a classifier using the training data. USING CART() classifier
    //var classifier = ee.Classifier.cart().train({
    var classifier = ee.Classifier.randomForest(100).train({
      features: training,
      classProperty: 'Land_Cover_Type_1',
    });

    // Apply the classifier to the original composite.
    var upsampled = landsatComposite.classify(classifier);

    // Draw the upsampled landcover image.
    Map.addLayer(upsampled, landcoverVisualization, 'Classified Image');

    // Show the training area.
    Map.addLayer(ee.Image().paint(geometry, 1, 2), null, 'Training region');

    var jsonCoordString = geometry.toGeoJSON();

    Export.image.toDrive({
      image: upsampled,
      description: '2011_classification',
      folder: "Syria",
      scale: 30,
      fileFormat: 'GeoTIFF',
      formatOptions: {
        cloudOptimized: true
      },
      region: jsonCoordString
      });


texts: Using QGIS for current SAR (Search And Rescue) operation?

texts:
A friend of mine, Christiaan Wilson is missing since 15 may during descending the Dhaulagiri mountain in the Himalayas. One searching team has already finished one operation. A second is on the way.

More information about the search: https://www.christiaanwilson.nl/english/

But now we are looking for other options instead of helicopters and ground teams. From my Hydrology background I know that there are a lot of possibilities with satellites. For that I have already contact with a teacher of mine, but with more people involved the process can be faster. Do you know someone who can help me further with this? Or can you help me further? It would really help me and the people who are working day and night to find Christiaan.

It would be really nice, because every minute is vital.


texts: Mixing geometry types in one PostGIS table

texts:
I'm faced with the following problem. I have to migrate from Oracle database to PostgreSQL+PostGIS. Currently all geometries of all types are stored in one table and each record contains a "lid" field which indicates features of the same layer.

What are the pros and cons of using such a method? Should I break up data into multiple tables if I don't need to use the database with third-party software? What about performance of spatial queries, will the indexes help me?


texts: Leaflet GeoJSON features not showing up in map

texts:
New to GeoServer and Leaflet...loving it so far, but I am weak in JavaScript. I have just built a simple web page that should load points from a GeoServer WFS service using GeoJson. However, even though the page loads without errors and I can see the GeoJSON returned from GeoServer using Firebug, the points are not showing up on the map. Can anyone look at my code and see what I am missing?

  <link rel="stylesheet" href="http://cdn.leafletjs.com/leaflet-0.5/leaflet.css" />
  <title>title name here</title>


texts: QGIS, CRSs and metric measurements of distance matrix and nearest neighbour?

texts:
I have thousands of points on a multiple layers and I need to get the distances between each point and every other point within a layer (so i can get min, max, mean and median distances) and the nearest neighbour values. I have used the distance matrix to achieve distances between each point, and the nearest neighbour tool to get the nearest neighbour value.

It all seemed to be working fine until I realised that the distance values are wrong. I assumed they were in km as the measuring tool displayed the distances of known distances correctly. However, when running it through the distance matrix and nearest neighbour, the values are massively different. In a similar vein to this question: Getting result in metric from QGIS Distance Matrix tool in QGIS?

Example:
Point 1 to Point 2 is 1.12km in reality
Manual measuring tool on QGIS says 1.12km
Distance matrix says 0.010682
Nearest neighbour mean observed distance says 0.010682

I have tried changing CRS and reprojecting and saving as a new shapefile, but it still doesn't work. I can get the measuring tool to display the correct measurements but as soon as i run it through the distance matrix or the nearest neighbour i get the same wrong values no matter what the CRS is.

I thought it might be in degrees, but how do I then convert it to km's so when I am writing up my results I can say that 'the mean value was Xkm' rather than some obscure number which doesn't readily mean anything to anyone?

Can anyone offer any advice?

Is there a projection I should be using (my points are plotted from long/lat)?

Is there a setting I need to change for the plugins?

Is there another tool which will do this for larger numbers of points? I

s there a formula I can use to change all the 'wrong' values (possibly degrees) into m/km's?

I don't know python, but is there a clever bit of code I can use to do this?

This worked for a while, but now it's reverted again. I managed to batch execute the reprojection, and then batch distance matrix and nearest neighbour and it seemed to work. For all except one layer. I've tried reprojecting just one file using save as manually, but it won't work. I've installed it on a second computer and tried it won't work. I went back to the original data, imported it into QGIS from csv, straight into ESPG:27700 with the project set to ESPG:27700 and it's still measuring in degrees. It's not the coord's, i've checked and know what the measurements should be. I then ran it with the second year of data ready to disregard one bad file, and it is refusing to work. I really don't understand the problem.

It says it's reprojected the file, but it isn't, or it is but something is making QGIS continue to measure in degrees. I didn't changed any settings from the set which worked, to the set which didn't.

Is this a glitch in the program?

Is there another free GIS program that can do distance matrix and nearest neighbour which might save me this continuous circling?

I have just taken an original unsplit file and reprojected the entire year as one file. I saved and then split it by date, before using distance matrix. It is still running and outputting files and this is likely to take some time but the first couple of outputs so far appear to be in metres. I believe the problems I am experiencing are related to a bug. I will try to report it as such, but for anyone experiencing similar problems all I can suggest is that you disregard any files which aren't working and start right from the very beginning with the data, reprojecting before any other analysis as it seems something lodges itself in the file and stops it working properly when you try to run distance related analysis at some point.


texts: Finding the average value of a gradient in a certain region in QGIS

texts:
In one QGIS project, I have a few TIFFs. One is a black-to-white gradient map showing average precipitation. Each pixel has a value that signifies how much yearly precipitation it gets. I have another TIFF in the same project with many climate zones, each color signifying a different zone.
Would it be possible to find the average precipitation in each zone, like the average value on the gradient of each region? I'm not very good with QGIS.
Also I don't know how to code at all. I don't know if this is possible to do without having to code, but if it's necessary, would it be possible to make the instructions for simple or at least understandable for a beginner?


texts: Multiple attribute layers with same geometry in geopackage with ogr2ogr

texts:
I have multiple PostGIS tables I want to export to a single gpkg file. The tables have identical schema, and identical geometries -- the numerical attributes differ. How do I create a gpkg file, using ogr2ogr with these multiple tables, without duplicating the identical geometries? I need to use ogr2ogr because this needs to be scripted.

My question is related to this existing question Saving duplicate geometry layers into a geopackage? but I don't see a solution using ogr2ogr.

I am using the -update flag per the accepted answer here: Saving duplicate geometry layers into a geopackage?

Is there a way to tell ogr2ogr / gpkg to not duplicate the geometries? Or are the only solutions to either have all attributes in a single feature table (many many columns), or have 1 feature table and many attribute tables, and then perform joins in an application when opening the gpkg?


texts: How to create new shapes based on partial existing shapes and georeferenced scan?

texts:
I am trying to generate a new shapefile in QGIS based on a georeferenced scanned map of colonial powers in Africa and the NaturalEarth country shapes. The overlay of the georeferenced TIFF looks quite ok, and most of the borders still exist today.

How can I combine
a) Existing parts of a shape into a new, larger shape?
b) Existing parts of a shape with some new borders I draw to make a new shapefile?


texts: Why doesn't the gml 3.1.1 schema validate with xmllint?

texts:
Running xmllint on the gml 3.1.1 schema fails, even though the release notes claim that this release was created to get the schema to validate using 'common tools'.  There's questions about this on the web, but no satisfying answers.

xmllint --schema http://schemas.opengis.net/gml/3.1.1/base/gml.xsd \
        --noout foo.xml

returns something like:

  element complexType: Schemas parser error :
  local list type: A
  type, derived by list or union, must have
  the simple ur-type definition as base type

This is not the same problem as the one described here, which was a Xerces C++ issue.

If there is some other xmllint issue here, what xml validator is known to work with GML?

I have this problem on a few different flavors of Ubuntu, CentOS, and RHEL 5.x

Update:

The 3.1.1 and 2.1.2 GML schemas validate using the following tools:

The W3C XML Schema Validator
Sun's msv:

running java -jar msv.jar -warning gml\base\gml.xsd

returned 1 warning about "anyType" in gmlBase.
Using the validator bundled in eclipse worked with no errors or warnings.


texts: How to update attributes from list of values stored in a Excel or CSV table?

texts:
I need a Python program regarding update attributes. I have 20 shapefiles and in each shapefile, attributes to be updated based another field

My requirement is to update the codes in CODE field as per the name (Example list is given right side) for 50000 features, like if name is AAA, then code will be 419, if name is ABC, then code is 420.

I have 850 different codes and I want to update all 850 at once in shapefile. I can update two or three values by using below code but I don't know how to update 850 codes.

import arcpy
from arcpy import env

env.workspace = r"C:\\data\test.shp"
env.workspace = fc
fclist = arcpy.ListFeatureClasses()

for fc in fclist :
   fields = ["name", "code"]

   with arcpy.da.UpdateCursor(fc, fields) as cursor:
      for row in cursor:
          if row[0] == "AAA"
              row[1] =  419
          elif:
          if .row[1] == "BBB"
              row[6] = 420
              cursor.updateRow(row)

print "codes updated"

If I want to update all 850 codes, need to write if...else statement 850 times  to update all codes. I have saved all names and codes in Excel. So instead of writing many if else statements, please help how to update all codes based on name by using Excel or CSV file.


texts:
I have MODIS NDVI TIFF files for a year
I want to first clip it with a irregular polygon (not by its extent only) of my study area and save it as a multi band raster in TIFF format
The TIFF image are already projected to the needed coordinate system and datum (WGS 84, UTM 45N)
I used the writeRaster() function for the rasterbrick object but it gives me an error.

If anyone has faced similar problems or could you suggest me better way doing this task.

    writeRaster(ndvit1, filename="ndvi_2015t1.tif", bandorder='BSQ', format="GTiff", overwrite=TRUE)
    Error in rgdal::putRasterData(x@file@transient, vv, band = i, offset = off) :
    Failure during raster IO
    In addition: Warning message:
    closing unused connection 4 (C:\Users\Ghimirebr\AppData\Local\Temp\RtmpwJschv\raster\r_tmp_2016-04-11_121550_3352_50703.gri)

## Clip a raster (brick) with the polygon (shape file) in its extent
## and mask value nan for the remaining area.

library(raster)
library(rgdal)

# Setting working directory:

setwd("C:/Working_Folder/R")
ndvi_path <- "ndvi_2015"

# Reading the shapefile
myShapeInR <- readOGR("C:/working_folder/Boundary",
                      "Chitwan_Nawal_Boundary")

# To view the summary of vector file
summary(myShapeInR)

# Getting the spatial extent of the shapefile
e <- extent(myShapeInR)

# Reading the raster you want to crop
all_ndvi <- list.files(ndvi_path,
                       full.names = TRUE,
                       pattern = ".tif$")
ndvi_stack <- stack(all_ndvi)
ndvi_brick <- brick(ndvi_stack)

# Converting from -1 to 1 (NDVI)
ndvit1 <- ndvi_brick*0.0001

# Saving the RASTER BRICK
writeRaster(ndvit1, filename="ndvi_2015t1.tif", bandorder='BSQ',
            format="GTiff", overwrite=TRUE)
# myraster <- raster(myraster.filename)

# Cropping the raster to the shapefile spatial extent
ndvit1.crop <- crop(raster, e, snap="out")

# Dummy raster with a spatial extension equal to the cropped raster,
# but full of NA values
crop <- setValues(ndvit1.crop, NA)

# Rasterize the catchment boundaries, with NA outside the catchment boundaries
myShapeInR.r <- rasterize(myShapeInR, crop)

# Putting NA values in all the raster cells outside the shapefile boundaries
ndvit1.masked <- mask(x=ndvit1.crop, mask=myShapeInR.r)


texts: How to get a group layer delete signal in QGIS

texts:
I'm developing a plugin into QGIS 2.18. I need to catch if the user deletes a group layer. I tried the signal layerWillBeRemoved from QgsMapLayerRegistry class but it sends only layers not group layer.

What should I do?


texts: Compiling Postgis 2.2.0 on Raspian error doesn't install postgis extension

texts:
Installed postgresql 9.4 following this tutorial using the method of already compiled pkg. All seems fine with postgres installation.

testdb=# SELECT name, setting FROM pg_settings where category='File Locations';
       name        |                 setting
-------------------+------------------------------------------
 config_file       | /etc/postgresql/9.4/main/postgresql.conf
 data_directory    | /var/lib/postgresql/9.4/main
 external_pid_file | /var/run/postgresql/9.4-main.pid
 hba_file          | /etc/postgresql/9.4/main/pg_hba.conf
 ident_file        | /etc/postgresql/9.4/main/pg_ident.conf
(5 rows)

Then tried to install postgis 2.2.0 folowing the steps at this trac page (im fairly new with *nix env's and this build from source stuff stills confusing me a lot so go gentle with me) .

So downloaded postgis-2.2.0.tar.gz in a /temp_postgis and untarred .
List item
Folowed steps for install from source in this trac for debian whezzypage
During the install it requested many librarys like geos , proj4 and others that i installed.

Now the problem when in try to create postgis extension it gives me an error and i can see the extension at the pgadmin as option to install extension on databases.

receiving the following error:

testdb=# CREATE EXTENSION postgis;
ERROR:  could not open extension control file "/usr/share/postgresql/9.4/extension/postgis.control": No such file or directory
testdb=#

Any clues what could have go wrong ? Regards.

EDIT 1:
After make check -

PATH is /usr/lib/postgresql/9.4/bin:/usr/local/sbin:/usr/local/bin:/usr/sbin:/usr/bin:/sbin:/bin
Checking for shp2pgsql ... found
Checking for pgsql2shp ... found
TMPDIR is /tmp/pgis_reg
psql: FATAL:  role "root" does not exist
Creating database 'postgis_reg'
createdb: could not connect to database template1: FATAL:  role "root" does not exist
psql: FATAL:  role "root" does not exist
Loading PostGIS into 'postgis_reg'
 failed (Error encountered loading /temp_postgis/postgis-2.2.0/regress/00-regress-install/share/contrib/postgis/postgis.sql: /tmp/pgis_reg/regress_log)
Makefile:265: recipe for target 'check' failed
make[1]: *** [check] Error 1
make[1]: Leaving directory '/temp_postgis/postgis-2.2.0/regress'
GNUmakefile:16: recipe for target 'check' failed
make: *** [check] Error 1

Now on the other side- i was able to run the postgis.sql and spatial_ref_sys.sql on pgadmin sql window, all ended successfully , i can now see postgis extension as option to install postgis via pgadmin and pslql - recreated a few other db's and installed postgis and seems to be fine . Now the question - if make check tests failed should i worry ? i belleive so :-)


texts: Looking for sample code or documentation for gpxpy waypoint extension data?

texts:
I am translating some of my Ruby scripts that manipulate GPX files into Python so I can use them from within QGIS but I can't find any detailed documentation.  The readme on github gives a few very basic examples and leaves it at that.  In particular I am trying to figure out how to access extension data for waypoints.

I wondered if anyone had made their own cheat sheet that they would care to share?


texts: QGIS Raster Calculator via python: no substraction and division possible

texts:
Subtractions and divisions are not possible when using the QGIS Raster Calculator from a python script within the QGIS python console: Error code:4 (Parsing Error?). Additions and Multiplications work well. PYthon 2.7, QGIS 2.14 & 16.

import qgis.core, processing, glob,
from qgis.core import QgsMapLayerRegistry, QgsRasterLayer
from qgis.analysis import QgsRasterCalculator, QgsRasterCalculatorEntry

InputPath = "D:/Heckmann.M/NigerABN/1_Projekte/14_Landsat/trialcbs/"
band2Path = glob.glob(InputPath+"*_sr_band2.tif")[0]
band5Path = glob.glob(InputPath+"*_sr_band5.tif")[0]
rband2 = QgsRasterLayer(band2Path, "name2")
rband5 = QgsRasterLayer(band5Path, "name5")
iface.addRasterLayer(band2Path, "Band2")
iface.addRasterLayer(band5Path, "Band5")
outputPathQ = "D:/Heckmann.M/NigerABN/1_Projekte/14_Landsat/trialcbs/QGIS_35_sub7.tif"

def calc_MNDWI(green, MIR, rasterFile, MNDWI_file):
    entries = []
    # Define band GREEN
    green = QgsRasterCalculatorEntry()
    green.ref = 'green@1'
    green.raster = rband2
    green.bandNumber = 1
    entries.append(green)

    MIR = QgsRasterCalculatorEntry()
    MIR.ref = 'MIR@1'
    MIR.raster = rband5
    MIR.bandNumber = 1
    entries.append(MIR)

    calc = QgsRasterCalculator('(MIR@1-green@1)', MNDWI_file, 'GTiff',  rasterFile.extent(), rasterFile.width(), rasterFile.height(), entries)
    print "Calculation Result", calc.processCalculation()
    iface.addRasterLayer(MNDWI_file,'MNDWI QGIS')

calc_MNDWI(band2Path,band5Path, rband2, outputPathQ)


texts:
How can i add OSM layer in QGIS?

The help says there is an OSM plugin i can use. But i don't see it the list of plugins.

Or maybe there is a WMS URL for OSM i can use?

  QGIS version 1.7.0-Wroclaw on Windows 7 64bit


texts: QGIS Small Square Dots Appearing on Postcode Sector Polygons

texts:
I dissolved UK postcode polygons at the sector level to create the postcode sector polygons I need. However, I seem to have more postcode polygons than I expected, about 800 additional postcodes that aren't in the national statistics data.

I added a binary indicator to all of the postcode sectors I got from the national statistics office in the UK (ONS), and visualised them. I got the following:

And also at a higher zoom, this:

I'm not able to select these small square dots and I don't seem able to find them in the attributes data. For example, when I search for polygons based on area and rank them, the sectors with the smallest areas are too big to be these dots, for example:

Pcd.Sector Area
M60 8   73
VGU02   129
PH2 1   149
VUB01   616
DL981   807
VSE04   847
M60 7   1101
OL133   1263
W1A 3   1264
EC4N1   1450
EC3A2   1537

I'm really not sure what these are, has anyone seen these before?


texts: Relying on order of layer categories in PyQGIS

texts:
Follow-up to this question, after input considering this other question:

To figure out which feature in a QGis vector layer, styled with a categorized renderer, is currently visible (i.e. enabled), I need to

a) Get the value of each category, e.g. with

print [cat.value() for cat in layer.rendererV2().categories()]

--> [u'', u'1', u'9']

b) Get the status of the checkbox, e.g. with

ltl = QgsProject.instance().layerTreeRoot().findLayer( layer.id() )
ltm = iface.layerTreeView().model()
legendNodes = ltm.layerLegendNodes( ltl )
print [node.data(Qt.CheckStateRole) for node in legendNodes]

--> [2, 0, 2]

So far, so good. Can I rely on the order of those two being the same or do I have to do some sort of cross-checking (and if so, any ideas how? The symbol instances that I can access through both paths don't seem to match...)


texts: labeling in qgis 2.0: if one letter in string is "X"

texts:
I have a shapefile and I want to label all shapes with a "B" or "A" in its name (string). All others should be ignored.

Like:
If [field] contains "B" draw label

Im using qgis 2.0.


texts: Showing arrowheads in line end using QGIS?

texts:
I would like to show arrowheads in the endings of every drawn linestring.
Is it possible in QGIS?

I have been trying to use a marker line but repeats all de arroheads all over the line, it is not restricted to the ends:


texts: How to find coordinates of a triangle with azimuth?

texts:
I am just starting out in GIS and it's been a long time since I took math in university. I am teaching myself through a textbook but there are a few questions I am stumped on.

Essentially I have a triangle. Say one point is A, the other is B and the third is C. I have the coordinates for point A and B. I have the azimuths from A to C and from B to C. I am trying to find out what the coordinates of point C.

Does anyone have any suggestions?

Note: VERY basic math, not a coordinate system or curvature of the Earth. Just the basic concept.