texts: Why does it matter which geometry I transform in a spatial query?

texts:
I am using PostGIS and have two tables. One containing points of start-locations in EPSG:4326 (about 3 million records in total) and another one containing county-polygons in UTM EPSG:25832 (about 30 records). I wanted to calculate the amount of start points per county. So I used this query:
SELECT
count(f.ride_uid) as amount,
l.name,
l.u_key
FROM schema_a.rides as f

LEFT JOIN
schema_b.counties as l
ON ST_Intersects(f.start_geom,ST_Transform(l.geometry,4326)) -- f.start_geom are points in 4326 and l.geometry are polygons in 25832 (UTM)

WHERE
f.start_date_year = 2021
and f.start_date_month in (1,2,3)

GROUP BY
l.name,
l.u_key

But it ran endlessly and I cancelled it after 2 hours, wondering whats wrong. Tried it again after checking the code a second time, but again no result after 15 minutes. Then trying
SELECT
count(f.ride_uid) as amount,
l.name,
l.u_key
FROM schema_a.rides as f

LEFT JOIN
schema_b.counties as l
ON ST_Intersects(ST_Transform(f.start_geom,25832),l.geometry) -- f.start_geom are points in 4326 and l.geometry are polygons in 25832 (UTM)

WHERE
f.start_date_year = 2021
and f.start_date_month in (1,2,3)

GROUP BY
l.name,
l.u_key

gives the correct result in about 1 second.
Why does it matter which geometry I transform? Why does ST_Intersects(ST_Transform(f.start_geom,25832),l.geometry) work so fast and ST_Intersects(f.start_geom,ST_Transform(l.geometry,4326)), well, not at all?


texts: How to create a Winkel-Tripel map of countries with different colors?

texts:
I need a good quality rasterized world map, Winkel-Tripel projection, in which all countries have different (random?) colors. No borders, text, grid or anything else on the map.

Is GIS the right way to go? I found a lot of map data on naturalearthdata.com. Can I create this map with QGIS? Do I need ArcGIS?

I'm completely new to GIS; it looks like an exciting field. I'd appreciate any help to get started.


texts: How can I add data from Points to Polygons (intersected or near) in ArcGIS 10

texts:
I'm using ArcGIS 10. I have a Point and a Polygon shapefile. The Points have a NAME field filled, and the Polygon the same but emptied a minority (350 out of 3000). I need to catch these names from the points only to polygons which have no data in this Field.

Is there any tool, or anything from ArcGIS which offer this solution?
Screenshot att:

Thanks!!


texts: Exporting data from ArcSDE gdb atribute table to file or personal geodatabase?

texts:
I'm currently learning GIS and I'm working on a project, I'm doing my own research and trying to accomplish the following with or without Python:

Export data from an ArcSDE geodatabase attribute table to a File or Personal Geodatabase OR
Export data from an ArcSDE geodatabase attribute table to a dbf table or excel table, then from here I can export it to a personal geodatabase .mdb


texts: Prevent ModelBuilder Tools From Running if Published to ArcGIS Server

texts:
I'm looking for a way to prevent an ArcGIS ModelBuilder tool from running if it's been published to ArcGIS Server. Or, even better, for a way to prevent the tool from being published to Server in the first place.

Ideally, I'd like to find a way to keep the tools from running without a license. This has been discussed in Copy protection via license file for ModelBuilder models? although no specific code is given. If licensing isn't possible, I'd like to keep these tools from being published to Server. Does anyone have ideas for Python coding (or a method directly in ModelBuilder) that could achieve this?


texts: Calculate walking time on existing paths (hiking) in ArcGIS

texts:
I would like to calculate the walking time on existing paths in ArcGIS (hiking trails).

I have a dem of my area (30px cell size) and I consider that the slope is the main factor which affect the walking speed. When the path is uphill, the speed is slower so the walking time is higher and when the path is downhill the speed is higher and the walking time is lower.

I have read about Tobler's hiking function which calculates a vertical factor table ( first_column slope_degree, second_colum time(in hours) to cross 1 meter)

TIME (HOURS) TO CROSS 1 METER =0.000166666*(EXP(3.5*(ABS(TAN(RADIANS(slope_deg))+0.05))))

but I can not understand how to use to in my case.
Is there any idea?


texts: Finding the children name of menubar in QGIS Python

texts:
I have try the code in Finding name of QGIS toolbar in Python? and Name of children of toolbar in QGIS python and how to disable? to show the children name of menuBar. But, get error. The children name of menuBar are : Project, Edit, View, Layer, Settings, Plugins, Vector, Raster, Database, Web, Processing, and Help.

I have input the code :

from PyQt4.QtGui import QToolBar, QDockWidget, QMenuBar
for x in iface.mainWindow().findChildren(QMenuBar):
    print x.objectName()

until here, it is work. The output is menubar. When I type the code :

for icon in iface.menubar().actions():
    print icon.objectName()

It gives error : AttributeError: 'QgisInterface' object has no attribute 'menubar'.

To test the name is true, I want to remove many of children name. For example, I want to remove Edit, View, Layer, and Settings. So, how the code in python?


texts:
I have installed pyMODIS along with Numpy and the necassary GDAL bindings.

How do I go about running the scripts such as modis_convert.py and modis_download_from_list.py.

Are these ran from the command prompt?

I have also attempted to run them from the Python IDLE but have got Syntax errors.


texts:
I have a satellite image as a PNG with known coordinates (bottom, left, top, right). I want to draw a shapefile on it (the limits of a county) in R.

