texts:
I am trying to compare the NDVI-values of different land cover classes over a certain time period and display it in a chart.
So in a first step, I'd like to extract the crop land cover and  use its geometry to calculate just the NDVI of the croplands
This is my code so far.

var area = ee.FeatureCollection('users/jonathanreith/KK')
Map.centerObject(area,9);
var lc  = ee.Image("ESA/GLOBCOVER_L4_200901_200912_V2_3").select('landcover');

var crop = lc.eq(30)

// Convert the crop landcover to vectors.
var crop = crop.reduceToVectors({
  geometry: area,
  geometryType: 'polygon',
  eightConnected: false,
  labelProperty: 'landcover',
  scale: 300
});

print(crop)
var cropu = crop.union(1)
print(cropu)

//chart
var NDVITimeSeries = ui.Chart.image.seriesByRegion({
    imageCollection: l8,
    regions: cropu,
    reducer: ee.Reducer.mean()
});

Unfortunately, crop consists of thousands of features, while union somehow erases all the information and just leaves the outer bounding box.

So how do I create an union feature with its geometries intact?
My Code
I already checked these Q&As:

Dissolving featurecollection in google earth engine?; and
Polygons vs MultiPolygons in Google Earth Engine


texts: I am using the corner label feature in Arcmap 10.  This allows you to set the "Numer of principal digits" and the "Number of base digits".  What do these properties refer to?


texts: Error in creating changes map in MOLUSCE

texts:
I'm using QGIS 2.8.9 Wein and the plug-in MOLUSCE(ver. 3. 0. 13).
Now I'm trying to create a change map, but get the following error.

Traceback (most recent call last):
  File "C:/Users/.../.qgis2/python/plugins\molusce\moluscedialog.py", line 479, in changeMapDone
    self.inputs["changeMap"].save(self.inputs["changeMapName"])
AttributeError: 'NoneType' object has no attribute 'save'

Python version:
2.7.5 (default, May 15 2013, 22:44:16) [MSC v.1500 64 bit (AMD64)]

I'm new to QGIS and not familiar with Python.

So I wonder that "NoneType" something means that the preprocessing was wrong... but not sure.

EDIT: Updated to QGIS 2.18, but the same issue persists.


texts:
Does anyone know where I can find a data model for land use?

I know about data models available on Esri Support site.

Are there any other resources?


texts:
I have a county polygon. In it, I have a field (column) that is either yes or no. There are 100 rows to the attribute table and 75% of the attributes are yes, 25% are No.
Is there a way I can symbolise this in ArcGIS Online? I was thinking of a pie chart, where the user would click on the county polygon and it would show 3/4 Yes (Green) and 1/4 No (Red). If it would be possible to click on that chart to bring up a pop-up, that would be even better.


texts: Connecting Oracle Database 11g in ArcMap?

texts:
I installed Oracle Database Enterprise Edition 11g(11.2.0.1) in my machine. And I installed ArcMap Desktop 10.3.1 on another machine. Whenever I am trying to create a direct connection of Oracle DB from ArcMap without ArcSDE installation, I am getting the following error like "Failed to connect to the specified server. Do you want to continue? Bad login user". Even I have enabled the port in Oracle DB machine. I was giving the oracle server is the DB IPAddress, Instance Name is "orcl" and port is 1521.

Do I need to install Oracle Client software where the Oracle DB is existing and, if not, where I am lacking in this concept?


texts: Error in batch process converting raster from .img to .tif in QGIS 3.4.8

texts:
I'm a long time GIS user but new to QGIS. I can convert single raster files in .img format to .tif files in QGIS 3.4.8 without any issues. However, I cannot batch process the same files. The Python error message is:

2019-05-29T16:58:13     WARNING    Traceback (most recent call last):
              File "C:/PROGRA~1/QGIS3~1.4/apps/qgis-ltr/./python/plugins\processing\gui\BatchAlgorithmDialog.py", line 88, in runAlgorithm
              parameters[param.name()] = wrapper.parameterValue()
             AttributeError: 'NoneType' object has no attribute 'parameterValue'

The parameters in the batch file window are the same as in the single file processing window. Can anybody point me to a solution?

The log in the processing window is empty. The general messages window contains the message below which only points me to the Python Error tab.

2019-05-30T06:57:52     WARNING    Python error : An error has occurred while executing Python code: See message log (Python Error) for more details.

When processing a single file (which works fine) the Log window displays this information:

Processing algorithm…
Algorithm 'Translate (convert format)' starting…
Input parameters:
{ 'COPY_SUBDATASETS' : False, 'DATA_TYPE' : 0, 'INPUT' : 'G:/Snow Damage Aerials/LRT Test/orthoc190426a_0493.img', 'NODATA' : None, 'OPTIONS' : '', 'OUTPUT' : 'G:/Snow Damage Aerials/LRT Test/orthoc190426a_0493.tif', 'TARGET_CRS' : None }

GDAL command:
gdal_translate -of GTiff "G:/Snow Damage Aerials/LRT Test/orthoc190426a_0493.img" "G:/Snow Damage Aerials/LRT Test/orthoc190426a_0493.tif"
GDAL command output:
Input file size is 4208, 5931

0...10...20...30...40...50...60...70...80...90...100 - done.

Execution completed in 9.20 seconds
Results:
{'OUTPUT': 'G:/Snow Damage Aerials/LRT Test/orthoc190426a_0493.tif'}

Loading resulting layers
Algorithm 'Translate (convert format)' finished


texts: I am using ArcGis 9.3 and I need to split a large number of circles in quadrants. The circles corresponde to a 200m buffer of a point layer and, for each point, I need to split the buffer in 8 quadrants. Is there any ArcGis tool or extension that does this for all points/buffers at the same time?


texts:
I try to add labels of a vector I added in openlayers with a stylefunction. The function works for stroke, but no text is displayed and I can't find a solution for it.

My code is the following:

var myfeature = new ol.layer.Vector({
  title: '&nbsp; myfeature',
  source: new ol.source.Vector({
      url: 'data/myfeature.geojson',
      format: new ol.format.GeoJSON()
  }),
  style: styleFunction
});

function styleFunction() {
  return [
  new ol.style.Style({
    stroke: new ol.style.Stroke({
      color: [0, 0, 0, 1.0],
      width: 1,
      lineDash: [1, 5, 3, 5]
    }),
    text: new ol.style.Text({
      font: '50px Calibri'
      }),
      text: myfeature.get('label')
    })
  ];
}

The fault is probably in text: myfeature.get('label'), but I can't say what I do wrong nor worked anything else I tried.


texts: How to create varying opacity within a layer in Mapserver?

texts:
Suppose you have a mapfile containing a background layer and a layer with a lot of overlapping linestrings.

Now, for the overlapping linestrings, it might be good to visualize their density, thus at places where a lot of lines do overlap they might be drawn more "solid" than at places where only a few lines are present.

Using QGIS, for example, I can do this by just turning the layer opacity to some small value (say five percent). In Mapserver, however, setting the opacity of the layer, the opacity will just make the whole layer see through, but will not add up for each single element.

How can I archive the desired effect?


texts: Creating spatial layer from MS SQL Server table?

texts:
I have a table in MSSQL that I am trying to use to create a spatial layer from using the QgsVectorLayer function to load the table as a layer, and then using Fiona to create a shapefile but I am having problems in working out whether I am going about this in the correct way!

My code is below - I am able to connect to my MSSQL table so have it loaded as a QgsVectorLayer, but I am struggling with instructing Fiona / Shapely to convert the table to a spatial layer by using the XREF and YREF fields from my table as the geometry.

Looking at my code, am I along the right lines here? It does not produce any file so I obviously need to adapt this!

import sys
from qgis.core import *
from PyQt4.QtCore import QVariant

from osgeo import ogr
import os

from shapely import wkb, geometry
import fiona
import ogrtools

# Initialize QGIS Application
QgsApplication.setPrefixPath("C:\\Program Files\\QGIS Pisa\\apps\\qgis", True)
app = QgsApplication([], True)
QgsApplication.initQgis()
print 'QGIS Application initialized'

# Add the path to Processing framework
sys.path.append('C:\\Program Files\\QGIS Pisa\\apps\\qgis\\python\\plugins')

# Import and initialize Processing framework
from processing.core.Processing import Processing
Processing.initialize()
import processing
print 'QGIS processing framework initialized'

driverName = "ESRI Shapefile"
driver = ogr.GetDriverByName(driverName)
GazProp_shp = "X:\\llpg_test.shp"
GGP_LLPG_shp = "X:\\GGP_LLPG.shp"
GGP_LLPG_test = "X:\\ggp_llpg_test.shp"

uri = QgsDataSourceURI()
# set host name, port, database name, username and password
uri.setConnection("myserver", "", "GGPNLPGRDC", "GGPAccessRDC", "GGPAccessRDC")
# set database schema, table name, geometry column and optionally
# subset (WHERE clause)
uri.setDataSource("dbo", "ADDRESSTABLE","")

vlayer = QgsVectorLayer(uri.uri(), "ADDRESSTABLE", "mssql")
if not vlayer.isValid():
    print "Layer failed to load!"
else:
    print "llpg_tan layer loaded successfully"

with fiona.open(vlayer, 'r') as input:
    input.schema['geometry'] = "Point"
    with fiona.open(GGP_LLPG_test, 'w', driverName, input.schema.copy(), 'EPSG:27700') as output:
        for elem in input:
            #geometry
            points=Point(input['XREF'], input['YREF'])
            output.write({'properties': elem['properties'],'geometry': mapping(points)})
        print 'Job done!'


texts: Applying Symbology to group layer using ArcPy?

texts:
I currently have around 200 features that all need have the same symbology applied to them. I have used the ApplySymbologyfromlayer_management however this takes much to long to process, at least a few hours and the labels arent brought across either.

import arcpy
mxd = arcpy.mapping.MapDocument("CURRENT")
inputlayers = ['ST13003-00108100-009','ST13003-00116584-073','ST13003-00115372-080','ST13003-00112948-011','ST13003-00111736-021','ST13003-00122544-083','ST13003-00099716-018']
symbologylayer = 'Layer'

for layer in inputlayers:
    arcpy.ApplySymbologyFromLayer_management(layer, symbologylayer)

This works however justs takes way to long.

What I want to do export all lines from ArcMap as a layer or layer package and apply the symbology across all of the features without doing this manually.

When I try this however this error occurs.

  Runtime error Trace back (most recent call last):   File "",
  line 1, in    File "c:\program files
  (x86)\arcgis\desktop10.1\arcpy\arcpy\management.py", line 5696, in
  ApplySymbologyFromLayer     raise e ExecuteError: ERROR 000967: The
  input layer type must be a feature layer, raster layer, or tin layer.

I looked at creating Representations, is that the solution or is there another way around this?


texts: How does a GPS unit ensure it has the same time reference as a GPS satellite?

texts:
I've been reading about GPS to find how it works, but one thing still bothers me.

To calculate distance from GPS unit/receiver to GPS satellite, it uses the difference between radio signal sending time and radio signal receiving time, which is broadcast by the satellites.

To do that calculation, of course, the GPS unit time has to be in-sync/the same with the GPS satellite time.

I read that GPS units first make a (false) distance calculation (or X,Y,Z coordinate) to three satellites, and then make the adjustment from a fourth satellite to make all the distance radii intersect at a single point. How exactly does a GPS unit do this?

I read that it achieves this by the means of pseudo-random code; how?


texts: How to sum up values of neighbouring rasters?

texts: I have a raster (25X25 meter) and I will make a new raster where each cell get the value of the sum of all the raster cells in a radius of 500m. I use QGIS. What is the best option to do this?


texts:
I am trying to plot some color code (only 3 colors by the moment) cars in my map instead the regular color code points. I am using torque cat to do it, but all the time the car is in black!, I need some help. This is my code:

