texts: Can QGIS analyse a personally created DEM to find transport links or place infrastructure?

texts:
I'm working on game development and I'm trying to use QGIS to help me show some features on a created landscape.

I have used a diamond square algorithm in excel to create randomly generated 'worlds', taken the data for one of these worlds and converted to a csv file to import into QGIS as a vector layer. From this I have created a DEM raster layer, calculated channel networks and assigned land use.

My next hope is to add transport links and infrastructure onto this map, and I wondered if this was possible to do in QGIS or whether the problem was too complex?


texts: Setting border-style to 'no pen' using QgsGraduatedSymbolRendererV2

texts:
I am developing a QGIS plugin that, after a series of inputs, will produce a thematic map using the QgsGraduatedSymbolRendererV2. Because the polygon are very small, the boarders make it such that when zoomed out the entire map almost looks completely black. I want to remove the boarders to do this but no documentation or form I have found have demonstrated how to do this. Below is what I have to create the thematic map, but somewhere in there I need to set the border style to 'no pen'. But how!?

  uri = QgsDataSourceURI()
  # set host name, port, database name, username and password
  uri.setConnection("localhost", "5432", "public", "postgres", "password")

  uri.setDataSource("susip", str(filename)+"_mce_results", "the_geom")
  layer = QgsVectorLayer(uri.uri(), str(filename)+"_mce_results", "postgres")

  fieldName = "score"
  fieldIndex = layer.fieldNameIndex( columnName )
  provider = layer.dataProvider()
  numberOfClasses = 10
  color1 = QColor("red")
  color2 = QColor("green")
  color3 = QColor("white")
  ramp = QgsVectorGradientColorRampV2(color1, color2)

  renderer = QgsGraduatedSymbolRendererV2.createRenderer(layer, columnName, numberOfClasses, QgsGraduatedSymbolRendererV2.Quantile, QgsSymbolV2.defaultSymbol(layer.geometryType()), ramp)

  layer.setRendererV2( renderer )
  QgsMapLayerRegistry.instance().addMapLayer( layer )


texts: Setting up TileCache to make use of Mapnik and OpenStreetMaps