My first thought was to convert the PNG (it's read as a three-layered matrix, for red, green and blue) to a raster. But creating a three-layered raster may be an overkill.

And simpler solution to do this?

I found the function rasterImage(), and I'm doing like that:

Map0 <- readPNG('myMap0.png')
cM <- c(xmin,ymin,xmax,ymax)
png('myMap.png',10240,8610)
par(mar=c(0,0,0,0))
plot(myShapeFile,border=NA,xlim=c(cM[1],cM[3]),ylim=c(cM[2],cM[4]))
rasterImage(Map0,cM[1],cM[2],cM[3],cM[4])
plot(myShapeFile,add=T,border='white',lwd=5)
dev.off()

It's almost fine, except for a white border around the image. Since I used the dimensions of the original image (10240x8610), this border is probably distorting the new image, because the painted area is smaller than the original. I had this problem before, but could never solve it.

I thought the problem was with png(), so I did the following:

par(mar=c(0,0,0,0),bg=NA)
plot(myShapeFile,border='white',xlim=c(cM[1],cM[3]),ylim=c(cM[2],cM[4]),lwd=5)
dev.print(width=10240,height=8610,'myBorder.png',dev=png)

So I could try to merge the two images as layers in Gimp, but I see that the border remains, causing the shapefile to be drawn in the wrong position (it aligns perfectly with myMap.png, but not with myMap0.png). In Google Earth I know that the coordinates in cM are precise, so, since I'm giving mar=0, and xlim and ylim should refer to the limits of the plot, why this margin keeps on disturbing the image created?

Following @mdsumner comment, I tried the following:

map0 <- brick('myMap0.png')
extent(map0) <- c(cM[1],cM[3],cM[2],cM[4])
png('myMap.png',10240,8610)
par(mar=c(0,0,0,0))
plotRGB(map0,maxpixels=10240*8610)
plot(myShapeFile,add=T,border='white',lwd=5)
dev.off()

Now it's almost right, except for the left and right white margins, each 27 pixels wide).


texts: Using GDAL's Rasterize (Vector to Raster) utility?

texts:
I am trying to convert a point layer to a raster file using QGIS python consle. I want the output raster with the extent of a raster file. The new raster should have the point values from the specified field of the input vector and a value of 0 elsewhere.

I am using GDAL's Rasterize (Vector to Raster) utility (https://docs.qgis.org/2.8/en/docs/user_manual/processing_algs/gdalogr/gdal_conversion/rasterize.html).
This documentation is for QGIS 2.8 version while I am running my script QGIS v. 3.6.3 and 3.4.8.

I am getting error:

'Incorrect parameter value for UNITS' for the last line of the code, that is for the output.

I have seen several existing posts on QGIS but I still couldn't solve this issue. Like for instance, when I use the script given here: Rasterizing shapefiles with GDAL and Python? I get an output raster with a 0 value for all pixels.

I also tried by giving the pixel size instead of extent by following this post - extent parameter pyqgis gdal processing rasterize. But then too I get the same error of 'Incorrect parameter value for UNITS'.

Using the code given in another post Python, creating layers and rasterizing polygons in GDAL
I get the error:

'AttributeError: 'NoneType' object has no attribute 'SetGeoTransform'.

I checked documentation for SetGeoTransform and I failed to figure out what went missing. With different codes I am getting different error.

I am attaching the code that I have used.

I just do not seem to understand why the 'Units' error is occurring. Both input data have same projection.

import os
from osgeo import gdal, ogr
import processing
#input = iface.activeLayer()
input = 'C:\Trials\points.shp'
#layer = QgsVectorLayer(input.source(),"point","ogr")
layer = QgsVectorLayer(input,"point", "ogr")

#Raster extent for the new layer
St_Area_Rst = 'C:\Trials\Extent.tif'
#collecting dimensions of the above raster
dim = gdal.Open(St_Area_Rst)
xoff,a,b,yoff,c,d = dim.GetGeoTransform()
band1= dim.GetRasterBand(1)
rows = dim.RasterYSize
cols = dim.RasterXSize
Total_pixel = (rows*cols)
#specifying extent of input layer
extent = layer.extent()
xmin = xoff
xmax = xoff + (a*cols)
ymin = yoff + (d*rows)
ymax = yoff
output = 'C:\Trials\points_rst.tif'
processing.run("gdal:rasterize",
               {"INPUT":layer,
               "FIELD":"Rstrz",
               "DIMENSIONS":0,
               "WIDTH":cols,
               "HEIGHT":rows,
               "RAST_EXT":"%f,%f,%f,%f"% (xmin, xmax, ymin, ymax),
               "TFW":1,
               "RTYPE":4,
               "NO_DATA":0,
               "COMPRESS":0,
               "JPEGCOMPRESSION":1,
               "ZLEVEL":1,
               "PREDICTOR":1,
               "TILED":False,
               "BIGTIFF":2,
               "EXTRA": '',
               "OUTPUT":output})


texts: Adding row number in QGIS 3 using Python

texts:
I want to add a row number in QGIS 3 using Python.
First of all I add a new field with the code below and it works.
Then I want to fill the table with the row number. In QGIS Field Calculator there is the function @row_number. That is exactly what I want. But I want to do it automatically with Python.
The code below doesn't work.
## join ID field
vlayer=##
myField1 = QgsField('ID', QVariant.Double)
vlayer.startEditing()
vlayer.dataProvider().addAttributes([myField1])
vlayer.updateFields()
idx=vlayer.dataProvider().fieldNameIndex('ID')
e = QgsExpression('$rownum')

##fill the attributes with ID-Number
for f in vlayer.getFeatures():
    f[idx]=e
    vlayer.updateFeature(f)

vlayer.commitChanges()

Is there an easy way to add a row number or an "ID" with Python?


texts: Creating a plugin for adding and editing attribute data through a GUI

texts:
I'm wondering if it's possible to create a plugin using python that uses a GUI to edit attribute data instead of the normal attribute table. I want to make a plugin that makes it easier for the less computerly abled to add point data to an already existing point layer, input basic information about that point (preferably through roll-down choices etc.) and then edit the point afterwards.

I'm completely new to python as well as QGIS but I want to try do as much as possible myself, but I'm not sure if QGIS allows for adding your own GUI etc.

What is possible and what is not?


texts: Access layer of raster stored in memory using ArcPy to create NDVI from subset of multi-band raster

texts:
I am trying to create a final output of an NDVI raster for a number of Landsat images run through a for loop in Python 3 (for ArcGIS Pro). I am able to produce the NDVI rasters when I save the intermediary Landsat subset from ExtractByMask, but receive the following error when I instead try to store the intermediary subset in memory:

nir_raster = arcpy.Raster(nir)

  RuntimeError:  ERROR 000732: Input Raster: Dataset in_memory\myMemoryRaster\Layer_4 does not exist or is not supported

Here is the code where I attempt to refer to layers of the raster in memory (in my case layers 4 and 3 refer to NIR and red bands, respectively in my particular multi-band raster).

import arcpy, os

arcpy.CheckOutExtension('spatial')
arcpy.env.overwriteOutput = True
arcpy.env.workspace = 'in_memory'

# Define area of interest for subset
aoi = 'P:/aoi/aoi.shp'
aoi_name = os.path.basename(aoi)

# Select Landsat Images
images= "P:/Landsat/L7_p39r37_01242016_b2-5_7.img; P:/Landsat/L8_p39r37_02172016_b2-7.img"
image_list = images.split(";") # used here because real script takes multi-value parameter for images

# Select output folder to which NDVI rasters will be saved to
output_workspace = 'P:/NDVI'

# 1. Generate NDVI raster
for i in image_list:

    image_name = os.path.basename(i)
    image_name_chunks = image_name.split('_')
    ndvi_name = 'ndvi_' + image_name_chunks[2] +'.img'
    ndvi = os.path.join(output_workspace, ndvi_name)

    # 1.1 Subset image
    outExtractByMask = arcpy.sa.ExtractByMask(in_raster = i, in_mask_data = aoi)
    output = "in_memory" + "\\" + "myMemoryRaster"
    outExtractByMask.save(output)

    # 1.2 Calculate NDVI on subset
    nir = output + '\Layer_4'
    nir_raster = arcpy.Raster(nir)

    red = output + '\Layer_3'
    red_raster = arcpy.Raster(red)

    numerator = arcpy.sa.Float(nir_raster - red_raster)
    denominator = arcpy.sa.Float(nir_raster + red_raster)

    ndvi_output = arcpy.sa.Divide(numerator, denominator)
    ndvi_output.save(ndvi)

arcpy.Delete_management('in_memory')

However, when I substitute the following code in sections 1.1 and 1.2 (i.e. not saving subset to memory) my code runs fine:

# 1.1 Subset image
subset_name = image_name[:-4] + '_' + aoi_name[:-8] + '_subset.img'

subset = os.path.join(output_workspace, subset_name)

subset_output = arcpy.sa.ExtractByMask(in_raster = i, in_mask_data = aoi)
subset_output.save(subset)

# 1.2 Calculate NDVI on subset

nir = subset + '\Layer_4'
nir_raster = arcpy.Raster(nir)
red = subset + '\Layer_3'
red_raster = arcpy.Raster(red)

I've visited: Access individual bands and use them in map algebra
and
What is the proper syntax and usage for ArcGIS "in_memory" workspace?

Update: adding multi-line version of user2856's Describe note below:

for child in arcpy.da.Describe(outExtractByMask)['children']:
    print(child['catalogPath'])


texts: How to directly add new fields at the right position with PyQGIS?

texts:
I have an existing Attribute Table.
I add new fields with PyQGIS but I would like to set this fields at the second and third position and, not at the end. Because by default, the new field is added at the last field position.

Example: I have fields A, B, C, D, E, F, G and I want to add the x and y at second and third position: A, x, y, B, C, D, E, F, G

Here is my pyqgis script to create new field:

for champ in champs :
    provider.addAttributes([QgsField(champ, QVariant.String)])
    inlayer.updateFields()
    indice = inlayer.fieldNameIndex(champ)

I know I can do it with refactor fields from processing but I would like to do it with PyQGIS.


texts: How to use QGIS Atlas with multiple layers/files per coverage?

texts:
I've got 30 contour datasets that all sit over the same extent.  Each contour set is in its own file. I also have 30 extent polygons for each contour set as a separate file.

Can I use Atlas to generate a series of maps - each displaying a 'pair' of datasets - ie. contours and extent for each pair over the same extent?

Everything I've read indicates that you can loop through features in a layer with Atlas, but not through multiple layers. Could I work around this by adding a 'layer display list' as an attribute of the coverage layer?

Maybe with a python script?

Update:

This is going back to a previous question from a few years ago that I asked and didn't find a solution for.  QGIS has advanced a few versions since then - so I'm hoping there is either an updated solution, or someone could provide guidance to develop a script.

Advanced warning - I'm a geologist - layers of rock are how I think! :-)

I have a project area for which I want to produce a series of maps (an atlas).  Each map will cover exactly the same extent, but will show properties for a different layer of rock.  In this example - contours of the depth to the top of 50 different rock layers.

The desired output would be 50 maps, same extent, different layer (or layer set) displayed

To add more complexity - I'd like to add 'sets' of layers.  Rock layer 1 also has unique shapefiles showing the extent polygon layer and also a label_point layer.

Current Method

I've tried merging all the contour layers into one shapefile and using this as the Atlas coverage layer. This method doesn't quite produce the desired result because:

I can't filter the combined layer to only show one layer (not one feature) at a time.
I can filter using this rule based style expression on the combined layer! (where unit is the rock layer name):

attribute($currentfeature , 'unit') =  attribute(@atlas_feature, 'unit')

Critical issue - I get one page per individual contour (ie. one per feature - so hundreds of pages per rock layer). I can't work out how to filter so I only get one page per rock layer (or if that is possible).
I can't merge the different 'sets' (points, polylines, polygons). Workaround - I could merge and filter based on geometry type - same as above.
Merging 50 input files into one master shapefile seems a bit redundant - surely one more code loop could negate the need for this? Workaround - use the processing framework to combine input files.


texts:
I am using gstat::krige to interpolate a data set of arund 2700 points. However I get an error message "Error in seq.default(zrng[1], zrng[2], length.out = cuts + 2):'from' must be a finite number" and I don't understand what does it means.
My confusion is that when I do cross validation in the same data set, it works fine. When I do kriging interpolation in less than 2000 points sampled from the same data set, it works fine. But when I do kriging interpolation in the entire data set, predicted values (var1.pred) and variance (var1.var) results is NA.

I can't localized what throws the error. Below you can see what I did and in the link are my files:https://www.dropbox.com/sh/ff2231r2qobg4ud/AAB8aYHlObhakLKdySvLZFsSa?dl=0.
In the same link you will find the output of this.

library(sp)
library(gstat)
t1 <- read.table("SelecPts.iter.txt")
g1 <- read.table("Dh9709.grid.txt")
coordinates(t1) <- ~ x+y
coordinates(g1) <- ~ x+y
gridded(g1) <- T
Var <- variogram(dh9709~x+y, t1, alpha=c(45,135))
m <- fit.variogram(Var, vgm("Exp"))

Cross-validation works

cv1<-krige.cv(dh9709 ~ x+y, t1, model = m, verbose=T))

Kriging with 2700 point doesn't work

kg1<-krige(dh9709 ~ x+y, locations = t1, newdata = g1, model = m))


texts:
I am using the GDAL csv driver and I can open a csv and specify geometry columns with ogrinfo using the -oo flags. However, when I want to specify these as part of another command-line tool, for example, gdal_rasterize, I can't seem to find the syntax to include them.

For instance, if I run

gdal_rasterize [--rasterize-options] myfile.csv --oo X_POSSIBLE_NAMES=xcol myoutputraster

I get an error on the --oo option.

Is the only way to do this using the command line to create a VRT file?


texts: gdallocationinfo used in javascript with node-gdal?

texts:
I cannot figure it out if it's possible to get pixel information from my rasters using gdallocationinfo in javascript. I found node-gdal, but after reading its documentation I'm not sure if gdallocationinfo is usable.

I basically need this:

gdallocationinfo mytiff.vrt lonlat.lon lonlat.lat -valonly -geoloc

So I added to my javascript the following as a test (after node-gdal was succesfully installed as explained in the link above):

var gdal = require("gdal");
var dataset = gdal.open("mytiff.vrt");
console.log("number of bands: " + dataset.bands.count());
console.log("width: " + dataset.rasterSize.x);
console.log("height: " + dataset.rasterSize.y);
console.log("geotransform: " + dataset.geoTransform);
console.log("srs: " + (ds.srs ? ds.srs.toWKT() : 'null'));

but I get:

Uncaught ReferenceError: require is not defined

It seems to me that node-gdal is not being recognized, do I need to add some kind of src path to my HTML file first to get node-gdal working?


texts:
Hi I had loaded a satellite imagery in PostgreSql using command
"raster2pgsql -s 4326 -I -C -M E:\20Mar\55M1_Geo.tif -F -t 100x100 public.55M1_Geo | psql -d rasterdata -p 5437 -U postgres" and committed succesfully.

Now I want to open it in QGIS3 using DB Manager and during add to map canvas, it gives warning error "WARNING   PG:  dbname=rasterdata host=localhost user=postgres port=5437 mode=2 schema=public column=rast table=55m1_geo is an invalid layer - not loaded"

Below the PostgreSql credential it shows "Can not open GDAL dataset PG: dbname=rasterdata host=localhost user=postgres password=pgadmin port=5437 mode=2 schema=public column=rast table=55m1_geo", "testlayer" Couldnt establish a database connection"

A python error was occur after click database credentials as below:

File "C:/PROGRA~1/QGIS3~1.0/apps/qgis/./python/plugins\db_manager\db_plugins__init__.py", line 76, in
                initDbPluginList()
              File "C:/PROGRA~1/QGIS3~1.0/apps/qgis/./python/plugins\db_manager\db_plugins__init__.py", line 47, in initDbPluginList
                exec(u"from .%s import plugin as mod" % name, globals())
              File "string", line 1, in
              File "C:/PROGRA~1/QGIS3~1.0/apps/qgis/./python\qgis\utils.py", line 664, in _import
                mod = _builtin_import(name, globals, locals, fromlist, level)
              File "", line 1009, in _handle_fromlist
              File "", line 205, in _call_with_frames_removed
              File "", line 961, in _find_and_load
              File "", line 950, in _find_and_load_unlocked
              File "", line 655, in _load_unlocked
              File "", line 674, in exec_module
              File "", line 781, in get_code
              File "", line 741, in source_to_code
              File "", line 205, in _call_with_frames_removed