/** torque_cat visualization */

Map {
  -torque-frame-count:64;
  -torque-animation-duration:1;
  -torque-time-attribute:"date";
  -torque-aggregation-function:"CDB_Math_Mode(torque_category)";
  -torque-resolution:2;
  -torque-data-aggregation:linear;
}

#untitled_table_1{
  comp-op: lighter;
  marker-fill-opacity: 0.9;
  marker-line-color: #FFF;
  marker-line-width: 0;
  marker-line-opacity: 1;
  marker-type: ellipse;
  marker-width: 6;
  marker-fill: #FF9900;
  marker-file: url(http://www.clipartbest.com/cliparts/nTX/AGX/nTXAGXyTB.svg);
}

#untitled_table_1[value=1] {
  marker-fill: #B81609;
}
#untitled_table_1[value=2] {
 marker-fill: #012700;
}
#untitled_table_1[value=3] {
  marker-fill: #7B00B4;
}


texts: Draw specific polygon and select points within

texts:
My goal is to select the points that are inside a rectangle.

How do I create a rectangle with 100 meter width and then only select and extract the points that fall inside?

I'm using ArcMap 10.5.


texts: Leaflet Routing Machine calculate routes in loop

texts:
I have an array with coordinates ([[54.73, 55.74], [55.66, 54.73]...]). I want to get route from each point to another. E.g. if I have 3 points a1, a2, a3, I want to calculate routes a1a2, a1a3, a2a1, a2a3 etc. I need this to make distance matrix later. Here is my loop where coords is array of coordinates

function routingLoop() {
            for (var i = 0; i < coords.length; i++){
                for (var j = 0; j < coords.length; j++){
                    var routing;
                    routing = L.routing.control({
                        waypoints:[
                            L.latLng(coords[i][0], coords[i][1]),
                            L.latLng(coords[j][0], coords[j][1]),
                        ],
                    }).addTo(Map);
                }
            }
        }

But I get this error Routing error:
Object { status: -3, message: "TypeError: this._altContainer is undefined" }. Any ideas what's wrong with that code?


texts: Assigning polygon layer attributes to point layer features in QGIS

texts:
I have a task which I’m sure is quite common, but being new to QGIS I’m hoping someone can point me in the right direction.  I have a point layer with PLANT_ID and SPECIES of individual plants, and three polygon layers: SOIL, DRAINAGE, and CANOPY.  For each individual plant, I would like to add three attributes, one for each of the three polygon layers.  The result would be something like this:

So for this example, Plant 1 was Species A, and it was located in a SOIL polygon with attribute Soil 1, and was located in a DRAINAGE polygon with attribute Drainage 4, and was located in a CANOPY polygon with attribute Canopy 4.  I suspect that this has something to do with a Spatial Join in MMQGS, but I've played around with that and was not able to get what is shown in the table above.

Beyond this, I was planning to apply some data mining techniques to create a decision tree like this:

If SOIL 1 and DRAINAGE 2 and CANOPY 3, then SPECIES D

Might there be a QGIS function that would create this type of tree?


texts: Tips for improving GeoServer WPS Performance

texts:
I have the following filter:

function getFilter(polygon) {
var filter =
'<ogc:Filter>' +
'<ogc:Intersects>' +
    '<ogc:PropertyName>the_geom</ogc:PropertyName>' +
    '<ogc:Function name="collectGeometries">' +
     '<ogc:Function name="queryCollection">' +
        '<ogc:Literal>ws:layer_2</ogc:Literal>' +
        '<ogc:Literal>the_geom</ogc:Literal>' +
        '<ogc:Literal>INTERSECTS(the_geom,' + polygon + ')</ogc:Literal>' +
      '</ogc:Function>' +
    '</ogc:Function>' +
'</ogc:Intersects>' +
'</ogc:Filter>';

return filter;

}

I've incorporated this filter in a gs:Aggreate WPSExecute request. However, processing takes a very long time, sometimes even reaching a minute. All operations are on my localhost.

Tips on how to improve performance of WPS requests will be highly appreciated. I'm using GeoServer 2.1.3, PostGIS 1.5, and OpenLayers 2.11.


texts: Elevation distribution of protected areas using QGIS

texts:
I am using QGIS.

Based on a DEM and a vector layer of polygons for protected areas in the Alps, I would like to calculate the areas of the respective protected area polygons in certain altitude classes (0-500m; 500-1000m;1000-1500m; etc.).

Is there a direct way of realizing such a calculation or do I need to convert the DEM to shapefiles according to the altitude classes and then clip/intereect those polygons (as described here https://gis.stackexchange.com/a/202435/152823)?

So far I have tried to create contours and now I want to try to clip this shp with an overlay with my protected area polygons.


texts: Is there a way to change the QGIS default CRS?

texts:
I'm working in a UTM projection and bringing in dxf layers from a different program.  The dxf files are all in the correct UTM coordinates already, but when I bring them in I need to change the CRS to UTM 17N instead of WGS 84.  Is there a way to change this default CRS?

Working in QGIS 2.6 Brighton


texts: Plugin for average distance from node to polyline in QGIS?

texts:
Is there a plugin available that will calculate this?

EDIT 1:

The red polyline is the average distance of the point to the black polyline. The red polyline is what I want to display. Alternatively, is there something which would calculate the radius of the circle so I can use a buffer?


texts: Openlayes - Loading 3 different KMLs via proxy; one appears and the others don't

texts:
I'm loading 3 KMLs using an ashx handler taken from http://code.google.com/p/iisproxy. The 3 KMLs I'm trying to load are Openlayers sundials ( http://openlayers.org/dev/examples/kml/sundials.kml ), Marine Traffic KML ( http://www.marinetraffic.com/ais/ge_marinetraffic.kml ), and flightwise ( http://mw1.google.com/mw-earth-vectordb/flightwise/flightwise_root.kml ).

The Sundial is working for me, but the other 2 don't.
I'm afraid the problem is because the other 2 KMLs have file links inside while sundials KML has data only.

Is there a way I can make the other two work on my IIS?

This is the code I used:

<!-- language: lang-js -->

function init() {

    var map = new OpenLayers.Map("map");
    OpenLayers.ProxyHost = "/proxy.ashx?url=";

    var osm = new OpenLayers.Layer.OSM();

    var flight = new OpenLayers.Layer.Vector("Flightwise", {
        strategies: [new OpenLayers.Strategy.Fixed()],
        protocol: new OpenLayers.Protocol.HTTP({
            url: "http://mw1.google.com/mw-earth-vectordb/flightwise/flightwise_root.kml",
            format: new OpenLayers.Format.KML({
                extractStyles: true,
                extractAttributes: true
            })
        })
    })

    var marine = new OpenLayers.Layer.Vector("Marine", {
        strategies: [new OpenLayers.Strategy.Fixed()],
        protocol: new OpenLayers.Protocol.HTTP({
            url: "http://www.marinetraffic.com/ais/ge_marinetraffic.kml",
            format: new OpenLayers.Format.KML({
                extractStyles: true,
                extractAttributes: true
            })
        })
    })

    var sundials = new OpenLayers.Layer.Vector("Sundials", {
        strategies: [new OpenLayers.Strategy.Fixed()],
        protocol: new OpenLayers.Protocol.HTTP({
            url: "http://openlayers.org/dev/examples/kml/sundials.kml",
            format: new OpenLayers.Format.KML({
                extractStyles: true,
                extractAttributes: true
            })
        })
    })

    map.addLayers([osm, flight, marine, sundials]);

    var controls = [
                    new OpenLayers.Control.Navigation(),
                    new OpenLayers.Control.LayerSwitcher(),
                    new OpenLayers.Control.ScaleLine(),
                    new OpenLayers.Control.MousePosition(),
                    /*new OpenLayers.Control.OverviewMap(),*/
                    new OpenLayers.Control.KeyboardDefaults()
                ];

    map.addControls(controls);

    map.zoomToMaxExtent();
}


texts: Looking for mobile web app map control that can handle layer redraw with good performance

texts:
I would like to develop mobile web app to track GPS coordinates in real time (refresh every 10 seconds) from Geoserver.

Currently I have desktop site I wrote in ExtJS 4 and it works pretty well (Openlayers + GeoJson from WFS) - the layer redraws are not noticeable and it gives a feeling that the points are moving (because there are no flickering when layer redraw).

I've tried doing the same using sencha touch 2 (found leaflet and Openlayers plugins on Github) . here is what I achieved thus far:

Tried Google map v3 API - the map renders fast but when drawing the geojson layer or kml layer , in order to update I need to delete the previous features and redraw which cause all the markers to flicker.
Openlayers - The map renders very slow when panning / zooming.
Leaflet - Also , the map renders very slow (and sometimes idle...) when trying to pan fast or zoom multiple times.

At this point I'm stuck with a nice UI but poor performance map control (or maybe I'm not using the Google API correctly)

This is how I implemented it using Google API:

var gMap = Ext.getCmp('map').getMap();
georssLayer.setMap(null);
georssLayer = new google.maps.KmlLayer({
   url: 'http://apps.utilis.co.il/geoserver2/wms/reflect?layers=ClearcutCB:DeviceInfo& format=rss&rnd=' + Math.random()});
georssLayer.setMap(gMap);

Will be more than happy having your advise on this one.


texts: Using dynamic label to give all results in capital letters?

texts:
I need this dynamic label to give me all the results in capital letters:

[dyn type="page" property="attribute" field="Predio" domainlookup="true"]


texts: Checking if vector coordinates are inside raster file using pyqgis?

texts:
I have some tifs files inside a folder and want to check if a shape coordinates are inside one of the tifs with pyqgis.

I tried this from Determining coordinates of corners of raster layer using PyQGIS? to get the raster corners so I can check if my shape is inside it :

from osgeo import gdal, osr

bag = gdal.Open('LC82220762017084LGN00_B2.TIF')  # replace it with your file # raster is projected
bag_gtrn = bag.GetGeoTransform()
bag_proj = bag.GetProjectionRef()
bag_srs = osr.SpatialReference(bag_proj)
geo_srs =bag_srs.CloneGeogCS()                 # new srs obj to go
from x,y -> φ,λ
transform = osr.CoordinateTransformation( bag_srs, geo_srs)

bag_bbox_cells = (
(0., 0.),
(0, bag.RasterYSize),
(bag.RasterXSize, bag.RasterYSize),
(bag.RasterXSize, 0),
)

geo_pts = []
for x, y in bag_bbox_cells:
    x2 = bag_gtrn[0] + bag_gtrn[1] * x + bag_gtrn[2] * y
    y2 = bag_gtrn[3] + bag_gtrn[4] * x + bag_gtrn[5] * y
    geo_pt = transform.TransformPoint(x2, y2)[:2]
    geo_pts.append(geo_pt)
    print x, y, '->', x2, y2, '->', geo_pt

But it returned

SyntaxError: Non-ASCII character '\xcf' in file intersect.py on line 8, but no encoding declared; see http://python.org/dev/peps/pep-0263/ for details

And tried this Determining if shapefile and raster overlap in Python using OGR/GDAL? to check if my shape overlaps the raster:

import ogr, gdal

raster = gdal.Open('LC82220762017084LGN00_B2.TIF')
vector = ogr.Open('Faz_Esperanca.shp')

# Get raster geometry
transform = raster.GetGeoTransform()
pixelWidth = transform[1]
pixelHeight = transform[5]
cols = raster.RasterXSize
rows = raster.RasterYSize

xLeft = transform[0]
yTop = transform[3]
xRight = xLeft+cols*pixelWidth
yBottom = yTop-rows*pixelHeight

ring = ogr.Geometry(ogr.wkbLinearRing)
ring.AddPoint(xLeft, yTop)
ring.AddPoint(xLeft, yBottom)
ring.AddPoint(xRight, yTop)
ring.AddPoint(xRight, yBottom)
ring.AddPoint(xLeft, yTop)
rasterGeometry = ogr.Geometry(ogr.wkbPolygon)
rasterGeometry.AddGeometry(ring)

layer = vector.GetLayer()
feature = layer.GetFeature(0)
vectorGeometry = feature.GetGeometryRef()

print rasterGeometry.Intersect(vectorGeometry)

Running this with my rasters and vector I get False returning every time even if the vector overlaps the raster or not.

Are there any function on Qgis that I not used?


texts: Error when saving geodataframe to geojson after a spatial join

texts:
I am writing a simple python code to perform sjoin on two geodataframe.
    gdf1 = gpd.read_file(polygon1.geojson)
    gdf2 = gpd.read_file(polygon2.geojson)
    gdf = gpd.sjoin(gdf1, gdf2, how="left", op='intersects')
    gdf.to_file(filename= 'myfile.geojson', driver='GeoJSON')

However, I encounter this error when writing the result to geojson.
I can avoid this error by writing to a shapefile then convert the shapefile to geojson. But this sounds too odd.
Is there anyone know the root cause of such an error?

C:\ProgramData\Anaconda3\lib\site-packages\geopandas\geodataframe.py
in to_file(self, filename, driver, schema, index, **kwargs)
531         from geopandas.io.file import to_file
532
--> 533         to_file(self, filename, driver, schema, index, **kwargs)
534
535     def to_crs(self, crs=None, epsg=None, inplace=False):
C:\ProgramData\Anaconda3\lib\site-packages\geopandas\io\file.py in
to_file(df, filename, driver, schema, index, **kwargs)
191             filename, "w", driver=driver, crs_wkt=crs_wkt, schema=schema, **kwargs
192         ) as colxn:
--> 193             colxn.writerecords(df.iterfeatures())
194
195
C:\ProgramData\Anaconda3\lib\site-packages\fiona\collection.py in
writerecords(self, records)
347         if self.mode not in ('a', 'w'):
348             raise IOError("collection not open for writing")
--> 349         self.session.writerecs(records, self)
350         self._len = self.session.get_length()
351         self._bounds = self.session.get_extent()
fiona/ogrext.pyx in fiona.ogrext.WritingSession.writerecs()
fiona/ogrext.pyx in fiona.ogrext.OGRFeatureBuilder.build()