texts:
First of all, I'd say that I'm completely newb at maps related stuff... I am trying to install TileCache and configure it to use Mapnik to render tiles from OpenStreetMaps. I managed to setup the Tilecache basic CGI conf, but now I do not know how can I continue with setting up Mapnik and my own tiles because I am not able to find any tutorial on the net explaining the steps for a complete TileCache setup for OSM and Mapnik (neither in TileCache's web it is explained), and, since it is the first time I deal with WMS and maps rendering and so on, I find difficult to understand the proccess... Could somebody help me, please?

BTW, I need this to run on Windows Vista/7

Thanks in advance!


texts:
I would like to be able to change the order of the layers in leaflet.

I tried implementing two leaflet plugins in order to change the layer order.

First I tried leaflet-control-orderlayers, which seems to work, but does not support layerGroups, here is my filed bug.

Next I tried mapbbcode, which does support layerGroups, but does not seem to work or be maintained.

Does anyone have any suggestions for how to implement this feature?


texts: Updating text elements not working when using APRX path vs "CURRENT"

texts:
I am beginner in ArcGIS Pro and I am trying to expand into using python with ArcGIS Pro.
I have a working script that will update a specified text element in all layouts in ArcGIS Pro, in my case it's a project name.

project_name = 'Test Project'
aprx = arcpy.mp.ArcGISProject("CURRENT")
for lyt in aprx.listLayouts():
    for elm in lyt.listElements("TEXT_ELEMENT", "Project Name"):
        elm.text = project_name
aprx.save()

My problem comes when I try to run the script outside of a running ArcGIS Pro instance and use the APRX path rather than "Current" in arcpy.mp.ArcGISProject()

import arcpy

project_number = "12345"
project_name = "Test Project"
project_number_group = project_number[:2]
newpath = "N:\\Projects\\" + project_number_group + "000\\" + project_number + "_" + project_name.replace(" ","_")
aprx = arcpy.mp.ArcGISProject(newpath + "\\APRX\\Project_Figures.aprx")
for lyt in aprx.listLayouts():
    for elm in lyt.listElements("TEXT_ELEMENT", "Project Name"):
        elm.text = project_name
aprx.save()

The odd thing is when I go line by line things appear to be normal, the aprx variable calls the correct project, the path appears correctly, it's just nothing actually updates when I run it. No errors spit out either, so I am at a loss currently.


texts: Is it possible to get the positions of all the stoplights for a given area?

texts:
I would like to know if there is a database of longitude and latitude values that locate all stoplights for a given area in the US like a city or a state? How would one go about obtaining this info?

Looking to gain the locations of stoplights (traffic signals) to do some traffic research.


texts: Best way to relate PostgreSQL table to layer

texts:
There is an online PostgreSQL database (non-spatial) with water monitoring data which is regularly being updated, in another side is a spatial layer of rivers (for example) which should be represented in ArcGIS Server.

PostgreSQL database is being updated by other users. There is also a web interface to access this data.

I need to relate (not join) layer to this database based on the code (ID) to display respective content (table records) in ArcGIS Server interface.


texts:
I am using OpenLayers to plot a point on a map; the point should be located in the UK; but keeps appearing in the middle of the ocean (0,0). I have check my projection and I think this is correct. What am I doing wrong?

Here is the code I use to define the map in OpenLayers3

var map = new ol.Map({
        target: 'map',
        layers: [raster, vector_layer,vectorLayer],
        view: new ol.View({
            center: ol.proj.transform([51.714090, -0.413285], 'EPSG:900913', 'EPSG:4326'),
            zoom: 2,
                projection: 'EPSG:900913',
            displayProjection: 'EPSG:4326',
        })
    });


texts:
I'm using OpenLayers 3 and i need to load a Geoserver group-layer

In OL2 i used a code like this and it worked:

var layer_group = new OpenLayers.Layer.WMS(
        "Sample Layer Group",
        "http://localhost:8081/geoserver/wms",
        {layers: "layer_group_name", transparent: true},
        {isBaseLayer: false, visibility: false}
    );
    map.addLayers([layer_group]);

I have no idea how OL3 loads a group.. simple layers he can load correctly


texts: Dissolving all overlapping, intersecting, touching polygons on single layer using QGIS?

texts:
I'm quite new to Qgis and have a problem which is easy to solve, I think. I've done research for hours but did not succeed.

I have a data set which contains thousands of polygons in one layer.
What I want to do:
Dissolve all polygons which overlap, intersect or touch with others. This means that there should not be any borders between overlapping, intersecting or touching polygons anymore. I don't have any helping attributes to achieve this.

Finally, I would like to export a new .shp file of the new data set.

I tried "dissolve all" but this makes just one big polygon, which I don't want. As a first step i thought to solve it with "Intersect...". But this works just with two layers. With the advanced digitizing tools I would have to select all of them manually to "merge selected features".

I'm using Qgis version 2.14.15


texts:
I have data that represents Airport floor plans stored in GeoJSON format.

Are there any standards, guidance, or good practices for Indoor Maps Styling?
Any references or resources would be useful.


texts: Distance between points along polyline ArcGIS for Desktop (pt. 2)?

texts:
I read over the solution in Distance between points along polyline ArcGIS for Desktop? and I have a slight variation of that problem:

I need to find the distance of each red point (transformers) from the yellow point (the main generator) along the distribution network.

Each transformer receives power from the main generator. I've included a screenshot of (a) large-scale view and (b) small-scale view (green box). Each line has a distance attribute in feet, which I've labeled. The lines do not inherently "flow" toward or away from the generator, directionality is arbitrary (see blue circled).

A line exists between each transformer, so we never need to find a "midpoint" distance per-se, transformers always intersect line endpoints/startpoints.

If there are multiple routes to reach the main generator, we need to find the the shortest route (I believe there is always only one conceivable route from transformer to generator in this data, but just in case...).

What is the best way to go about this?

Modify the "chainage" code mentioned above to process multiple line segments en route to a specific point an unknown distance away? Challenge is how would script know which path to take when the network forks en route to the generator?
Create a script to perform some sort of "Select by Location" ==> select network lines that intersect generator, then select network lines that intersect selected network lines recursively, causing the whole selection to propagate upon each iteration until the whole thing intersects the transformer of interest, then whittle the selection down to the route between the transformer of interest and generator (somehow), and sum the distances of lines along the route?
Read up on Linear Referencing and apply it here. I have certainly never used Linear Referencing in my career so far, but if M-enabling the data will solve my problems, I will gladly do due diligence.
Pretend our data is in Google Maps and use the Google Maps Distance Matrix API via Python wrapper. The Origin is the generator. The List of Destinations would be all the transformers. Get distance for each route. This option is obviously facetious, but could we use Network Analyst to do some of the logic above?


texts: Limitations when running the Union tool with large datasets?

texts:
I'm running ArcGIS 10.3 with an Advanced License on a Windows 7 64x system with 20 GB of RAM and a 4th Gen Intel i7 Quadcore-CPU.

I have a problem with the normally very efficient Union tool. What I'd like to achieve is a simple union of two feature classes in order to fill holes in my main feature class, say "flat areas" (blue). I therefore created a large surrounding polygon "union polygon" (green). The below images explain the situation (overview and zoomed in):

My blue feature class is large (60 MB), but not huge. It was created combining different features, then I ran the dissolve and buffer tools on it.

The problem now is that when I try to unify the two feature classes, ArcCatalog returns "Error 999999: Error executing Function. Invalid topology". I have read this answer about processing large data and followed as many of the instructions as possible:

No network drive
SSD
dissolved polgons

So finally my questions:

Does anyone have an idea what could cause the problem? Hardware issue?
Any further tipps or recommendations about the Union tool?
Or even better - an alternative in order to fill the holes in the "flat areas" feature class?

Update:

I've meanwhile found this post: Invalid Topology [Topoengine error] from ArcGIS for Desktop? which seems to be about the same problem as I'm experiencing, however there isn't a solution either. I'm currently running the Repair Geometry tool.

Also I've tried to run identity instead of union since that would also work for my task. The same error showed again.


texts: Rasterizing polygon vector grid cell. Define raster cell position?

texts:
I am doing many variables calculation over a vector 1x1 km² grid layer in QGIS. Each cell has many attributes, different spatial indexes.

I'd like to rasterize the layer by certain attributes but make raster grid cell overlap the vector 1x1 km² grid objects exactly.

So beyond defining cell size in rasterize algorithms I should define its position. Is there any way to do so?

Thanks.
I was already trying with use layer extent but raster layer is anyway displaced from vector one.

My input layer is EPSG:25831 but raster output is automatically projected to EPSG:3043. The displacement could be caused by that? Any way to fix it?


texts: By far the most of my data is taken up a column of long strings pointing to the source of the data for each point. I need to keep that. However, it would be neat if CartoDB could use this information (i.e. build the URL) without me providing the full URL, only the relevant ending (and save on the size of the table by not repeating the common part as a long string). Is this doable in their info window?


texts:
I have a tile server (MapCache) with a personalized map.

For client side's rendering I use OpenLayers.

I'd like my personalized cache to be accessible only on my web site,
i.e. I don't want another website to use the mapcache on my MapCache server.

Is there any solution ?


texts: How to define layer order in Openlayers?

texts:
I have a marker layer, two vector layers and two raster layers.
The marker layer gets covered by the two raster layers.

Is there a way to move the marker layer always on top?

Edited: The order of Z-indexes works well in the vector layer, but when I add another WMS raster layer on the map with a lower Z-index, the raster layer still covers the vector layer.

Edited: Could anyone point out the difference of the rendering order in terms of "LayerIndex" and "Z-Index" of a Layer? I think the layer index deal with the overlay order and the Zindex deals with the feature order in a vector layer, am I right?


texts: QGIS HTML labels only for objects visible on map

texts:
I'm printing QR-code that is generated on the fly in QGIS Print Composer. The URL to get the QR-code ist stored in a PostgreSQL database in the attribute-field "url".
<img src = "[% "url" %]"
    style="width: auto; height: 100%; display: block; margin: 0px 0px 0px auto;" />

In Print Composer this works fine as long as the map, I use is part of an Atlas. When the map switches to another Atlas-feature, the generated QR-code is updated.
My problem is, that I want to use a map without the Atlas functionality.
How can I get the URL only for the feature (supposing that only one feature is displayed) that is displayed on the map in Print Composer?
I was looking trough the expressions from QGIS, but did not find a suitable function to catch the data of a feature on a defined layer within the extent of a a map in Print Composer.


texts: When I create a stream segmentation using Archydro,“Fdr” and “Str” are the inputs for the Flow Direction Grid and the Stream Grid respectively. The output is the Link Grid. How is defined considering both of Fdr and Str simultaneously?


texts: QGIS Map Composer goes blank, outputs corrupt PDF

texts:
I am using 64bit 2.10.1-Pisa on Windows 10.

When I set up my composer to print, everything looks fine until I try to print or export. The composer goes completely blank and outputs a corrupted PDF using the Export to PDF function. The PDF file won't even open.

I had the same result exporting to PNG.

When I try to use the PDF document printer instead of exporting, QGIS stops responding indefinitely.


texts: (When) should I use a GIS-enabled database?

texts:
I am a long-time programmer, new to GIS. I am trying to get a feel for the database aspects of GIS and I understand that there are some databases which are specialized for GIS use. Basically I am trying to understand whether to use such a database or to stick with the hugely popular, well established, widely supported (and free) MySql.

The sort of applications I am liable to code would things like:

fleet management (land or ocean based vehicles)
employee tracking
inventory control (at a granularity of a physical position accurate to with a metre or so)
errm, that's about it, really

I would normally expect to be tracking no more than 1 few hundred (maximum, possibly a few thousand) items. Sometimes the items will be located in the same (large) building, series of building, city, country or world-wide, depending on the application.

I will occasionally be visually representing them either on a custom floor plan of a building or, more likely, in Google Earth or the like (more of which in another question).

It seems to me that I am as well off with MySql and adding columns for lat/long or other positional data, but that might just be because I know MySql.

Is there any reason why I should be looking at a more specialized database?


texts: Populate buildings polygons with address from TIGER address range

texts:
We have OSM shapefile of buildings polygons connected to OSM shpaefile of roads. And then from TIGER we have their well known address range table connected to roads again.

So via roads we can connect between buildings polygons (OSM) & address ranges (TIGER). {All this because in TIGER we don't have shapefile of buildings polygons}.

Now, Is it good idea to automatically assign each polygon with an address? For example If we have range of 17, and there are 17 polygons - With the direction data - Isn't it should be a trivial task? Is it true for any case of matched number of polygons and numbers in the range?


texts: Compiling resources in QGIS 2.99 (Windows) for creating plugin?

texts:
I tried out to create, for Windows, the Where Am I? QGIS 2.99 plugin (similar to Gary Sherman's book) by using Plugin Builder 3. In spite of being created, it was generated an error that it can be visualized in following image:

and this is an extract of produced text:

  TypeError: keyPressEvent(self, QKeyEvent): first argument of unbound
  method must have type 'QDialog'

  Traceback (most recent call last):
  File
  "C:/Users/Usuario/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\pluginbuilder3\plugin_builder_dialog.py",
  line 169, in keyPressEvent
      QDialog.keyPressEvent(event) TypeError: keyPressEvent(self, QKeyEvent): first argument of unbound method must have type 'QDialog'

However, I think that it could be not important because I found out that is was always produced by pressing Shift key.

As I don't have compilation resources in Windows, plugin was compilated (pyqt5-dev-tools) in my Debian GNU/Linux with make (without any error):

pyrcc5 -o resources.py resources.qrc

Afterward, it was newly exported to Windows at folder:

C:\Users\Usuario\AppData\Roaming\QGIS\QGIS3\profiles\default\python\plugins

When I try to activate plugin, I got that it is not compatible for this version of QGIS; as it can be observed at following image:

When I open folder of each installed plugin (reloader, builder, zoom to Belgium), at top they always have a __pycache__ folder; as it can be observed at following image:

However, created plugin doesn't:

I also corroborated that resources.py was adequately compilated because it has this header:

# -*- coding: utf-8 -*-

# Resource object code
#
# Created: Mon Nov 27 15:50:04 2017
#      by: The Resource Compiler for PyQt (Qt v5.3.2)
#
# WARNING! All changes made in this file will be lost!

from PyQt5 import QtCore
.
.
.

and similar structure in all plugins (however my compiler version for PyQt  is older; v5.3.2 vs v5.7.1).

My question is why was not adequately produced __pycache__ folder with Plugin Builder 3 or plugins created by user need to be placed in another folder in Windows environment?

Editing Note:

After creating Where Am I? plugin, I edited Qt5 objects in Qt5Designer from Linux and I included required lines to show location of a point click on the map; as it can be observed at following image. Plugin works as expected.


texts:
When I load a bunch of layers into QGIS from PostGIS, they seem to load in reverse alphabetical order.

Is there any tools or plugins that allow you to quickly sort the QGIS TOC in a different order?  For example I can think of a couple of options one might want:

Sort layers by geometry, where points are on top, then lines, polys, & rasters.
Sort Layers by name (alphabetical, reverse alphabetical).
Sort first by geometry, then by name.  This would be a combination of #1 & #2.  Points layers are sorted to the top (same as #1), but then points layers are sorted alphabetically (point-a, point-b, point-c, etc).

When loading a large number of layers, this would help to rapidly build up the right ordering for map production.


texts: OpenLayers 4.2 - Add a source.vector from String, not url

texts:
I'm using OpenLayers 4.2.
Currently I add a source vector to a layer vector like this :

var vectorSource = new ol.source.Vector({
    url: 'http://example/carto.geojson',
    format: new ol.format.GeoJSON()
});

var vectorLayer = new ol.layer.Vector({
        source: vectorSource
    });

I'd like to load my source from a string.

(url = "'{"type":"FeatureCollection",
"features":[{"type":"Feature",
  "id":0,
  "geometry":{"type":"Point",
    "coordinates":[-1.68117,
    48.3717]},
  "properties":{"titre":"Nutus",
  "datedebut":"1 avril 2018",
  "datefin":"12 avril 2018",
  "lieu":"ETANG AUX MOINES"}},
{"type":"Feature",
  "id":1,
  "geometry":{"type":"Point",
    "coordinates":[-4.50069,
    48.4083]},
  "properties":{"titre":"Nutus",
  "datedebut":"1 avril 2018",
  "datefin":"12 avril 2018",
  "lieu":"BREST "}}]}"';)

Is it possible and how?


texts:
I want to use this service for exporting images based on longitude/latitude coordinates. This API accepts esriMeter as parameter for bounding box, so I have to convert long/lat to esriMeter. How I can achieve that in Python?
For example, converting [49.958961724999995, 40.42190548000002, 49.95907198100002, 40.421990584000014] to esriMeter format.


texts: Is it easy to digitize in ArcGIS using the Wacom tablet?

texts: I do just enough digitizing (and touching up/ drawing when post-processing) to perhaps warrant the purchase of a Wacom tablet. Just wondering how well it plays with ArcGIS. What are your experiences? Is there a learning curve. Can the tablet and pen replace the mouse all together?


texts: Using pgRouting final disconnected nodes from given node

texts:
I am new to using the pgRouting library and haven been given a standard edge table
 id, source, target, cost, reverse_cost

I am comfortable with shortest path algorithms, etc, but this has me stumped
I am trying to find all the nodes where no graph exits, or is disconnected, from a given node and I am not sure what tools I can use to achieve this.


texts:
I have MODIS NDVI images, which needs to be masked with VI_Quality. I am using GRASS software and executing GRASS commands using shell scripts. I need to perform raster calculation (multiplication) for NDVI images with the VI_Quality. I want to multiply ith NDVI image with the ith image of VI_Quality. I have used the following codes :

for file in $NDVIFILES;
do
r.mapcalc "$file = $NDVIFILES * $QAFILES;
done

i=1
for file in $NDVIFILES;
do
r.mapcalc "file $((i++)) = NDVIFILES $((i++)) * QAFILES $((i++))";
done

I am unable to do this raster multiplication in a loop. Please let me know the correct syntax to do this operation.

Thanks


texts: How do I crop a TIFF down to a specific region?


texts: Error in rasterizing shapefile by using pygeoprocessing

texts:
While using pygeoprocessing rasterize() function to convert the shapefile into raster format I got the following error:

0warning 1: failed to find field wt on layer pointshp, skipping.

My script is as follows:

def vec2rast():

    source = ogr.Open("G:\Data\shashankdata\pointshp\pointshp.shp")
    layer = source.GetLayer()
    schema = {}
    ldefn = layer.GetLayerDefn()
    for n in range(ldefn.GetFieldCount()):
        fdefn = ldefn.GetFieldDefn(n)
        schema[n] = fdefn.name
    print schema
    index = int(input("Enter index number corresponding to field"))
    pg.geoprocessing.rasterize("G:\Data\shashankdata\pointshp\pointshp.shp",
                               "G:\Data\shashankdata\LSMDATASET\MOSIAC\MOSIAC.tif", burn_values=[255],
                               option_list=["ATTRIBUTE= %s" % schema[index]], layer_index=0)


texts:
Using Leaflet 1.0.3, I'm trying to create a legend with L.control() and Leaflet plugin MakiMarkers. Can't seem to get the Maki Marker icon to show up with description. Pretty sure it's my lack of proficiency with HTML and innerHTML. Any help would be great, I've spent way to much time on this.

//Maki Markers,  markers are declared globall,

var greenTree = L.MakiMarkers.icon({
    icon: "park",
    color: "3F9110",
    size: "s"
});

var redTree = L.MakiMarkers.icon({
    icon: "park",
    color: "F04441",
    size: "s"
});

var blueTree = L.MakiMarkers.icon({
    icon: "park",
    color: "442DB5",
    size: "s"
});

//    Legend Control

var legend = L.control({position: 'topright'});

legend.onAdd = function (map) {
    var div = L.DomUtil.create('div', 'legend');
div.innerHTML +=  '' + '
BCT Owned Land
'
div.innerHTML +=  ''  + '
Conservation Restriction on Private Land
'
div.innerHTML +=  ''   + '
Conservation Restriction on Town Land
'

    return div;
};

legend.addTo(map);


texts: How to get selected Features/Geometry from Windows Clipboard in ArcMap/C#?

texts:
I have a custom tool in ArcMap 10. The Workflow is the following:

Select Features in a Layer with the "Select Features by Rectangle" tool
Copy the Features to the Clipboard using STRG+c
Select the custom tool
The tool checks the Clipboard. In case of finding a copied Feature in the Clipboard the Workflow changes

I found an example on http://kiwigis.blogspot.com/2011/07/how-to-create-drop-target-for-layers-in.html where similar stuff is implemented using drag and drop.

Following and adapting the example and trying to copy layers in the TOC in ArcMap to the clipboard using STRG+c is working for me:

//Get the Data from the Clipboard
IDataObject clipBoardDataObject = Clipboard.GetDataObject();

//The format of the object data is "ESRI Layers" (this can be checked using clipBoardDataObject.GetFormats();)
MemoryStream geomStream = clipBoardDataObject.GetData("ESRI Layers") as MemoryStream;
byte[] bytes = geomStream.ToArray();

IMemoryBlobStreamVariant memoryBlobStreamVariant = new MemoryBlobStreamClass();
memoryBlobStreamVariant.ImportFromVariant(bytes);
IMemoryBlobStream2 memoryBlobStream = memoryBlobStreamVariant as IMemoryBlobStream2;
IStream stream = memoryBlobStream as IStream;
IObjectStream objectStream = new ObjectStreamClass();
objectStream.Stream = stream;

byte pv;
uint cb = sizeof(int);
uint pcbRead;

objectStream.RemoteRead(out pv, cb, out pcbRead);
int count = Convert.ToInt32(pv);

//Guid you find in the ESRI ILayer Interface
Guid iLayerGuid = new Guid("34C20002-4D3C-11D0-92D8-00805F7C28B0");

for (int i = 0; i < count; i++)
{
  object o = objectStream.LoadObject(ref iLayerGuid, null);
  ILayer layer = o as ILayer;
}

I can cast the object to ILayer, this is working!

Trying to achieve the same with copied Features is not working for me, at least I don't really know what to do.

//Get the Data from the Clipboard
IDataObject clipBoardDataObject = Clipboard.GetDataObject();

MemoryStream geomStream = clipBoardDataObject.GetData("ESRI Geometry List") as MemoryStream;
byte[] bytes = geomStream.ToArray();

IMemoryBlobStreamVariant memoryBlobStreamVariant = new MemoryBlobStreamClass();
memoryBlobStreamVariant.ImportFromVariant(bytes);
IMemoryBlobStream2 memoryBlobStream = memoryBlobStreamVariant as IMemoryBlobStream2;
IStream stream = memoryBlobStream as IStream;
IObjectStream objectStream = new ObjectStreamClass();
objectStream.Stream = stream;

byte pv;
uint cb = sizeof(int);
uint pcbRead;

objectStream.RemoteRead(out pv, cb, out pcbRead);
int count = Convert.ToInt32(pv);

//What to do next? What object do I need to unpack to?

You see here that the object format is "ESRI Geometry List" (of which I can't find any information online). So how can I "unpack" this object? IFeatureSelection, IFeature, IGeometry wont work.

Any Ideas?


texts: Using GeoPandas to find crossing azimuths

texts:
I have got a series of UTM coordinates and azimuths like that:
Point one: Z23, 603484.00m E,7781714.00m S, Azimuth 222°
Point two: Z23, 600593.00m E,7781310.00m S, Azimuth 137°
Those azimuths (should) cross in someplace. I'd like to get this crossing point coordinates.
Is it possible to do that using GeoPandas? If I don't, how can I get this using another Python tool?


texts:
I would like to have a tooltip for each of my layers based on their display which is done by a checkbox.

I started a function with a condition depending on the layer but I have the error "layer is not defined".
How can I declare a variable that retrieves the name of my layers to be able to filter for the tooltips next?

function displayTooltip(evt) {
    var pixel = evt.pixel;
    var map = $('#map').data('map');
    var feature = map.forEachFeatureAtPixel(pixel,
     function(feature, layer) {
       return feature;
        });

    var tooltip = document.getElementById('tooltip');

    var overlay = new ol.Overlay({
       element: tooltip,
        offset: [0, 0],
        positioning: 'bottom-left'

    });

    map.addOverlay(overlay);

    tooltip.style.display = feature ? '' : 'none';

    if (layer ===  missionsEAU | missionsEAU2) {
        overlay.setPosition(evt.coordinate);
        tooltip.innerHTML = feature.getProperties().DOCUMENT_QUOTE;

    } else if(layer ===  scot_admin){
        overlay.setPosition(evt.coordinate);
        tooltip.innerHTML = feature.getProperties().SCOT_2016;

    } else if(layer ===  epci_admin){
        overlay.setPosition(evt.coordinate);
        tooltip.innerHTML = feature.getProperties().EPCI_LATEST;

    } else if(layer ===  commune_admin){
        overlay.setPosition(evt.coordinate);
        tooltip.innerHTML = feature.getProperties().NOM_COMMUNE_LATEST;
    }
};


texts: When I open a QGIS file, QGIS 1.8 resets the scale factors for all labels.  Is there a way to save the file so that it knows to also save the label scales?


texts: I converted a KML file into a PostGIS db using ogr2ogr. How do I publish all the tables as seperate layers?

texts: As the title says I converted/imported a KML file into my PostGIS database using the ogr2ogr command. Now I am trying to publish each table as a seperate layer onto a GeoServer. I know about cURL, but as far as I can see I would have to have a new command for each table in the db. I have 124 tables I need to publish so doing that one by one would be labor intensive. Is there anyway straight-forward way to streamline this process?


texts: Accessing ArcGIS 9.3 geodatabase from 9.2?

texts:
If I'm using ArcGIS 9.2 is it possible to access a feature class from a geodatabase that was created by ArcGIS 9.3?

If not, are there any other ways to access the data?


texts:
I try setup a QGIS debug environment under Windows10 64bit.
I followed the instrictions given at https://github.com/qgis/QGIS/blob/master/INSTALL#L423

I also had a look at
https://www.shaeffer.co/compiling-qgis-on-windows/
and
https://www.cursosgis.com/compilar-y-debuggear-qgis-en-windows/

When I try to start debugging in Visual Studio I get following error message:

Could not load qgis_app.dll
module not found
Check D:\Projects\QGIS\ms-windows\osgeo4w\build-qgis-test-x86_64\output\bin\RelWithDebInfo\qgis.env for correct environment paths.

The qgis_app.dll exists in the QGIS\ms-windows\osgeo4w\build-qgis-test-x86_64\output\bin\RelWithDebInfo folder.

During the build there is a crash of the crssync.exe which is called during the build:

Problemsignatur:
  Problemereignisname:  APPCRASH
  Anwendungsname:   crssync.exe
  Anwendungsversion:    0.0.0.0
  Anwendungszeitstempel:    5ddd3d2f
  Fehlermodulname:  qgis_core.dll
  Fehlermodulversion:   6.3.9600.19464
  Fehlermodulzeitstempel:   5d6727f2
  Ausnahmecode: c0000135
  Ausnahmeoffset:   00000000000ecf30
  Betriebsystemversion: 6.3.9600.2.0.0.272.7
  Gebietsschema-ID: 1031
  Zusatzinformation 1:  ac05
  Zusatzinformation 2:  ac0507478d1c5bd693cfc4fe3987e900
  Zusatzinformation 3:  ac05
  Zusatzinformation 4:  ac0507478d1c5bd693cfc4fe3987e900

What can I do to narrow down the source of the error and get Visual Developer to build the ALL_Build successfully?

The error messages vary. often the module syncsrsdb and qgis_core fail but sometimes also other modules cause trouble:
One example for error messages:

    129>------ Build started: Project: qgis_authcryptotest, Configuration: RelWithDebInfo x64 ------
127>       CustomBuild:
127>         Running crssync
127>     1>
127>C:\Program Files (x86)\MSBuild\Microsoft.Cpp\v4.0\V140\Microsoft.CppCommon.targets(171,5): error MSB6006: "cmd.exe" exited with code -1073741515.
127>
127>Build FAILED.

C:\Program Files (x86)\Windows Kits\10\Lib\10.0.14393.0\um\x86\Version.Lib : warning LNK4272: library machine type 'X86' conflicts with target machine type 'x64'

C:\Projects\QGIS\ms-windows\osgeo4w\build-qgis-test-x86_64\output\bin\RelWithDebInfo\qgis_core.dll : fatal error LNK1120: 3 unresolved externals

One message looks like a problem caused by a mixture between 32bit and 64 bit components

("C:\Program Files (x86)\Windows Kits\10\Lib\10.0.14393.0\um\x86\Version.Lib : warning LNK4272: library machine type 'X86' conflicts with target machine type 'x64'")

but I could not find any special 64 bit versions of the Windows SDK


texts:
I looked into using getFlags() for this as described here,
Get road type (clazz) in osm2po
but for some reason I'm getting 65551 instead of 15, and 65557 instead of 21. Seems like if I unset the 16th bit, I get the right value (rrs.getFlags() & ~(1 << 16)). Is there some significance to that 16th bit, and would there be some other bits set? Do we just take the lower 16 bits? I'm using osm2po 5.0.0.


texts: Exporting only certain columns to CSV file in ArcGIS for Desktop?

texts:
I've written a python script using arcpy which outputs a polygon feature class into a File Geodatabase.  I have added a function to export the attributes to a separate CSV file.  I am using the code which I have found in this post which works perfectly.  However, that code exports every column in the feature class.  I only want to export the fields that do not have the following names: OBJECTID, Shape, or Shape_Length.

My CSV file generates successfully, and it correctly does not include the OBJECTID or Shape_Length fields.  However, the Shape field is written to the file.  An example value that is written to that field is:

<geoprocessing describe geometry object object at 0x28CB90A0>

I have added a line to print the field names as it iterates through them and surprisingly, Shape is not printed.  It's as if ArcGIS is hiding it or gives it a different name.

The code for my function is below:

def exportToTable():
    """
        Exports the final outputs to a CSV File.
    """

    # Create path to CSV File (note the varialbe outputPath is declared elsewhere).
    CSVFile = outputPath+'\\FinalOutput.csv'
    arcpy.AddMessage("Created CSV File: %s" %CSVFile)

    # Get all fields in FinalOutput feature class and remove unwanted fields.
    fields = arcpy.ListFields('FinalOutput')
    for field in fields:
        arcpy.AddMessage("Field.name is:"+field.name) #not printing 'Shape' field name
        if field.name in (['OBJECTID', 'Shape', 'Shape_Length']):
            fields.remove(field)

    i = 1
    f=open(CSVFile, 'w')
    for field in fields:
        #--write the wanted field names to the output file
        if i < len(fields):
            f.write('%s,' % field.name)
            i += 1
        else:
            f.write('%s\n' % field.name)

    # Use a search cursor to iterate through the rows of the table and write them to the CSV file.
    rows = arcpy.SearchCursor('FinalOutput')
    for row in rows:
        i = 1
        for field in fields:
            if i < len(fields):
                f.write('%s,' % row.getValue(field.name))
                i += 1
            else:
                f.write('%s\n' % row.getValue(field.name))
    del rows
    f.close()

Does anyone know what's going on here?

I modified my code to follow the advice of @sgrieve and it was still writing the Shape field.  If I add a line to print the field names as it iterates through them, it lists all the fields except the Shape field, yet it still writes to the CSV.  It also added the X and Y coordinates of the polygon as two new columns and the columns no longer aligned with the column names.

I've modified the line where @sgrieve declares the fields to be the following:

fields = [f.name for f in arcpy.ListFields('FinalCadastre') if f.type <> 'Geometry']

The new code works fine, but I'm still unsure what the problem was.  Does anyone know what was going on?  What's the deal with the Shape field?


texts: How do I automatically export a release of the EPSG registry?

texts:
I am trying to use the EPSG Registry API to export a release
of the registry for offline use but am receiving a registry error.

After following the Developer Guide I am HTTP POSTing the following XML request to the API URL:

<?xml version="1.0" encoding="UTF-8"?>
<GetRecords
    xmlns="http://www.opengis.net/cat/csw"
    xmlns:ogc="http://www.opengis.net/ogc"
    startPosition="1" maxRecords="100"
    outputFormat="application/xml; charset=UTF-8"
    resultType="results">
  <!--
    Retrieve the ReleaseObject releated to the current Version History record.
    Then use the ReleaseObject identifier to obtain the actual compressed
    GML Dictionary file containing the EPSG Dataset via a getRepositoryItem
    request.
  -->
  <Query typeNames="Association=a RegistryObject=release">
    <ElementSetName typeNames="release">full</ElementSetName>
    <Constraint version="1.1.0">
      <ogc:Filter>
        <ogc:And>
          <!-- find the correct association -->
          <ogc:PropertyIsEqualTo>
            <ogc:PropertyName>$a/targetObject</ogc:PropertyName>
            <ogc:Literal>urn:ogc:def:version-history:EPSG::6.15</ogc:Literal>
          </ogc:PropertyIsEqualTo>
          <ogc:PropertyIsEqualTo>
            <ogc:PropertyName>$a/associationType</ogc:PropertyName>
            <ogc:Literal>
              urn:x-ogp:def:AssociationType:EPSG:ReleaseFor</ogc:Literal>
          </ogc:PropertyIsEqualTo>
          <!-- now relate it to the ReleaseObject -->
          <ogc:PropertyIsEqualTo>
            <ogc:PropertyName>$a/sourceObject</ogc:PropertyName>
            <ogc:PropertyName>$release/id</ogc:PropertyName>
          </ogc:PropertyIsEqualTo>
       </ogc:And>
      </ogc:Filter>
    </Constraint>
  </Query>
</GetRecords>

The XML is copied from section 5.3 'How To... Export a Release' of the guide, and the URL to send GetRecords requests to is suggested by section 4.1.1 'GetRecordById'. The XML error response I get is as follows:

<?xml version="1.0"?>
<ows:ExceptionReport xmlns:ows="http://www.opengis.net/ows" version="2.0">
  <ows:Exception exceptionCode="InvalidRequest">
    <ows:ExceptionText>XML Schema validation error</ows:ExceptionText>
    <ows:ExceptionText>
Error:Error
 Line:    2
 Column:  193
 Message: cvc-elt.1: Cannot find the declaration of element 'GetRecords'.
    </ows:ExceptionText>
  </ows:Exception>
</ows:ExceptionReport>

I am using the curl command line program to run the query as follows:

curl -H 'Content-Type: appliation/xml' --data @request.xml 'http://www.epsg-registry.org/indicio/query'

Can you see where I'm going wrong? If not is there any other way to
automatically download the registry? I am aware that I can login to the Registry HTML interface and download the registry from there but I
would like to automate the process.


texts: Creating Mosaic with multiband layers using QGIS?

texts:
I have twice scenes of landsat and I would like to create a mosaic with multiband layers with them. I have tried to do it with Qgis. However, the geoalgorithm which qgis has not able me to do it. Here is the message that I received before ran saga merge tool.

How can I do a physic mosaic with multiband raster layers using qgis?


texts: Erase area from GeoTIFF covered by area in another GeoTIFF

texts:
I have 2 GeoTIFFs, A and B and I want to produce GeoTIFF C:

Both A and B have 4 bands (RGBA) and the same projection
A covers large area and has lower resolution (bigger pixel size) than B
B covers smaller area, all within area of A
B has nodata areas

I would like to produce C which would be the same as A but with all the areas that are in B as non-nodata (ie. has data) to be set to fully transparent (band 4 and optionally all other bands to zero).
In other words I want to produce C = A - B where the "hole" created by B would be RGBA=0,0,0,0.
I need to do it on Linux (GDAL, QGIS, ...). Note that both files are pretty big (A - 56GB, B - 295GB).
Any ideas? As a last resort I would maybe use gdal_calc but I am not sure how would it handle different pixel sizes.


texts: Commenting Out Definition Queries in ArcGIS Desktop?

texts:
I use Definition Queries a lot in my daily GIS life, but one of the many things I've found annoying about ArcMap is the inability to comment out Definition Queries. I'd like to have this option as sometimes I want to turn on/off the Definition Query or I want to have multiple Queries on a single Shapefile, but only one active at a time.

So far to get around this, I usually just have a notepad document in the background and copy and paste as necessary. But the ideal solution would be to use the definition query tab in the layer properties to replace this document. I might have multiple different queries on a layer, but all commented out.

As I understand, Definition Queries are written in SQL. To comment in SQL I believe you either start a line with: -- or encompass a statement in brackets: { }. I've tried both of these:

Original: "DATA_AQUISITION" LIKE 'MISSING XSECTION'
Dashes:   --"DATA_AQUISITION" LIKE 'MISSING XSECTION'
Brackets: {"DATA_AQUISITION" LIKE 'MISSING XSECTION'}

The last two statements return a SQL expression error and won't draw in ArcMap. I'm not looking for a programming solution, as that would just be too much work.


texts: How to do "random" distribution of certain classes to certain amount of polygons with respect to spatial constraints

texts:
I have a polygon layer.

Now, I would like to assign/allocate/classify these Polygon into 13 different classes respecting spatial constraint, in my case: neighboring polygons (sharing line or nodes) are not allowed to have the some class.

How can I do this in ArcPy?


texts: Adding item description from PostgreSQL/PostGIS data?

texts:
Is it possible to have metadata or an item description with a dataset from a PostgreSQL/PostGIS database without using SDE?

I now see an XML-error when trying to view the item description and there is no possibility to add or edit it.


texts:
I want to use OSM data to make network dataset in ArcGIS 10.2.
However, this is not installed successfully. It keeps telling me I don't have ArcGIS 10.2 that is exactly what I am using now.
I think others have same issues.
Anyone help me to resolve this problem?
Or is there any way to make network dataset without ArcGIS Editor 10.2 for OSM?

Further information for this issue, I am using 10.2.2 desktop with 64 bit. And this is concurrent user version supported by my school for the researcher.
I can use all of functions and extensions with this privilege.
Also, I followed the direction indicated below:
http://www.arcgis.com/home/item.html?id=16970017f81349548d0a9eead0ebba39

As you can see comments in the website, other people have same issues like me.


texts: Does OSM data contain direction of travel?

texts:
I was looking for Direction of Travel (DOT) information inside OpenStreetMap shapefiles, but only I've just found "ONEWAY = 1 or 0", but nothing about its direction.

If I work on an OSM editor, I can see its arrows which declare its direction. But as you can see, once downloaded as Shapefile type, there is no DOT information about.


texts:
I have LiDAR point cloud from 26 X 1 ha forest plots. I am trying to calculate roughness index for different vegetation layers in each plot with following code using lidR package.

files <- list.files(path= "FilesDirectory/", pattern= "*.laz", full.names = TRUE, recursive = FALSE)
O = lapply(files, function(x) {

  las = readLAS(x, select = "xyzc", filter = "-drop_z_below 0")

  las_l1 = lasfilter(las, las$Z <= 1.0) #first vegetation layer <= 1 m
  las_l2 <- lasfilter(las, las$Z > 1.0 & las$Z <=8.0) # 8.0 >= second layer >1
  las_l3 <- lasfilter(las, las$Z > 8.0) # and third layer > 8.0

  chm_l1 = grid_canopy(las_l1, 0.2, p2r())
  chm_l2 = grid_canopy(las_l2, 0.2, p2r())
  chm_l3 = grid_canopy(las_l3, 0.2, p2r())

  roughness_l1 <- rumple_index(chm_l1)
  roughness_l2 <- rumple_index(chm_l2)
  roughness_l3 <- rumple_index(chm_l3)

  return(data.frame(file=x, roughness_l1=roughness_l1,roughness_l2=roughness_l2,
                    roughness_l3=roughness_l3))
})

However, this code is returning error

Error: The point cloud contains 0 point

I think the cause of error could be that in some files there may not be any point on the third layer. Therefore, I want to include if condition that if the the third layer doesn't have any point print ("no points on third layer") otherwise calculate chm_l3 for the third layer with following code.

if (length(las_l3)==0) {"no points on third layer"} else {chm_l3 = grid_canopy(las_l3, 0.2, p2r())}

But it didn't help.


texts: QGIS 2.18 Add partial file path to composer

texts:
How to add the current project path in Print Composer? doesn't quite do what I am hoping for as it gives the whole filepath (which the client doesn't need/want to see) - is there any way of just getting the last few folders?

At the  moment, using
[%concat ( @project_path,  @project_title )%] I get this:

  /media/jonathan/3C5VCG82K32FDB/Data/GIS/Projects/0020_Client/0019_Mynameforclientproject/Projects/QGISProjectfilename.qgs

when I just want

  /0020_Client/0019_Mynameforclientproject/Projects/QGISProjectfilename.qgs

Is this possible?


texts: Looking for a clip showing how cool GIS is

texts:
Next year I will be teaching a GIS course since this is not a compulsory course I have to make the students register. I thought of producing an e-mail brochure with a few lines of text and a clip showing what is GIS and how cool it can be.

Does anyone have an idea to such clip?


texts: Retrieve features from WFS using geotools

texts:
I am new to GI systems. I am trying to retrieve features from a WFS using geotools:

this is my snippet :

String url = "http://geodaten.metropoleruhr.de/poi/poi?request=GetFeature&service=wfs&version=1.1.0";
Map<String, Object> connectionParameters = new HashMap<String, Object>();
connectionParameters.put("WFSDataStoreFactory:GET_CAPABILITIES_URL", url );

//connection
DataStore data = DataStoreFinder.getDataStore(connectionParameters);

but I get NullPointerException.
Do I miss something ?


texts: Removing feature class or table from TOC using C#?

texts:
What I want to achieve is to remove specific layers from the TOC either based on their array place [n] or their name (which would be ideal in case their place in the TOC changed).

The only thing I can find is the ClearLayers() method to clear all of the layers from the TOC.

Also I am aiming to remove a featurelayer and a table.


texts: Algorithm for estimating position of many objects and differentiating those?

texts:
Let's assume a crowd-sourced tree location use case where a research team asks hikers to GPS the location of a specific tree type in case they encounter it.

This will lead to a large number of longitude-latitude points each identifying the measured (noisy) GPS location of a tree. The researchers would now like to compute from those measurements a map of those trees.

That means they have to differentiate location point estimates for close trees and after assigning the estimates to a specific tree calculate the location of that tree.

Now trees don't move and GPS measurements can be pretty inacurate. What clustering algorithm would be the canoncial choice for this application?

(If you have an R package::function to suggest please tell me)


texts: How to get EPSG:4326 coordinates on Drag and drop in OpenLayers

texts:
I followed the example on how to implement a drag and drop interaction with openLayers.
I now want to read the coordinate of every features like the following
this.dragAndDropInteraction.on('addfeatures', (event: DragAndDropEvent) => {
      const features: Feature<Geometry>[] = event.features

      features.forEach((feature, index) => {
        const properties = feature.getProperties()
        const coordinates = properties?.geometry

        console.log(coordinates.flatCoordinates)
      })
    })

But I get [799785.5152789096,...] and I wanted to have it with the EPSG:4326 format.
I did try to change the projections when I initialized the map like this
 const dragAndDropInteraction = new DragAndDrop({
    formatConstructors: [<any>GPX, KML],
    projections: 'EPSG:4326',
  })

this.map = new Map({
  projections: 'EPSG:4326',
  interactions: defaultInteractions().extend([this.dragAndDropInteraction]),
  target: target,
  layers: [
    new TileLayer({
      source: new OSM(),
    }),
  ],
  view: this.view,
})

But I keep get the wrong coordinate.
Is it normal?
I'm importing a KML file
Here is his content
<?xml version="1.0" encoding="UTF-8"?>
<kml xmlns="http://www.opengis.net/kml/2.2"
  xmlns:gx="http://www.google.com/kml/ext/2.2"
  xmlns:kml="http://www.opengis.net/kml/2.2"
  xmlns:atom="http://www.w3.org/2005/Atom">
  <Document>
    <name>Traces du passé</name>
    <open>1</open>
    <Style id="default0">
    </Style>
    <StyleMap id="default">
      <Pair>
        <key>normal</key>
        <styleUrl>#default0</styleUrl>
      </Pair>
      <Pair>
        <key>highlight</key>
        <styleUrl>#hl</styleUrl>
      </Pair>
    </StyleMap>
    <Style id="hl">
      <IconStyle>
        <scale>1.2</scale>
      </IconStyle>
    </Style>
    <Folder id="kml_ft_Sentier">
      <name>Sentier</name>
      <Placemark>
        <name>31</name>
        <description>La porte de Fully</description>
        <LookAt>
          <longitude>7.184550738826426</longitude>
          <latitude>46.17214547339596</latitude>
          <altitude>0</altitude>
          <heading>-0.0005252280984056049</heading>
          <tilt>60.00131024332176</tilt>
          <range>168.3472321103629</range>
          <gx:altitudeMode>relativeToSeaFloor</gx:altitudeMode>
        </LookAt>
        <styleUrl>#default</styleUrl>
        <Point>
          <gx:drawOrder>1</gx:drawOrder>
          <coordinates>7.184330180845047,46.17181721801122,0</coordinates>
        </Point>
      </Placemark>
    </Document>
  </kml>


texts:
I have made a map using the Leaflet Draw plugin which allows the user to download their drawn items. These drawnItems are exported as GeoJSON using the following code from here:

document.getElementById('export').onclick = function(e) {
// Extract GeoJSON from featureGroup
var data = featureGroup.toGeoJSON();

// Stringify the GeoJSON
var convertedData = 'text/json;charset=utf-8,' +
encodeURIComponent(JSON.stringify(data));

// Create export
document.getElementById('export').setAttribute('href', 'data:' + convertedData);
document.getElementById('export').setAttribute('download','data.geojson');
}

This works perfectly, but it would be even more ideal if the GeoJSON was converted to .kml before being exported. I am aware of the toKml plugin but I am struggling to get it to work (I am still quite new to all of this). Where would I add:

var kml = tokml(geojsonObject);

And I am sure I would have to edit the //create export part of the GeoJSON code


texts: How to display charts next to features in the map?

texts:
This question is meant to be primarily on QGIS software but all other software solutions are welcome as far as it is possible to achieve it as easy as possible, preferably out of the box. Though the limitation is static map (to be printed) so no interactive web map solutions applicable here.

Is there a way how to create charts (tool / plugin in QGIS) for time series data / table for points (can be also lines / polygons) and display it next to features?

I am aware of simple diagram / histogram tool. I also got an idea to load charts as SVG but that means I need to export a lot of charts first and then link them as symbols manually. Is there some more profound out of the box solution?

What I tried so far looks like this:

What I would like to achieve looks like this:


texts: Raster centroid and relationships with geoTIFF

texts:
I have a list (s) containing information on the probable locations of many animals in South America. For example, this is the type of stored information and what it looks like when plotted for the first individual.

Example:

> s[1]
[[1]]
class       : RasterLayer
dimensions  : 418, 313, 130834  (nrow, ncol, ncell)
resolution  : 0.16666, 0.16666  (x, y)
extent      : -86.333, -34.16842, -55.91633, 13.74755  (xmin, xmax, ymin, ymax)
coord. ref. : NA
data source : in memory
names       : layer
values      : 0, 1  (min, max)
> plot(s[[1]])

Note: the green areas are all "likely" locations and the grey areas are "unlikely" locations.

I also have spatial data in geoTIFF format from NASA.

I would like to determine the centroid of the likely areas for each individual, and identify the value for the GeoTIFF spatial data that matches with the centroid (i.e., likely a different value for each individual).

UPDATED
To determine the centroid, I've tried the following:

indi1<-s[[1]] #restrict list to 1st individual
colMeans(xyFromCell(indi1, which(indi1[]==1)))
        x         y
-52.62926 -18.16871

However, even though this just uses coordinates where the value was "likely" (1's), you can see that the centroid falls in an area that is "unlikely".

par(mfrow=c(1,2))
plot(indi1)
points(colMeans(xyFromCell(indi1, which(indi1[]==1)))[1],
       colMeans(xyFromCell(indi1, which(indi1[]==1)))[2],
       col="red")
plot(indi1, xlim=c(-55,-50), ylim=c(-20,-15))
points(colMeans(xyFromCell(indi1, which(indi1[]==1)))[1],
       colMeans(xyFromCell(indi1, which(indi1[]==1)))[2],
       col="red")

So I still need a solution for restricting the final centroid value to an area that is "likely".  And, ultimately, associating this location with a value in the GeoTIFF.

Note: this question is related to another question I have here.


texts: Summing two columns in ArcGIS Desktop, one of them has Null values

texts:
I need to sum values of 2 columns in field calculator and populate a new column with the result, one of the values always in Null. I already had found a useful answer In here, but it doesn't help with the last couple of columns, that contain a date in DD.MM.YYYY (string) format- geoprocessing report shows successful operation, but the new column stays empty. Is there a small tweak i could use for the code, to make it work?

def stack(Y, Y_1):
  itemList = [Y, Y_1]
  myList = [Y for Y in itemList if (Y != None)]
  return sum(myList)

Field value= stack( !Y!, !Y_1!)


texts:
I am trying to display a netCDF file that I obtained from NOAA in ArcMap. However, it keeps crashing ArcMap. I have used precipitation netCDFs from NWS before with no issues.

The netCDF is just a sample they sent me of many netCDF files they will start making available as part of a new lightning data release. So I would like to write a script that corrects the netCDF files.

They told me there are "some issues with data formatting (e.g., unsigned integers for event lat and lon) and that it won't be fixed for some time."

I used the following code to open the netCDF to look at the data. I only included the coordinates output and not the data variables and attributes.

xr.open_dataset('C:\\lightning\\sample.nc'

<xarray.Dataset>
Dimensions:                                 (number_of_events: 5879,
number_of_field_of_view_bounds: 2, number_of_flashes: 687, number_of_groups:
3454, number_of_time_bounds: 2, number_of_wavelength_bounds: 2)
Coordinates:
    event_id                                (number_of_events) int32 -940849349 ...
    event_time_offset                       (number_of_events) datetime64[ns] 2017-09-08T19:30:19.917000 ...
    event_lat                               (number_of_events) float64 -1.687 ...
    event_lon                               (number_of_events) float64 -188.5 ...
    event_parent_group_id                   (number_of_events) int32 1324574251 ...
    group_id                                (number_of_groups) int32 1324574251 ...
    group_time_offset                       (number_of_groups) datetime64[ns] 2017-09-08T19:30:19.834000 ...
    group_lat                               (number_of_groups) float32 -1.68695 ...
    group_lon                               (number_of_groups) float32 -55.3391 ...
    group_parent_flash_id                   (number_of_groups) int16 12590 ...
    flash_id                                (number_of_flashes) int16 12590 ...
    flash_time_offset_of_first_event        (number_of_flashes) datetime64[ns] 2017-09-08T19:30:19.834000 ...
    flash_time_offset_of_last_event         (number_of_flashes) datetime64[ns] 2017-09-08T19:30:19.834000 ...
    flash_lat                               (number_of_flashes) float32 -1.68695 ...
    flash_lon                               (number_of_flashes) float32 -55.3391 ...
    product_time                            datetime64[ns] 2017-09-08T19:30:20
    lightning_wavelength                    float32 777.4
    group_time_threshold                    float32 0.0
    flash_time_threshold                    float32 3.33
    lat_field_of_view                       float32 0.0
    lon_field_of_view                       float32 -89.5

I then used np.savetxt to save event_lat, event_lon, group_lat, and group_lon to CSV files so I could ananlyze the coordinates.

What I found was that group_lat values ranged from 42 to -51 and group_lon values range from -44 to -142. These are correct because the file is supposed to display data from North and South America and parts of adjacent oceans.

But, event_lat values ranged from -0.03 to -133 and event_lon values ranged from -90.2 to -222.29. Looking at event_lat CSV, there is a chunk of data that ranged from -.03 to -51 (this should be data for southern hemisphere) and a chunk of data ranging from -90.5 to -132 (this should be the northern hemisphere data). However, that range should actually be 0.5 to 42 latitude.

I'm not sure what's going on with the event_lon data. -90 to -142 coordinates look similar between event_lon and group_lon, so I think that is fine. However, event_lon also has data ranging from -177 to -222 and no data greater than -90. The group_lon data ranges from -44 to -142.5.

Is there a way I can fix these files with python?


texts: How to handle Update error: Transaction does not have authorization for mytable3:mytable3.5?

texts:
I am building collaborative mapping app using openlayers 5.3.0, geoserver 2.14.0 and postgresql 10.4.
I want the user to be able to update features using wfs-t but before updating, the user needs to lock the feature so that no other users can update this feature before the first transaction committed to the db (actually don't know if this should be handled manually or geoserver will handle it itself)...

This is the lock feature request:

$.ajax('http://localhost:8080/geoserver/wfs', {
    type: 'POST',
    data: {
      service: 'WFS',
      version: '1.1.0',
      request: 'LockFeature',
      typename: 'mytable3',
      cql_filter: 'id=' + f.values_.id
    }
  });

And after this request Iam trying to update the feature using this request:

$.ajax('http://localhost:8080/geoserver/TOTO1/mytable3/wfs', {
  type: 'POST',
  dataType: 'xml',
  processData: false,
  contentType: 'text/xml',
  data: str,
}).done(function () {
  WFSsource.clear();
});

the update feature request without the lock works properly but when I use the lock the update feature request return with:

    <ows:ExceptionText>Update error: Transaction does not have authorization for mytable3:mytable3.5</ows:ExceptionText>


texts: What resampling technique should be used to reproject an altitude raster?

texts: Does anyone know whether I should choose: nearest, bilinear, cubic or majority when projecting the WorldClim altitude raster layer?


texts: When is it worthwhile to reproject in QGIS?

texts:
My QGIS project has various raster and (complex) vector layers from different sources, in different CRSes. Using reprojection on the fly, I don't have to reproject them. But should I, especially for the layers I've decided are part of my base map set, for performance or other reasons?

In particular, is the performance overhead associated with ongoing on-the-fly reprojection of a shading (DEM-derived) raster layer, or a complex and dense vector layer such as contours or vectorized waterbodies, worth manually reprojecting and saving in the displayed CRS?

Editing to add: my goal is a usable, responsive map that will continue to be improved with adding new information, rather than analysis of embedded GIS data.


texts: Importing KML and style information into Geoserver

texts:
So I'm working on a project that's been in development for over a year, but I only recently got involved. My task involves converting roughly 2,500 KML files into another format, importing these files into a Geoserver, and then copying over the associated style information for each file. I realize that Geoserver allows you to import KML files directly, however, our files never seem to render correctly, and my boss would like to move away from KML format for philosophical reasons.

Using ogr2ogr, I've tried to convert these files into ESRI Shapefiles, however the conversion process only worked on roughly a third of the files and erased some of the metadata called for in the style information. The files that didn't convert code for polygons and points, which I've read can be an error for other programs like QGIS.

I've also tried to put all of the KML files into a PostGIS database and convert them like that. However, that seems to delete all of the data leaving essentially blank tables for import.

TL;DR I need to find a way to convert 2,500 KML files into a format usable by Geoserver, and will work with XML style information. I'm a complete GIS noob and this is my first time posting, so I'm sorry if this seems a bit ridiculous.


texts: How to select feature first and then start the zoom animation

texts:
I'm trying to select a feature first and then start the zoom animation but openlayers (v6) first start the zoom animation and when it's done the feature gets selected, what's going on here?

Here is my code:

function zoom(){
    selectSingleClick.getFeatures().push(featToZoom);

    var vectorSource = new ol.source.Vector();
    vectorSource.addFeature(featToZoom);

    map.getView().fit(vectorSource.getExtent(), {
            duration: 5000
    });
}


texts: Dynamic label generated by python in composer template?

texts:
I just recently started to play with python and items in print composer. Given item ID it is possible to do many great things with the item through code, one of them is changing text of text label.

My template contain text label for date, coordinate reference system / projection and I would add project file path if possible to update it automatically.

Is it somehow possible to save python code within the template so whenever the template is loaded / exported it will update items like date / CRS / project file path?


texts: Why do raster band values have different values than raster statistics in QGIS?

texts:
I am working on a raster file (1 Band).

According to the raster metadata, my minimum value is 0 and maximum is 7458. This seems to be confirmed by the Style options that return some category values in the value range mentioned above. In addition, these numbers make total sense with the nature of my data.

However, when I use the "Identify Features" tool, I get some values completely out of the range (these values look like this one generally: 4.5761380462188e-315).

I don't understand why and it's a problem because I would like sometimes to get the value of a single pixel. Do you have an idea why I have this discrepancy between the Identify Feature tool and the raster statistics?

I use QGIS 2.2.0-Valmiera on Windows 7.


texts: Is it possible to use ESRI:54009 (World_Mollweid) proj with rasterIo?

texts:
In a python environment I try to project from EPSG:4326 to ESRI:54009 to compute surfaces.
The ESRI description doesn't seem to be recognise by rasterio lib :
from rasterio.crs import CRS

CRS.from_string('ESRI:54009')

And I get the following error

CRSError: The WKT could not be parsed. OGR Error code 5

instead if I use pyproj :
from pyproj import CRS

CRS.from_string('ESRI:54009')

returns
<Projected CRS: ESRI:54009>
Name: World_Mollweide
Axis Info [cartesian]:
- E[east]: Easting (metre)
- N[north]: Northing (metre)
Area of Use:
- name: World
- bounds: (-180.0, -90.0, 180.0, 90.0)
Coordinate Operation:
- name: World_Mollweide
- method: Mollweide
Datum: World Geodetic System 1984
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich

Unfortunately the rasterio calculate_default_transform requires a rasterio CRS. Is there a way to use the World_Mollweid in this lib ? (as it's incompatible with gdal python binding, I cannot use a simple gdal.Warp)


texts: Matching buffer in Openlayers and same buffer in PostGIS query

texts:
I am using a buffer tool to create a buffer on my base map in Openlayers. In addition to creating and visualizing this buffer, I pass the center point and radius to a query. The query collects all features of a certain table which fall into this buffer, and shows them on the map.

The problem is that in these two cases buffers are generated separately and are not exactly the same. This causes the query to collect less (not all) features inside the buffer. Below is an image to clarify the situation. You can see the buffer created in Openlayers in Orange color, as well as features inside it that are collected from PostGIS in black and gold colors (different features).

Its easy to see how much difference is between these two buffers and the one that is generated by the query is more like an ellipsoid rather than circle.

How can I solve this?

The code for generating buffer in Openlayers:

var sides = 40;
var new_geom = OpenLayers.Geometry.Polygon.createGeodesicPolygon(centroid, radius, sides, 45, projection);

and I use the following WHERE clause inside my SELECT query:

WHERE ST_DWithin(ST_Transform(h.geom, 26918), (SELECT ST_Transform(ST_GeomFromText($center_buffer, 900913), 26918)), $radius_buffer)

more information: h is the alias of my table where I want to collect features that fall inside the circular buffer. and center_buffer and radius_buffer are the variables pass by OpenLayers (the same used for creating the buffer in OpenLayers).

As mentioned by a site user below, the city data is around Hamburg and therefore I changed my SRID in my query from 26918 to 26931. With this, the problem of strange ellipsoid shape of my features inside buffer is solved but there is still a big difference between the buffer and the data retrieved from PostGIS. See image below:


texts: On the fly projection not working in QuantumGIS

texts: I am trying to put a point layer of cities on top of a shapefile of counties in QuantumGIS.  I have on the 'On the Fly' turned on and both files are in the same projection system.  However, the point file is not showing up on the county file and both have different coordinate systems for some reason.  How to I get the layers aligned?


texts: How to extract values from intersection of different size buffers with raster?

texts:
I have a polygon shape with buffers. However, the buffers are of different size from 4m radius to 15m. In total I have 7 different buffers.

I have a landcover raster and I would like to extract values from that raster for each of my buffers. I would like to get either the proportion of each cover type within each buffer or the major cover type for each buffer.

I am using ArcGIS 10. I tried using GME and intersect polygon with raster and almost all my fields are treated as noData... Both my raster and my polygon files have the same coordinate system and projection.

I was about to use CON() for each landcover, then neighborhood and then extract values to point but there might be a faster way.

Does anyone have an idea?


texts: QGIS suddenly very slow generating atlas on PC?

texts:
I have been using happily QGIS for about 6 months, in the last week it has suddenly got much slower when generating an atlas using print composer.

I am running QGIS 2.10.1 on a PC laptop (Windows 7, core i7 2.4GHz with 16GB RAM).

The projects I've been working on have about 10-15 shapefiles (100-200MB each) and around 10 print composers.

Rendering in the main screen doesn't seem to have been affected but as soon as I try and generate an atlas from one of the print composers (even for only one atlas feature) everything is greyed out and 'not responding' for at least 20 minutes.
This occurred even when I set up a simple project with 5 shapefiles and one print composer with one atlas feature.

I've tried uninstalling and re-installing qgis but clearly didn't manage a clean re-install as the recent file list is still populated (also, it made no difference to performance). I find it hard to believe my projects can be stretching the limits of the software or my computer to this extent. The one thing I did recently change was I installed the mmqgis plugin but uninstlling had no effect on the reduced performance.

Has anyone had a similar problem with the print composer? and have suggestions as to how to fix it?


texts: How to generate zonal stats for multiple rasters using rasterstats Python

texts:
I'm trying to generate the zonal stats for multiple rasters without using ArcPy.

I've successfully obtained results using rasterstats for one image and one polygon. I  would like to use the same polygon but many rasters. My code is below and not working. Does anyone know what I'm doing wrong?

Lista_raster = [
   "NDVI_2019-05-05.tiff", "NDVI_2019-05-25.tiff", "NDVI_2019-06-04.tiff",
   "NDVI_2019-06-09.tiff", "NDVI_2019-06-14.tiff", "NDVI_2019-06-19.tiff",
   "NDVI_2019-06-24.tiff", "NDV_2019-05-05.tiff"
]
for raster in Lista_raster:
    stats = zonal_stats("Pivo_4_wgs.shp", Lista_raster,
            stats=['min', 'max', 'median', 'majority', 'std'],
            geojson_out=True)
result = {"type": "FeatureCollection","features": stats}
import json
with open('lista.geojson', 'w') as outfile:
   json.dump(result, outfile)


texts: Automating Creation of Network Dataset Using Polyline Input and ModelBuilder?

texts:
Is there a way to automate the creation of a Network Dataset in Model Builder using a line shapefile as an input? (ArcGIS 10.4)

The only documentation I have been able to find on creating a network dataset involves right clicking on a feature dataset > New > Network Dataset.

I have tried using the "Build Network" stand-alone tool, but it will not accept either a polyline or polylineM as an input. It seems like the "Build Network" tool simply re-builds any pre-existing network dataset.


texts:
I have a point featureclass with a point for every 25cm - in Total 400,000 Points.
I would like to reduce this resolution to one point every Meter. How could I instruct ArcGIS Advanced(ArcInfo) to do this? I can´t see a thinning tool or generalization tool for Point feature classes and don´t want to create a raster from the data yet to resample.


texts: What is the difference between spatial heterogeneity and spatial nonstationarity?

texts:
Until recently I though both terms refer to the same concept. Most scientific literature seems to use these terms also interchangibly.

However, here the authors say:

  Spatial heterogeneity refers to the uneven distribution of a variable’s values across space

  Spatial nonstationarity refers to variations in the relationship between an outcome variable and a set of predictor variables across space

So, who is right? Is there a difference between spatial heterogeneity and spatial nonstationarity? Is spatial nonstationarity a special case of spatial heterogeneity?


texts: Controlling grids assigned to multiband QGIS Virtual Raster

texts:
I see this has been asked 4 years ago but I see no useful answer for me

Multiband virtual raster with gdalbuildvrt: Control how bands are sorted

I do not understand how to control the RGB assignation order and what the end result is.

Using Virtual Raster to build an RGB and I cannot control which grid becomes which band. I cannot work out which grid ended up where because statistics in the 3 band VRT do not resemble those of the grids. I cannot tell going in because the Multiple Selection windows often lists a different order to the Layers window, which can be different to the alphabetical order.

Which grid goes where in an RGB is incredibly important in radiometrics but QGIS appears to be random in this case.

Please tell me I am wrong and how to do it.


texts:
I want to work on QGIS using GRASS HRU. it needs QGIS software to boot first, but it is not booting. I tried reinstalling the software also. I can open GRASS but its not giving any result.
I installed QGIS. GRASS plugin is already there. it is taking input but not giving any output. in manual it was written that when you will open QGIS it will boot. but what I have is not booting. I have QGIS in ILMS software.  `what can I do to resolve this problem?


texts: Calculating raster area by converting to points to account for latitude changes?

texts:
I am looking at a large geographic area (e.g., continental Africa). I generated a raster of modeled habitat (using MaxEnt) defined in WGS 1984 depicting the hypothetical distribution of an organism at 2.5 Arc-minutes. I am looking to calculate a number of seemingly simple area statistics (using ArcGIS 10.3) but I am having trouble accounting for changes in cell sizes.

For example, I have an expert delineated range polygon shapefile with defined presence zones for the species defined in World Cylindrical Equal Area. I want to know how many pixels (and how much area in km2) of my raster are found in each of the zones, but my raster spans continental Africa. I've run zonal statistics but direct pixel counts provide inaccurate estimates of area size. I've tried projecting my raster into World Cylindrical Equal Area but subsequent error messages have deduced this is not possible.

I've thought about converting the raster pixels to a pointfile, projecting the pointfile to World Cylindrical and then counting the number of pixels per zone as defined in the expert delineated shapefile. If the pointfile and range polygon are both projected into World Cylindrical Equal Area projection would I be able to then 'sum' the number of points found in each zone delineated in the expert range shapefile and then multiply the number of points by the cell-size of the original raster (2.5 Arc-minutes)? Or am I losing information here?

I've looked around on the forum and found similar questions (particularly this one) but I'm a bit lost on how to actually generate a reliable estimate of area spanning large changes in latitude.

I'm not entirely new to GIS but I am a bit rusty.

It seems like this would be a relatively straightforward procedure as most articles I read on MaxEnt outputs report pixel percents and area coverage, but studies spanning large geographic regions must be accommodating for the change in cell size somehow.

I receive the following errors using the 'Project Raster' tool from WGS 1984 to World Cylindrical Equal Area.

ERROR 999999: Error executing function.
FDO error: -2147217370 [binary_ProjectRaster]
No spatial reference exists.
The table was not found. [binary_ProjectRaster]
Failed to execute (Project Raster).
Failed at Sat Jun 18 15:18:20 2016 (Elapsed Time: 54.82 seconds)

Screenshots of layer projection and tool interface


texts: How to make a WPF Window be a well-behaved child window of ArcMap?

texts:
I am pretty familiar with Windows Forms in ArcGIS but WPF, which I'm just starting to tinker with, seems to be a different animal. I've created an ArcMap add-in with a button that displays a Window with some controls on it and would like to display it as a modal dialog and have it behave like other dialogs in ArcMap. That is, it should:

Initially display centered in the middle of the ArcMap window
Block the ArcMap main thread from continuing until the window is closed
Not allow the ArcMap application's window to take the focus
Not have its own task bar button or show up in the alt-tab list
Always draw on top of the ArcMap window

I've gotten 1 and 2 taken care of, but the rest I haven't figured out yet.

Anyone already tackled this?

Related: How stable is WPF in Arcmap?


texts: Do I still have to use a template DB for PostGIS?

texts:
createdb mydb -T template_postgis

Is this -T template_postgis still necessary in new versions of PostGIS?

I never saw this db being different, it always has those 3 sql files you run there. Is this still necessary?


texts: I am looking for a solution that addresses this problem, but I have not been successful in my searches. After vectoring points on a raster, they are minimally aligned, as shown in the appended illustration. Then you get an initial line, the blue line, and you want this line to be adjusted according to solution 1 or solution 2, ignoring points more distant. I do not know if I could do this with PyQGIS or PostGIS. Has anyone had to provide similar solution can you point me a north?


texts:
Can QGIS read layers in OGR's virtual vector format? It seems to work with VRT raster files, but it reports my layer is an invalid data source for vectors.

http://www.gdal.org/drv_vrt.html

I'm running QGIS 2.6 on Windows with GDAL 1.11.0 (which supports virtual vectors).

I've created a virtual layer as follows, which concatenates 4 shapefile layers.

The reason for wanting to do this is that the data is too big to fit into a single shapefile. I could use another format, but it would be more convenient if I could just use a virtual raster. I'm sure there a plenty of valid reasons for wanting to do this.

<OGRVRTDataSource>
    <OGRVRTUnionLayer name="unionLayer">
        <OGRVRTLayer name="OSMM_SPNE">
            <SrcDataSource>OSMM_SPNE.shp</SrcDataSource>
        </OGRVRTLayer>
        <OGRVRTLayer name="OSMM_SPNW">
            <SrcDataSource>OSMM_SPNW.shp</SrcDataSource>
        </OGRVRTLayer>
        <OGRVRTLayer name="OSMM_SPSE">
            <SrcDataSource>OSMM_SPSE.shp</SrcDataSource>
        </OGRVRTLayer>
        <OGRVRTLayer name="OSMM_SPSW">
            <SrcDataSource>OSMM_SPSW.shp</SrcDataSource>
        </OGRVRTLayer>
    </OGRVRTUnionLayer>
</OGRVRTDataSource>

The ogrinfo tool can read the layer correctly:

$ ogrinfo -al -so virtual.xml
INFO: Open of `virtual.xml'
      using driver `VRT' successful.

Layer name: unionLayer
Geometry: Polygon
Feature Count: 1469256
Extent: (398123.820000, 196969.360000) - (501938.600000, 260658.500000)
Layer SRS WKT:
PROJCS["British National Grid",
    GEOGCS["Ordnance Survey Great Brit",
        DATUM["Ordnance Survey Great Brit",
            SPHEROID["Airy 1830",6377563.396,299.3249646]],
        PRIMEM["Greenwich",0.0],
        UNIT["Degree",0.0174532925199433]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["False_Easting",400000.0],
    PARAMETER["False_Northing",-100000.0],
    PARAMETER["Central_Meridian",-2.0],
    PARAMETER["Scale_Factor",0.9996012717],
    PARAMETER["Latitude_Of_Origin",49.0],
    UNIT["Meter",1.0]]
TOID: String (20.0)
FEATCODE: Integer (9.0)
VERSION: Integer (9.0)
VERDATE: Date (10.0)
THEME: String (80.0)
CALCAREA: Real (19.11)
CHANGE: String (80.0)
DESCGROUP: String (150.0)
DESCTERM: String (150.0)
MAKE: String (20.0)
PHYSLEVEL: Integer (9.0)
PHYSPRES: String (20.0)
BROKEN: Integer (4.0)
LOADDATE: Date (10.0)
AREA: Real (19.11)
LEN: Real (19.11)
SYMBOLCODE: Integer (4.0)
SHAPE_Leng: Real (19.11)
SHAPE_Area: Real (19.11)


texts:
I am trying to plot a simple single band raster in R using ggplot. The only problem with my plot is that the colors are inverted. How can I fix that?
### For plotting I am using ggplot2 and ggpspatial packages
## Map 4.13
# Convert the Band 5 raster to a dataframe in string format
B5 = raster("~B5.TIF")

B5_df = as.data.frame(B5, xy = T)
names(B5_df)[3] = "Values"
# Values represents Band 5 column values in the data-frame
ggplot() +
  layer_spatial(points, aes(col = value)) +
  geom_raster(data = B5_df, aes(x = x, y = y, alpha = Values)) +
  annotation_scale(location = "bl") +
  annotation_north_arrow(location = "tr", which_north= "true") +
  labs(title = paste("Map 4.13")) +
  theme(plot.title = element_text(hjust = 0.5))

My plot

Actual Plot

The raster can be downloaded from here


texts:
I'm trying to mosaic several rasters with gdal.BuildVRT in Python 2.7 (Windows system), but code runs with logical problems.

When I try to merge 3 Landsat images datasets, resulting image is a tile (Float32) that include 'nan' values in raster areas of two datasets (see next image).

When I do a mosaic with GdalBuildVRT in Python, it shows me the following problem:

The code used is the following:

from osgeo import gdal
import osr, numpy as np

def save_raster ( output_name, dataset, driver ,NaN_Value):
    """
    A function to save a 1-band raster using GDAL to the file indicated
    by ``output_name``. It requires a GDAL-accesible dataset to collect
    the projection and geotransform.
    """
    # Open the reference dataset
    g = ( dataset )
    # Get the Geotransform vector
    geo_transform = g.GetGeoTransform ()
    x_size = g.RasterXSize # Raster xsize
    y_size = g.RasterYSize # Raster ysize
    srs = g.GetProjectionRef () # Projection
    raster_data = g.ReadAsArray()
    NaN_rast = NaN_Value
    # raster_data[raster_data == NaN_rast] = 'NaN'
    raster_data[raster_data == NaN_rast] = np.NaN
    # Need a driver object. By default, we use GeoTIFF
    driver = gdal.GetDriverByName ( driver )
    dataset_out = driver.Create ( output_name, x_size, y_size, 1, \
            gdal.GDT_Float32 )
    dataset_out.SetGeoTransform ( geo_transform )
    dataset_out.SetProjection ( srs )
    dataset_out.GetRasterBand ( 1 ).WriteArray ( \
            raster_data.astype(np.float32) )

TILE1 = "C:\Users\italo\Documents\GeoPRADS\LANDSAT\NUEVAS IMGS LANDSAT\Input ETa\TILE1.TIF"
TILE2 = "C:\Users\italo\Documents\GeoPRADS\LANDSAT\NUEVAS IMGS LANDSAT\Input ETa\TILE2.TIF"
TILE3 = "C:\Users\italo\Documents\GeoPRADS\LANDSAT\NUEVAS IMGS LANDSAT\Input ETa\TILE3.TIF"

in_subruta = TILE1+'-'+TILE2+'-'+TILE3 #string
subruta = in_subruta.split('-') #to list
gdal.BuildVRTOptions(VRTNodata= 0)
b1 = gdal.BuildVRT("dummy", subruta,VRTNodata= "NaN" )
save_raster(pathout+'/mosaic.tif',b1,"GTiff",0)

VRTNodata take the 'nan' data as a value.

So, where is the problem in the code?


texts: Can't properly map Lat/Long data on a googlemap

texts:
At the very beginning I should mention that I am a novice in spatial data handling. I am in it because I am recently doing a project where certain mining information needed to be embedded on a global map. So my question may be naive. Please bear with me.

Here is what I want to do. I have information about mines (from USGS site) which I want to embed on a global map. To do this, I am using R's RgoogleMap package along with googleVis package. I am using the function gvisGeoMap() from RgoogleMap.

The first few lines of my data:

     Deposite.or.District.Name       Country   Latitude          Longitude
4                   Adiounedj         Mali      20.1833333    0.066666666667
5                 Agate Mountain     Namibia     -18.45     -12.0833333333
6                   Aley             Canada      -56.45      23.7333333333

For gvisGeoMap(), I needed to add two extra columns: locationVar, and numVar.
locationvar of latitude and longitude data needs to be inserted in a special format. locationvar="Latitude:Longitude" and numvar is currently only 1s

But I am not getting the exact rendering when I run the following command in R:

g1 <- gvisGeoMap(cleaned_lat_long_df,
                 locationvar = "LatLong",
                 numvar = "numVar",
                 hovervar="Deposite.or.District.Name",
                 options= list(dataMode = "markers",
                               width = 400, height = 300))

plot(g1)

Here is what I am getting.

Can anybody help me with it ? As I said, I am a novice with these type of data and visualization.


texts: Set Default File Format on ESRI Print Widget

texts:
I'm using the Print widget in the ArcGIS JavaScript API (4.4) and would like to set the File format to default to PDF instead of PNG8.  I've read through the documentation and don't see a straightforward way of setting that parameter.

Does someone know of a way to do this?

I'd like to minimize the number of steps my users have to take in order to get a document they can print successfully - PDFs are the best option in my case.


texts: Mapping Sales Data With ArcMap's Timeslider?

texts:
I have some made-up data describing "Brenner" (a made-up brand) beer sales per bar, restaurant, supermarket (i.e. "Account"). The volume of beer is measured in cases. And, the time range is 4 weeks. I'm trying to create an animation, using the Timeslider tool in ArcMap 10.3, showing how the volume of Brenner beer cases changes per account week-to-week.

My inspiration is this: http://tipstrategies.com/geography-of-jobs/

So far, I've created a graduated symbol map for one week. This shows, basically, little dots for accounts with low sales and bigger dots for accounts with higher sales. What I want is to have the graduated symbols change over time - to reflect the sales per account week-to-week. From my trials, I can't seem to find a way to combine the time variable with another spatial-analysis tool (e.g. graduated symbols). It seems to be one or the other.

Does anyone have experience in this regard? Can you pass on some of your wisdom?


texts: Calculating the length of two shapefiles in QGIS

texts:
I created the yellow area as a shapefile by hand. Then I calculated the perimeter. I also created another shapefile (blue) manually and calculated the length, see image 1.
The next step was to subtract the perimeter of the yellow area with the length of the blue line and add the remaining green length (total green length - blue line) to calculate the length of the purple line, which was created manually, see image 2.
Does anyone know if there is another way to calculate the purple length? I tried using Vector > Intersection to automatically generate the yellow area. Unfortunately, this did not work (the layer was empty). After that, I used Vector > Union. But I have no access to the yellow area, see image 3.
I am at a loss at the moment.

I discarded the approach with Union and Intersection and instead applied split the line for both shapefiles. Then examined the attribute tables and created a query using a Virtual Layer.
SplitContourLine returns five values where only one is relevant (row 1), see image 4.
I solved the problem by taking out all 0 and greater 5000 values, see image 5.
Now I want to access this inner line without using this condition c.ContourLen < 5000 to be more flexible when using this method for other examples.


texts: How to merge Shapefiles from multiple folders into Geodatabase Featureclass?

texts:
I have a question regarding the merge of a lot of folders, which has a lot subfolders, which has the shapefiles stored. I would like to merge all those in one geodatabase.
So I thought to use a Python script.

import os, arcpy.da

print os.getcwd()

for dirname, dirnames, filenames in os.walk('.'):
    for subdirname in dirnames:
        current_dir =  os.path.join(dirname, subdirname)
        arcpy.env.workspace = current_dir
        fcList = arcpy.ListFeatureClasses(".shp")
        Merge = r"D:\Test\Merge_multiple.gdb\Totals"
        for fc in fcList:
            print fc
            arcpy.Merge_management([fc,destination],Merge)
            arcpy.Delete_management(desination)
            arcpy.Rename_management(Merge, destination)
        break

So, I started the script and after a minute it stops and unfortunately, no result. Because there are several different names in the shapefiles, I thought to use only the suffix .shp.

Could you help me?

[Edit]
I get the path of the working directory in the python window. After this message, the window closes without any errors.

[Edit2]
Based on PolyGeo's advice, I've run the script also in the IDLE environment and this is the result:


texts: Using the same mask on two rasters, but I get different extents in R

texts:
I loaded a shapefile in R, and cropped and masked it with the same mask, but I get different extents. Any suggestions?

Code below:

library(raster)
library(rgdal)
library(GISTools)
library(sp)
library(maptools)

##Loading the first mask
Mask <- readOGR("C:\\Users\\rameshv\\Downloads\\Climate Stability\\SPV_Mask\\mask.shp")
proj4string(Mask) <- CRS("+init=epsg:4326")

#Loading the Rasters For Present and LGM and clipping them to the right extent
bio5 <- raster("C:\\Users\\rameshv\\Downloads\\Climate Stability\\Data_LGM_Present\\Present\\1_RawData\\bio_5")
proj4string(bio5) <- CRS("+init=epsg:4326")
lg5 <- raster("C:\\Users\\rameshv\\Downloads\\Climate Stability\\Data_LGM_Present\\LGM\\1_RawData\\cclgmbi5.tif")
proj4string(lg5) <- CRS("+init=epsg:4326")

##Cropping by using the Crop and Mask functions
cr<- crop(bio5, extent(Mask))
bio5 <- mask(cr, Mask)

cr2<- crop(lg5, extent(Mask))
lg5<- mask(cr2, Mask)

> bio5
 class       : RasterLayer
 dimensions  : 339, 246, 83394  (nrow, ncol, ncell)
 resolution  : 0.04166667, 0.04166667  (x, y)
 extent      : 72.45835, 82.70835, 8.083337, 22.20834  (xmin, xmax, ymin, ymax)
 coord. ref. : +init=epsg:4326 +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
 data source : in memory
 names       : bio_5
 values      : 207, 432  (min, max)
 attributes  :
    ID COUNT
 from: -86     1
 to  : 489    38

 > lg5
 class       : RasterLayer
 dimensions  : 339, 246, 83394  (nrow, ncol, ncell)
 resolution  : 0.04166667, 0.04166667  (x, y)
 extent      : 72.45833, 82.70833, 8.083333, 22.20833  (xmin, xmax, ymin, ymax)
 coord. ref. : +init=epsg:4326 +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
 data source : in memory
 names       : cclgmbi5
 values      : 187, 392  (min, max)

EDIT

Following a projection, I am able to align extents of the rasters, however, if I coerce them to a vector, I am getting differences in number of elements. Any thoughts?

new.lg5 <- projectRaster(lg5,bio5)

x <- as(bio5,"SpatialPixelsDataFrame")
y <- as(new.lg5, "SpatialPixelsDataFrame")

> str(x)
 Formal class 'SpatialPixelsDataFrame' [package "sp"] with 7 slots
  ..@ data       :'data.frame': 54202 obs. of  1 variable:

> str(y)
  Formal class 'SpatialPixelsDataFrame' [package "sp"] with 7 slots
  ..@ data       :'data.frame': 54747 obs. of  1 variable:


texts: Alaska Series E projection in ArcGIS for Desktop?

texts:
I am trying to digitize a scanned map that appears to be in the Alaska Series E projection. I have also seen this referred to as Modified Transverse Mercator. While this projection is listed as supported by ArcGIS, I can't seem to find it when changing the coordinate system in Data Frame properties. Has anyone worked with the Alaska Series E projection and how do you get it to project properly in ArcMap? I am running ArcGIS 10.3 and I have an Advanced license.

ArcGIS entry for Alaska Series E:
http://help.arcgis.com/EN/arcgisdesktop/10.0/help/index.html#/Alaska_Series_E/003r0000001m000000/


texts: Adding point shapefile layer to polygon shapefile window in R?

texts:
I'm new to Spatstat for R and I'm trying to add a point layer to a window I've created from a polygon. I see questions like this have been asked before but I still can't get what I need from the responses given.  So far I have:

##Run libraries:
library(“spatstat”)
library(“maptools”)
library(“sp”)
library(“rgdal”)

##set the working directory
setwd(“Z:/CP/spatstat/data/”)

##read in the shapefiles
Points <- readOGR(getwd(), "MyPointShapefile")
StudyArea <- readOGR(getwd(), "MyPolygonShapefile")

##create the window from polygon shp
w <- as.owin(StudyArea)

How do I add the two together so the points are projected in the window?


texts: Arcpy.da.UpdateCursor, Error - Sequence size must match size of row

texts:
I am trying to use an UpdateCursor fill in a field but I keep getting an error, can anyone tell me what I am doing wrong? My code is below:

import arcpy

fc= "C:\Users\Keri\Desktop\GEP_690\Project Data\FINAL_PROJECT.gdb\NYC_Demographics_2010"
with arcpy.da.UpdateCursor(fc, ["Total1","White","Pct_White"]) as cursor:
    for field in cursor:
        if field[0]==0:
            field[2]=0
        else:
            field[2] = field[1] / field[0]
        cursor.updateRow([field])
print "Complete."

del row

del cursor

This is the error I keep getting:

Traceback (most recent call last):
  File "C:\Python27\ArcGIS10.2\Lib\site-packages\pythonwin\pywin\framework\scriptutils.py", line 326, in RunScript
    exec codeObject in __main__.__dict__
  File "C:\GEP_662\Lab9.py", line 15, in <module>
    cursor.updateRow([field])
TypeError: sequence size must match size of the row


texts: Create a mosaic from images with no spatial reference system

texts: I have 10 raster images (.tif format) with no spatial reference system (SRS). Is it possible to mosaic them by taking into account that they do not have a SRS?


texts: Differences between ST_ClosestPoint and ST_LineInterpolatePoint?

texts:
As a bit of context : my exact goal is taking the set of points (some GPS data) and correcting their position so that they intersect some lines (OSM roads data), to have them nicely sitting on the roads. I will then use those corrected positions to compute several things.

While working on projecting a set of points onto a set of lines, I have stumbled upon two functions : ST_ClosestPoint and ST_LineInterpolatePoint (used with ST_LineLocatePoint).

Which I used this way :

Create table projected_points as (
    Select
        p.id,
        p.geom as original_pt,
        ST_ClosestPoint(l.geom, p.geom) as closest_pt,
        ST_LineInterpolatePoint(l.geom,
                                ST_LineLocatePoint(l.geom, p.geom)
        ) as interpolated_pt,
    From points p, lines l
);

I've tried out both functions, and found that both functions almost never return points that actually intersect the lines. Which in itself made me wonder what was going on, so I used an ST_Distance to try and understand what went wrong.

Using the following query :

 Select
        ST_Distance(p.closest_pt, l.geom),
        ST_Distance(p.interpolated_pt, l.geom),
    From projected_points p, line l

On the left is what I get with ST_ClosestPoint. On the right, what I get using ST_LineInterpolatePoint. In both cases, having a distance of "0" doesn't necessarily mean that the point intersects the line. I also notice that the distances that were computed were all the same digits.

But it seems like ST_LineInterpolatePoint works better. On the downside, it requires me to use ST_LineLocatePoint first, which makes my code very heavy.

What makes the two functions different, yet so similar ?


texts: How to create buffer zone in metres for OSGB1936 layer?

texts:
I am creating buffer zones in metres on top of UK Ordnance Survey map files that are in OSGB1936 which uses degrees.

How do I draw the buffers in meters?


texts: Upload raster to PostGIS problem: ERROR: unterminated quoted string at or near

texts:
I used this command:
raster2pgsql -s 4326 -I -C -M T31TEJ_20210216T104019_B01.tif -F -t 100x100 public.testraster | PGPASSWORD=1122 psql -d ter -h localhost -U dtn -p 5432
On MacOS Bigsur 11.2.3 it works very good. But, when I deploy to Ubuntu 20.04 It gives the error: ERROR:  unterminated quoted string at or near "'010000 and free(): invalid pointer.
Can you help me?
I configured the PostGIS like on Mac.

Here is the extensions which was installed.


texts:
I have downloaded the following OSM bbox file: https://overpass-api.de/api/map?bbox=-105.32,39.95,-105.2,40.08
And want to use it to find all points of interests. My teacher have given me the following code in Python but I cannot make it work:
import pandas as pd
import gzip
import xml.etree.ElementTree as ET

## Functions
def pick_from_dict(d, keys):
    return {key: d[key] for key in keys}
def merge_dicts(d1, d2):
    return dict(set(d1.items() | set(d2.items())))
def get_amenity(node):
    return {'amenity': next(filter(lambda child: child.attrib.get('k') == 'amenity', node)).attrib['v']}
def has_amenity_child(node):
    return any(child.attrib.get('k') == 'amenity' for child in node)
def get_pois(root):
    return pd.DataFrame(
        map(lambda node: merge_dicts(get_amenity(node), pick_from_dict(node.attrib, ['lat', 'lon', 'timestamp'])),
            filter(has_amenity_child, root.findall('node')))
    ).assign(timestamp=lambda d: pd.to_datetime(d.timestamp))\
    .sort_values(by='amenity')\
    .assign(lat=lambda df: df.lat.astype(float), lon=lambda df: df.lon.astype(float))

## Load data
df_pois = get_pois(ET.parse(gzip.open('map.xml.gz', 'rb')).getroot())
df_pois.head()

From the code I can see that he is using an xml.gz file. However when I download the bbox file the file format is imply "File". How can I convert the file and do I need to?


texts: How to make a model with more than 2 input layers for buffering and merging?

texts:
Would it be possible to add more inputs to Buffer and Merge? From what I can see there is a limit of two input layers per both. As far as I know this is possible in ArcGIS's Model Builder but from the looks of it, I can't find it.

My goal is to make a model for pre-processing 5 shapefiles, the goal is to output a single one.

What are my options in this case?


texts: Setting Custom Cluster Icon with Picture Marker Symbols using ArcGIS API for JavaScript

texts:
Problem:
Unsure how to set a custom icon for the cluster symbol, not the individual points or even the label of the cluster but the symbology for the cluster icon itself.
Context:
A Feature Layer is being rendered by unique value with picture marker symbols for each value. When clustering is applied it picks one of the PictureMarkerSymbols to use as the symbology for the cluster point itself. I would like to have a completely different symbol for the cluster itself.
Is there a way to customise the cluster icon?
TLDR: Esri Leaflet seems to accomplish exactly what I am describing (see link and gif):
https://esri.github.io/esri-leaflet/examples/styling-clusters.html
Is the cluster symbol being different from the feature symbol possible to do in the ESRI JS API?


texts: Merge multiple polygons in separate clusters all in the same layer

texts: I need to merge each cluster of adjacent parcels into one parcel.  No attributes is needed as this will be used in a spatial join. Please refer to the image as an additional aid to my explanation.