After that I had tried with "iface.addRasterLayer("PG: dbname=rasterdata host=localhost user=postgres password=pgadmin port=5437 mode=2 schema=public column=rast table=55m1_geo", "testlayer")" but get the error "2018-05-23T08:48:11 CRITICAL    Invalid Layer : GDAL provider Cannot get GDAL raster band:
            Raster layer Provider is not valid (provider: gdal, URI: PG: dbname=rasterdata host=localhost user=postgres password=pgadmin port=5437 mode=2 schema=public column=rast table=55m1_geo"

The same way I had loaded SRTM DEM. It display perfectly there is no issue in that. What the problem with satellite images. Any one help.


texts: Displaying composite raster in colors other than RGB in ArcMap?

texts:
I have a composite raster file with three bands in ArcMap 10.1. I am trying to display them in colors different than standard RGB. For example, I want band one to be blue, band two be brown and band three be green.

How can I do this task?


texts: Making curved labels for polygons in QGIS?

texts:
I have been able to find answers on how to do this in older versions of QGIS, but can't seem to come across anything that deals with the latest installment.

I need a way to make labels curve to irregularly shaped polygons; in this case rivers.

I have been assuming that this option would be listed under "Placement" in the "Labels" options but have been unable to find it anywhere.

I know it was present up until QGIS 2.2 and think it must be included in 2.6.1. I'm just not seeing it.

I am aware of the Easy custom labeling plugin curved option, but that doesn't seem to work as it should.


texts: How to choose a derived attribute as a parameter in FME?

texts:
Is it possible to have an FME parameter read attributes from a feature class that has been defined in a previous attribute?

So, in other words, I when selecting the Prompt and Run option, the user chooses a feature class.  In the below example, it's a feature class from a geodatabase, but could just as easily be any format.  Once they choose the feature class, I want a dropdown to appear showing the list of attributes that are in that feature class.

This would be similar to the functionality in an ArcGIS tool that has a derived parameter.

I was wondering if a python parameter could read the attributes and then return a list for the user to choose.

Is this even possible in FME?


texts: Zoom to link button in Popup, killing the Map

texts:
I am using the traditional Popup of JavaScript API and it seems to not work properly in my application. When you click on the map, the Popup comes fine but when you click on the Zoom to link button, the map zooms in and then if you move the map, the map disappears.

It only returns when you refresh the app. Not sure what is causing it?


texts: Clipping several raster layers to the same exact size QGIS

texts:
I'm new to QGIS, sorry if this is a newbie question... I need to clip about 40 different raster layers from one vector layer (that has the boundaries of 4 neighboring countries) but I need them all to have the same exact size, which didn't happen when I tried clipping the rasters one by one... How can I do this?

I thought about using the first cut raster to clip all others but I can't clip a raster with another raster... or can I?

Thanks for any help


texts: Python Gdal ERROR 6: GEOS support not enabled

texts:
I installed GDAL from here (https://pypi.python.org/pypi/GDAL/1.9.1)

Before attempting the python GDAL I installed the prerequisite gdal libraries using Kyngchaos's description some time ago (http://www.kyngchaos.com/software/frameworks). I have verified that gdal is running on at least the minimum required version:

gdalinfo --version
GDAL 1.11.1, released 2014/09/24

However, when I attempt to start python and access the gdal libraries I get errors on simple operations:

from osgeo import ogr
point = """{"type":"Point","coordinates":[-90.02001032734194,35.127405562848075]}"""
geom = ogr.CreateGeometryFromJson(point)
geom.Centroid()
ERROR 6: GEOS support not enabled.
<osgeo.ogr.Geometry; proxy of <Swig Object of type 'OGRGeometryShadow *' at 0x10344d9f0> >

I have googled this error and looked at how it has been mentioned on this site but no solutions seem to pertain to my situation. The Kyngchaos GDAL/OGR installs are built with GEOS libraries, so I'm very confused as to why python GDAL/OGR is saying GEOS is not enabled.

Any ideas on how I can fix or go about debugging this (running on Mac OSX)?

Edit 1:

I have verified that gdal 1.11.1 is being loaded with osgeo by forcing this specific version to top of path.

import sys
oldpaths = sys.path
sys.path = ['/Library/Frameworks/GDAL.framework/Versions/1.11/Python']
sys.path.extend(oldpaths)

from osgeo import gdal
print gdal.__version__
1.11.1

however I still get ERROR 6: GEOS support not enabled. on trying to take the centroid of a geom


texts: Geofencing: Query if a POINT is within a POLYGON

texts:
I am trying to use the following query to search if a point is within a defined POLYGON (aka a geofence I want to define):

https://spatial.virtualearth.net/REST/v1/data/[DataSourceID]/[DataSourceName]/[EntityName]?SpatialFilter=intersects('POINT(9.874434 57.003766)')&$format=json&key=[BingMapsKey]

My Data Source Schema looks like this:

<FourthCoffeeShops>
   <EntityID>1000</EntityID>
   <AddressLine>1 Microsoft Way</AddressLine>
   <Locality>Redmond</Locality>
   <AdminDistrict>WA</AdminDistrict>
   <PostalCode>98052</PostalCode>
   <CountryRegion>United States</CountryRegion>
   <Phone>303-555-0188</Phone>
   <Manager>Alan Steiner</Manager>
   <Latitude>57.003766</Latitude>
   <Longitude>9.874434</Longitude>
   <Confidence>Medium</Confidence>
   <GeoFence>POLYGON((9.86445 57.13876,9.89266 57.13876, 9.89266 56.94234,9.86445, 56.94234,9.86445 57.13876))</Fence>
</FourthCoffeeShops>

Essentially what I am trying to do is pass in Latitude and Longitude values to the POINT parameter in my query, and check if that location is within my GeoFence entity(defined POLYGON) or not. Please help me with creating a correct Data Source schema to run this kind of query.

I am using Bing SDS to hold all of this information. I was trying to follow this tutorial http://blogs.msdn.com/b/bingdevcenter/archive/2014/04/03/geo-fencing-with-bing-spatial-data-services-and-azure-mobile-services.aspx
But cannot get my DataSource schema and query quite right.


texts:
I am having an heavy TIFF map around 2Go and I want to index it in ES in order to geolocate data of this file. With GDAL, I could polygonize only a part of it in GeoJSON, but seems that big polygons are not handled.
Do you know a better way to upload a map and its data into Elasticsearch?


texts:
The OpenLayers plugin for QGIS allows one to add background maps from various sources to the data frame. These sources include OSM, Google, Bing, MapQuest, and Apple. However, Nokia's HERE Maps, is not among them.

Does anyone know how to add this, or is there another option available to use HERE Maps as a background in QGIS?

As far as I know, there is no straightforward way to add HERE Maps as a background in ArcMap, but it is possible in FME Data Inspector, as long as you have a HERE Developer account.


texts: How to build a geographic database of GPS logs?

texts:
On our aerial wildlife surveys we typically generate dozens of GPX files, and usually convert most of them into shapefiles for display in a GIS application. GPX files are a terrific way of keeping all (or almost all) of your GPS data in one place – they are xml files that store tracklog, waypoint & route information in one convenient format. They store multiple tracks, even old saved tracks from a Garmin.

I'd like to store all this track data in a spatialite database, in parallel with other tables that have point observational data from the same survey.

ogr2ogr will let me convert the track to a spatialite database:

ogr2ogr -append -f "SQLite" -dsco SPATIALITE=yes # add to table with Spatialite filetype
  -dsco INIT_WITH_EPSG=yes -t_srs epsg:4326 # using WGS84
  SURVEYDATA.sqlite GPXFILE.gpx tracks -nln tracktable; # adding tracks to 'tracktable'

However, this imports each saved track as a line – losing the time / location information for each track point!

Anyone know of any clever ways of building a tracklog database that will preserve that information?

EDIT:

Tracklogs are not just a line - they are collections of sequential points, each with an elevation and time attribute. Each point belongs to a track segment, which in turn belong to a given track, and points within each track segment were taken without interruption and may be assumed to be linked. For example:

<trk><name>ACTIVE LOG</name>
<trkseg>
<trkpt lat="-13.471403" lon="31.382982"><ele>467.818725585938</ele><time>2009-09-09T03:53:38Z</time></trkpt>
<trkpt lat="-13.471403" lon="31.382982"><ele>468.780029296875</ele><time>2009-09-09T03:53:50Z</time></trkpt>
<trkpt lat="-13.471403" lon="31.382982"><ele>465.896118164063</ele><time>2009-09-09T03:54:43Z</time></trkpt>
</trkseg>

Thus, to save a tracklog you need to save this collection of points with their associated data.

One solution might be to load these tracks as points, with additional columns for track name, segment, elevation and time; I don't think ogr will convert the tracks as points, though.

Solution

Thanks to @scruss who pointed out gpx2spatialite, which I installed (OS X 10.9) with sudo easy_install gpx2spatialite. This set of tools is designed for the 'drawinglife' project (and thus has some 'city' information that is probably unnecessary for most of us); the main gpx2spatialite page doesn't explain much, but the drawinglife wiki is better. Still, lots of undocumented features - like removing duplicates!

Create a spatialite database for GPX files which will hold track segments, tracks, trackpoints and waypoints. It also has 'users' set up (which seems odd, since sqlite is not a user-access system, but was meant for drawinglife) designed to label the various tracks:

gpx2spatialite create_db MyGPXArchive.db

Load a whole folder of GPX files, including all subfolders, checking and removing duplicates along the way:

gpx2spatialite -s -d MyGPXArchive.db -u SL25 Selous2013/

Options:

-s= don't check if the trackpoints are in the city database (speeds it up a LOT)
-d MyGPXArchive.db = database to use;
-u SL25 = add the user SL25 (just a label I used for the set of tracks collected);
Selous2013 = the folder (with subfolders) that has the tracks I'm interested in.

Bugs:

If you have saved tracks on an older Garmin that discards time stamp info, it'll crash the import (remove those saved tracks using a utility like GPSU). If you have gpsbabel, filter the file first with something like:

gpsbabel -w -r -t -i gpx -f [INPUT.gpx] -x track,start=20000101 -o gpx -F [OUTPUT.gpx]


texts: Webmapping js library with best support to WFS-T

texts: I've tried to use WFS-T in ol3 but it was really a hard task, it seems that WFS-T is badly undocumented feature and it's much complicated to implement, and i'm looking for recommandations concerning the easiest manner to deal with WFS-T, and what js library to work with.


texts: Calculating volume of buildings using QGIS?

texts:
I am working with QGIS 2 or 3

I have DSM and DTM already so I can calculate DSM - DTM:

Plus, I have shapefile with buildings' footprints.

How to calculate volume of every single building using DSM-DTM layer and add this value in m^3 to feature table? (Like with Zonal stats?)


texts: Adding Attachments in ArcGIS 10.0 (problem - can’t find the instrument Enable Attachments)

texts:
I want to create attachments to shp themes in ArcGIS (doc\pdf files, web-links). In the Help of ArcGIS I found the description of this operation:

Adding attachments to features

Add Attachments Tools (Data Management)

But before attachments can be added using this tool, they must first be enabled using the Enable Attachments tool.

So, I have to find Enable Attachments tool in the Data Management Tools.
But I don't have this instrument in my ArcGIS Data Management Tools, so can’t add attachments (I tried to do it without use Enable Attachments tool, but in this case I just don't have "Open Attachment Manager button" which is needed for add attachment)).

Could you advise me, how I can find the instrument “Enable Attachments tool”?
Maybe I have to turn on something in the “properties”? Or this instrument has another name?


texts:
I have a Garmin eTrex 30 that I am using with QGIS 2.0.1 and Mac OS 10.9. I want to load paths from QGIS onto my Garmin so that I can view them on my Garmin map while in the field.

I have a shapefile with vector lines describing old logging roads, which I exported from QGIS as a .gpx file in UTF-8 format and WGS84 projection. I tried skipping attribute creation, also tried putting FORCE_GPX_TRACK=YES in the Layer OGR creating field. I put the .gpx files onto my Garmin and every time I try to view them they do not show up on the screen.