texts: Add new field with the name based on another field

texts:
I would like to know how can I add a new field in an attribute table with a name based on the value of another field.  I mean, if I had a field with values 1,2,3 and 4, I need to add four new fields like A1, A2, A3 and A4 (A is aleatory, only for giving a name to new field).  I did it in ModelBuilder but I would like to do the same in a script for python.

Can anyone help me with this task?

I'm using ArcGIS 10.2.1 and Python 2.7.5


texts: What resampling method is used in the raster 'Save as...' dialog?

texts:
I received a raster that was resampled using the Save as dialog in QGIS. Now I need to describe the resampling process in a report. But I could not figure out the method used in the resampling. Is it Nearest Neighbor? Bicubic?

This question has been asked before (Resampling method used in QGIS Raster save as?) but the answer does not specify the method used.

I also tried looking at the QGis source code in GitHub, but unfortunately, I could not find where the magic takes place (I'm not a programmer)


texts: Unable to save point features to QGIS on MySQL database

texts:
I created a table with point geometry in MySQL. Subsequently added a vector layer from MySQL table. The table could be accessed from QGIS. However on pressing the add button the system threw warning message that no CRS was defined hence WGS84 has been used. I then changed the CRS to EPSG24378 as was desired. Now I select this vector layer and toggle to editing mode to add points. The MySQL table is shown when I click at the desired location. However when I attempt to save then QGIS trows error as
Could not commit changes to layer MySQL:xxxxxxxxxxx
OGR error creating feature-2:
MySQL error message: SRID value is out of range in 'st_geomfromtext'
Description: INSERT INTO tbl_loc(loc_grid, Name) VALUES (ST_Geomfromtext('POINT( 3383449.16856926 973171.230534863)',-2), Garud)

Obviously SRID of -2 is an error. How do I correct this?
QGIS VER 3.4.5
MySQL VER 8.0.3


texts: Google Earth Engine: How to calculate regular day interval means in image collections, maintaining "system:time.start" field?

texts:
I've merged three Landsat collections in Google Earth Engine, now I need to calculate median values at regular time intervals, let's say 32 days, because I will create a collection to join with MODIS. However, when I do that, I lose the "system:time.start" information from Landsat Merged collections. Does anyone knows how to calculate this regular medians based on DOY? I want to filter the Landsat_merge collection to 32 days regular interval!

This is the code so far:
Code to access MODIS Landsat Joined collections


texts: Has anyone played around with the Feature Service URL in Survey 123 Connect? There appears to be a feature service URL section under Data in the form Settings. I don't see any documentation on this function. Has anyone tried using it?


texts: Installing OpenStreetMap Tile Server with Custom Projection Systems?

texts:
I am trying to install the OSM tile server on a local machine with a custom projection system like epsg:3573. For the server installation I followed this tutorial (http://switch2osm.org/serving-tiles/building-a-tile-server-from-packages/) and it works fine with the default settings. For making changes, I played around with osm.xml and changed the srs information in there but still no luck.

Can anyone give me some hints on customizing the OSM tile server with a custom projection system?


texts: Black and white aerial image classification with ecognition fails to execute

texts: I am trying to classify an image in ecognition, it gives a good number of well shaped polygons, but when I try to execute the classification, it only shows one class; I've done the same process for color images, that's why I believe the error has to do with the B&W aerial photo.


texts: GeoWebCache tiles not showing up in folder

texts:
I am using GeoServer 2.5.1 and am attempting to use the integrated GeoWebCache. I have installed tomcat and am using the web app process for the GeoServer. I managed to add all the tiles and got the JAI working.

When I had tested out the first time of installing GeoServer and the tiling on my local host it worked in seeding the layers and creating cached tiles. The only problem is I don't exactly remember what I did and all the things that I remember what was don't I did it.

I then installed the Tomcat and GeoServer on our developer server so it is now using an IP address instead of a local host, but for some reason all the things I tried with the localhost does not seem to work on the developer one.

I also noticed when I either add:

<context-param>
    <param-name>GEOWEBCACHE_CACHE_DIR</parm-name>
    <parm-value> C:\temp</param-value>
</context-param>

then it seems to make pink tiles on my map.

and when I create a folder after setting up the environmental variables at the system variables in the C:\temp, when I create a folder in my C drive that matches the environmental variables cause I thought you might have to create the folder for it to work but it just seems to make my GeoServer work disappear and makes it as if has installed all over again but when I delete the folder then I get my data back.

Can any one help with this problem. I have looked through all the questions that vaguely match what I'm trying to do but to much avail I'm still having problems and the GeoServer and GeoWebCache documentation is not making too much sense.

When I click on the side of the GeoServer under service Capabilities when I click the WMS section it goes to an HTTP 404 page on either one of those but when I click on my localhost it work fine and gives me an XML page.


texts: Exporting/Printing to PDF, raster basemap has stripe(s) or band(s)

texts:
Now that my map is pretty and perfect, I'm discovering that I cannot export it properly. I have tried Exporting the Map to PDF, PNG, TIFF, etc. I have tried various settings of image quality and DPI. I have tried various compression methods. Gave up on exporting.

Tried using the Adobe Printer to print a PDF instead. That gives me worse results, with more bands or other artifacts. I've played with some of the Adobe Printer settings with no positive result.

In every type of export/print, I have at least one very obvious horizontal band.

The basemap is not transparent at all (though I wish it was). It is sourced from NCOneMap, a public North Carolina mapping portal/GIS service.

Now that I've spent my morning trying to export a single map (shakes fist at ESRI), does anyone else have some ideas?


texts:
I have a grid over a rasterlayer and a second a line .file.

I want to place the grid paralell to the line shape or dock it on the shape file.

I tried it with turning the grid manually, but it takes a huge amount of time to fit it correctly.

Is there a way to turn and place this grid?
Maybe there is also a way to genereate the gird directly along the line?

Best, mike


texts: QGIS rule-based point style responding to number of co-located objects?

texts:
I am mapping the geolocated location information for 5.5k Twitter profiles as point data. Many of them have the same location ie "London, UK" meaning that they return identical coordinates. On the map these overlaid features appear the same as a single feature, meaning density is hidden.

I would like to create a symbol style that aggregates all the features with the same co-ordindates as a single symbol sized according to how many features are located there. Can anyone help with this?

Multiple overlaid features, next to a single feature, look the same


texts: Reading binary rainfall data (*.trm) in ArcGIS for Desktop?

texts:
I have some *.trm binary data that is merged Tropical Rainfall Measuring Mission (TRMM) rainfall data.

I have to read it through ArcGIS for Desktop to extract rainfall data from this.


texts: Seeking bulk geocoding service for Brazil that allows permanent lat/long data storage

texts:
I have to geocode a huge dataset, with about 35 million addresses from Brazil. I was not able to find a good solution, even considering paid geocoders. So, I'd like to know if there is any geocoder service that has good accuracy and matching rate and, at the same time, allows bulk geocoding and also allows the storage of the lat/long results in a permanent database.
My initial idea was to do the geocoding using R, which has packages for many geocoders, like OSM/Nominatim (geocode_OSM), Google (ggmap) and HERE (HereR). Another possibility was to use QGIS, which also has geocoding plugins with the same 3 geocoders.  So, I did a short benchmark with these geocoders and noticed that it would not be possible to use any of them:

OSM/Nominatim (https://nominatim.org/): It's free, but does not allow bulk geocoding. Its accuracy and matching rate is also below my requirements.
Google Geocoding API: Has a good accuracy and matching rate. However, terms of use do not allow geocoding results storage. According to https://cloud.google.com/maps-platform/terms/maps-service-terms: "Customer can temporarily cache latitude (lat) and longitude (lng) values from the Geocoding API for up to 30 consecutive calendar days, after which Customer must delete the cached latitude and longitude values." Because of this restriction, I understand that it would not be allowed to do what I have to do using Google API;
HERE Geocoding API: Also has a good accuracy and matching rate, but with my test dataset both, accuracy and matching rate, were a little bit bellow than Googles API results. However, HERE API has similar terms of use restrictions than Google, not allowing results storage (https://knowledge.here.com/csm_kb?id=public_kb_csm_details&number=KB0016408).

I noticed that most of the API geocoders have similar restrictions than the observed in Google and HERE APIs. That's probably because they are web services created for apps, but not for bulk geocoding.
So, I’d like to know if there are any other bulk geocoding options (free or paid), without lat/long results storage restrictions.


texts: Configuring GeoServer on AWS to serve Cloud Optimized GeoTIFFs (COGs) from S3

texts:
I want to use S3 buckets as data stores for ImageMosaics on GeoServer in AWS.  I'm using the Support for S3 GeoTIFF and the ImageMosaic module. I'm using geoserver 2.16.
First of all, is this possible?
The S3 Support module documentation is a bit thin, but it states that you can refer to images using a convention from the module configuration file.
Is it possible to use that convention for each granule in the image mosaic?
Also, the configuration file for the S3 module requires a username and password.
What account is this referring to?
If I have enabled programmatic access for the EC2 instance to read and write to the S3 bucket, why would I need a user name and password?


texts: Creating layer depicting area closest to specific polygons in ArcMap?

texts:
I have a set of oddly shaped and vastly different sized polygons spread across a large area. There are three 'hub' polygons in this feature dataset. I need to figure out which of the hub polygons is closest for all of the other polygons in the feature dataset. distance just needs to measure from the edge of one polygon to another, not the central point. Optimally I need to create an output which maps the boundaries between the three hubs to show at what point the closest hub polygon swaps.

I've tried some different methods with the proximity tools but can't seem to get it working correctly.


texts: Merge Linestring that intersects without making them MultiLinestring

texts:
I am working on PostGIS 2.5.3 and PostgreSQL 11.2

As seen in the example below, I have 3 LineString that I would like to Merge to form one single Linestring. The squares represent the vertices, so we see there isn't any vertex at the intersection of the green and the blue LineString. Currently the only option I found is to use ST_Union(), which produce a MultiLineString since it simply merges the geometries without any special consideration.

Do you know a way to have the proper result, here a single LineString which crosses itself (as I would draw if I were to trace it manually) ?

EDIT :
Improving from obchardon's answer, here is an example to work with :

SELECT ST_AsText(ST_LineMerge(ST_Union(geom)))
FROM
-- Three dummy LineString for the example
(SELECT (ST_Dump(ST_GeomFromText('MULTILINESTRING((0 0,3 3),(3 3,1 3), (1 3, 3 2))'))).geom) AS T


texts: Setting Flag in layout to get selected vector features to appear selected in output PDF using PyQGIS

texts:
It would seem that I am using the wrong properties setting if I had to guess.  Anyhow, I have a layout template in a project that I load using python.  Within the project canvas, I can select a set of polygon features and zoom the layout map to those features.  When I export the layout as a pdf the selected polygons don't show up as selected.
Is there a setting (as in the title above) that can turn layout rendering on for the exported pdf?
Here is the line I've used:
settings = iface.mapCanvas().mapSettings()
settings.setFlag(QgsMapSettings.Flag.DrawSelection,True)

This does nothing...am I actually setting the flag for the map canvas which isn't picked up by the layout map item?


texts: Temporarily disable snapping with Arcobjects add-in tool

texts:
I've built an Arcobjects add-in tool that creates point features OnMouseDown. I set up the tool's snapping environment with ISnappingEnvironment using code from the Arcobjects Help page Working with the ArcGIS snapping environment and it works great.

When using tools in ArcMap (e.g. Measurement tool, Create features), users have the option to temporarily disable snapping by holding down the spacebar. However, when using the add-in tool I created, holding down the space does not temporarily disable snapping.

How can I plug that spacebar functionality into my programmatic snapping environment?

I've tried checking if snapping is enabled via ISnappingEnvironment.Enabled and then UnInitializing the snapping feedback, but it seems that holding the spacebar does not change ISnappingEnvironment.Enabled to false.

The snapping code I am using is from "Working with the ArcGIS snapping environment", but here it is:

ISnappingEnvironment pSnappingEnv;
IPointSnapper pSnapper;
ISnappingFeedback pSnappingFeedback;
IPoint pPosition;

protected override void OnActivate()
    {
        InitializeSnapping();
    }

protected override bool OnDeactivate()
    {
        if (pSnappingFeedback != null)
        {
             pSnappingFeedback.UnInitialize();
        }
    }
private void InitializeSnapping()
    {
        // Get the snap environment and initialize the feedback
        UID snapUID = new UID();
        snapUID.Value = "{E07B4C52-C894-4558-B8D4-D4050018D1DA}";

        pSnappingEnv = ArcMap.Application.FindExtensionByCLSID(snapUID) as
        ISnappingEnvironment;
        pSnapper = pSnappingEnv.PointSnapper;

        pSnappingFeedback = new SnappingFeedbackClass();
        pSnappingFeedback.Initialize(ArcMap.Application, pSnappingEnv, true);
    }

protected override void OnRefresh(int hDC)
    {
        if (pSnappingFeedback != null)
        {
            pSnappingFeedback.Refresh(hDC);
        }
    }

protected override void OnMouseMove(MouseEventArgs arg)
    {
        //Get the current position in map units... This is the cursor location before snapping
        pPosition = ArcMap.Document.ActiveView.ScreenDisplay.DisplayTransformation.ToMapPoint(arg.X, arg.Y);

        if (pPosition != null)
        {
            ISnappingResult snapResult = null;
            //Try to snap the current position
            snapResult = pSnapper.Snap(pPosition);
            pSnappingFeedback.Update(snapResult, 0);

            if (snapResult != null)
            {
                pPosition = snapResult.Location; //Snapping occurred. Set the current position to the snapped location
            }
            else
            {
                IPoint point = ArcMap.Document.ActiveView.ScreenDisplay.DisplayTransformation.ToMapPoint(arg.X, arg.Y);
            }
        }
    }


texts: Simplify a road network to obtain the intersections

texts:
In QGIS 2.6.1 in Ubuntu, I want to measure the density of block in this city (FYI my final aim is to measure the TOD level). I have the road network from French national IGN dataset BDTopo.
I have measured the number of intersections using the QGIS tool "Intersection" by computing the intersections of the network.

As shown on the image of the intersection window tool, I selected the same layer twice 'CHEMIN' and used the idea of each edge 'ID'.

The result is shown here. My problem is that some of the generated intersections are not real intersection. They are just two edges with different characteristics that connect, like the one selected here, which have two different physical width (shown in the map). The field 'ID' was used for the intersection, I tried with a different field, the name of the street and the result is roughly the same, there are too many intersection, more than the ones I want to obtain.

Now I want to simplify my network and avoid having a change of edge ID each time the width of the street changes or when there is a bridge, as is currently the case in the dataset. I would like to have less edges and merge edges when they belong to the same street basically.


texts: Most appropriate weights for interpolating temperature for a point from raster data

texts:
I have raster data for minimum and maximum temperature from NOAA's CPC Global Daily Temperature data set found here https://psl.noaa.gov/data/gridded/data.cpc.globaltemp.html
I also have survey data, which includes GPS/date information. When I merge the temperature data from the raster file into the survey data set, some of the temperature data is missing for coastal cities. NOAA only interpolates land squares, and not water squares, so some squares on the coast are mostly water and not interpolated.
My solution is to take an average of the temperatures of the adjacent non-missing squares weighted by the distance from their center to the GPS coordinates of the survey site.
What sort of weight is typical to use for interpolating temperature data? I'm currently using inverse-squared, but wondering if a linear approximation is better for temperature?

Update in case anyone else ever sees this post...
Since I didn't hear back, I did some regression analysis for different types of weights. I imputed temperature for survey sites that fell into raster squares I had data for as the average of the 8 adjacent squares using inverse distance weights, inverse-squared weights and weighting the closest square 100%. I then regressed the non-missing data on the imputed temperature data for each model. I also varied the regressions by starting with all available data, and then gradually dropping observations with more non-missing surrounding raster squares. I found that the linear weights worked best for all data (mostly inland squares) while the "closest square" weight worked best for observations with 4 or fewer non-missing adjacent squares (presumably the coastal squares). The inverse-square distance weight were a close second at both extremes of the analysis, so I decided to go with that.


texts:
I was applying the vegetation index
($n1_ttc0003r_16(2) - $n1_ttc0003r_16(4)) / ($n1_ttc0003r_16(2) + $n1_ttc0003r_16(4))
to a multispectral image, and the raster resulting on this operation had all its values equal to -1.#IND. Can anyone help me understanding what's happening? I think this means that some pixel values of the output raster are the result of a division by zero. How can I solve this?


texts: Event listeners after each plotting of a point of a Polygon Vector in OpenLayers

texts:
I am quite new to OpenLayers. I have a polygonLayer for users to draw polygon vectors. The plotting feature is working fine.

polygonLayer = new OpenLayers.Layer.Vector("Polygon Layer");
map.addLayers([ polygonLayer ]);
var polygonControl = new OpenLayers.Control.DrawFeature(polygonLayer, OpenLayers.Handler.Polygon);
map.addControl(polygonControl);
polygonControl.activate();
var editControl = new OpenLayers.Control.ModifyFeature(polygonLayer);
map.addControl(editControl);
editControl.activate(RESIZE);

What I want to do is to have event listeners after each plotting of a point. Supposedly, I want to show popups to assist the user in between each point of plotting. How could I achieve it?


texts: Water mask from NDVI from 8 day NDVI in Earth Engine

texts:
With the following code, I am trying to make a water mask from the Landsat 8's LANDSAT/LC08/C01/T1_8DAY_NDVI dataset.

The bit I am stuck on is how to create an if statement that will change the values of the pixels to 1 if greater than 0 and to -1 if less than zero.

The goal is to have the water in blue and everything is in white, with the ability to investigate the rivers visual size and shape at times of interest.

I have tried using ee.Algorithm.If but I am missing something in the syntax or it is not the best method?

https://code.earthengine.google.com/deab597e9a22de0bdde2876c1ab9f1ff

//Water NDVI Mask
var NDVI_dataset = ee.ImageCollection('LANDSAT/LC08/C01/T1_8DAY_NDVI')
                  .filterDate('2016-10-01', '2016-10-15')
                                    ;
var NDVI_Visualisation = NDVI_dataset.select('NDVI');

var water = ee.Algorithms.If(NDVI_Visualisation.set(>=0, 1, -1));

//Display the NDVI
Map.addLayer(NDVI_Visualisation,
 {bands: 'NDVI', min: -1, max: 1, gamma: 1.6, opacity: 1,},
 'NDVI Visualisation', false);

 var VIS_WATER_MASK = {
  min: -1,
  max: 1,
  palette: ['blue', 'white']
};

Map.addLayer({
  eeObject: NDVI_Visualisation,
  visParams: VIS_WATER_MASK,
  name: 'Water',
  shown: true
});


texts:
I am having difficulty with layer order. I have 4 layers but zIndex is not working for L.tileLayer. Here's a part of  my code: https://pastecode.xyz/view/d726ccb1

 var overlay_Numerydziaek1 = L.tileLayer.wms('http://mapy.geoportal.gov.pl/wss/service/pub/guest/G2_GO_WMS/MapServer/WMSServer', {
    layers: 'NumeryDzialek',
    opacity: 1.0,
    format: 'image/png8',
    transparent: false,
    continuousWorld : true,
    zIndex:110,
}).addTo(map);

var overlay_Dziaki2 = L.tileLayer.wms('http://mapy.geoportal.gov.pl/wss/service/pub/guest/G2_GO_WMS/MapServer/WMSServer', {
    layers: 'Dzialki',
    format: 'image/png',
    uppercase: true,
    transparent: false,
    continuousWorld : true,
    tiled: true,
    info_format: 'text/html',
    opacity: 1,
    zIndex:100,
    identify: false,
});

map.addLayer(overlay_Dziaki2);
overlay_Dziaki2.bringToFront()
var overlay_googlemapsstelita_1 = L.tileLayer('http://mt0.google.com/vt/lyrs=s&hl=en&x={x}&y={y}&z={z}', {
    opacity: 1.0
    zIndex:50,
});
var overlay_OpenStreetMap_1 = L.tileLayer('http://a.tile.openstreetmap.org/{z}/{x}/{y}.png', {
    opacity: 1.0,
    zIndex:10,
});
map.addLayer(overlay_OpenStreetMap_1);
overlay_OpenStreetMap_1.bringToBack()

With this code var overlay_googlemapsstelita_1 appears on top
and the second is var overlay_OpenStreetMap_1.
Then var overlay_Numerydziaek1 and var overlay_Dziaki2.
I want overlay_Numerydziaek1 to appear on top.
How can I do this ?
For other non L.tileLayer, zIndex is working fine.


texts:
I'm planning on building a transportation related GIS App.

In the app, users will be allowed to add their own links to the network. For example, let's say a user wants to add a "driveway" in front of their house or a trail that's not on the OSM map. Is it possible to do so? Add to that, these additions should only be visible to the user that added those links and to whomever else they give access too.


texts:
I have a little bit of dilema here.

I have a hierachical collection of classExtensions that I very much would like to keep them as they are now:

Timestamper (master concrete class Extension)
1.1) BaseCadastralExtension (derives from Timestamper)
Derived classes (from BaseCadastralExtension)

But I want to implement custom IObjectInspectors for each of the derived classes. I can't keep them deriving from BaseCadastralExtensions, because the IObjectInspector needs to derive from UserControl (am I right?). So, I cannot do this with multiple inheritance :P

Is it possible to implement an IObjectInspector without trashing the class structure? Or do I need to rewrite? I wonder that composition in this matter (baseCadastralExtension to have a "funcionalities" member to take care of all the events, etc) and implement for each a custom IObjectInspector would do better?

I'm not sure if I'm being clear, but I would like suggestions from you expert ESRI devs.


texts: How do I upload a file geodatabase to FME Server 2013?

texts:
I periodically load an SDE geodatabase with about 1 GB of data that is delivered to me in an ArcGIS file geodatabase (.gdb) format. In order to run that job on FME Server I need to upload the .gdb to FME Server. I tried to upload it in various ways but since the .gdb is a directory containing lots of files in Windows, the upload dialog box does not accept the directory, instead it wants a file. I also tried to zip the .gdb, and then I got a success but no data was read from the .gdb. This is an excerpt from the log file:

2014-05-02 12:29:24|   1.9|  1.6|WARN  |Generic Reader: Unable to guess a valid reader name from dataset 'D:\FME_TEMP\fmeservertst_10_0.gdb.zip_1399026562740_14984'. A NULL reader will be used instead.
2014-05-02 12:29:24|   1.9|  0.0|INFORM|Generic Reader: No format specific settings have been found for the NULL Reader -- default values as documented in the Reader/Writer manual will be used

Can anyone put me on the right track? I use FME Server 2013 and ArcGIS 10.0.


texts: Finding route of fixed length, starting and ending at same point using ArcGIS Network Analyst?

texts:
This is a problem I would like to solve with Network Analyst in ArcGIS 10.3x.  The goal is to find a running route for a road race, using a slope attribute (which I have assigned to road segments) as the cost to find the flattest route meeting the given length and given start/end (same point).  I have previously used Network Analyst to find the shortest/least cost route from point A to point B, but I'm not sure how to accomplish the problem I stated above.

Does anyone have an idea of how to do this?


texts:
How do I draw a polygon in a shapefile that I've created from joined points of 2 .csv files?

Which specific tool do I use?


texts:
I have proj4 version 4.9.2 installed in a non-standard location within a Linux environment (/p/bin/proj4). I know that proj_def.dat did not go out with version 4.9.1, but was of the understanding that it was included in version 4.9.2. I have verified that it works, I have added /p/bin/proj4/share/proj to my PATH. In addition, I downloaded the datum shifting files for North America to proj/conus. When I run R and load rgdal, I receive the message that Path to PROJ.4 shared files is autodetected. This is in contrast to gdal, which shows the correct path to my installation of gdal. I provided the following argument when I installed rgdal "--with-proj-share=/p/bin/proj 4/share/proj".

However, when I run raster(), I receive the following message:

  rgdal::checkCRSArgs:: no proj_defs.dat in PROJ.4 shared files

I have verified that there is a proj_def.dat in proj/ folder.

Is their an argument that I can provide when loading rgdal that will point to the proj_defs.dat files.


texts: Using linguistic phonetic alphabet in QGIS?

texts:
I want to know if there is a way to use the International Phonetic Alphabet in QGIS?

I am working with a group of linguistics that want to map dialectic varibility in a language (Nahuatl) and I am not able to see the specific characters of a word, usually they appear as interrogations characters.

I tried loading the data in UTF-8 and system.

I am working with windows 10 64 bits in Spanish.


texts: Algorithm to get the nearest line segment (link) to a point

texts: What is the fastest algorithm to get the nearest line segment to a point (that is not on any line segment) in a GIS database?


texts: Publishing Mapbox styles using V2.1 and above

texts:
I am trying to integrate mapbox maps/styles in HTML code. I succeeded in doing so using:

<script src='https://api.tiles.mapbox.com/mapbox-gl-js/v0.14.1/mapbox-gl.js'></script>
<link href='https://api.tiles.mapbox.com/mapbox-gl-js/v0.14.1/mapbox-gl.css' rel='stylesheet' />

mapboxgl.accessToken = 'pk.eyJ1IjoibHB0ZXNzaWVyMTMiLCJhIjoiMEhrZUZaTSJ9.xU3ECSdq6K3LUkUuE8mZpQ';
var map = new mapboxgl.Map({
    container: 'map',
    style: 'mapbox://styles/lptessier13/cin5q5ilo00skctnfgeyrglfc',
    center: [6.08, 46.22],
    zoom: 11
});

Unfortunately, I cannot publish my maps/styles using the following:

<script src='https://api.tiles.mapbox.com/mapbox.js/v2.1.5/mapbox.js'></script>
<link href='https://api.tiles.mapbox.com/mapbox.js/v2.1.5/mapbox.css' rel='stylesheet' />

L.mapbox.accessToken = 'pk.eyJ1IjoibHB0ZXNzaWVyMTMiLCJhIjoiMEhrZUZaTSJ9.xU3ECSdq6K3LUkUuE8mZpQ';
    var map = L.mapbox.map('map','mapbox://styles/lptessier13/cin5q5ilo00skctnfgeyrglfc')
        .setView([6.08, 46.22], 11);
    map.scrollWheelZoom.disable();

I am not yet well versed in coding, but I thought that the problem may be residing in the API (mapbox-gl-js/v0.14.1 VS mapbox.js/v2.1.5)

If necessary, I could add the whole HTML file which is not very long.


texts: Combining adjacent polygons sharing similar trait into single polygon using ArcGIS Desktop?

texts:
I have a large dataset of habitat classes which I have further simplified down into my own unique classes.

However, it resulted in there being habitats adjacent to each other with the same forest classification system.

Is there any way to join these two (or however many there are touching that are classed as the same type) adjacent polygons into a single larger polygon?

When I use the tool Dissolve, it takes all of the habitats of the same class for the entire map and creates a single polygon.


texts:
I have a map of fracture distributions over an area and I would like to calculate the density distribution of fractures that intersect a moving window with a specified radius.

Is there a way to calculate this in QGIS, similar to the line density calculation in ArcGIS Spatial Analysis?


texts: Free Open Source Geocoders using Local Address Data

texts:
We publish a geocoding service using ArcGIS server but are interested in moving to a free and open source alternative but can't seem to find anything that allows us to use local address data. Are there any geolocators that are both free open source and have the ability to utilize local data? The alternatives that I've evaluated (including those mentioned here) are either not free or don't support the use of local data:

TIGER geocoder for PostGIS: built primarily for TIGER features and address ranges. Can't really find any documentation on how to modify it to incorporate local data instead of TIGER data
Geocommons: Free and open source, but doesn't accept local data
Geocode.xyz: Same issues as Geocommons
Gisgraphy: Only accepts OSM and GeoNames data
Mapzen/Pelias: Free up to 25K requests/month; doesn't support local data
Photon: Same issues as Mapzen


texts:
I installed a mapserver on my own server with mapnik, apache and mod_tile.

I build an application with WorldWind and my map server. My problem is that there are blank tiles instead of blue tiles to represent the sea.

Is there a simple way to remplace these tiles? Or do I have to change files like generate_tiles.py for example?


texts: Appending polygon feature classes with already repaired geometry => WARNING 000442: self intersections from Check Geometry?

texts:
I have an ArcGIS 10.2.2 for Desktop process which:

Uses the Data Interoperability extension and the Feature Class To Feature Class tool to convert 88 MapInfo TAB files (with the same schema, representing different geographic areas) to 88 file geodatabase feature classes.
Performs a Repair Geometry on each of those 88 feature classes
Performs a Create Feature Class, Project Define (using one of the 88 feature classes as the spatial reference source) and then Append (with NOTEST) on those 88 feature classes to get a merged feature class.  Note: I did not use the Merge tool because when I tried to do that only about 6 feature classes were merged before the tool reported that it had completed successfully.

The process appeared to work fine to create the merged feature class of 2.5 million polygons.  However, when I ran Check Geometry on this feature class it started to report self intersections before appearing to hang and then failing with an ERROR 999999 five minutes later..

Executing: CheckGeometry C:\Temp\DCDB\Staging.gdb\DCDB C:\Temp\DCDB\RWC_CheckGeom.dbf
Start Time: Fri Jun 03 09:42:08 2016
WARNING 000442: self intersections at 883 in C:\Temp\DCDB\Staging.gdb\DCDB
WARNING 000442: self intersections at 1479 in C:\Temp\DCDB\Staging.gdb\DCDB
WARNING 000442: self intersections at 6003 in C:\Temp\DCDB\Staging.gdb\DCDB
WARNING 000442: self intersections at 6004 in C:\Temp\DCDB\Staging.gdb\DCDB
WARNING 000442: self intersections at 6015 in C:\Temp\DCDB\Staging.gdb\DCDB
WARNING 000442: self intersections at 6016 in C:\Temp\DCDB\Staging.gdb\DCDB
WARNING 000442: self intersections at 6417 in C:\Temp\DCDB\Staging.gdb\DCDB
WARNING 000442: self intersections at 7285 in C:\Temp\DCDB\Staging.gdb\DCDB
WARNING 000442: self intersections at 11813 in C:\Temp\DCDB\Staging.gdb\DCDB
WARNING 000442: self intersections at 12127 in C:\Temp\DCDB\Staging.gdb\DCDB
WARNING 000442: self intersections at 14861 in C:\Temp\DCDB\Staging.gdb\DCDB
WARNING 000442: self intersections at 16328 in C:\Temp\DCDB\Staging.gdb\DCDB
WARNING 000442: self intersections at 16847 in C:\Temp\DCDB\Staging.gdb\DCDB
WARNING 000442: self intersections at 16893 in C:\Temp\DCDB\Staging.gdb\DCDB
WARNING 000442: self intersections at 21008 in C:\Temp\DCDB\Staging.gdb\DCDB
WARNING 000442: self intersections at 22254 in C:\Temp\DCDB\Staging.gdb\DCDB
WARNING 000442: self intersections at 25064 in C:\Temp\DCDB\Staging.gdb\DCDB
WARNING 000442: self intersections at 26671 in C:\Temp\DCDB\Staging.gdb\DCDB
WARNING 000442: self intersections at 26872 in C:\Temp\DCDB\Staging.gdb\DCDB
WARNING 000442: self intersections at 27513 in C:\Temp\DCDB\Staging.gdb\DCDB
WARNING 000442: self intersections at 27521 in C:\Temp\DCDB\Staging.gdb\DCDB
WARNING 000442: self intersections at 27540 in C:\Temp\DCDB\Staging.gdb\DCDB
WARNING 000442: self intersections at 27579 in C:\Temp\DCDB\Staging.gdb\DCDB
WARNING 000442: self intersections at 27600 in C:\Temp\DCDB\Staging.gdb\DCDB
WARNING 000442: self intersections at 27603 in C:\Temp\DCDB\Staging.gdb\DCDB
ERROR 999999: Error executing function.
Failed to execute (CheckGeometry).
Failed at Fri Jun 03 09:47:19 2016 (Elapsed Time: 5 minutes 10 seconds)

Should appending polygon feature classes with already repaired geometry be expected to produce a feature class that can pass the Check Geometry test?


texts: Accurate GPS apps for smartphones available?

texts:
I'm working on a project where we need to make accurate GPS measurements and the required precision would be in centimeter level. I've been doing some searching and not yet found any solid solutions so I make the question:

"Are there any free accurate GPS apps available for smartphones?"


texts: Splitting lines at multiple points in ArcGIS

texts:
I have two shapefiles:

First shapefile consists of lines features (approx 1500); it represents the public transport lines
Second shapefile consists of point features (approx 3000); it represents the public transport stops on the lines

Now, I want to do a analysis on it and I want to split each line feature into half according to the points. For example I have one line with starting point A and ending point B and 3 point features (lets say a,b,c) on it. I want to split the line and calculate the distance as follows:

Split the line in two parts one is from A-a and a-B
Split the same line two part one is from A-b and b-B
Split the same line two part one is from A-c and c-B

This is what I want to do for the entire network but doing it manually will take a lot of time. How can I  identify the tool which can do it fast and without any error?


texts: Add raster attribute table to raster using arcpy

texts:
I have downloaded US soil raster data by state from https://gdg.sc.egov.usda.gov/GDGOrder.aspx at 10m resolution.

I am now trying to add a modified raster attribute table to tif that was exported from the geodatabase using arcpy, but have been unable to do so. The existing tif has an attribute table with an ID that matches the IDs in the modified attribute table.

I have tried to (1) make the raster a layer (not sure this is the right approach), (2) add a table using AddJoin_management, and (3) export the results. Obviously, I am missing something.

import arcpy
from arcpy import env

env.workspace = "C:/Data"
inFeatures = "InRaster.tif" ## This is tif raster exported from geodatabase
layerName = "raster_layer"
rat = "rat.dbf" ## attribute table for joining
arcpy.MakeRaster_management(inFeatures,layerName) ## This doesn't work, with or without extension.

tmp = arcpy.AddJoin_management(layerName,"Value",rat,"ID") # This works, but is a result and I am not sure how to export it.


texts:
I am currently working with GOES-R on AWS. I've downloaded recent data, which is in the netCDF file format. Specifically I am using band 3 (M3C03).

When loading it into GDAL, it is taking upwards of 5 minutes.

goes_dir = "goes"
goes_files = glob.glob(goes_dir + '/' + str(doy) + '/18/*M3C03*.nc')
goes_files[:3]

['goes/257/18/OR_ABI-L1b-RadF-M3C03_G16_s20172571815379_e20172571826146_c20172571826188.nc',
 'goes/257/18/OR_ABI-L1b-RadF-M3C03_G16_s20172571845379_e20172571856146_c20172571856196.nc',
 'goes/257/18/OR_ABI-L1b-RadF-M3C03_G16_s20172571830379_e20172571841146_c20172571841193.nc']

go = gdal.Open('NETCDF:' + goes_files[1] + ':Rad')
import time
start_time = time.time()
goes_arr = go.GetRasterBand(1).ReadAsArray()
print("--- %s seconds ---" % (time.time() - start_time))
# --- 277.96904397 seconds ---

What can I do to speed it up? The files are only 50 MB maximum.


texts: QGIS Marker Lines (every N map units) shifting as I pan the map

texts:
I'm drawing a map of the West Highland Way, a long-distance hiking trail in Scotland.

I want to place markers every 1 km and every 1 mile on either side of the line.

I'm using QGIS 2.18.0

Successfully downloaded the trail from OSM using the QuickOSM plugin - really easy in this case, as it's a single feature.

To tidy up the geometry, I used the Join Multiple Lines plugin to make sure all the lines point in the same direction.

I normally use QChainage plugin for this but thought I'd try Marker Lines style, using "mark every ... map units" (I'm using EPSG:27000 for this line, OTF is also on EPSG:27700)

I've added a marker line style to confirm that the lines are joined end-to-end in a consistent direction.

Here's a composite image showing the problem. Note how the kilometer markers have shifted after panning.

Can anyone explain this?

I'm not sure if a bug, or a consequence of how things are rendered (perhaps QGIS doesn't render the whole line, and is starting the 'chainage' somewhere closer to the start of the canvas extent?)

It doesn't seem to happen in Print Composer (possibly because the whole feature is rendered in one go?)


texts: Why are updates made in MXD by my tool not visible on some machines?

texts:
I've developed a tool that makes updates in many mxd files.
When it runs in my machine, everything works fine, and I can open the mxd and see the updates made by the tool.
When the tool runs in any other machine, the updates are made, but when open the mxd to see the updates, they can't be seen. When I open in my machine, I can see the updates.
Both machines has the same ArcMap configurations, ans service packs.
Any Ideas?


texts: arcgis 10.2 .shp with no projection info but .xml file


texts:
I'm trying to optimize a route SQL function and playing with some ST_* functions I found around.

But even found a good one, I don't want to simply copy/paste what someone wrote. I want to understand what is happening in underground. Yes, I know: RFM. but even reading this I don't know why the function combinations can make one better then other.

Can someone help me by describing what each function is doing? Other optimizations will be wellcome!

$1 = source
$2 = target
$3 = k
$4 = directed

Option 1: No optimization.

    SELECT
        *
    FROM
        pgr_ksp(
         'SELECT id, source, target, cost, reverse_cost FROM osm_2po_4pgr',$1, $2, $3,
directed:=$4)

Option 2 (faster): Explain ST_Expand(ST_Extent(geom_way),4).

SELECT
    *
FROM
    pgr_ksp(
     'SELECT id, source, target, cost, reverse_cost FROM osm_2po_4pgr as r,
             (SELECT ST_Expand(ST_Extent(geom_way),4) as box FROM osm_2po_4pgr as l1
            WHERE l1.source = ' || $1 || ' OR l1.target = ' || $2 || ') as box
            WHERE r.geom_way && box.box',$1, $2, $3, directed:=$4
    )

Option 3 (even faster): Explain ST_Buffer(ST_Envelope(ST_Collect(geom_way)), 4).

SELECT
    *
FROM
    pgr_ksp(
     'SELECT id, source, target, cost, reverse_cost FROM osm_2po_4pgr as r,
             (SELECT st_buffer(st_envelope(st_collect(geom_way)), 4) as box FROM osm_2po_4pgr as l1
            WHERE l1.source = ' || $1 || ' OR l1.target = ' || $2 || ') as box
            WHERE r.geom_way && box.box',$1, $2, $3, directed:=$4
    )


texts:
I used ArcCatalog 10.2.2 to connect to enterprise database (SQL SERVER 2012) I've created before.

Some feature classes load just fine. But others show this error when I try to view data table

This feature class works fine in personal DB, this error only happened in Enterprise DB.

Thing I've tried:

Re-install SQL Native Client 12
Export error feature class and re-import to DB

Anyone has a solution for this?


texts: Merge connected lines IF they share adjacent polygon(s) (arcgis/arcpy)

texts:
I have lots of lines that I need to merge IF and only if they share the same polygon(s). All lines are made from the polygons (copy-paste into line) so they all snap.

In this little example line 499, 451, 447, 452 and 448 could be merged since they all share polygon 136 and none polygon on the outer side of the line.
Similarly line 442, 443 and 427 could be merged since they all share polygon 133 on the inner side and polygon 136 on the outer side.

My guess so far is to make points on the middle of each line (Features->Feature Vertices to Point and point type [MID]). And then do a Spatial Join where I use this new point layer as Target and my polygonlayer as Join [Join-One-To-Many]. Then I get a new feature that tells me the line-to-polygon relation; line = ORIG_FID and polygon JOIN_FID.

Line 427 is adjacent to polygon 133 and 136.
Line 443 is adjacent to polygon 133 and 136.
Line 447 is adjacent only to polygon 136.

And this is where I stop.
I now need a script that checks which polygon(s) each line is adjacent to and if two or more lines share exactly the same polygon or polygons then merge.

UPDATE:
I figured out that if I could make an unique ID on each link based on adjacent polygon(s) [Pol_ID], I could use this new link-ID to merge links with the same new_ID value. Here I made it simple and added the Pol_ID value(s) from adjacent polygon(s); see black values in figure.

New_ID 136 = Pol_ID value 136.
New_ID 269 = Pol_ID values 136+133.

So far so well, and I easily could merge all links with the same new_ID.

But then a new problem arose: the lines also have to be connected. Therefore the two lines with value 133 (blue arrows) shall not be merged since they are not connected to each other.

The missing part is now to: check if links with same new_ID value is connected, then merge the connected ones.

The red values represent the unique value [X-coordinate + Y-coordinate] to the start-/endpoint of the line. Thought I could use this value to check if lines were connected (share start/endpoint), but haven't figured it out yet..


texts: Inconsistencies in st_contains function in R sf package

texts:
I am trying to identify, from a set of overlapping buffers, the smallest of their intersections. In other words (easier to understand if you look at the picture below), identify amonst all the intersections which are not containing another, smaller, one (the yellow one, below). I was very excited to find the st_contains() function in R sf a while ago, but it doesn't work as I thought.
Here is a working example.

> library(sf)
#creating a set of three points
> dfpt<-cbind.data.frame(X=c(328900,328400,331200),Y=c(403000,403400,402300),ID=1:3)
> sfpt<-st_as_sf(dfpt,coords=c('X','Y'))
# creating buffers of 2000m
> sfb<-st_buffer(sfpt,dist=2000)

# intersecting the buffers
> sfint<-st_intersection(sfb,sfb)

# I am interested in two buffers that are embedded
> sfsel<-subset(sfint,(ID==3&ID.1==1)|(ID==3&ID.1==2)) # whether it's ID or ID.1 might change

#check what it looks like
> plot(st_geometry(sfb),col=NA,border='black')
> plot(st_geometry(sfsel[1,]),col='blue',add=TRUE)
> plot(st_geometry(sfsel[2,]),col='yellow',add=TRUE)

# checking which of the three buffers contain which of the two polygons
> st_contains(sfb,sfsel)
Sparse geometry binary predicate list of length 3, where the predicate was `contains'
 1: 2
 2: (empty)
 3: (empty)

(the expected answer should be 1:2, 2:1,2,3:1,2, which I got in other, similar cases)

# checking if the two polygons contain each other (the first should contain the second!)
> st_contains(sfsel,sfsel)
Sparse geometry binary predicate list of length 2, where the predicate was `contains'
 1: 1
 2: 2

(the expected answer should be 1:1,2, 2:2)

What happens there? I think it might be the approximation of the polygons arcs that lead R to believe the border of the two polygons (blue and yellow above) are not perfectly aligned. I tried to change the scale (by increasing all distances by 1000) and it did improve the function in some cases, but not many. What can I do?

(For info)

> R.version
               _
platform       x86_64-w64-mingw32
arch           x86_64
os             mingw32
system         x86_64, mingw32
status
major          3
minor          6.1
year           2019
month          07
day            05
svn rev        76782
language       R
version.string R version 3.6.1 (2019-07-05)
nickname       Action of the Toes


texts:
I have an Excel file (with .xlsx extension) containing some non-spatial data that I want to use in combination with some other spatial data.

I generally know how to load an external source (Loading Layers from the PyQGIS Cookbook), but I haven't found any reference for files with .xlsx extension, so I tried to import it in this way:

uri = 'path_to_the_file/test.xlsx'
layer = QgsVectorLayer(uri, 'test', 'ogr')

but when I try to read something, for example all the elements from the first column:

for item in layer:
    print item[0]

I get this error:

'QgsVectorLayer' object is not iterable

What I'm doing wrong? If I previously convert my file to the .csv format, then there is not problem with the reading of the items in PyQGIS, but I would like to know if there is a direct way in order to avoid this passage.


texts: Why does the method partCount in ArcGIS always give back a value of 1?

texts:
Does anyone knows why method "partCount" in arcgis version 10.0 always gives me back value 1, even when I have object with interniorring or multipart object?

partnum = 0
partcount = feat.partCount
while partnum < partcount:
    outFile.write(str(inRow.GetValue(inDesc.OIDFieldName)) + " " + str(partnum) + "\n")
    gp.AddMessage(str(inRow.GetValue(inDesc.OIDFieldName)) + " \\ " + str(partnum))
    part = feat.getpart(partnum)
    part.reset()
    pnt = part.next()
    pnt_count = 0
    while pnt:
        outLine = str(pnt_count) + " " + str(pnt.x) + " " + str(pnt.y) + " " + str(pnt.z) + " " + str(pnt.m) + "\n"
        outFile.write(outLine.replace(".", sepchar))
        pnt = part.next()
        pnt_count += 1
        if not pnt: # ?
            pnt = part.next()
            if pnt:
                outFile.write("InteriorRing\n")


texts: Corruption of export using Merge feature in QGIS 3

texts:
Using QGIS 3.2.0 - Bonn
I'm getting a nasty corruption when trying to use the Merge feature to pull together a selection of .map (OziExplorer) files and merge into a single .bmp file. They load in perfectly into QGIS -

Then I select them all as input layers, change the data type to 'Byte' I assume that's the right one for .bmp - as 'Float32' default does not work?

Then I get an .bmp export of 4GB, but this odd corruption occurs:

If I switch to .tif - export works perfect 5.42GB

Am I just using the incorrect 'Output data type' ?


texts: Using function in pgPointCloud extension to PostgreSQL with PostGIS

texts:
I am using the pgPointCloud extension to PostgreSQL and have created a table of blocks/patches of points in a point cloud. Now, I am just trying to run a few of the functions in the pgPointCloud/postgis extensions on the blocks, such as PC_BoundingDiagonalGeometry(pcpatch) and PC_EnvelopeGeometry(pcpatch). However, both of these functions are failing with the same error message. For example, running

SELECT ST_AsText(PC_BoundingDiagonalGeometry(pa)) FROM patches;

returns

ERROR:  function pc_boundingdiagonalgeometry(pcpatch) does not exist
LINE 1: SELECT ST_AsText(PC_BoundingDiagonalGeometry(pa)) FROM patch...
                         ^ HINT:  No function matches the given name and argument types.
You might need to add explicit type casts.
SQL state: 42883
Character: 18

I do not understand what is causing the problem. I have already created the POSTGIS, POINTCLOUD, and POINTCLOUD_POSTGIS extensions in my database.I have also created the patches table accordingly as described on this page:
https://github.com/pgpointcloud/pointcloud


texts: Setting seed value for random selection in QGIS?

texts: I am using QGIS 2.18.11 and need to select 250 points from the existing points in a map. I know that I can randomly select points through Vector> Research Tools> Random Selection but I don't know how to make it "reproducible" (replicable). I don't know much about Python.


texts: Seeking online reference available for Fiona and Shapely libraries?

texts:
Is there an online reference available for the Fiona and Shapely libraries?

I'm thinking a real reference, beside the tutorial-style documentation:

http://toblerity.github.io/shapely/manual.html

http://toblerity.github.io/fiona/manual.html


texts:
After reading this Esri article on how to distribute Python Toolboxes, I'm wondering how redistribution works?  I have a few toolboxes that are in constant development, but would be useful to my co-workers now.  If I distribute them as the article suggests and they show up in the System Toolbox in ArcMap can I use the same process to update them?  Or is there some extra code that will need to remove the previous version of toolboxes before I distribute the updated versions.

I'm sure I just need to be pointed to a help file where it is all laid out, but I haven't found it yet.

Thanks


texts: Can't perform Helmert transformation in Proj.4

texts:
I'm trying to convert a point from one SRS to another using proj4 strings and GDAL. My origin SRS uses Helmert transformation to convert point to EPSG:2229 so I tried to use pipeline operator (https://proj4.org/operations/pipeline.html) along with helmert transformation (https://proj4.org/operations/transformations/helmert.html#examples) and got something like this

+proj=pipeline
+step +proj=helmert +convention=coordinate_frame +x=0 +y=0 +theta=0 +scale=1
+step +proj=lcc +lat_1=35.46666666666667 +lat_2=34.03333333333333 +lat_0=33.5 +lon_0=-118 +x_0=2000000.0001016 +y_0=500000.0001016001 +datum=NAD83 +units=us-ft +no_defs

Here is my code in python (simplified case):

    proj4_4326 = "+proj=longlat +datum=WGS84 +no_defs"
    proj4_2229 = "+proj=lcc +lat_1=35.46666666666667 +lat_2=34.03333333333333 +lat_0=33.5 +lon_0=-118 +x_0=2000000.0001016 +y_0=500000.0001016001 +datum=NAD83 +units=us-ft +no_defs"
    proj4_helmert = "+proj=helmert +convention=coordinate_frame +x=0 +y=0 +theta=0 +scale=1"
    srs_4326 = osr.SpatialReference()
    srs_4326.ImportFromProj4(proj4_4326)
    srs_2229 = osr.SpatialReference()
    srs_2229.ImportFromProj4("+proj=pipeline +step {} +step {}".format(proj4_helmert, proj4_2229))
    global_coords_calculated = convert_coordinates_cs(6427213.13, 1984970.215, srs_2229, srs_4326)
    print global_coords_calculated

But it ends up with error TypeError: in method 'Geometry_Transform', argument 2 of type 'OSRCoordinateTransformationShadow *. I've tried to google it and found this https://stackoverflow.com/questions/27572190/gdal-reprojection-error-in-method-geometry-transform-argument-2-of-type-osr
Seems that it means that some library is missing, but I've installed pyproj along with GDAL.
I am using virtualenv, python2 and macOS.


texts: How to create 4m x 4m survey grid in QGIS

texts:
I am trying to build a 4 km2 grid comprised of 4m x 4m squares in QGIS using the vector grid tool. However, in the parameters section of the vector grid window, I am unable to enter values less than 0.0001.

Is this a built-in limitation? Is there another way I should be trying to build such a grid?


texts:
Is there a way to know what version of pgrouting I have?
My OS is ubuntu 12.04.
I want to know my pgrouting version because the 1.05 version has some issues with shooting start and the directed graphs and the last stable version which has this algorithm is 1.03.


texts:
I have an openlayers vector layer with features added from a text file. I am now trying to access these features via their attributes (if layer.features.attributes.some_attribute equals to something hide the specific feature). However, I don't seem to be able to access the features at all. I cannot even count them. I think I need some basic advice on how this accessing is done.

My original question is Hiding specific icons/markers on text generated vector layer, and I'm trying to accomplish something like Openlayers - Vector Layer - Combining hide/show features using style property.

I thought I had all the clues I needed but must be making some basic coding error...


texts: I'd like to produce a simple walking map with some points of interest. Ideally, I'd like to use the rendered OSM data, as you can get from the page. Using the QGIS OSM plugin, I can only download a fraction of my area, and it is tricky to download a few PNGs at a high enough zoom level and stitch them together from OSM's "Export" tool. I only need to do this a few times, so I don't really see much point in downloading all the XML/Shapefile data, and re-rendering with Mapnik. Is there an easy way of getting hold of the already-rendered raster data fro a given region at a given zoom level in a format that I can easily use with QGIS/GDAL/OGR?


texts: QGIS TimeManager plugin: how to define the settings from the python console?

texts: rather than defining interactively the 'Settings' of the QGIS Time Manager plugin, I would like to have the possibility to set them using the python console.  This would save a lot of time in case one wants to animate a long list of rasters.


texts: Difference between SnakeGrid and latitude and longitude

texts:
I am not a geography subject expert and maybe my this comparison make no sense as i don't have enough knowledge (one of the reason to ask this question here).

I have a project where I am doing conversion from unity co-ordinate system to GPS(Lat/Lon) using a library. The conversion working fine but now I got the environment model which is based on SnakeGrid  and the same conversion is not accurate. So my question is  how is a SnakeGrid based environment different from a latitude and longitude based coordinate system?


texts: How to draw into selection cache in ArcMap

texts:
I have a tool,something like Edit Tool in ArcMap,when I double click on mouse,and it selects a feature (polyline for example) and draw a sketch on it too,but if I change zoom,the sketch disapears because I draw sketch on the scrren.I know if I can draw sketch on into selection cache (bitmap) it will stay till selection changes,so it doesn't disapear with changing zoom.
here are the steps i pass to draw the sketch with current situation:

I take cache count into cachecount variable using screen display object
I take a handle to device countext for cacheIndex=cachecount-1
start drawing using 2 parameters in steps 1 and 2
4.vdraw polyline using screen display DrayPolyline
FinishDrawing

.............................................

HRESULT hr;

//Get handle to memory device context

HDC HandleDC;

short cacheCount;

this->m_ipSDisp->get_CacheCount(&cacheCount);

short cacheIdx=cacheCount-1;

this->m_ipSDisp->get_CacheMemDC(cacheIdx,(OLE_HANDLE*)&HandleDC);

if(FAILED(hr=this->m_ipSDisp->StartDrawing((OLE_HANDLE)HandleDC,cacheIdx)))
return hr;

//Draw feature here

this->m_ipSDisp->FinishDrawing();

.............................................

Because ArcMap has at least 3 caches,1st for feature,2nd for annotation and graphics and 3rd for selection,i supposed the last feature index must be for selection cach,if I choose another index less than this,it doesn't even draw any thing.I checked this number and it always returns 2 for cachecount,I think because I don't have any graphics or annotation on map,it returns 2(at least I wish it is a true consumption).

but unfortunately it draws polyline on the screen because drawn polyline disapears at the first refresh of map.
Is there something wrong with these procedure? May you tell me how I can draw a sketch on the selected feature just like Edit Tool in ArcMap?


texts: Display the position of a feature/polygon borders while editing it

texts:
I have a feature displayed on the map which can be modified. And I wish to be able to see the exact coordinate of the feature when hovering on it with the little indicator. I tried using

ol.control.MousePosition

like this: https://jsfiddle.net/Almia/qpqnrsnb/

But it’s not precise enough. When you circle around the dot it gives a different coordinate (which I suppose it only gives the exact data of where the mouse is)

What allows me to get the exact coordinate of the feature/polygon borders when hovering on it with the mouse?


texts: QGIS: Merging Rasters (Digital Surface Models)

texts:
I got several Raster Layers as Digital Surface Models (DSM) in QGIS. I now want to merge them by using Raster -> Miscellaneous -> Merge, but if I do, the resulting merged raster layer is just blank space, as seen in the picture below. Does anyone know what's the problem here? Changing the data type doesn't d anything, I used float32. I also tried merging only two, also didn't work...
UPDATE I got it working by using Raster -> Miscellaneous -> Build Virtual Raster but I'd still like to know what's the problem here. Does the virtual raster have any disadvantages vs the merged one?


texts: Is there a way to calculate the geometric mean of a raster in ArcGIS?

texts:
For the parameters whose values are log-normally distributed, we should calculate the geometric mean rather than arithmetic mean. In the Focal Statistics, there is option to calculate the arithmetic mean but I am interested to calculate the the geometric mean of my raster file. Is there any way to calculate the geometric mean of a raster in ArcGIS?
I want to calculate geometric mean of the surrounding neighborhood. The same way focal statistics calculate the arithmetic mean with defined map unit.


texts: Merging information of two maps using ArcGIS Desktop?

texts:
I have a raster file with 68M pixels and a shapefile with 83 polygons covering the same area.

I would like to find, for each pixel, which polygon it is located. If there is more than one polygon per pixel, I would like to choose one based on the predominant area, but if it is very "machine expensive", it could be a random choice as well.

I have both layers in a PostGIS DB, and I tried to do it using PostGIS ST_Intersects but it will take ages to do that, even in my very fast workstation and using geo-indexes (I tested with a small part of the raster data).

How can I do this using ArcGIS Desktop or python?


texts: Why the WGS84 in Trimble Business Center software is different than real?

texts:
Why is the WGS84 in Trimble Business Center software different than another software than Leica Infinity or the public WGS84?

In the picture, the flattening(1/f) is 298.25722293287 in TBCv3.5 software. But the public WGS84, the flattening(1/f) is 298.257223563.

It will have error 1cm when switched from Lat, Long, Heigh to X, Y, Z when using the TBC software


texts:
I have two tables:

CREATE TABLE periphereies
(
  gid serial NOT NULL,
  per character varying(60),
  geom geometry(MultiPolygon),
  CONSTRAINT periphereies_pkey PRIMARY KEY (gid)
)
WITH (
  OIDS=FALSE
);

CREATE TABLE archaia_theatra
(
  gid serial NOT NULL,
  "ΘΕΑΤΡΟ" character varying(254),
  "X" numeric,
  "Y" numeric,
  geom geometry(Point),
  CONSTRAINT archaia_theatra_pkey PRIMARY KEY (gid)
)
WITH (
  OIDS=FALSE
);

First of all i am trying to write a query that return's me all the periphereies record's that contain at least one archaia_theatra.geom point record..
Then i want to select the periphereies which contains the maximum number of archaia_theatra (compared with all the other perifereies records)

I am trying this:

CREATE OR REPLACE VIEW per_max_theatra AS
SELECT count(p.per) FROM periphereies AS p JOIN archaia_theatra AS a ON
ST_Contains(p.geom, a.geom) ORDER BY count(p.per) DESC LIMIT 1

but it's not working

Any help?

Thx in advance


texts: Increasing a specific landcover class by a % in ArcGIS

texts:
Is it possible to increase a specific landcover class by a % figure in ArcGIS. For example, I am making a flood model in HEC-HMS and how rainfall events are affected by an increase in urban landcover, the gov projections for the major city in my catchment area states an increase in population by 13%.
Would there be a way to increase the urban landcover by 13%?


texts: Changing geometry types from kml points Z to just points in QGIS

texts: How do I change geometry types from pointZ to point in QGIS from a kmz file?


texts: Change shapefile field names using ogr AlterFieldDefn

texts:
I need to change field names of a shapefile to upper case using ogr. I want to do this automatically, versus calling out each field and replacing the name with upper case.  When I try to implement using layer AlterFieldlDefn, I get a TypeError. Here is the code I am using:

Shapefile = r"C:\test.shp"
driver = ogr.GetDriverByName('ESRI Shapefile')
dataSource = driver.Open(Shapefile, 1)

if dataSource is None:
    print 'Could not open %s' % (Shapefile)
else:
    print 'Opened %s' % (Shapefile)
    layer = dataSource.GetLayer()
    layer_defn = layer.GetLayerDefn()
    for i in range(layer_defn.GetFieldCount()):
        print layer_defn.GetFieldDefn(i).GetName()
        uppercase_name = layer_defn.GetFieldDefn(i).GetName().upper()
        layer.AlterFieldDefn(i, uppercase_name, 1)

The error I get is this:

  TypeError: in method 'Layer_AlterFieldDefn', argument 3 of type
  'OGRFieldDefnShadow

How do I correctly uses AlterFieldDefn to change the field names to upper case?  I do not have ArcGIS 10.2, so AlterField_management is not an option.


texts: What labeling settings do I need when I have a polygon layer with several areas that have touching borders. I want labels that follow the border of each object INSIDE itself even if that means that the label is upside down. This is a often used style seen in maps.


texts: Unable to transform feature classes having GCS_WGS_1984 to the current data frame layer (NAD_1983_UTM_Zone_10N)

texts: This is my first time posting here, so I apologize if my question seems too basic. I am designing a map of various locations that are of interest to my organization. I've obtained all of the layers from the City of Burnaby OpenData website. Upon adding them into ArcMap 10.3, I noticed that the feature classes were operating on the GCS_WGS_1984 datum. Given a city like Burnaby (British Columbia, Canada), I reasoned that the current datum distorts the map far too much. I thus changed the data frame to NAD_1983_UTM_Zone_10N, but the respective feature classes still operated on the GCS_WGS_1984 datum. Did I miss a few steps in between the transformation?


texts: Need to preserve file structure of project file?

texts: I understand that we have to keep the file structure unchanged for our data files (shapefiles, etc.), so our project files can find them. But I had no idea that a similar thing may apply to the locations of our QGIS project files. I recently "reorganized" some of my computer files (not including my data files), and found that upon opening a project, I had to tell it where to find all the data files. Is that correct, or am I missing something?


texts: Point layer with random size and rotation

texts:
I have a point layer representing trees. This layer uses a SVG symbol.
I set the size with rand(3, 5) in map unit (I'm in meters). And the rotation with rand(0, 359). With this, my map looks a bit more "natural"!
My problem is that the random is re-calculated each time I zoom in and out. I want that every new point uses the SVG with a size between 3 and 5 meters, and a rotation between 0 and 359. But once it's done, they don't move anymore.
Is it possible with the geometry generator? If yes, how?
The reason is if I send a map to my customers, and they ask me just to change the title for example, I don't want all my symbol to be also modified!


texts: Qgis action adding layers using relative paths

texts:
in QGIS 2 windows
I have some index layers I use for loading raster and vector tiles through an Python action. Currently I am using absolute paths which works fine but I would like to set it up so other teams can use the tiles. I have been trying to get the path of the active layer (the index layer) and use that but had no luck.

I would also like to be able to load the new layers in a specific group to keep projects tidy - can this be done?


texts: How to vectorize creating a Shapely Polygon in Pandas

texts:
I have a GeoPandas DataFrame with a polygon in each row. I want to add a column with coordinates of a bounding box of each polygon. I can do it this way:
def create_bbox(row):
    xmin, ymin, xmax, ymax = row.geometry.bounds
    return Polygon.from_bounds(xmin, ymin, xmax, ymax)

osm_buildings['bbox'] = osm_buildings.apply(lambda row: create_bbox(row), axis=1)

However, because of the size of the dataset, I need to speed this process up. I want to use vectorization. What I've tried is this:
osm_buildings['bbox'] = Polygon.from_bounds(
     osm_buildings.geometry.bounds.minx,
     osm_buildings.geometry.bounds.miny,
     osm_buildings.geometry.bounds.maxx,
     osm_buildings.geometry.bounds.maxy)

However, I get
*** ValueError: The truth value of a Series is ambiguous. Use a.empty, a.bool(), a.item(),
a.any() or a.all().

I don't understand very well what "truth value" is the error referring to and how can I fix it.
Question: How can I make this operation faster than using apply()? I am asking specifically about vectorization but if there is something else to speed it up I'm interested as well.


texts: Is Corine Land Cover 2018 available from WFS server?

texts:
I need to work with CLC 2018 across Europe from mobile devices with little storage capacity. It would be ideal if I could access it through a WFS server and extract just the areas and/or classes I need. I have found a few national websites with part of the dataset as WFS service but not the entire European dataset.
How can I do it?


texts: Checking for overlapping geometries in one table using PostGIS?

texts:
I have a table with geometries where some geometries overlaps with other geometries in the same table.

I already found an answer for this question here:
Check whether table has overlapping polygons, in PostGIS?

But when I try to create a view of the query as result I do not get any geometries for the query only a table, this also helps a lot but with geometries it would help even more.

The query is

SELECT *
FROM schema.table1 a
INNER JOIN schema.table1 b ON
(a.geom && b.geom AND ST_Overlaps(a.geom, b.geom))
WHERE a.id != b.id;

With this query I get the table and if I substitute st_overlaps with st_intersection  I get the query error that the intersection must be type boolean and not type geometry.


texts: Adding a table nulls all values in one column

texts:
I've read (and had to deal with) values being changed to "NULL" when joining data from a table to a shapefile. Now I'm having an even more basic problem: for some tables, when I add the table, one of the number/integer columns fills with "Null". Here's what I've done so far:

My csv file has three columns, one is the County FIPS Code and the other two contain some numeric demographic information (average age and average income).
I created a .csvt file that states: "String","Integer","Integer"
Both files are saved in the same directory.
I've tried adding the table as a vector file (no geometry), as well as a text delimited file. In both cases, the second integer column fills with "Null".
I've closed an re-opened QGIS several times between each attempt.

Any suggestions for what I can try next?

Image of attribute table here.

Image of original data table in Excel:


texts:
I have a polygon feature which returns Bad envelope when run through the Check Geometry tool.

According to the help:

  Bad envelope : The envelope does not match the coordinate extent of
  the geometry.

I'm not sure what this means. I tried googling feature envelopes but couldn't find anything discussing.

I ran a Feature Envelope To Polygon on the bad feature and The polygon created fit the polygon perfectly:

I checked the original polygon's extent min/max x/ys and they match the envelope polygon exactly.

>>> print inExt.YMin
13693868.258
>>> print envExt.YMin
13693868.258
>>> print inExt.YMax
13694050.7098
>>> print envExt.YMax
13694050.7098
>>> print inExt.XMin
1814653.01959
>>> print envExt.XMin
1814653.01959
>>> print inExt.XMax
1815376.21679
>>> print envExt.XMax
1815376.21679

Question: What does it mean for a feature to have a bad envelope?

Note: I ran a Repair Geometry on the feature. The tool indicates the envelope has been fixed, but after this is done a Check Geometry on the same feature returns Bad envelope still.


texts: Rasterizing shapefile with attribute value as pixel value labels with GDAL in Python

texts:
I am new to GDAL and geom shapefiles.
I have a 3-band .img file (with other associated files) as raster images.
And also a  shapefile which includes 10000+ features (multi-polygons), the properties of each includes 'shape_leng', and 'BM' (corresponding to the geographic class of this area, e.g. water/forest).

I want to rasterize the shapefile with property 'BM' as labeled pixel value (single/three channels) to put it together with my raster image so that I can visualize the segmentation or further train an image segmentation network using this as a training dataset.
How to implement it using Python? Or is there any tutorial/notebook that I can learn from?


texts: CartoCSS detect odd vs even for use in labeling conditional

texts:
I liked the part in the advanced label styling about randomly asigning a 0 or 1 to control label placement and then using as a conditional

#labels {
  text-name: "[name]";
  text-face-name: "OpenSans Regular";
  text-placement-type: simple;
  text-placements: "E,NE,SE,W,NW,SW";
  [dir=1] { text-placements: "W,NW,SW,E,NE,SE"; }
}

Rather than create a new variable, could I just use my FID and detect whether or not it is odd/even? Is there a way to detect this property in the CartoCSS?


texts: Grouping overlapping points by UTC Time, sum and average?

texts:
I have data from airplane calibrations where the plane flies in a crossing pattern; i.e. going E->W and W->E, then N->S and S->N.

I need to sum and average the values for both flights in each direction (i.e. when the plane goes E->W, and then W->E) and get the absolute difference.

I'm able to assign 'NS' or 'EW' using SAGA's Separate Points By Direction algorithm and then using a field calculator. However, I don't know how to differentiate the E->W flight from the W->E flight. Any ideas?

The only idea I have is that there is a UTCTime attribute; so we can group EW data into two separate groups by time somehow, but the flights are always at different hours of the day on successive surveys.

I'd like to be able to do this using a graphical modeller so no 3rd party plugins if possible.


texts: Rasterize using irregular points as centroids for each cell

texts:
I have a dataframe (xyz), which I loaded in ArcGIS using add xyz-files, saved it as a new shapefile. The distance between the points is irregular. I want to rasterize the pointshape with each point as centroid in the equivalent cell.
I thought the conversion-tool: point to raster would do this automatically. But it does not give me the result I was looking for.
Example result:
I defined the cell size as: 9950,9950
CRS is: WGS_1984_UTM_32N
ArcGIS version 10.6