However they DO show up when I open Garmin BaseCamp on my computer. Any ideas what I should try?

I've looked here and here, and other places as well, but I can't seem to find anything that solves my issue.


texts:
I'm trying to add a new field to an existing shapefile and filling it afterwards. I started doing this with the field calculator but as it involves many consecutive steps I thought it might be better to do it with the Python console.
As I never used the Python console before my code below is based on this video:
https://www.youtube.com/watch?v=NgvdgjyY03s
# define path to shape-file to work with
uri = 'D:/myfilepath/Windzonen_Deutschland.shp'

#adds layer to project
layer = QgsVectorLayer(fn, 'Windzonen_Deutschland_test', 'ogr')
pv = layer.dataProvider()

# create Field
pv.addAttributes([QgsField('Windzonen', QVariant.Int)])

layer.updateFields()

# create expressions
expression1 = QgsExpression('''case when "NAME_0"='Germany' then '1' end''')

# excute expressions for created fields
context = QgsExpressionContext()
context.appendScopes(QgsExpressionContextUtils.globalProjectLayerScopes(layer))

with edit(layer):
    for f in layer.getFeatures():
        context.setFeature(f)
        f['Windzonen'] = expression1.evaluate(context)
        layer.updateFeature(f)

If I run the code from the console as a block I receive this error:
Traceback (most recent call last):
  File "C:\PROGRA~1\QGIS3~1.10\apps\Python37\lib\code.py", line 90, in runcode
    exec(code, self.locals)
  File "<input>", line 1, in <module>
  File "<string>", line 20, in <module>
  File "C:/PROGRA~1/QGIS3~1.10/apps/qgis-ltr/./python\qgis\core\additions\edit.py", line 38, in __enter__
    assert self.layer.startEditing()
AssertionError

If I run in line by line I get the following error for this line of code pv.addAttributes([QgsField('Windzonen', QVariant.Int)]):
False

I use QGIS 3.10.5.A Coruna
What am I doing wrong? I tried to add layer.startEditing() to my code but I did not solve the poblem and I think because of with edit(layer) this should not be necessary...
Any hints or links?


texts:
I have feature classes of property parcels that are polygons. My client maintains data about the parcels in a google doc spreadsheet (online) and updates the information regularly. I am required to update this information as it relates to the polygons in ArcMap 10.3. Right now I am doing a simple join to an Excel spreadsheet, which I update by copying the online data into. The problem is that I really need to do a one to many since there needs to be more than one record pertaining to each individual polygon.

I can convert the Excel file to a geodatabase table to then do a relate, but then how do I easily update the gdb table?
Right now the only way I can see to do this is update the Excel spreadsheet each week, then re-convert the Excel file to a gdb table, overwriting the previous week, and use a relate within the map doc only; given that a gdb relationship class would be voided if I write over the data table.

The other caveat is that I then export the polygons with the joined data to a KML that is used in meeting. I found that the data in a "relate" does not export but the data in a "join" does export to the KML.

Am I understanding all of the ramifications? Could I use some other method of table (i.e. access table) whereby I can copy/paste new data over old data and not lose links? It almost seems that my only option is to make duplicate polygons and add a suffix -1, -2 for the multiple records. (they are linked by the parcel number)


texts: Get the list of layers in a store from geoserver

texts:
Is there a way to get the all the layers present in a store in geoserver?

My JS file contains code for each layer present in a store. So when I publish a new layer in geoserver I have to manually include the code of the new layer in JS. But what I want is to put a for loop in JS with the layer names present in a store, like

for(i<=all layers present in store)
{
   var wmslayer2 = new OpenLayers.Layer.WMS("Layer2","http://localhost:8080/geoserver/wms",
   {'layers': '**cite:layer[i]**', 'format':'image/png', 'transparent':'true'},
   {'opacity': 1.0, 'isBaseLayer': false, 'visibility': true});
}


texts: Raster Analysis - DEM and Population data

texts:
I'm trying to undertake some analysis on several raster data sets but am having some trouble. I have a DEM of my study area and a raster of population data for my study area. I am now trying to calculate the population numbers below a certain elevation using the elevation data derived from the DEM.

I know I have to use the raster calculator but am not savvy whatsoever with the use of python commands. Are there any sources where i can learn to better use python for this purpose or direct information?


texts: Reclassification Issues - Conditional Statements to Change Raster Values

texts:
I am looking for help with the statement to put into Raster Calculator that would turn some incorrect raster values (discrete values from 1385-1399) into a "1" without changing the majority of raster cells with correctly assigned values (1-14).

Background: When I run the Reclassify tool (Spatial Analyst) on my DEM in ArcGIS, which has an elevation range of 1385-2629 masl, it never works perfectly. The first class is supposed to be <1400, but instead I get classes for 1385, 1386, 1387, 1388, 1389, 1 (1390-1399), 2 (1400-1499) etc. I work with large, high rez datasets (100 km2 area, 1m cell size) and reclassification issues always seem to arise where at least one class is created containing only the minimum or maximum value(s) regardless. I read that this is a known issue with earlier versions of Arc but I am using 10.4.


texts: Reverse Geocoding - calculate every time, or store result in database?

texts:
Several questions, such as this one, discuss whether to store only lat/long in the database & calculate the street address every time it is needed, but they seem to miss a few points & seem to be opinion, rather than authoritative.

Hence this question, seeking a final, definitive clarification, with justification.

1) reverse geocode every time
pro - smaller database size (why should everyone duplicate Google's database?)
cons - might occasionally fail whereas storing the result will not

2) store address in the database
pro - less likely to hit a Google imposed limit on the number of requests
con - could be a large amount of data. A few hundred vehicles storing a new location every 5 minutes or so, 16 hours per day... probably more coding

If I take approach 2), I plan a new MySql table for addresses, indexed by lat/long and storing the street address. Every time a vehicle reports position, look for a matching address. If none, calculate & add.

In addition: later, when a client requests movement history, if the address is blank, attempt to calculate & store if found.

I am restricted to using a MySQL database.


texts:
I have a spring boot application which is making request to GeoServer to get a layer as a vector tile instead of the front end application making the request directly. The problem is I am not getting the response as a .pbf file instead I am getting it as a file. My code is below:

#controller

@ResponseBody
@GetMapping(value = "/tile-layer",produces= {"application/x-protobuf"})
public Object getWmsVectorLayerMap(@RequestParam(name = "SRS", required = true) String crs,
        @RequestParam(name = "LAYERS", required = true) String layers) {
    try {
    //  logger.info("making wms request...");
        String layer = Utilities.getLayers(layers).replace(" ", "%20");
        String url = wmsTileLayerUrl + "layers=" + layer +"&srs=" + crs;
        HttpHeaders header = new HttpHeaders();
        String auth = geoServerUsername + ":" + geoServerPassword;
        byte[] authentication = auth.getBytes();
        byte[] base64Authentication = Base64Utils.encode(authentication);
        String baseCredential = new String(base64Authentication);
        header.add(HttpHeaders.AUTHORIZATION, "Basic " + baseCredential);
        header.setAccept(Collections.singletonList(MediaType.ALL));
        header.setContentType(MediaType.parseMediaType("application/x-protobuf"));
        HttpEntity<String> request = new HttpEntity<String>(header);
        //restTemplate.getForEntity(url, Object.class);
        ResponseEntity<Object> rsp = restTemplate.exchange(new URI(url), HttpMethod.GET, request, Object.class);
        return rsp.getBody();
    } catch (Exception e) {
        logger.error("Error while making wms request because: "+e.getMessage());
        e.printStackTrace();
        return null;
    }
}

Is there anything I am doing wrong? How do I solve it please?


texts: Similar to layer stack feature in ERDAS what functionality is available in GRASS. I am aware of r.composite, but it can accept only 3 layers. To obtain a stack of more than 3 layers what can be done?


texts: How to create a buffer for catchment area of an historical settlement

texts:
I'm looking for a way to create a buffer from a specific point that represents the catchment area in a selected time in QGIS. I've found many tools that allow me to do it using driving time or public transport. I need to do the same only using walking time in a natural environment using the geophysical information of a DTM as criteria.
More in detail, I'm trying to reconstruct the catchment area of some different roman period settlements and to do that I'd like to create some isochrones of different walking times.
I've searched all over the web but I have not been able to find a solution to this specific problem. I've spent time trying to use SAGA algorithm of Accumulated costs but somehow I'm not able to run this tool.
I hope someone could help me or suggest me a proper tutorial.


texts: Need a simple python program to convert EPSG/SRID 102002 (LCC) to WGS84

texts:
I'm going to convert some coordinates from Canada LCC (EPSG/SRID 102002) into WGS84 (EPSG/SRID 4326).
I found a python package for this purpose which is pyproj. However, seems it doesn't support 102002.
since I get this result:
Proj("+init=EPSG:102002")

CRSError: Invalid projection: +init=epsg:102002 +type=crs: (Internal Proj Error: proj_create: crs not found)

Transformer.from_crs(102002, 4326)

Invalid projection: epsg:102002: (Internal Proj Error: proj_create: crs not found)

I know there is a way to define a new coordinate with pyproj but it needs expertise in GIS which I don't have. Building a Coordinate Reference System in pyproj
Would you please help me to create that coordinate reference system or give me a simple python code to convert those coordinate to each other?
I found this but it doesn't work fine:
>>> import pyproj
>>> p = pyproj.Proj("+proj=lcc +lat_1=50 +lat_2=70 +lat_0=40 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs")
>>> p(7678236.84856942, 1280719.40571689,inverse=True)
-27.952924124739766 15.502198604012452

but I expect a pair of number near (45.4330538,-75.722731)


texts:
I'm working on an application that renders GeoJSON into OpenLayers map.

The idea is about not having PNG images downloaded from external services (OSM or Bing) but having geo-information in GeoJSON file and rendering them in real-time. From my research this technology seems to be called 'vector tile'. I had a look at some of the guys' websites such as
http://openstreetmap.us/~migurski/vector-datasource/
and
http://mike.teczno.com/notes/postgreslessness-mapnik-vectiles.html (same guy).

His work seems to be what I want to achieve but I do not really get how it works.

From my understanding, since I use OpenLayers, a client retrieves map data in GeoJSON format and transfers them into vectors using OpenLayers' GeoJSON parser. But it's just the terms used in this technology that I do not get.

In his blog he described that there will be a server that provides data that is called 'MapNik Vector Tile' format. What's the difference between vector tile and GeoJSON? does it simply mean tiled GeoJSON that corresponds to the bounding box? Or is it a new image type?

If you have a look at his example, you can observe the labels on roads are aligned with the strokes of the roads. I tried to do the same on OpenLayers by changing Renderer code in OpenLayers but I'm only able to rotate labels but not make them bend them along the roads (E.g., when roads have corners, map displays labels in 'L' shape. ** + very well rendered!

It looks like Mike (the author of the blog) uses some different map besides OpenLayers but I cannot figure out how he did it. There are so many questions to ask but I'd like to understand the above first.


texts: Dissolve large dataset in ArcGIS and QGIS

texts: I have more than 80k point data. I have created few buffer from 1 to 5 km of them. Now I want to get the covered area by 1 km to 5 km buffered range. For this I have tried to dissolve those individual buffered file. In ArcGIS 10 it took 6 hours and then replied an error. Also in progress bar it was showing 33% completed. I think this 33% progress is some kind of bug as it stay 6 hours. In QGIS it took 8 hours to complete 61% then stopped. Is there any other way or tools to dissolve those large dataset in a reasonable time?


texts: Using QGIS coordinate transformation in user code

texts:
I am now working with coordinate transformations between different CRSs, e.g. latlong to UTM.

After comparisons, I found out that the coordinate capture plugin in QGIS is quite accurate for the results.

Is there a way I can use the methods QGIS uses to do transformation in my own project?

I have checked that QGIS is using the proj4 library to do the transformation. I have implemented the proj4 transformation, but got slightly different results.

pj_svy21 = pj_init_plus("+proj=tmerc +lat_0=1.366666666666667 +lon_0=103.8333333333333 +k=1 +x_0=28001.642 +y_0=38744.572 +ellps=WGS84 +units=m +no_defs")

pj_latlong = pj_init_plus("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")

double test1 = 37621.665;
double test2 = 12056.915;

projResult = pj_transform(pj_svy21, pj_latlong, 1, 1, &test1, &test2, NULL);

My result (latlong):
1.1253115030392806 103.9197679941208889

QGIS results:
1.35650722603 103.690059679 (check on map more accurate)

Can anyone help me on this?

Need to include iostream and proj_api.h:

int main ()

{

projPJ pj_svy21, pj_latlong

pj_svy21 = pj_init_plus("+proj=tmerc +lat_0=1.366666666666667 +lon_0=103.8333333333333 +k=1 +x_0=28001.642 +y_0=38744.572 +ellps=WGS84 +units=m +no_defs")

pj_latlong = pj_init_plus("+proj=latlong +ellps=WGS84 +datum=WGS84 +no_defs")

double test1 = 12056.915;

double test2 = 37621.665;

projResult = pj_transform(pj_svy21, pj_latlong, 1, 0, &test1, &test2, NULL );

std::cout << "latlon: " << test1 * RAD_TO_DEG << ' ' << test2 * RAD_TO_DEG << std::endl;

return 0;

}


texts: Can I use a ArcGIS 10.1 book for ArcGIS 10.3?

texts:
To learn ArcGIS, I'm going to purchase a home license of ArcGIS 10.3. But because 10.3 just came out, there are no books out yet. Will a 10.1 book suffice, or are there too many changes?

Looking at this book due to its good reviews.


texts: Keeping QGIS from writing all CRS to Spatialite database?

texts:
I am wanting to move from shapefiles to spatialite database files, but am new to QGIS and spatialite.

How do I keep QGIS from creating 4924 records in the spatial_ref_sys table and 4873 records in the spatial_ref_sys_aux table when it creates a spatialite database?

In effect a completely empty spatialite database created by QGIS is 5,460 KB in size. And it appears that the bulk of that is in the two spatial_ref_sys tables.

All I am needing for this data is epsg:4326 WGS 84.

I found an article that mentioned deleting all records from the table except what was needed, but when I try that I get a foreign key error and they will not delete.

DELETE FROM "spatial_ref_sys" where srid!= 4326;

Returns

FOREIGN KEY constraint failed: DELETE FROM "spatial_ref_sys" where srid!= 4326;

Ideally, I would like to know if QGIS has an option to only include certain CRS in the database at time of creation.

Or at worst, how to actually delete the unwanted records.


texts: Using PostgreSQL functions in a QGIS Plugin

texts:
I'm planning to write a QGIS plugin which calculates isochrones.

I have created an tested a network in PostgreSQL and now I want to create a plugin. The Idea is to use a point layer as input and calculate isochrones using the network and nodes I've already created

In PostgreSQL I used the pgDrivingDistance function and alphashape for the calculation of the isochrones.

Is there a way of using the pg function in a QGIS Plugin or are there any other Ideas how to implement a function with the same function as pgDrivingDistance?


texts: Hillshading for small, flat areas in QGIS?

texts:
I have a grid of approximately 172 points with elevation heights for one mound on a flat land. The surface of the mound is 200x200 meters, and the height towards the center is 3 meters.

I did TIN interpolation and created a DEM file which looks excellent. But, when I try to apply hillshade, I get mostly noticeable quadrants which is understandable because of the small surface and the not so big density of the points.

How can I improve the hillshade of my DEM? I experimented with the z-factor, azimuth etc, but still no luck.


texts: FME Server for Linux and Excel reader/writer

texts:
I have an FME model with Excel reader and writer.
In FME Desktop help it is stated: "the format is available only on Windows."

So, as I understand it will not be possible to run this model on FME Server for linux?
Or there is any possibility to do it?

PS we have no FME Server yet. Right now we try to choose the right version FME Server for linux or for windows. Linux is prefered


texts: Got "ImportError: No module named osgeo" in gdal_calc but not in main script

texts:
I could import gdal without any problem in my main Python script with:

from osgeo import gdal

I could also use gdalwarp, gdal_translate with no issue. But when it comes to gdal_calc, I got the error:

Traceback (most recent call last):
File "C:\Program Files (x86)\GDAL\gdal_calc.py", line 51, in <module>
    from osgeo import gdal
ImportError: No module named osgeo

I have tried to find out if anyone had the same issue but found nothing.

(In case you need more info: I installed gdal following the instructions from https://sandbox.idre.ucla.edu/sandbox/tutorials/installing-gdal-for-windows)


texts: Calculate Geometry for Area is 0 in Arcmap

texts: I had a Nation wide data set was in WGS 1984, and then used the project tool  to NAD 1927 Contiguous USA Albers so that I could use the calculate geometry tool in the attribute table to get the Area in square miles, but when I do, I just get 0's when using the new projected file. It isn't disabled, but it is just giving 0 for outputs.


texts: Drawing bus routes on OpenStreetMap with mapnik

texts:
I have successfully set up PostgreSQL with PostGIS and MapNik to serve map tiles. The database is populated with data from OpenStreetMap project. I can see that bus routes are in the data-set, but I have hard time displaying them on the tiles.

Can anybody provide some advice how to do it?


texts: How I can sum 360 raster files in a folder with the Arcgis raster calculator?

texts: In a folder I have 360 precipitation raster files and I want to sum them all to the raster calculator


texts: One SelectFeature control for every Layer?

texts:
I would like to have one SelectFeature control for every layer so that I can have different representational behaviour for every layer.
Means in one layer I want to highlight the features when the mouse hovers in the next not.

When I do this, only the most upper layer is producing onSelect-Events!


texts: ogr2ogr MIF to KML fails on writing kml in linux

texts:
I have multiple MIF files which I would like to translate to KML format.

ogr2ogr -f "KML" geo-mif geo-kml
FAILURE:
Unable to open datasource `geo-kml' with the following drivers.

It seems that ogr2ogr has issues writing KML.
Do I need to install something additional?
Other formats like Mapinfo file works like a charm e.g.:  ogr2ogr -f "Mapinfo file" geo-kml geo-mif works perfectly


texts: Change direction of polylines from DEM data

texts: I am trying to define the direction of my polylines (which represent gutters in a drainage system). Does anyone know if it is possible to define the direction based on elevation data extracted from a DEM? I would like my polylines to be oriented such that they begin at the upstream point and end at the downstream point. I know you can use the "flip" gp tool to change the polyline direction, but this is tedious/impossible to do for a large data set.


texts: Using Mapbox studio style and exported mbtiles to create web map?

texts:
I used Mapbox Studio to create a style using a local tile source. I exported the source to get vector mbtiles.

Now how do I use these mbtiles and the style (which is still in Mapbox Studio) to actually host my slippy maps?

I can't host my data with Mapbox.

There are plenty of options to read mbtiles and serve maps but I can't figure out how to style those maps on the fly using my CartoCSS from Mapbox Studio.


texts: What dem data are Flex Viewers using to create elevation profiles?

texts:
I asked a similar question a year or two ago here:
Can Elevation Profile data in Arc's Flex Viewer 3.0 be exported? which appears to not have an answer.

Another question regarding the Flex Viewer (3.3 now) is where does it grab the dem data for these static elevation profile images it creates?
Is this dem data a WMS that can be utilized in ArcGIS desktop ?


texts: Intersection between an angle and a polygon (postgreSQL + postGIS or Python2.7 + shapely + psycopg2)

texts:
I have a point geom feature (in wkb) in a PostgreSQL database table (with postGIS extension).
There is also an azimuth value (degrees), as float, in an other field of the same table.

How to perform, as simple as it would be to cut a part of a cake, an intersection between a x[m] buffer around the point and a 20° angle around the azimuth (±10° on both side, the azimuth should be the bisecting line of the angle)?

From now I've calculated the buffer as an other geom (as a polygon also in wkb) which is stored in an other field of the table.

I wonder what would be the faster and more CPU effective calculation to intersect this circle with the angle?

The only postGIS function I have found yet is ST_azimuth:
http://www.postgis.org/docs/ST_Azimuth.html
but it may be rather "complicated" to implement as it would need to calculate two "virtual secondary points" that are "far enough" from the buffer and apart respectively from -10 and +10° from the existing azimuth and then make the difference between these two new azimuth angles...
And after, I really can't figure out how to intersect an angle, which is some kind of "not clearly defined" geometry, with the circular buffer.

An other point; at the end I would need to implement this in a Python script (I'm basically using shapely and psycopg2 modules yet).
But for now it could be done in an SQL query directly on the postgreSQL DB.

--

http://toblerity.org/shapely/manual.html

http://initd.org/psycopg/


texts: Database of labeled hyperspectral imagery

texts: Are there any datasets of hyperspectral imagery that can be used for AI training and classification, comparable to ImageNet? The Universitat de Valencia has one listed on their site, HyperLabelMe, but it seems to have been taken down.


texts: Call attribute window with ArcObjects C#

texts:
I want to programatically call AttributeWindow using C# ARCGIS SDK.

UID uidDockWin = new UIDClass();
uidDockWin.Value = "esriEditor.CreateFeatureDockWin";
dockWindow = ArcMap.DockableWindowManager.GetDockableWindow(uidDockWin);
dockWindow.Show();

above code works well and call CreateFeatures Window. I tried following for AttributeWindow,

uidDockWin.Value = "esriEditor.AttributeWindow"; AND
uidDockWin.Value = "{44276914-98C1-11D1-8464-0000F875B9C6}";

however my program crashes when i click on Toolbar Add-in button where above code executes in OnClick method.


texts: Exporting PDF with Hyperlinked Attributes from ArcMap?

texts:
I'm using ArcGIS 10.1 and would like to export a map to a pdf, preserving feature's attributes. I would like to be able to identify, for example, a point feature from the map on pdf with its underlying features (essentially hyperlink to a document or URL). I know one can export to such a pdf in ArcMap, with an advanced option of exporting features with attributes. I use Adobe Acrobat Pro v.10, and cannot seem to be able to find a way to bring these attributes up.

Is this possible?


texts: QGIS: Categorized styling - by multiple columns

texts: I´ve got some polygons that I want to colourize by multiple columns. Can this be done or will I have to edit the attributes? Some polygons only have data-attributes in 1 one column, where as other have data in both column 1,2 and 3. I´m thinking of the concept like the pie chart  function in the diagrams-function in layer properties. In this case I want my polygons splitted in colors by the values in the columns.


texts: How to Geocode an address where the reference data has an address range in a single field

texts:
I would like use something like the "US Address - One Range"(Commonly used locator styles) where there is a from address and to address. However it requires side of street information and my range of addresses are in a single field.

My reference data looks like this:

+----------------+---------+
|     Field      |  Value  |
+----------------+---------+
| Shape          | Polygon |
| Street_Number  | 45-55   |
| Street_Name    | Foo     |
+----------------+---------+

So overall I would like to Geocode multiple addresses to a single point by creating a Address Locator that takes a range of house numbers.

I have looked at this Customizing ArcGIS 10 locators (An Esri Geocoding Technical Paper)

It looks promising but much more complicated than what I am looking for. Is this the only way to Geocode my reference data?


texts: Making centroids for vector grid in QGIS?

texts:
I need to add polygon centroids to vector grid but I can't...
I have output grid like polygon, but I can't choose layer to make centroids

I need it for R studio and pollution model of city.

QGIS 2,4,0


texts: Updating SDE point feature classes with SQL SERVER 2008 R2

texts:
I am contacting you in hope that you can help me resolve an issue using SQL Server 2008 R2 to push updates to a SDE point feature class.

We are using a C# app that uses ArcGIS Runtime to edit the attribute table and create new address records with their cooresponding points.   The problem is we can see the updates in ArcMap, the new records are showing up, and when you click on the record the place where the point should be is highlighted, however it doesn't graphically display the point in ArcMap.

I have been working with ESRI and they can't seem to figure out the problem.  So I thought I would reach out and see if you have any suggestions.


texts: How to set up a default value for key/value (hstore) control in QGIS 3.x?

texts: It seems that PostgreSQL sintax doesn´t work. Any suggestions?


texts: Displaying null locations using CARTO Builder?

texts:
I'm using CARTO to display locations of conference. Some of the conferences were "teleconference" so they don't have an actual location.

Is it possible to display them on the map with the annotation "teleconference" somewhere?

And is it possible to do so in a time-series animated map?


texts: PostGIS: Lines from location touches but not intersecting polygons

texts:
I have OSM .shp file of buildings. I created vertices along each line of the building in that way: vertices (red) & mid points (green). The mid points made by ST_Line_Interpolate_Point() after breaking the buildings to lines.

I have single location (blue star) from which I want to make lines to all vertices, such that the lines not intersects the buildings (blue lines OK, black lines NOT OK).

I started that way:

WITH cross_lines AS (
  SELECT DISTINCT ST_MakeLine(vertices.geom, ST_GeomFromText('POINT(lon lat)', 4326)) AS lines_geom
  FROM vertices, buildings  --vertices are green+red
)
SELECT DISTINCT cross_lines.lines_geom
FROM cross_lines, buildings
WHERE ST_Touches(cross_lines.lines_geom, buildings.geom);

The result is as follows:

I use ST_Touches and not ST_Intersects because from some reason it gives me best results. It gives me almost all the lines but-
It ignores the mid points made by ST_Line_Interpolate_Point() (green).

Is it possible ST_Line_Interpolate_Point not locating points exactly on the lines of the buildings?


texts: Overlaying a (say) 1Km x 1Km to a CartoDB or Here map

texts: Is there a simple way to add a (say) 1Km x 1Km grid layer (or overlay) to a cartoDB map or Here map so that as when zooms in or out the grid lines remain accurate. (NB. An absolute origination/reference point, while desirable, is not necessary, a consistently implemented grid is sufficient.)


texts: Replace null features with features from other field

texts:
I have a polygon shapefile with two fields X and Y. I'd like to create a new field Z reporting the feature/value of X when Y is missing, and the feature of Y otherwise. The function in QGIS Field Calculator would be if(Y IS NULL, X, Y). I'd like to find a way to write it in a PyQGIS code.
On PyQGIS I've attempted:
my_layer = QgsVectorLayer(my_dir, "my_layer", "ogr")

my_layer.dataProvider().addAttributes([QgsField("Z", QVariant.String)])
my_layer.updateFields()

with edit(my_layer):
    for f in my_layer.getFeatures():
        if f["Y"] == NULL:
            my_layer.changeAttributeValue(f.id(), my_layer.fields().lookupField("Z"), f["X"])
            my_layer.updateFeature(f)
        else:
            my_layer.changeAttributeValue(f.id(), my_layer.fields().lookupField("Z"), f["Y"])
            my_layer.updateFeature(f)

but it does not work. Any idea on how to proceed?


texts: Having title which includes dynamic text of layer name?

texts:
What I'm currently working on: approximately 20 historical topographic maps of a specific area spanning 70 years.

What I'd like to do: instead of saving each historical topo as a different file, I would like to just enable each layer (in this case rasters .tif) and print or export to pdf.

For example the particular text in question would say 1938 HISTORICAL TOPOGRAPHIC MAP or 1940 HISTORICAL TOPOGRAPHIC MAP which corresponds to the layer name.

I don't think Data Driven pages applies here since I'm using the same extent (please correct me if I'm wrong).

It would appear that maybe a python script is the best way to go and I found this one here by https://gis.stackexchange.com/users/1608/artwork21.

mxd = arcpy.mapping.MapDocument("CURRENT")
df = arcpy.mapping.ListDataFrames(mxd, "Layers")[0] # default Layers dataframe
for lyr in arcpy.mapping.ListLayers(mxd, "*", df):
     if lyr.visible == True:
          for elm in arcpy.mapping.ListLayoutElements(mxd, "TEXT_ELEMENT"):
               if elm.text == 'Text': # whatever your text element is named here
                    elm.text = lyr.name
                    break

However, it's been years since I used python and I'm not sure how to implement this code in ArcPy.


texts: How do I perform a Union (as per ArcGIS and QGIS) in MapInfo

texts:
I am running MapInfo 15.0 and need to perform a Union of two polygon layers (i.e. combine the two and preserve all boundaries and associated attributes).  I know there is some ambiguity over the term Union but I essentially want to do what Esri considers a Union, only in MapInfo.  More information on Esri's take on Union can be found at:

http://webhelp.esri.com/arcgisdesktop/9.2/index.cfm?TopicName=union_(analysis)

Any help will be great!


texts:
Recently I am using Geoserver grammatically through its REST api. When I request for a map it gives me a map image.

My request code is like this:

public string GetMap(string name)
            {
                WebRequest req = WebRequest.Create("http://localhost:8888/geoserver/wms/reflect?layers=" + name);

                using (WebResponse resp = req.GetResponse())
                {
                    FileStream imageStream;
                    using (imageStream = new FileStream(@"E:\image.png", FileMode.OpenOrCreate))
                    {
                        var responseStream = resp.GetResponseStream();
                        if (responseStream != null) responseStream.CopyTo(imageStream);
                    }
                }
                return @"E:\image.png";
            }

It gives me an image like this:

But i need the whole control of Geoserver map like zooming or others:

How can i do that?


texts: How to link a shapefile with points to a shapefile with lines?

texts: I have a shapefile with points (nodes) and a shapefile with lines (pipeline) which I have extracted from a water distribution network. I have been able to transform the points into the right projection and the layer looks good, but the lines I dont know how to project or illustrate. In the attribute table there is the name of the start node and end node. Can I link the shapefile with lines to the points?


texts: Creating custom icons within Leaflet Marker (only icon/ keep marker style)

texts:
I'm looking to adjust the actual icon within a marker rather than replacing the entire marker with a custom image.

If I'm trying to use a custom .svg file for the 'icon' property, how can i force leaflet (through L.Icon or a similar method) to replace what the white icon is while preserving the actual tear drop shape of the marker?

I just want to use my own icon image instead of calling the MakiMarker api. I am open to adjusting to another library/method compared to L.MakiMarker.

Below is where I want to pass a custom .svg icon to the icon property.


texts: How to use the tool "Polygon to Raster" with arcpy when the polygon is a row from a Shapefile?

texts:
I tried to iterate over the rows from a Shapefile with the search cursor. Every row is a polygon and I want to make a raster out of this polygon so I tried this with arcpy.PolygonToRaster_conversion. But it is not working.  What did I wrong?
This is what I have so far:

import os
import arcpy

arcpy.env.workspace = r"C:\Jul\birds_bear\BOTW.gdb" #path to workspace
arcpy.env.overwriteOutput= True
inputshp = r"C:\Jul\birds_bear\BOTW.gdb\kleiner100"
outraster= r"C:\Jul\birds_bear\BOTW.gdb"
rows = arcpy.SearchCursor(inputshp)
count= 0
for row in rows:
    print row
    valField = "VALUE"
    outraster = os.path.join(outraster, "outRaster" + str(count))
    arcpy.PolygonToRaster_conversion(row, valField, outraster,"","","0,0833333")
    count = count + 1

The shapefile has 16.400 rows. At the end I have to sum the overlapping raster with, I hope it will work for 16400 rasters, see my original question here . I can work with 10.0 and 10.2.


texts: Is it possible to calculate shadow areas of buildings or simulate shadows of buildings in a city using e-cognition software, with the DSM of the buildings and the sun angle and azimuth known?....If not, could you please tell me the software that can do this work.


texts: Using PROJ.4 library to transform from local coordinate system coordinates to global coordinate system using ground control points?

texts:
I have a point cloud whose coordinates are with respect to a local coordinate system. I also have ground control points with GPS values. Can I convert these local coordinates to a global coordinate system using PROJ.4 or any other library?

Any code in Python for the above stated problem would be a great help.


texts: Find the vertices of the edge of the polygon where line intersects using Shapely

texts:
I have this function to get the intersection point of line and a polygon.
How to find the vertices of the polygon(red dots marked in the image) edge where the line intersects using shapely?
I have managed to find the nearest point using shapely nearestPoint function.

# get the intersection point.
def checkIntersection2(polyX, polyY, linePtX, linePtY):
    poly = LineString([(x, y) for x, y in zip(polyX, polyY)])
    line = LineString([(x, y) for x, y in zip(linePtX, linePtY)])
    intPoints = poly.intersection(line)
    return intPoints


texts: ArcGIS Server will not start on new amazon EC2 windows instance

texts:
We have had ArcGIS server running on an EC2 instance without problems for the last 2 years. It is currently at version 10.2. I am now trying to change to a larger instance. I followed Amazons steps to create a new instance

Stop existing instance
Make a new image from the instance
Launch the new image specifying the new larger instance size

I then tried reauthorizing the software which completed successfully (on a side note should this even be necessary since the images are supposed to be identical?)
Everything appeared to go smoothly. I can log into the new instance and all data and files are there. The problem is ArcGIS server will not start. I tried manually starting the service from the service console but it stops after a few seconds. It doesn't look like the tomcat web server gets launched. service_error.log shows the following

INFO: getting data at url: http://169.254.169.254/latest/meta-data/local-ipv4
com.esri.arcgis.discovery.admin.AdminException: Could not create directory path '\SITEHOST\config-store'. This is required to maintain the Site configuration. Please check that the location is valid and that the ArcGIS Server account has permissions to the location.

SITEHOST is the name of the server machine in ArcGIS server. Typing \SITEHOST\config-store into a browser shows page not found, but in the working instance it shows a directory with the expected config files. This is why I think the problem is with tomcat starting. Does anybody know what could be causing this?


texts: Can we include any URL path for the Shapefile while creating a feature store in Deegree GIS?

texts:
I have used the Documentation and the examples of Deegree GIS application and I am able to give the path for my Shapefile in my local system, while creating a feature store just like below.
<ShapeFeatureStore xmlns="http://www.deegree.org/datasource/feature/shape" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" configVersion="3.4.0" xsi:schemaLocation="http://www.deegree.org/datasource/feature/shape http://schemas.deegree.org/datasource/feature/shape/3.4.0/shape.xsd">
  <StorageCRS>EPSG:26912</StorageCRS>
  <File>../../data/utah/vector/SGID93_LOCATION_UDOTMap_CityLocations</File>
</ShapeFeatureStore>

I am having a Shapefile in my Amazon S3 Cloud system, which can be accessed via the URL. Is there any way to include that URL anywhere in the code above, while creating a feature store ?


texts: Transparent NoData through WMS in Leaflet

texts:
I'm trying to create a Leaflet map with an OpenStreetMap as a basemap and a layer from a GeoServer WMS.

I use:

var wmsLayer = L.tileLayer.wms('http://localhost:8080/geoserver/workspace/wms', {
        layers: 'workspace:imgfile.tif',
        transparent: true
        }).addTo(map);

to include the layer.

Unfortunately the no data values in and around the map are shown as white.

The NoData in the underlying .tif are set to 255. The .tif shows as intended in QGIS.

I included <ColorMapEntry color="#FFFFFF" quantity="255" opacity="0"/> in the geoserver-style (SLD) used. I also tried adding an alpha channel with gdalwarp -dstalpha but it doesn't work.

I have observed that the area outside the extent of the .tif is also white. So maybe the L.tileLayer.wms() is not understood as an overlay but Leaflet is just showing this layer only?

Using the same WMS source in OpenLayers (ol.source.TileWMS()) shows the data correctly (transparent at NoData and outside extent) which indicates further that the data is OK and the issue is to be resolved within leaflet (unlike this related solution).

Links for reference: 1 2 3


texts: Dividing stream network of specified lengths in ArcGIS for Desktop?

texts: we are trying to find hydro power potential of river basin by using ArcGis and SWAT ,we have created stream network ,but unable to divide stream network at 500mts interval, unable to find the tool


texts: Is there String Size Restriction in ArcGIS Python Script Tool parameter?

texts:
I'm trying to post a list of unique ID's (as a string) to a Python script tool for processing by ArcGIS Desktop 10.0.

The list may be 50k+ id's long (~450kb).

e.g. "237434,237435,237436,237423,237424,237425,237426"

When defining the input param for my tool, I have tried using the String, and 'Any Value' data types (different version, but described here).

The problem is: the parameter string is being truncated at 4285 Id's.
I have been unable to find any info on the size restrictions of Strings passed in as tool params?

The number of Id's will be increasing over time, so I need to find a more robust way of dealing with large strings being passed in.

The obvious choice is to slice the string into digestable chunks - however, the way that multi value params are defined it appears that they must be preconfigured (I must know, or have a fixed number of chunks).

I want a dynamic solution; any suggestions?

I am able to get the Id's if I read them from a static text file, but the Id list will be a dynamic list being sent from client application for processing.

So basically I need to find a window in the tool params to get the list in.


texts: Global variables in QGIS mixed on some computers

texts:
I create and set some global variables in my QGIS plugin as below:

QgsExpressionContextUtils.setGlobalVariable('variable_name', 'some_value')

The name is unique and specific, so different from any QGIS global variable name and not a regular word. There are 5 variables set in the plugin, some on the init and some later. Values passed are always strings. The problem is that on some computers values are assigned to wrong variable names. Eg.:

it should be:
'name_one': 'variable_one'
'name_two': 'variable_two'
'name_three': 'variable_three'
'name_four': 'variable_four'

it is:
'name_one': 'variable_three'
'name_two': empty
'name_three': 'variable_one'
'name_four': 'variable_three'
# totally mixed up

Has anyone experienced it? I'm using different subversions of QGIS 2.18 on Windows 7 and Ubuntu 16.04.3. Specifically:

Does not work:

Ubuntu 16.04 x64, QGIS 2.18.5
Ubuntu 16.04 x64, QGIS 2.18.13
Windows 7 x86 (Virtualbox on [4]), QGIS 2.18.18

Works

Ubuntu 16.04 x64, QGIS 2.18.17
Windows 7 x86 (Virtualbox guest on [4]), QGIS 2.18.16
Windows 7 x64, QGIS 2.18.9 32-bit

EDIT:
So far I got it working (on OSs that it hadn't before) by declaring and setting to an empty string all of the global variables, that are going to be used later, and it was done at the very beginning of the plugin code. Without purging QGIS config it kind of 'repaired itself' this way. Drawback 1) - it requires knowing in advance what globals exactly will you have. Drawback 2) - I still have no clue what happened here and while it works now, I wouldn't be surprised if it magically stops one day.


texts: How to make traffic time contour map in ArcMap?

texts: I want to make travel time contour map which shows the area people can go from an origin (e.g. downtown, airport) within certain time interval. I have a roadway linear shapefile with the information of travel time and some example maps found online, but lack of the idea how to make the map. Thanks for your help!


texts: gdal_translate using /vsistdin/ and /vsistdout/ to pipe commands

texts:
I am trying to pipe some of my commands with gdal_translate's output. These commands should read from standard input, so I'm writing to /vsistdout/. However, I am getting some errors. To be a bit more generic, I tried some other commands to illustrate the issues I face:

$ gdal_translate -of JPEG2000 Lenna.png /vsistdout/ | gdal_translate -of PNG /vsistdin/ Lenna01.png
ERROR 6: Read or update mode not supported on /vsistdout
ERROR 3: Unable to create file /vsistdout/.

ERROR 4: `/vsistdin/' not recognised as a supported file format.

GDALOpen failed - 4
`/vsistdin/' not recognised as a supported file format.

However, the following works:

$ gdal_translate -of JPEG Lenna.png /vsistdout/ | gdal_translate -of PNG /vsistdin/ Lenna0.png
Input file size is 512, 512
0...10...20...30...40...50...60...70...80...90...100 - done.

(In my tests, I am using my image format, which gets the same error as JPEG2000).
Is it that gdal_translate has a per-format support of the /vsistdout/ (and /vsistdin/) parameter?


texts:
I have two polygon layers - Parcels (in yellow below) and Sites (in Red hatch) - and I am trying to select all polygons from layer 1 that have any part of layer 2 in them, ignoring where edges touch.

In this example I want Parcels B, F, M, N, Q, R, W, X to be selected.  I do not want Parcels A, C, etc. selected even though their edges touch the Sites

So far the only option under Select by Location that I've found to mostly work is "Intersect" with a negative Search Distance.

This negative Search Distance isn't ideal as I need to be sure I'm finding all features (although with 5mm I'm pretty close).  In my testing I found that setting a smaller Search Distance like -0.2 will still select the additional unwanted polygons.  I have checked, there is no overlap here.

Ultimately I am wanting to select all Parcel features that do not have a site, which will be incorporated into a much larger arcpy script.

Is there a geoprocessing tool that ignores touching edges when selecting polygons?

import arcpy
mxd = arcpy.mapping.MapDocument("CURRENT")
inFeat = "TestParcel"
selFeat = "TestSite"
arcpy.SelectLayerByLocation_management(inFeat, "INTERSECT", selFeat, "-0.5 Centimeters", "NEW_SELECTION", "INVERT")

I have tested the other selection Relationships like "Contains" however these don't select all the features I'd expect or select other features I don't want.  The negative Search Distance on Intersect is the best result I've found so far.

Note: My polygons (in both layers) are varying sizes and shapes, and aren't always centred/aligned with each other.  The screenshots here are just of my test data.


texts:
I'm trying to open a NetCDF file from the CHIRPS dataset with yearly precipitation in QGIS 2.14 (OsX) but it hasn't worked at all. I've tried using both the NetCDF browser plugin and opening the NetCDF as raster file (common way of doing it) but according to QGIS all the raster values in the dimension (years) are 0 (no rain). I thought it was an error in the NetCDF but I opened it in ArcMap and it shows all the different values of rain per year pretty well.

Does any of you had this same experience in QGIS with NetCDFs from other sources than CHIRPS? I wonder if it is a bug in GDAL or some very specific characteristic of the type of NetCDF file that the CHIRPS folks are exporting.

Any help or comments are welcome!


texts: Geoserver extremely slow performance on Raspberry pi model 3b

texts:
I have a postGIS database and tomcat8 with Geoserver installed in my raspberry pi model 3b running Ubuntu MATE.

The problem is, geoserver runs really slow, and I mean really slow. It takes up to 15 minutes to load just the presentation view in 127.0.0.1/geoserver.

This is not really an inconvenient anyway, since it's already configured (took forever), but since this is a service I have to use online for a map-based application, it's certainly not good.

We (my team) are working with leaflet, angular and spring, complementary to postgres and the geoserver for the layers. It takes literally 5 to 10 minutes to load the map every time a user moves through it or just zooms in/out.

I have no idea why, connection is good 100mbps up/down, the rest of the tomcat apps are running allright, it doesn't even have any load, it's 2 users at the same time max.


texts:
I have point layer of land cadastre centroids, which have the cadastre number in an attribute table in  ArcMap 10.6 and I am trying to export the layer to AutoCAD 3D Civil so that every point has a unique number as a symbol in AutoCAD.
So far all I managed to do is import a point layers with points as symbology in to AutoCAD.
I have to so some sort of conversion from a table with xy cordinates and edit symbology in the same time.
Does anyone has an idea if that is even possible and how I can do it?


texts: How to determine Distance between points using ArcGIS for Desktop?

texts:
I am trying to figure out the distance between points based on 1000+ points using ArcMap 10.2.

I have transit stops and am trying to figure out how far from stop 1 to stop 2, stop 2 to stop 3, stop 3 to stop 4, stop 4 to stop 5 etc.

Is there a way to do this in ArcMap through the use of a tool?

I'm not good at programming so writing a script is out. Any ideas?


texts: How to make analyzeGPS distanceGPS group by trip?

texts:
I have a dataset with id, lat, lon and time, where each trip has a unique id.

When I try to use mutate with the distanceGPS function in the AnalyzeGPS package and group by trip using this code:

gps%>%group_by(id) %>%
  arrange(time) %>%mutate(d <- c(0,
  distanceGPS(lat1 = lat[1:(n()-1)],
              lon1 = lon[1:(n()-1)],
              lat2 = lat[2:n()],
              lon2 = lon[2:n()]
              )
              )
  )

I get this error:
Error Colum ...<-NULL must be lenght 1(group size), not 3

This works:

d <- distanceGPS(lat1 = gps$lat[1:(length(gps$lat) - 1)],
                 lon1 = gps$lon[1:(length(gps$lon) -1)],
                 lat2 = gps$lat[2:length(gps$lat)],
                 lon2 = gps$lon[2:length(gps$lon)])
trips$d <- c(0, d)

... but it just has one trip with multiple points.


texts: How to digitize from a map with additional detail views?

texts:
I have a map with a river that includes islands within it.  Some islands have a blown up version (detail views) outside of the river showing land ownership details.

I would like to make polygons (shapefiles) of the large islands on the outside of the river and then copy these polygons and shrink them to fit within the lines of their matching islands within the river and then safe the result.  I suspect that the copy and paste part is possible but can the "shrink or resize" part be done?

I want this so that I can do distance queries along the river.  Is this possible and how would I do it?  Thanks


texts: Renaming attributes/fields in shapefile attribute table using QGIS?

texts:
I would like to rename some fields in my attribute table:

Like rename XRF_N3in something else. It seems quite a easy issue but I don't find any information about how can I do that... I've QGIS 1.8 in Mac Mountain Lion.


texts: Calculating distances between boundaries of polygons with Python

texts:
I have two sets of shape-files which have to be worked on within a Python Script. Those two Shape Files might differ from each other (same thematic content, but 2 different epochs).

If there occurred a change in geometry over the years within an object I am supposed to detect it, but only if the boundary changed more than x m in at least 1 place.

I am working completely with open source libraries within my scripts.

Does anyone have an idea how to solve that problem?


texts: How to monitor geoserver map performance using jmeter?

texts:
I want to monitor the geoserver map performance using jmeter which is been load balanced in Apache http server, now I want to know how should I achieve this and what all changes I have to do in jmeter.

Can anyone tell us the steps for configuring jmeter and how i can see the performance while using geoserver as well as by calling in openlayers.

These are the results when I call an application and how do I evaluate performance by this results.


texts:
I'd like to get a list of all metadata from GeoNetwork. According to this manual https://geonetwork-opensource.org/manuals/2.10.4/eng/developer/xml_services/metadata_xml_search_retrieve.html#search-metadata-xml-search I'm using this curl command
curl -X POST -u user:pass http://geonetwork_server/geonetwork/srv/eng/xml.search

I want to do it using Python, so I'm using requests to send POST request to the same url as above. I am also sending this string as 'data' parameter of Python's POST request:
<?xml version="1.0" encoding="UTF-8"?><request/>

The response comes with this element, among others:
<response from="1" to="100" selected="0" maxPageSize="100">

I'd like to get list of all metadata. The above request returns only 100 GeoNetwork records. How to get whole list, or another page of this list?
Is this link from above the official documentation for latest GeoNetwork version? I am using GN 3.10.


texts:
Is it possible to define a custom marker as symbolizer for a point in QGIS?

I'm looking into using a variation on the "Vector Field Marker" with a couple of arrows pointing at dynamic angles (the angles being read as point attributes).

I imagine it's doable in C++ (within the QGIS code) but I'm looking for a scripting solution (python).

I'm adding a representation of what I'm trying to achieve - for two points A and B with attributes x,y,z.


texts: Add an attribute to polygon layer based on if it intersects with another layer

texts: I have a layer of buildings and want to give it each polygon an attribute of 1 or 0, based on if it intersects a flood extent. Basically, a flooded or not flooded attribute.


texts: MapServer *.map / WMS show only a part of the whole shape

texts:
I have a problem with a MapServer configuration: test.map // It's installed under: C:\MS4W\ms4w\apps\mapserver_wms. The content of the test.map is follwing:

Map
  NAME "MAPSERVER_QUICKSTART"
  STATUS ON
  EXTENT -137 29 -53 88
  UNITS DD
  SHAPEPATH "data"
  #SIZE 800 600

  IMAGETYPE PNG24

  PROJECTION
    "init=epsg:4326"
  END

  WEB
    METADATA
      ows_title "MapServer Quickstart"
      ows_enable_request "*"
      ows_srs "EPSG:4326 EPSG:25832 EPSG:25833"
    END
  END

  LAYER
    NAME "Countries"
    STATUS ON
    TYPE POLYGON
    DATA "ne_10m_admin_0_countries"
    CLASS
      STYLE
        COLOR 246 241 223
        OUTLINECOLOR 0 0 0
      END
    END
  END

END

Now I would like to use this whole world map as a WMS Serice with the URL:

http://localhost:81/cgi-bin/mapserv.exe?map=C:/MS4W/ms4w/apps/mapserver_wms/test.map&SERVICE=WMS&REQUEST=GetCapabilities&VERSION=1.3.0

I put this WMS URL into QGIS and connect: This is working, but it show ONLY a small part of the whole area (world)! See the diffrence between blue and yellow (blue = imported the shape directly to QGIS // yellow = same shapes as WMS but only a part of the whole). How can I give out the whole world map in the WMS?


texts:
I am super new to QGIS so I may be doing something really dumb, but it doesn't seem to work to create a heatmap in QGIS Desktop 2.8.2 (just downloaded this yesterday so I think it is the most current stable release). I load my points in from a SQL server geometry column, then click Raster > Heatmap > Heatmap. Then I change none of the defaults on the "Heatmap Plugin" dialog except where to save the raster and hit "OK".  The resulting raster looks like it uses very few of the points in the "Input point layer". Is there any way to get this to use all the point when constructing the heatmap?

Here are some pictures showing this process:

After loading all my points:

Heatmap Dialog Box:

After hitting OK on the Heatmap dialog box above and it renders (a progress bar dialog that isn't in the picture pops up for a minute after hitting "OK"):

This is what the rendered tiff looks like when I open it in Photoshop:


texts: Linking Google sheet as a database for GIS

texts:
I’m working on a database update, at the moment this database is simply an excel spreadsheet not formatted to be used in a GIS project. The database has to be updated and modified over time. Some fields will be added along the way, some entries too. These operations will be done by different people not familiar with GIS. In the end the database have to be associated with polygons in existing shapefiles.

I would like to limit the number of manipulations between the modifiable database and the GIS software. So exporting a new .CSV every time a modification is done, is not what I’m looking for.

I would like to know how (or in what “software”) to setup the database so it would be easy to modify it for non-GIS users (like in Google Spreadsheet) and use it at the same time in QGIS or ArcGIS to link its data to existing shapefiles.

Google Spreadsheet is easily modifiable and shareable, if there was a way to link it with QGIS or ArcGIS so it would constantly update, it would be perfect. But if anyone now another way to proceed, I’m listening!


texts: Snap point to point with St_Snap Postgis

texts:
I write a request that allows to snap a point to another.
Here is the request:

 SELECT
        f.gid as gid,
        ST_Snap (f.Geom, g.Geom, min, 5) as geom
 FROM
        boiten as f,
        (SELECT ST_Collect (Geom) as Geom
        FROM ft_chambre) as g

It works but the snapped points missed the closest position

Update:

For the CRS it's 2154 ,I want to snap the point to the red ones

I modified my code to this:

SELECT
   f.gid as gid,
   ST_Snap(f.Geom, g.Geom, st_distance(f.Geom, g.Geom)*1.2) as geom
FROM
   boiten as f,
   (SELECT ST_Collect(Geom) as Geom
   FROM ft_chambre) as g

But like you see there are points witch snapped to the wrong place or has snapped to more than one.

For exemple
 This point has snapped to the wrong place where the other point has already snapped to


texts: How to translate John Snow data sets to map coordinates

texts:
In the HistData package for R (https://r-forge.r-project.org/R/?group_id=574) I have the data sets related to John Snow's
map of the cholera outbreak in London, 1854.  I believe they are
authoritative, having been carefully digitized under supervision
of Walter Tobler. Some details on these data sets are described by
John Mackenzie, at http://www1.udel.edu/johnmack/frec480/cholera/cholera2.html.

Unfortunately, the coordinates of deaths, pumps,
and streets use an arbitrary coordinate system, not map coordinates
suitable for other GIS applications or mapping software in R
(spatial packages, ggmap, etc.)

In http://freakonometrics.hypotheses.org/19213 Arthur Charpentier
uses ggmap with a version of the John Snow data from
http://www.rtwilson.com/downloads/SnowGIS_v2.zip.  The Cholera_Deaths.shp file, however lists only 489 deaths, not the 578 I have recorded in HistData::Snow.deaths.

One idea is to find the relations between the means and
standard deviations of the (x,y) coordinates and rescale linearly,
but maybe there is a better way?

Here's what I've tried so far

> data(Snow.deaths, package="HistData")
> D <- Snow.deaths[,2:3]
> colMeans(D)
       x        y
13.03312 11.69721
> var(D)
          x         y
x 3.8150987 0.3802654
y 0.3802654 2.7213828

Read the Cholera_deaths file

> folder <- "C:/Dropbox/R/data/Snow/SnowGIS_v2/SnowGIS"
> library(maptools)
> deaths <- readShapePoints(file.path(folder, "Cholera_Deaths"))
> head(deaths@coords)
  coords.x1 coords.x2
0  529308.7  181031.4
1  529312.2  181025.2
2  529314.4  181020.3
3  529317.4  181014.3
4  529320.7  181007.9
5  529336.7  181006.0
> # deaths has only 250 observations; 489 deaths
> sum(deaths@data$Count)
[1] 489

 > # try to relate to Snow.deaths
> X <- deaths@coords
> colnames(X) <- c("x", "y")
>
> XX <- data.frame(X, Freq=deaths@data$Count)
> XX <- vcdExtra::expand.dft(XX)
>
> colMeans(XX)
       x        y
529414.8 181031.9
> var(XX)
          x        y
x 10813.816 1521.693
y  1521.693 6227.924
>

OK, then I try to re-scale D to have the same means and standard
deviations as XX,  but something doesn't work correctly here-- the column means of Dscaled should have turned out to be equal to those of XX:

> # scale D to have the same means and standard deviations as XX
> Dscaled <- scale(D, center=TRUE, scale=TRUE)
> Dscaled <- scale(Dscaled, center=colMeans(XX), scale=sqrt(diag(var(XX))))
> colMeans(Dscaled)
        x         y
-5091.040 -2293.947
>

EDIT: It might be helpful in this problem to see Snow's map as drawn by the new function, SnowMap(axis.labels=TRUE) now in the development version of HistData (rev 102) on R-Forge.  The axis labels show the origin of the coordinate system in the lower left corner as they are in my data Snow.* data sets.


texts: QGIS Save location and information to Excel file

texts:
I have georeferenced an old map in QGIS and want to save the location of a place and then be able to type information about that captured location (owner, acreage etc) in a table.

I am an historical geographer struggling to learn QGIS so apologies if I have missed an already answered query - I did try and check.


texts: How to calculate the optimal zoom-level to display two or more points on a map

texts: We want to display several markers on a static map and want to calculate the optimal zoom-level like Google Maps does. We already calculated the bounding rectangle and the centerpoint of the map but now we have a hard time to calculate the correct zoomlevel to display the whole bounding rectangle. Can someone please point us into the right direction?


texts: Is it possible to change my Vector Strategy? Open Layers

texts:
I added a layer previously on my code, thats it's source scope.
const vectorSource = new ol.source.Vector({
    format: new ol.format.GeoJSON(),
    strategy: ol.loadingstrategy.bbox,
    loader: ()=> { ... }
});

I'd like to know if there is any way to change it's strategy using OpenLayers functions, like getSource().setStrategy() (i couldn't find it).
I also tried to create a new source but with ol.loadingstrategy.all as a strategy and after this used setSource() into the layer, but didn't worked as expected.
OpenLayers updated my feature's source but kept the old one working as well(somewhere in the memory).
Using a new source made my request double, now it requests the entire layer and partially using BBOX, how can i get rid of the old one? I think changing the strategy may be better, but i can't figure out how.
Deleting the entire layer and recreating it, is something that i don't want to do!


texts: Obtaining coordinates when computing distance of a polygon to a line

texts:
I have a shapefile of a county.  I am interested in computing its distance to a line, which I've been doing using Near.  I want to know the coordinates on the line that produces the minimum distance.  How do I obtain the coordinates?

I use ArcMap in ArcGIS 10.1