texts:
I have two raster files as numpy arrays. I'm combining the rasters into one file. The rasters overlap on several pixel but have different values for the the overlapping pixels (one raster has nodata value for pixels, the other has the actual elevation data point). How do I combine the two arrays such that I only take the actual elevation data point instead of the nodatavalue?

file one: [-999,-999,-999,-999]
file two: [7,4,3,2]

desired output for pixels after combination: [7,4,3,2]

My code that I've been using produces this image. http://i.imgur.com/2nFOBMH.jpg

from osgeo import gdal
from gdalconst import GA_ReadOnly
f1 = r'C:\Users\User\Data\NHDPlusCA\NHDPlus18\NEDSnapshot\Ned18a\elev_cm'
f2 = r'C:\Users\User\Data\NHDPlusCA\NHDPlus18\NEDSnapshot\Ned18b\elev_cm'

xmins = []
ymins = []
xmaxs = []
ymaxs = []

files = [f1, f2]

for f in files:
    raster = gdal.Open(f, GA_ReadOnly)
    projection = raster.GetProjection()
    cols = raster.RasterXSize
    rows = raster.RasterYSize
    f_xmin, f_pwidth, f_xskew, f_ymax, f_yskew, f_pheight = raster.GetGeoTransform()
    xmins.append(f_xmin)
    xmaxs.append(f_xmin + (f_pwidth*cols))
    ymaxs.append(f_ymax)
    ymins.append(f_ymax + (f_pheight*rows))
    del raster

#extents
xmax = max(xmaxs)
ymax = max(ymaxs)
xmin = min(xmins)
ymin = min(ymins)
newwidth = round((xmax - (xmin))/f_pwidth)
newheight = round((ymax -(ymin))/abs(f_pheight))
print('xmins', xmins)
print('ymins', ymins)
print('xmaxs', xmaxs)
print('ymaxs', ymaxs)
print('largest x', xmax)
print('smallest x', xmin)
print('largest y', ymax)
print('smallest y', ymin)
print('new raster width', newwidth)
print('new raster height', newheight)

driver = gdal.GetDriverByName('GTiff')
new_ds = driver.Create(
    r'C:\Users\User\Desktop\out.tif',
    newwidth,  # xsize
    newheight,  # ysize
    1,    # number of bands
    gdal.GDT_Int32,  # The datatype of the pixel values
    options=[  # Format-specific creation options.
        'TILED=YES',
        'BIGTIFF=IF_SAFER',
        'BLOCKXSIZE=256',  # must be a power of 2
        'BLOCKYSIZE=256',# also power of 2, need not match BLOCKXSIZE
        'COMPRESS=LZW'
    ])

# fill the new raster with nodata values
outband = new_ds.GetRasterBand(1)
outband.SetNoDataValue(-2147483647)
outband.Fill(-2147483647)
geotransform = [xmin, f_pwidth, 0, ymax, 0, f_pheight]
new_ds.SetGeoTransform(geotransform)
new_ds.SetProjection(projection)
for f in files:
    raster = gdal.Open(f, GA_ReadOnly)
    f_xmin, f_pwidth, f_xskew, f_ymax, f_yskew, f_pheight = raster.GetGeoTransform() #geotransform for the old files being iterated over
    cols = raster.RasterXSize
    rows = raster.RasterYSize
    xoffset = int((f_xmin - xmin) / f_pwidth)
    yoffset = int((f_ymax - ymax) / f_pheight)
    band = raster.GetRasterBand(1)
    data = band.ReadAsArray(0,0,cols,rows)
    outband.WriteArray(data,xoffset,yoffset)
    outband.FlushCache()

# When all references to new_ds are unset, the dataset is closed and flushed to disk
new_ds = None


texts: Identify partial intersections between a polygon layer and a line layer using PyQGIS

texts:
I have a drains layer (polygon) around which I have created buffers of varying distances. I want to find the intersection between the different buffer layers and pipes (line) layer and accordingly assign score in a field. I have used the following code:

pathName = os.path.dirname(unicode(drainLayer.dataProvider().dataSourceUri()))
waterLayer.startEditing()
for filename in os.listdir(pathName):
    if filename.endswith('1.shp') and ('buf' in filename):
        #name = ntpath.basename(filename).split('.')[0]
        buf = QgsVectorLayer(pathName,filename,"ogr")
        for f1 in waterLayer.getFeatures():
            for f2 in buf.getFeatures():
                if f1.geometry().intersects(f2.geometry()):
                    f1['drainScore']=3
                    waterLayer.updateFeature(f1)
waterLayer.commitChanges()
for filename in os.listdir(pathName):
    if filename.endswith('2.shp') and ('buf' in filename):
        #name = ntpath.basename(filename).split('.')[0]
        buf2 = QgsVectorLayer(pathName,filename,"ogr")
        for f1 in waterLayer.getFeatures():
            for f2 in buf2.getFeatures():
                if (f1.geometry().crosses(f2.geometry())):
                    f1['drainScore']=2
                    waterLayer.updateFeature(f1)
                else:
                    f1['drainScore'] =1
                    waterLayer.updateFeature(f1)

Wherever the intersections are complete it gives the correct output. But wherever there are partial intersections, it does not identify them. I have used crosses,touches and intersects method in the second for loop.

I have also tried processing algorithms like qgis:intersection, qgis:joinattributesbylocation (with crosses as PREDICATE) but none of them recognize partial intersections. When I use the vector tools of intersection and join attributes by location from QGIS 3 desktop, it gives the correct output, i.e., it recognizes partial as well as complete intersection. I have attached an image of the kind of partial intersection I want the code to recognize.


texts:
I'm dealing with several thousand polylines of Geologic data with mixed degrees of curvature.  I'm attempting to quantify just "how curved" these lines are and was curious if there was a means to do so in ArcGIS.

I've also exported by Dataset to Matlab in an attempt to find a solution there.

I am seeking a repeatable and robust means for calculating curvature.


texts: How to include `rasterizeWholeImage` when exporting QGIS layout?

texts:
I use something like the following to export a layout to pdf:
projectInstance = QgsProject.instance()
layoutmanager = projectInstance.layoutManager()
layout = layoutmanager.layoutByName("test")
exporter = QgsLayoutExporter(layout)
exporter.exportToPdf("C://test.pdf", QgsLayoutExporter.PdfExportSettings())

But I want to use the rasterizeWholeImage parameter when exporting otherwise it defaults to using vectors. How can I set this parameter when setting up QgsLayoutExporter.PdfExportSettings()?


texts: empty file geodatabase error (nowhere feature class)using arcpy

texts:
i want to import a python(2.7) script to arcmap 10.1 my script  has 2 inputs variables one for workspace and another raster dataset and i want o convert input raster as feature class polygon .
i want to create a file geodatabase and import in her a export feature class,run this script and finaly show me complete message(not error)...
but created the basis is empty,nowhere feature class.my feature class=Catchment
and i want to import in gdb(my file geodatabase) but if i run this script show me complete message but the gdb(my file geodatabase) is empty

import arcpy
import sys
import os
from arcpy import env as e
from arcpy.sa import *
arcpy.ImportToolbox("C:\Program Files (x86)\ArcGIS\Desktop10.1\ArcToolbox\Toolboxes\Arc Hydro Tools.tbx", "archydrotools")
try:
 e.workspace = arcpy.GetParameterAsText(0) # import workspace
 Cat = arcpy.GetParameterAsText(1)         # import raster dataset
 ws = e.workspace
 gdb = "teleia"                            # import name of geodatabase
 # Local variables:
 arcpy.CheckOutExtension("spatial")
 arcpy.env.overwriteOutput = True
 arcpy.CreateFileGDB_management(ws,gdb)
 out = ws + os.sep + gdb + os.sep
 # Process: Catchment Polygon Processing
 arcpy.CatchmentPolyProcessing_archydrotools(Cat, Catchment)
 Catchment =out + "Catchment"
except:
 arcpy.GetMessages()

any idea ?


texts: Merging two shapefiles only for those polygons that have same boundaries or within tolerance distance using ArcMap

texts:
I have two shapefiles for district boundaries of India in 2001 and 2020. Because boundary of districts changed over time since districts merge and split, I want to export only those districts (they are polygons in shapefile) that have the same boundaries in both years (or are located within the tolerance distance, e.g. 10km).
I am using ArcMap 10.4.1.


texts: What are the FOSS equivalents to these ArcGIS products?

texts:
I'm a long-term user of ESRI software, and I now have a need to use free and open source software.

I've been reading about OpenLayers, PostgreSQL, PostGIS, GeoServer and MapServer but I can't find a good overview stating what each product does, why it's needed, and how they all fit together.

For example, http://www.osgeo.org/ and http://freegis.org/ both list a bunch of products but don't give enough information for me to decide which ones are relevant.

In ArcGIS, I would use the following:

file geodatabase or ArcSDE geodatabase for data storage
ArcMap desktop to edit the data and compile the map document
ArcGIS Server to create web services
ArcGIS Server JavaScript API to create the end-user maps from the web services

What are the equivalent FOSS products?


texts: Loop for exporting shapefiles from an existing shapefile

texts:
I have  a pnt shp with 50 pnts.

I want to export each one of it in a separate shapefile with arcpy writting a loop.
My problem is that i want the output name i create has a different name, e.g. for 1st pnt i want the name to be pnt1_lyr for the 2nd pnt2_lyr etc. I wrote that code but i always have an error. Can someone correct me;

import arcpy
arcpy.env.workspace = "c:/Wing/Sk"
for i in range(0, 49):
    arcpy.MakeFeatureLayer_management("spatially50.shp", "pnt[i]_lyr", "FID = (i)")


texts: PYQGIS processing.runalg stops recognising shapefiles

texts:
I am running a processing algorithm qgis:dissolve with pyqgis outside QGIS in Ubuntu 16.04 with Qgis 2.18.4

I wrote a python script which takes shapefiles from a folder and runs this with shpfile as the input and generate a shapefile in another folder.

import sys
import glob
import subprocess
import time
import qgis
from qgis.core import *

def dissolve():
    for shpfile in sorted(glob.glob("/home/user/test/*.shp")):
        filename = shpfile.split("/")[-1]
        f_name = filename.split(".")[0]
        #filename looks like input_001
        num = int(f_name.split("_")[1])
        print filename
        processing.runalg("qgis:dissolve",shpfile,True,"ID","./shp_dissolved_test/"+filename)
        if num % 50 == 0:
            time.sleep(1)

def main():
    dissolve()

if __name__ == "__main__":
    sys.path.append("/usr/share/qgis/python/plugins")

    # Initialize QGIS Application
    QgsApplication.setPrefixPath("/usr", True)
    app = QgsApplication([], True)
    QgsApplication.initQgis()

    import processing
    from processing.core.Processing import Processing
    Processing.initialize()
    Processing.updateAlgsList()

    start_time = time.time()

    main()

    secs = time.time() - start_time
    minutes = secs/60
    print "---" + str(minutes) + "-- min"

I have around 600 shapefiles. But the proceesing fails after 200 or so shapefiles and gives error number 4

Unable to execute algorithm
Wrong parameter value: /home/user/test/input_197.shp
input_197
ERROR 4: Unable to open /home/user/test/input_197.shx or /home/user/test/input_197.SHX.
ERROR 4: Unable to open /home/user/test/input_197.shx or /home/user/test/input_197.SHX.
ERROR 4: `/home/user/test/input_197.shp' not recognised as a supported file format.

But after I rerun it again from that file (input_197) it starts to process again.

Do i need to give some wait time between each processing or is it something else?

Note: After this error I added a time.sleep(1) after every 50 shapefile processing but it still fails.


texts: Solution to draw objects on layers and share maps

texts:
I'm looking for an online solution to build custom maps on top of OSM/Bing/Google layers and share them via export or link. I need to find an address, draw and drop lines, polygons, circles, add and remove markers. I'd like to change color and opacity of those objects.

I tried :

Leaflet with Leaflet.draw, Leaflet.StyleEditor and Leaflet.print | it works but it shows limits when printing big maps to PDF (even the controls and boxes are printed) plus I can't save it to edit later.
share.mapbbcode.org and it's really good except I can't draw circle, see my contact with the author
Ilya Zverev (the author of mapbbcode.org) points me to umap which is great but have the same limitation with circles

I'm disappointed as I see it's close to fits my needs but I don't know how to make it.
Has someone else been successful at drawing lines, polygons, circles (and changing the color and opacity) and saving/sharing maps ?


texts: How to use non-native spatial data from MS SQL Server with ogr2ogr?

texts:
I have some tables in MS SQL Server (2012) with WKT string columns.
I want to export these data to a MapInfo tab file.

Currently I am using the following command:

ogr2ogr -f "MapInfo File" testsql.tab "MSSQL:Server=etc." -sql "select geography::STGeomFromText(FACE,4326),NAME from TABLE" -a_srs "EPSG:4326"

Where FACE is a string column containing WKT geometries.

This gives some problems, because the coordinates are switched. This is supposedly due to a change between SQL Server 2008 and 2012. Outputting it to WKB works better.

I was wondering if it is possible that ogr2ogr uses the WKT column directly without converting it to the MS SQL Spatial type first?


texts: Calculate the area covered by spectral end members in Google Earth Engine?

texts:
3
down vote
favorite
I am doing a project on calculating the growth of urban areas in different cities located in different parts of India. I have been able to classify the area into urban, vegetation and water. Now the problem is I cannot calculate the area covered by each pixel. I searched on Google and the available methods uses band values, what I need is to find the area covered by urban/water/vegetation based on spectral end member values that is been defined. If anyone has idea on how to perform it please share and if possible suggest an idea to export it into a graph. What I need is to find the area covered by the different spectral endmembers (urban,veg,water) values which has been defined.

Here is a link to the script:

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

I'm fairly new to Earth Engine so my knowledge is quite limited, please excuse


texts:
I want to try and split a GeoJSON FeatureCollection into the individual features. So that the name of the file I write to is the "ED_ID" field. So far I can only read the JSON.

I am new to python and programming in general.

import json
with open('ded.geojson') as f:
data = json.load(f)

My data looks like this I have pulled out most of the coordinate data and this is one feature the file contains multiple features which I wish to separate.

 {"type":"FeatureCollection","features":[{"type":"Feature","properties":{"OBJECTID_12":3001,"ED_ID":"197023","ED_ENGLISH":"BREEDOGE","ED_GAEILGE":" ","COUNTY":"ROSCOMMON","CONTAE":"Ros Comáin","PROVINCE":"Connacht","CENTROID_X":"573408","CENTROID_Y":"796269","OBJECTID":927,"ED_ID_1":"197023","FIRST_COUNTY":"ROSCOMMON","FIRST_CONTAE":"Ros Comáin","FIRST_PROVINCE":"Connacht","FIRST_ED_ENGLISH":"BREEDOGE","FIRST_ED_GAEILGE":" ","MaxSimpTol":100,"MinSimpTol":100,"GUID":"2AE196291AE513A3E055000000000001","ObjectID_1":927,"Shape__Area":0.00260221958729151,"Shape__Length":0.361585766792224},"geometry":{"type":"Polygon","coordinates":[[[-8.40526490835231,53.926763221445],[-8.39132968800989,53.9252835371864],[-8.38738628086417,53.9256611572713],[-8.38443774016629,53.927271794708],[-8.37396157235263,53.9250928131867]]]}}


texts: ArcGIS 10.2 Query layer ID field using double / floating field?

texts: I am using ArcGIS for Desktop 10.2.   Am I correct in understanding that ArcGIS will only let you select an integer or text fields as the "key" ID when adding a query layer?    Our Oracle enterprise DB mistakenly used NUMBER(22), Data Type: FLOAT, as unique identifiers even though the IDs are all whole numbers; however, I seem unable to select the 'number' fields as keys.  It will allow me to select integers and text (I guess you can combine a number with text to make a unique identifier).     I've tried CAST-ing them as integers in the query but the queries are crashing ArcGIS and I'm trying to figure out whether there is any way around the CAST statement to see if that is the problem.


texts: How to properly uninstall OSGeo4W+GRASS+QGIS on Windows7 64bit

texts:
Over the years I have tried QGIS and other mapping systems.  Lately I wanted to make a go at using the newer versions of QGIS but wanted the ability to run the new one, 3.4 or 3.6 and to keep the old LTR 2.18 running parallel.  I read that is possible through OSGeo4W.

Something went wrong in a DL or install and now I have versions of GRASS growing all through my Win7 64 bit system. (Growing Grass was a joke..) but GRASS is install all over in different versions, as is QGIS.

One helpful post suggested just deleting the normal C:\OSGeo4W directory.

I have tried just about everything and for my trouble, I have folders and files all around and they do not work or they error. (Although I must have deleted something else because now I can't produce an error, just nothing happens.

How can I clean this up?
In searching it looks like there is no uninstall facility listed.
ultimately I want to start clean with a new install.


texts: Can you create feature classes with the same name in different feature dataset

texts:
I've searched but can't seem to find any info on this which is surprising as this should be an easy one to answer.

I'm using ArcGIS 9.3 and trying to create a feature class called 'cadastre' in multiple feature datasets within the same file geodatabase.

ArcCatalog is giving me a 'The specified name already exists.  Please enter a different name' error.

It's been a while since I used ESRI software but I would have imagined this was possible since they are in seperate feature datasets.

Is this expected funcitonality of the geodatabase?


texts: Global name not defined, is this a syntax error?

texts:
I have been working with classes, and I have fixed several mistakes related to indentation, but in this case I have been checking all the possible errors. It might be a concept error. If someone can tell me what is my mistake, it will be very helpful.

From this code, I get the error "NameError": global name 'vector' is not defined

import qgis
from qgis.core import *
from qgis.gui import *
from qgis.utils import iface

class Varias_Funciones:

         vector = QgsVectorLayer("ROUTE","Municipios","ogr")

     #LA CLASE TENDRA LA CAPACIDAD DE ENTRAR EN EL INTERFAZ DE QGIS
         def __init__(self, iface):
             self.iface = iface

     #CARGA DE LA CAPA DE MUNICIPIOS (CAPA VECTORIAL)
         def capa_vectorial(self):
             QgsMapLayerRegistry.instance().addMapLayer(vector)

     #CARGA DE LA CAPA DE RASTER MTN50
         def capa_raster(self):
             raster = QgsRasterLayer("ROUTE", "MTN50")
             QgsMapLayerRegistry.instance().addMapLayers([raster])

     #CARGA DE SISTEMAS DE COORDEANDAS
         def coordenadas(self):
             crs = QgsCoordinateReferenceSystem(4258, QgsCoordinateReferenceSystem.EpsgCrsId)
             vector.setCrs(crs)

     #ABRIR LAS PROPIEDADES
         def abrir_propiedades(self):
             self.iface.showLayerProperties(vector)

comprobar_funciones = Varias_Funciones(qgis.utils.iface)
comprobar_funciones.capa_vectorial()
comprobar_funciones.capa_raster()
comprobar_funciones.coordenadas()
comprobar_funciones.abrir_propiedades()


texts: Is it possible to change the color model of an entire QGIS project from RGB to CYMK? I know that it's possible to change the color models of individual symbols/features (through the style manager), but I can't seem to find an option for the whole project. While my workflow usually sees me polishing my maps in Inkscape/Illustrator (where it's easy to switch between RGB/CMYK), it would be nice to have the option in QGIS as well.


texts: Remapping layers in QGIS project to other folders?

texts:
I created a QGIS project containing several layers which are on my PC. I saved the project to my local drive and then I saved it to the network and copied all the relevent project layers into a folder except the osmap shp files which are on a separate network drive.

I then hunted for a way to remap the QGIS project layers so it points to the network copies and not my local drive.

In AutoCAD this is easy, but I can't find a way to do it in QGIS.

I then tried just deleting the layers from the project and then re-adding them from the new location, but I got a minidump crash every time.

I ended up saving the print composers as templates, exporting the layer styles. Then I created a new project, using the files on the network, and created new print composers from the saved templates.

Is there an easier way to make this project available to other users?

I've seen some info about mapping bad layers, but I couldn't find a dialogue box to do that.

If nothing else, everything on the network will eventually end up on the archive drive and we need to be able to keep projects together.

For example, in CAD if I open a project and someone has moved a file into a new folder, it will tell me that it can't find the external reference. I then open the x-refs window, scroll to the saved path, click browse, then find the file in its new folder and click on it. That saves the changed path.

Based on responses so far, I now know that I need to know how to set a path to relative AND how to remap a path to a different folder.
One response gives a plugin for the latter, so I will try that.
I am close to understanding the former.

The answer seems to be as follows -
Save project. Then use the save button again which creates a project home folder. Use Qpackage to copy the shp files to the new destination. Than open that new project and save to create a new project home folder.
Then use change data source to remap all the layers to the project home folder (as QGIS does NOT save anything as relative at all. It saves all paths as absolute.)


texts: "Find by attribute" plugin UnicodeEncodeError

texts:
This is the error message, when launching the plugin. I hope this is the right place, where to report this kind of message.

  An error has occured while executing Python code:

  Traceback (most recent call last):
    File "/home/saabrs/.qgis//python/plugins/findbyattr/Findbyattr.py", line 66, in run
      tempname = str(name).rstrip('01234567890')
  UnicodeEncodeError: 'ascii' codec can't encode character u'\u0161' in position 5: ordinal not in range(128)

  Python version:
  2.7.3 (default, Apr 20 2012, 23:06:40)
  [GCC 4.6.3]

  QGIS version:
  1.8.0-Lisboa Lisboa, exported

  Python path: ['/usr/share/qgis/python', '/home/saabrs/.qgis//python', '/home/saabrs/.qgis//python/plugins', '/usr/share/qgis/python/plugins', '/usr/lib/python2.7', '/usr/lib/python2.7/plat-linux2', '/usr/lib/python2.7/lib-tk', '/usr/lib/python2.7/lib-old', '/usr/lib/python2.7/lib-dynload', '/usr/local/lib/python2.7/dist-packages', '/usr/lib/python2.7/dist-packages', '/usr/lib/python2.7/dist-packages/PIL', '/usr/lib/python2.7/dist-packages/gst-0.10', '/usr/lib/python2.7/dist-packages/gtk-2.0', '/usr/lib/pymodules/python2.7', '/usr/lib/python2.7/dist-packages/ubuntu-sso-client', '/usr/lib/python2.7/dist-packages/ubuntuone-client', '/usr/lib/python2.7/dist-packages/ubuntuone-control-panel', '/usr/lib/python2.7/dist-packages/ubuntuone-couch', '/usr/lib/python2.7/dist-packages/ubuntuone-installer', '/usr/lib/python2.7/dist-packages/ubuntuone-storage-protocol', '/usr/lib/python2.7/dist-packages/wx-2.8-gtk2-unicode', '~/.qgis/python', '/home/saabrs/.qgis/python/plugins/mmqgis/forms', '/home/saabrs/.qgis/python/plugins/QuickMultiAttributeEdit/forms', '/usr/share/qgis/python/plugins/fTools/tools', '/home/saabrs/Darbvirsma/mag/gis']


texts: QGIS error: could not perform 'gdal2xyz' conversion

texts:
GDAL command:
python3 -m gdal2xyz -band 1 C:/Users/Alexander/AppData/Local/Temp/processing_ZNucZn/b3f4565e0b4a4bf5b18081162ec805b0/OUTPUT.tif C:/Users/Alexander/AppData/Local/Temp/processing_ZNucZn/11f4e9949aa64c91b4193e8c006ae683/OUTPUT.csv
GDAL command output:
Traceback (most recent call last):

File "C:\PROGRA~1\QGIS 3.16\apps\Python37\lib\site-packages\osgeo\__init__.py", line 18, in swig_import_helper
return importlib.import_module(mname)
File "C:\PROGRA~1\QGIS 3.16\apps\Python37\lib\importlib\__init__.py", line 127, in import_module
return _bootstrap._gcd_import(name[level:], package, level)
File "<frozen importlib._bootstrap>", line 1006, in _gcd_import
File "<frozen importlib._bootstrap>", line 983, in _find_and_load
File "<frozen importlib._bootstrap>", line 967, in _find_and_load_unlocked
File "<frozen importlib._bootstrap>", line 670, in _load_unlocked
File "<frozen importlib._bootstrap>", line 583, in module_from_spec
File "<frozen importlib._bootstrap_external>", line 1043, in create_module
File "<frozen importlib._bootstrap>", line 219, in _call_with_frames_removed
ImportError: DLL load failed: The specified procedure could not be found.

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
File "C:\PROGRA~1\QGIS 3.16\apps\Python37\lib\runpy.py", line 193, in _run_module_as_main
"__main__", mod_spec)
File "C:\PROGRA~1\QGIS 3.16\apps\Python37\lib\runpy.py", line 85, in _run_code
exec(code, run_globals)
File "C:\PROGRA~1\QGIS 3.16\apps\Python37\Scripts\gdal2xyz.py", line 34, in <module>
from osgeo import gdal
File "C:\PROGRA~1\QGIS 3.16\apps\Python37\lib\site-packages\osgeo\__init__.py", line 32, in <module>
_gdal = swig_import_helper()
File "C:\PROGRA~1\QGIS 3.16\apps\Python37\lib\site-packages\osgeo\__init__.py", line 31, in swig_import_helper
return importlib.import_module('_gdal')
File "C:\PROGRA~1\QGIS 3.16\apps\Python37\lib\importlib\__init__.py", line 127, in import_module
return _bootstrap._gcd_import(name[level:], package, level)
ModuleNotFoundError: No module named '_gdal'

Execution completed in 0.16 seconds
Results:
{'OUTPUT': 'TEMPORARY_OUTPUT'}

Loading resulting layers
Algorithm 'gdal2xyz' finished


texts:
I'm rather new to python so my code might be horrid but I had a question about executing scripts from ArcMap and performance involving several if statements.

The main objective is to create a toolbox with several different possible executions depending on the users input.

Essentially I have two options that I see!?! --- note there might be errors as I am simply writing an example code that represents what I am doing. At the moment, iterating through only 5 (6000x6000cells) rasters my code runs in hours that in my mind should not take more than one.

import arcpy
arcpy.checkoutextension('spatial')

def main(inputfolder,outputfolder,x,y):
    inputfolder = arcpy.env.workspace
    RasterList = arcpy.ListRasters()

    for RasterImage in RasterList:
         try:
             if 'True' in x:
                  sa.Slope(.....)
             else:
                  if 'True' in y:
                       sa.Reclassify(....)
                  else:
                       RasterImage * 2
         except Exception e:
              arcpy.GetMessages(e)

or is defining separate functions better? i.e.

def main(inputfolder,outputfolder,x,y):
    inputfolder = arcpy.env.workspace
    arcpy.ListRasters()

    if 'True' in x:
        x(inputfolder, outputfolder, RasterList)
    elif 'True' in y:
        y(inputfolder, outputfolder, RasterList)
    else:
         for RasterImage in RasterList:
             try:
                 RasterImage *2

def x(inputfolder, outputfolder, RasterList)
    for RasterImage in RasterList:
         ....

def y(inputfolder, outputfolder, RasterList)
     .....

Any help would be most appreciated!

UPDATE

Solved my problem. Thanks for all the responses and help.

The problem was not with my code but rather that a few of the datasets that I had been given were corrupt. Using a try: except: I was able to pass through those corrupt files (write an error log of those files) and later delete them.


texts: Add a new field in the attribute table from derived feature

texts:
Is it possible to create a new field in the attribute table with the same value of the a specific derived feature?
For example: I would like to add "feature id" -from derived feature- to "ID_der" -attribute table-


texts: Creating artificial coordinate system to plot garden experiment?

texts:
I am trying to map plant death in a garden experiment. The location of the garden doesn't matter, I just care about each plant in relation to the others. They were all planted on a grid system. My data looks like this:

Column Row Plant State
1      1   3333  Alive
1      2   3334  Alive
1      3   3335  Dead
2      1   3336  Dead
2      2   3337  Alive
2      3   3338  Alive
3      1   3339  Alive
3      2   3340  Alive
3      3   3341  Dead

I want each plant to be a point on the grid with the column functioning as the X and the row as the Y coordinates. I'll then use symbology to change the point colors to reflect the state of the plant.

My questions:

Is this at all reasonable/feasible to do with GIS?
If it is, how do I create an artificial coordinate system that would work for my purposes?
If not, what program would you recommend (that would ideally take an .xlxs or .xls file)?

I can use either ArcGIS or QGIS so please give recommendations for which would be easier.


texts: Creating good quality picture marker symbols with ArcMap?

texts:
I have created a  nice clean symbol in Illustrator. After exporting to .emf format and appending to objects in ArcMap the symbol looks terrible (see Attachment).

Has can I import good vectors into ArcMap?

ArcMap Symbol

Illustrator Symbol


texts: How to get spatial reference of a shapefile using python?

texts: I always use arcpy.GetParameterAsText to get the spatial reference from the user. There should be a way to decide what is the spatial reference of a shapefile but I can't find it. Something like Describe....


texts: I'm a complete newbie at this so please excuse my lack of knowledge. I need to be able to place hundreds of thousands of points on a map. I'm using OpenLayers with Google Maps base layer and I just installed OpenGeo Suite. I also need to be able to add/delete points and manipulate them real-time. Placing those points as a vector or marker layer doesn't work so I think I need to have the server generate tiles? Where do I start with this? I've scoured the documentation and this website and can't find where to begin with generating tiles and doing it dynamically with my data. Thanks!


texts: Intersecting two shape problem using geopandas

texts:
I'm trying to intersect two polygons based on the following link,
intersecting two shapefiles from Python or command line

But problem occured that shows a "ValueError: Geometry column cannot contain mutiple geometry types when writing to file."

I'm new on these free open source GIS tool manipulation, can anyone give me some idea?
My source code is like this:

    from shapely.geometry import shape,Polygon,MultiPolygon,mapping
    import geopandas as gpd

    g1 = gpd.GeoDataFrame.from_file("./origin_test.shp")
    g2 = gpd.GeoDataFrame.from_file("./ref_test.shp")

    data=[]
    for index, orig in g1.iterrows():
        for index2, ref in g2.iterrows():
            if ref['geometry'].intersects(orig['geometry']):
               owdspd=orig['wdspd']
               data.append({'geometry':ref['geometry'].intersection(orig['geometry']),'wdspd':owdspd})

    data = data.set_geometry()
    df = gpd.GeoDataFrame(data,columns=['geometry','wdspd'])
    df.to_file('./intersection.shp')


texts:
I have a lot of mosaic datasets:

Zone8GK100k, Zone9GK100k, ... , Zone7GK25k, ... , Zone7Gk50K, ... , Zone7GK200K.

As I know, 7GK, 8GK... are utm zone numbers.

1) What does 25K, 50k, 100k, 200k mean?

2) What is the relation between these numbers and the ArcGis scales (1:1 000, 1:10 000, 1:24 000, 1:100 000, 1:250 000 ...)?

Update:

Thanks. I have layer Zone9GK200k. So, I set: "In beyond 1:200 000". What minimum scale(out beyond) for this layer do I set?


texts: Using geojson hosted on GitHub from Leaflet: it doesn't work

texts:
I'm trying to use a geojson hosted on GitHub (rif. https://raw.githubusercontent.com/cesaregerbino/DatiGeo/master/ScuoleTorino.geojson), from Leaflet.

here you are my code ....

    My first Leaflet map

<meta name="viewport" content="width=device-width, initial-scale=1.0">

<link rel="stylesheet" href="http://localhost/Leaflet/072/leaflet.css" />

The page doesn't show the markers of my geojson file: it seems that no return come from the request. Here you are the situation

If I put the same geojson file under my apache it works fine.

The new situation ...

I've tried to use the same geojson form QGIS desktopanche it works too.

The image ....

Any suggestions?


texts:
In geoserver is possible to generate a JSON or XML file with the information necessary to generate a legend in client?

I want to customize further the legend beyond the capabilities of GetLegendGraphic


texts:
I am trying to automate a process using python and arcpy. I have a scenario where a user will open an mxd and edit a polygon layer drawing a new polygon.
They will then want to use just this new polygon to undertake a select by location against every other layer (except the polygon itself obviously) in the table of contents.

I need help enabling my code to export one spreadsheet/csv of every layer that is within the search distance (the individual rows of data associated with each feature are not required, simply a yes or no). Below is the code I have so far.

#Imports the modules
import arcpy
#Sets the Map Document
mxd = arcpy.mapping.MapDocument("CURRENT")
#Sets the dataframe
df = arcpy.mapping.ListDataFrames(mxd)[0]
#Overwrites the Output
arcpy.env.overwriteOutput = True
#Output location for the new polygon feature class that will be used for the select by location.
Search_Area = "U:\\Test.gdb\\Search_Area"
#Selects the polygon with the MAX ObjectID from the Convey feature class.
arcpy.Select_analysis(Convey_L, Search_Area, "\"OBJECTID\" = (SELECT MAX( \"OBJECTID\" ) FROM Convey)")
#Defines the layers as everything in the table of contents.
layers = arcpy.mapping.ListLayers(mxd)
#Loops through the layers in the MXD and select everything that is within 200m of the Relevant Area (Search Area)
for layer in layers:
    if layer.name != Convey_L:
        arcpy.SelectLayerByLocation_management(layer, "INTERSECT", Convey_L, "200 meters", "ADD_TO_SELECTION", "NOT_INVERT")
for layer in layers:
     desc = arcpy.Describe(layer)
     if desc.fidSet:
         print layer.name, "has {0} features selected:".format(len(desc.fidSet.split(';')))
         for row in arcpy.da.SearchCursor(layer, "*"):
             print row
for lyr in arcpy.mapping.ListLayers(mxd,"",df):
    arcpy.SelectLayerByAttribute_management(lyr, "CLEAR_SELECTION")


texts: How to Save a Raster or Polygon to ArcMap Table of Contents in C#?

texts:
Just as the title says, I want to save the Raster to the table of contents.

I have a Raster displaying on the ArcMap IDisplay, but I can't figure out how to save it to the Table of Contents in ArcMap.  The save to TOC will come after a button is pressed in ArcMap.

Looking for code here.  C# of VB.NET will be fine!

Thanks all!


texts: javascript function to search all OpenLayers vector features covering a pixel

texts:
I have an OpenLayers map that draws features in a vector layer. The features are selectable and have a popup on select. Unfortunately in a lot of cases the features overlap so it can be impossible to select some features. I think that what I need to do to fix this is change my select control so that it uses a click handler and searches the map for features at this point. What sort of function do i need to write? are there any examples of this being implemented before?

This is how the features are drawn:

var vector_Layer = new OpenLayers.Layer.Vector();

function GetFeaturesFromKMLString (strKML) {
    var format = new OpenLayers.Format.KML({
    'internalProjection': new OpenLayers.Projection("EPSG:900913"),
    'extranalProjection': new OpenLayers.Projection("EPSG:4326")
    });
    return format.read(strKML);
};

vector_Layer.addFeatures(GetFeaturesFromKMLString('$newkml'));

And this is how Layers are currently selected:

var select = new OpenLayers.Control.SelectFeature(vector_Layer, {clickout: true});
        vector_Layer.events.on({
            "featureselected": onFeatureSelect,
            "featureunselected": onFeatureUnselect});

        map.addControl(select);
        select.activate();

        select.handlers['feature'].stopDown = false;
        select.handlers['feature'].stopUp = false;

here is the click event listener I was planning on using:

OpenLayers.Control.Click = OpenLayers.Class(OpenLayers.Control, {
            defaultHandlerOptions: {
                'single': true,
                'double': false,
                'pixelTolerance': 0,
                'stopSingle': false,
                'stopDouble': false
            },

            initialize: function(options) {
                this.handlerOptions = OpenLayers.Util.extend(
                    {}, this.defaultHandlerOptions
                );
                OpenLayers.Control.prototype.initialize.apply(
                    this, arguments
                );
                this.handler = new OpenLayers.Handler.Click(
                    this, {
                        'click': this.onClick,
                        'dblclick': this.onDblclick
                    }, this.handlerOptions
                );
            },

            onClick: function(evt) {
                //function that seachers for and selectas features at this point
            },


texts: Duplicating layer and applying filter using PyQGIS?

texts:
Is there a way, using Python, to duplicate a layer in QGIS as many times as the single occurrences of a variable contained in a field and with a filter applied on the duplicated layer based on the variable itself?

An example: I have layer0 with a field0 that contains var1, var2, var3, etc. I would like the layer0 duplicated to layer1, layer2, layer3 etc. so that layer1 has a filter where field0 = var1, layer2 has a filter where field0 = var2, etc.

The idea is similar to the LayersByField plugin but instead of creating a new layer for each occurrence of a variable I wish to just duplicate the "master" layer (so to avoid unnecessary files on disk)


texts:
I have a huge number of shape files that I receive on a recurring basis. I want to be able to compare the new version to the current version, and then only update the db with those that are different.

Currently, I am using arcpy.FeatureCompare_management to check Geometry & Attributes. However, this seems to be quite timely. I was considering if a byte-by-byte comparison of the .shp and .dbf files would work, but I run into discrepencies? This is what I have tried so far:

# These refer to the FQP of the Geometry and Attribute files for the same data-source
old_geo_file = r'\myfiles\q1\myshape.shp'
new_geo_file = r'\myfiles\q2\myshape.shp'

# These are only needed by filecmp. ArcPy grabs related files.
old_attr_file = r'\myfiles\q1\myshape.dbf'
new_attr_file = r'\myfiles\q2\myshape.dbf'

arcpy_geo = arcpy.FeatureCompare_management(new_geo_file, old_geo_file, 'MY_SORT_COL', 'GEOMETRY_ONLY')
arcpy_attr = arcpy.FeatureCompare_management(new_geo_file, old_geo_file, 'MY_SORT_COL', 'ATTRIBUTES_ONLY')

print arcpy_geo
print arcpy_attr

print('*' * 10)

# Clear the filecmp cache to ensure comparison
filecmp._cache = {}

filecmp_geo = filecmp.cmp(new_geo_file, old_geo_file, shallow=False)
filecmp_attr = filecmp.cmp(new_attr_file, old_attr_file, shallow=False)

print filecmp_geo
print filecmp_attr

This is the output:

True
False
**********
False
False

Can someone explain why a byte comparison would come out False, but ArcPy still considers them geometrically the same?


texts:
I have  a folder which contains (monthly) nc files from January 1980 to December 2015 with hourly data for different compounds (such as PM2.5, O3, NO2, NO, etc).The time is in UTC as follows:

int time(time) ;
               time:units = "seconds since 1979-12-02 00:00:00 UTC" ;
               time:long_name = "time" ;
               time:axis = "T" ;
               time:calendar = "standard" ;
               time:standard_name = "time" ;

My question is: Can I shift the time by + (or -) 1h, 2h,3h, etc, depending on the longitude, so I can get the local time? For example, for 15 degrees< longitude< 30 degrees the local time would be LT= UTC+1h. Any ideas?


texts: PyQt.QtSql - get database connection fail informations

texts:
Python 3 - PyQt4 - Qgis 2.14

I need to display message when my connection is success or fail with PyQt4.QtSql.

How retrieve message after db.connect() ?

Here is my code :

def testConnect(self):
        msg = QMessageBox()
        msg.setIcon(QMessageBox.Information)
        db = QSqlDatabase.addDatabase("QPSQL")
        db.setDatabaseName(str(dlg.txt_db.text()))
        db.setUserName(str(dlg.txt_user.text()))
        db.setPassword(str(dlg.txt_psw.text()))
        db.setPort(int(dlg.txt_port.text())) # int entry
        #db.open() - suggest to remove it
        if not db.open():
            msg.setText("Connection failed !")
            retval = msg.exec_()
        else:
            msg.setText("Connection success !")
            retval = msg.exec_()

Where is information to change :

  msg.setText("Connection failed !")

by real failed connection information from PyQ4.QtSql and debug ?


texts: Password does not work on custom install of CartoDB

texts:
I got an installation of CartoDB working from a recent git pull; however, I cannot login.  I am using the domain:

development.localhost.lan

and username 'development'

The password is the one that I gave to create_dev_user.  I checked that there is an entry in the database for this user but whatever I do I cannot get it to login.

Any suggestions for what to do to get it to work?


texts: Vector/Processing ToolBar not available QGIS 3.4/ Mac OSx 10.13.6

texts: When I installed QGIS, I performed the bug fix as directed (going into the Environment section, adding a variable, etc), and restarted both the program and eventually my computer, but there is still no Vector or Processing toolbox available. I've also uninstalled and reinstalled the program.  Any other ideas? We have to use this version for the class I'm taking.


texts: Querying Additional Styles from Geoserver

texts:
I'm using GeoServer 2.1.2, and I've got a layer to which I've applied added additional WMS styles. However, I'm having trouble working out how to query the layer and retrieve the list of available styles (for use in my web app).

Any ideas?

I haven't been able to find anything in the GeoServer docs, and I'm still fairly new to the whole WMS/Geoserver thing so I'm not sure what I should be looking for on Google.


texts:
I'm using QGIS 2.0.1 and want to send a GPX file to a GPS device, Garmin Astro 320.

When using the GPS Tool in QGIS and click the OK button I get the message

  Can not start GPSBabel!

What's wrong?


texts:
Just wanted to know if it is possible to define the direction of taken photo?
I took a lot of pictures with QFIELD (at my Phon) and I would like to know which direction has each Photo when I was taking it.
I only have x and y as attribute in a Point.shp.


texts: I'm working on a custom addin with C# + VS2010 + Arc10.1 Pre Release and I don't seem to find the way to make custom icons work for a button that the extension holds and for the extension itself. I pointed the Config.esriaddin appropiate attributes to the files which are both png with alpha transparencies (the icon is 16x16 and the extension img is 64x64). Is there any rule about creating images for icons I'm missing? Why would this not work? Thanks in advance.


texts: How can I make the background of the Leaflet map transparent.  I see opacity controls from the layers, but not the map background.  Should the be changed just with plain CSS?


texts: Have you all solve importing Z point from dxf to QGIS?

texts:
I have read all the answer about importing dxf to QGIS, but I still not understand yet.

Does anyone have a video link (youtube or another) that show me step by step to export dxf with Z point to QGIS?


texts: Where to download Mars DEMs and landcover?

texts: I'm looking on the websites of the Mars Global Surveyor and the NASA sites, but all of the data is in jpeg or img form. Does anyone know where to find useful raster data sets for elevation and land cover (soil types)? I'm sure the information is sparse, since the entire planet has not been covered.


texts: Strange behavior of of ST_Intersects in PostGIS


texts:
I have a couple of Landsat 7 images and they all have striping as a result of the SLC failure. Does anyone have any suggestions as to how to destripe the image using ENVI? I plan to use these images for NDVI analysis.

When I downloaded the data from GloVis, it also came with a Gap Mask folder...is that somehow related to the SLC-Off?


texts: Collecting data in related tables using Qfield?

texts:
recently I have started to use ArcGIS collector to collect data in the field. In collector I am able to collect a point feature or other feature class and then add entries to the related table that I have set up using a 1 to many relationship class. This allows me to collect tree entries while timber cruising or collecting data for CFI plots fairly easily.

That being said I have started to explore using QField for this work because I do not want to be stuck in the ArcGIS Online ecosystem and I like the functionality of the QGIS's Attribute form. The problem that I am running into, is not being able to collect the data on the related table that I have set up. In QGIS I am able to create a point, enter the data for that point and then add entries in the related table using the form that I have set up. Is this functionality present within QField? If it is, how can I set it up so that I can collect this related data?


texts: How to view .dwg group layers in ArcMap?

texts:
I'm trying to view 3 .dwg group layers in ArcMap 10.1.  Whenever I add them to the map, I get an error message that says they are missing a spatial reference.  That's fine.  I'm just trying to view them, but I can't seem to figure out how.

Is there an extension that needs to be enabled?  The .dwg files contain contour data, if that makes a difference.


texts: How can I make it that the font size of the label changes when I zoom in and out? I have a polygon and I want the label get bigger when I zoom in and smaller when I zoom out (just like the polygon does). I tried to use map points but it didn't work.


texts: Convert raster to a CSV values and lat long based on raster value in ArcGIS

texts:
I am looking to convert a raster into a csv file (or tab delimited) based on values in the raster (Value > 0).

The csv would have three values, latitude of center of cell, longitude of center of cell, and the raster value where the raster value is greater than 0.

The data in unprojected in geographic coordinate system WGS 1984. It is a GTiff.

I know I can go to point file and then add the latitude and longitude and then go to CSV but I wonder if a direct method exists as I have 130 million points.

Working on going Ras > Point (file Geodatabase) > Query Points > CSV in ArcGIS but it is slow.

I feel I need to build in the query at the Ras > Point section to speed things up as the entire raster has 130 billion cells but only 130 million have values. Rest are 0 or NoData.

Looking for ArcGIS solution.


texts:
I am looking to project all my data into EPSG::2959 (NAD83 CSRS UTM 18N) as I am working out of Ontario, Canada. Every time I start a project I set the project CRS to EPSG:2959, and then begin to add data. I typically use OSM from "QuickMapServices" (QGIS plugin) as a basemap for all projects, and change the CRS of OSM to EPSG:2959.

At this point I begin adding data, which typically changes the project CRS to EPSG::3857 or EPSG::4326. I typically change the CRS of both the project (again) and layers at this point. I often 'Save As' the data to the correct projection.

My issue now comes when I use the measure tool after everything is in EPSG::2959 (NAD83 CSRS UTM 18N). The distance measured is always incorrect.

How can I fix this to accurately measure distances in EPSG::2959?

I will take any advice or change any steps in order to fix this issue.

This may have been asked before but it's causing me an incredible amount of trouble and I can't seem to find a solution.


texts: Getting nearest point but with a filter conditional in QGIS 3.4?

texts: I have two point layers in CSV files, each with about 500 points/rows. Each point in both layers has a year column (attribute). I know how to use Distance Matrix to find the nearest point. BUT: For each point in Layer 1, I want to find the nearest point from Layer 2 that has one of the three previous years to the point in Layer 1. For example, if point x in Layer 1 has 1989 in the Year column, then it should search for nearest point only out of all points in Layer 2 that have 1986, 1987, 1988 in the Year column, and return the nearest one.


texts:
I`d like to find out if and how it is possible to visualize LoD2 data (source:https://www.opengeodata.nrw.de/produkte/geobasis/3d-gm/) in QGIS 3.4/Madeira.
So far I managed to load them into the project, but unlike other examples and tutorials I´ve seen, they come as lines instead of polygons.  After enabling 3D for the layer and opening the 3D-View, it only "kind-of" worked  with a second set of lines flying above the footprints - see screenshot (3D-View without background-map for a clearer image). So - unless there are issues with the data - I suppose I missed something...

I also tried to convert them into polygons - first as original data and then after saving them as Esri-shapes, but the result was always that even the "polygons" came as "lines":

Does anybody know if these data can be visualized with QGIS and if/how I have to prepare them or which steps in QGIS are necessary?


texts: Calculating a new attribute for a polygon based on its corners

texts:
i have a pointlayer with a real number as attribute named amount.
I want to calculate and visualize the difference of amount between neighbourhood points. So first I made a delaunay triangulation where the edges represent the neighbourhood relationship.

Now my problem is to calculate the mentioned
difference for each triangle edge (e.g. as a new attribute for this edge) of the delaunay triangels.

Because I didn't find out how to split each polygon (triangle) to 3 seperate lines (edges), it would be also okay for me to make a defined calculation¹ for each polygon (triangle) based by the attribut amount of the 3 points/nodes/corners around.

I hope you understand what I want to do and can help me to solve this in QGIS

¹ Calculation: New triangle-polygon attribute = (NodeA-NodeB)-(NodeA-NodeC)-(NodeB-NodeC)

where NodeA is amount of node A etc...


texts: Calculate Feature Area in Shapefile, which crs to use?

texts: For a global shapefile I would like to calculate the area per feature. If I understand correctly it is best practice to reproject your shapefile to a projected coordinate system. In the past I've used World Eckert IV. Are there any other projections that are suitable or even better for global area calculation, ideally with an epsg code?


texts: SQL Error in displaying long/lat from .csv

texts:
I tried to display long/lat (°) as points from a .csv but an

ArcMap drawing Error occurs:

  One or more layers failed to draw:

  2016-07-29_GPS_Collar09743_20160729125522.csv Events:  An invalid SQL
  statement was used.

My Data look like:

And this is the Display XY Data Window

What am I doing wrong if i want it to draw in WGS 84 ?

After the correction suggested by JGH it looks like this:

And if I use the Identify tool:


texts: How to display a thumbnail of a geotagged photo in QGIS?

texts:
I am using Qgis 3.10 on windows. I have created a shapefile of geotagged photos using the 'Import photo's' plugin. This gives me a map with an icon at the location of each photo. The attribute table has a filepath to the location of the photos and a relative filepath.

I want to replace that icon with a thumbnail of the photo. I have found some questions about this but none of the answers tell me fully how to do it step by step.
I have changed the symbology to 'Raster image marker'. I think I need to use 'Data defined override' to point to the image location somehow, but can't work out how.

I need a simple method that can be used by lots of people in an organisation that only occassionaly use GIS.


texts: Plugins to guide inexperienced users: updating shapefiles and managing data

texts: I am attempting to find a way to block adding features while certain attributes are not yet added. This would be useful in order to guide inexperienced clients and employees with the process of adding entities and updating shapefiles while managing data. For example, add feature in edition mode disables 'Okay' button until certain attributes are filled. Any suggestions on either a plug-in that exists or where to start for programming one?


texts: Reprojecting raster file with unknown datum using ArcGIS Desktop?

texts:
I am new to GIS. I had a grib format data file and it's coordinate system shown below. I want to project it to WGS84, however, the datum of original file is unknown. I know I need to provide the datum info while projecting.

How can I solve this using ArcGIS Desktop?


texts: Spatial join based on polygon centroid in QGIS

texts: I have 2 cleaned shapefiles, one containing building footprints and one land parcel shapes, both have attributes. I want to do a spatial join when data from the land parcel is merged into the footprint layer. However, as you can see, some footprints transcend across several parcels. As the data in the parcels are roughly similarly, I would like to do the join based on the building polygon's centroid, so that each footprint has data from 1 parcel shape. I know ArcGIS is able to do this with spatial join tool, anyone knows if QGIS can do this too? I looked at "join by location" but it doesn't seem to have the centroid function.


texts: Loss of geometry in shapefile after saving edits

texts:
I'm using QGIS 2.18.17

When I'm editing (e.g. removing objects) from a shapefile in QGIS and click to save the edits, all geometries are lost. All non-deleted objects are still present, but they do no longer contain a geometry (no objects visible in the layer and, when copying the objects from the attribute table in an Excel-document, the column "wkt_geom" is filled with NULL for all objects).

At first, I thought the problem only occurred with a point layer (shapefile), but today it also occurred with a polygon layer (shapefile).

The same thing occurred when I wanted to save a selection from a point layer as a new shapefile. The shapefile is created and the dbf-table is complete, but no geometries are present.

As I cannot imagine it concerns a bug in the used version of QGIS (as it is such a basic functionality), I assume the problem arises when saving the file with extension .shp (which stores the geometry) of the shapefile.

I'm working at a local government and the environment which I'm working on is a thin client (access through a Citrix viewer to a remote server) where a desktop environment is created which is identical for all users within the local government. All data is also stored on the same remote servers.

Anyone who had a similar problem in the past and found a way to solve it?

Edit: further investigation indicated that the cause eventually could be found in the way QGIS handles the network folders (on a Windows 2012 R2 Server), as the problem does not occur when writing on the local disk. Could it be something with some difference between relative and absolute paths? I already tried to switch this setting (without closing QGIS) but without result.


texts:
I need to run a script tool from ArcCatalog. It needs to access a feature class in an ArcSDE database so I can link the script to the layer as a source layer. I have connected the folder the script's toolbox is in, but the script still can't find the database. How can I connect the script to the database?

Normally I just paste the path and filename into the script tool's dialog, but that's not working. I have a script for shapefiles and a script for FC's that does the same thing. I will add a screenshot of this to the original post. The script for FC's connects to the database, but the one for shp's does not. We switched our server to a new company (MS Azure) a few months ago.

I just created the tool for shapefiles from the script for the feature classes today. I only changed one line that has nothing to do with accessing the source layer. So the other tool existed before the conversion to the new server, but this one was created after the conversion.

The script tool for the FC's has not given me this problem and this is the first time I've had to use the script for shp's.

Tool for FC's:

Tool for shp's:


texts:
Comparing results of zonal statistics in QGIS and R, I noticed that the results they both yield do somehow differ. Manually checking, I also found that both QGIS and R seem to partly produce erroneous results. I wonder whether there is a robust solution to calculating mean values of raster cells within polygons.

My use case:
I have a Raster Layer with cell size = 1 x 1 km² and a Polygon Layer with quadratic polygons of size = 250 x 250 m². Both are projected in EPSG:5677.

I used raster -> zonal stats with the "mean" option in QGIS ("QGIS.Result") and tried

R.Result1<-extract(Raster, Polygons, method="simple", fun=mean)

as well as

R.Result2<-extract(Raster, Polygons, method="simple", fun=mean, weights=TRUE, normalizeWeights=TRUE)

in R with the raster package.
The results look like this:

ID   QGIS.Result  R.Result1    R.Result2
1   -63.9505964  -64.0000000  -64.60000000
2   -67.0000000  -67.0000000  -67.00000000
3   -67.0000000  -67.0000000  -67.00000000
4   -62.3765259  -64.0000000  -61.85714286
...  ...          ...          ...
19  -99.9498123  -97.5999908  -99.19999186
20  -92.7999878  -92.7999878  -92.79998779
21  -92.7999878  -92.7999878  -92.79998779
22  -102.5493486 -92.2999954  -99.19999186
23  -97.5452805  -89.3999939  -92.79998779
24  -97.5452805  -89.3999939  -92.79998779
25  -87.4895458  -87.0000000  -87.33333333
26  -86.0000000  -86.0000000  -86.00000000
27  -86.0000000  -86.0000000  -86.00000000
...  ...          ...          ...

Results differ whenever a Polygon intersects with more than one raster cell - note in the example below that:

all three methods yield the same result for polygons 20,21,26,27 which all do intersect with only one raster cell and different results for the other polygons which do intersect with more than one raster cell.
the QGIS result gives a smaller value for polygon 22 than for polygon 19 although it should be the other way round. R.Result2 gives the same result for both polygons. Only R.Result1 performs better here.
polygon 25 intersects with Raster values 88 and 86. R.Result1 is the mean of both: 87. However, as polygon 25 contains a larger amount of value 88 than of value 86, the result should be biased towards 88 (like in QGIS.Result and R.Result2).

Thus, all three methods seem not to provide what I need: a robust, coverage-based mean of the raster layer for each polygon. Is there a way to obtain such (preferably in R)?

The used data can be found here:
https://www.dropbox.com/s/bctjidc4dz2j396/testdata.7z?dl=0


texts: OpenLayers 4.2 - Changing the interaction geometry type from polygon to multipolygon

texts:
I would like to make a selection with DragBox interaction and get the geometry type to Multipolygon.

Right now I'm getting a Polygon type back :

var source = new ol.source.Vector({wrapX: false});

var map = new ol.Map({
    target: 'map',
    layers: [
        new ol.layer.Tile({
            source: new ol.source.OSM()
        }),
        new ol.layer.Vector({
            source: source
        })
    ],
    view: new ol.View({
        center: ol.proj.fromLonLat([-3, 48]),
        zoom: 7
    }),
});

var select = new ol.interaction.Select();

var dragBox = new ol.interaction.DragBox({
    condition: ol.events.condition.platformModifierKeyOnly });

dragBox.on('boxend', function() {
    var geom = dragBox.getGeometry();
    var format = new ol.format.WKT();
    var wktRepresentation  = format.writeGeometry(geom);

    console.log("wktRepresentation: " + wktRepresentation);
});

Can someone help me?


texts: ARCGIS10: updatelayer results in error - LayerObject: Get attribute renderer does not exist

texts:
I'm quite new in using python in Arcgis 10.
I wanted to use the symbology of a geostatistical layer in the TOC in another geostatistical layer in the TOC. Therefore, one should use arcpy.mapping.updatelayer.
I saved the geostatistical layer as Layer file (TensorStyle.lyr).
The code to copy the symbology looks like

import arcpy

mxd = arcpy.mapping.MapDocument("CURRENT")

df = arcpy.mapping.ListDataFrames(mxd)[0]

updateLayer = arcpy.mapping.ListLayers(mxd, "*", df)

sourceLayer = arcpy.mapping.Layer("TensorStyle.lyr")

for lyr in updateLayer:

    arcpy.mapping.UpdateLayer(df, lyr, sourceLayer)

If I run the script I get the message:

  Runtime error : LayerObject: Get
  attribute renderer does not exist

Could you please help me.
Thank you in advance,

Ronny


texts: Find the most populous city within regular hexagon

texts: I have a map of europe as hexagonal grid. For each hexagon I have center(lat and long) and constant size. Is there a way to determine the most populous city within each hexagon? It may not be completely accurate so each city should be just point with a coordinates and no real area.


texts: Seeking alternatives to Business Analyst, AnySite/MapInfo, and Sites To Do Business?

texts:
This is a very general question and I hope it is ok. I've been "Googling" and can't find much. I work at a shop that currently uses MapInfo + AnySite to do real estate site analysis. I generate demographic reports on US zip codes and by radii (usually 1, 2, and 3 miles).Se

The workflow is either:

Select zip codes -> run "hot spot" analysis (which is an automated
z-score comparison within the software) and export table of data for
selected zip codes (in a preformatted table).
Create point for site -> draw desired radii -> generate report based
on these radii.

As AnySite and MapInfo cost us a hefty sum (even more than ArcGIS would), I'm tasked with finding an alternative. So far I have found ArcMap and Sites To Do Business (which I really hope we avoid).

Do you know of any others or of search terms I could use?

I'm stuck beyond the two I have already found.


texts: Creating offline basemap with small size using ArcGIS SDK for Android?

texts:
I want to develop an app for mobile devices .I'm using ArcGIS SDK for android and i want to have an offline basemap in the app . One solution is use to tpk file.but the size of tpk is too large for the app ( 70 megabyte) .

Is there any way to create an offline map with small size ?

Is it possible to use vector base map with high performance ?

Is there any API that i can to use offline map with small size ?


texts: how to creat a histogram from two raster layer?

texts:
I have two raster layer ( raster layer 1 and raster layer 2) that was made by kernel density. now, i want to create a graph that vertical axis shows the density of each cell in raster layer 1 and horizontal axis shows the density of each cell in raster layer 2.

How can i do that with arcgis?


texts: Get all related attribute values from features in list (QGIS)

texts:
According to this question: Get Attributes from 1:n relation in QGIS as Expression to use in Atlas I try to bulid waysigns for a hiking-GIS using the QGIS Atlas. The user should be able, to choose the destinations on the way sign.

I need to get only the attribute values from the relation "dist", which are written in a list ("wanted_dest"). As you can see at the screenshot below, at the time the way signs are labeled with all their destinations using this relation_Aggregate-expression:

relation_Aggregate(
relation:='dist',
aggregate:='concatenate',
expression:=concat("destination_id", ': ',round("network_cost"/1000,1), ' km  '), concatenator:=',')

But I don't want them all - I just wan't the destinations, of every list item (Place D, Place E and Place F) so I need to modify the aggregate or use a virtual layer. But I can't find a working solution.

Here is the used project for testing: http://www.map-site.de/tmp/WaySigns.zip


texts: OpenLayers adding circle feature with same vectorSource increases opacity of all

texts:
Background

Specs

OpenLayers 4.4.1
OSM

I'm fairly new to OpenLayers and have never used vectors before (primarily because I found out that I was using OpenLayers version 1, and had to relearn everything).

My application adds circles to a map relating to a position with a specific radius indicating position accuracy.

In its operation, multiple circles are added to the map at different times.

This is my code for loading the map:

var map = new ol.Map({
    layers: [
      new ol.layer.Tile({
        source: new ol.source.OSM()
      })
    ],
    target: 'mapdiv',
    controls: ol.control.defaults({
      attributionOptions: /** @type {olx.control.AttributionOptions} */ ({
        collapsible: false
      })
    }),
    view: new ol.View({
      //center: [0, 0],
      zoom: 16
    })
  });

  //this is where all map 'features' (circles) are stored
  var vectorSource = new ol.source.Vector({
   projection: 'EPSG:4326'
  });

As you can see, I load the 'vector source' right after the map as I understood that it holds all 'vectors' which are displayed on the map so long as you specify it as the 'source'.

This is the code I use to generate the circle (source) (I tweaked it at getPointResolution because the OP made a mistake):

  //code from https://stackoverflow.com/a/28299599
  function addCircle(map, vectorSource, radius) {
    var view = map.getView();
    var projection = view.getProjection();
    var resolutionAtEquator = view.getResolution();
    var center = view.getCenter();
    var pointResolution = ol.proj.getPointResolution(projection, resolutionAtEquator, center);
    var resolutionFactor = resolutionAtEquator/pointResolution;
    var radius = (radius / ol.proj.METERS_PER_UNIT.m) * resolutionFactor;

    var circle = new ol.geom.Circle(center, radius);
    var circleFeature = new ol.Feature(circle);

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

    map.addLayer(vectorLayer);
  }

Problem

Loading one circle goes normally, adds a blue stroked, opaque circle at the specified location with specified radius.

Loading a second circle appears more opaque than the last. Moving the map to the previous circle, it is also more opaque.

With each added circle, the apparent opacity increases for all displayed circles.

Running vectorLayer.getOpacity() in every circle generation results in 1, when clearly the circle is translucent, becoming increasingly opaque with every new circle.

Summary

Looking around, it appears that often it is the case that the developer is reloading the same circle over and over until many are stacked on top of one another. It almost seems like this is the case for me too, except I've triple-checked that I'm only running addCircle() once and the circle is in a different position than the last.

Is it possible that OpenLayers is redrawing all previous circles with every new circle?

Maybe this isn't related to getOpacity but has to do with the color as an rgba() combination...

Question

I want every circle to remain the same after drawing new circles. The default opacity and color is fine.

Am I doing something wrong?

UPDATE: Here's a fiddle as an example - https://jsfiddle.net/f5zrLt20/5/


texts: Find one centroid instead of multiple centroids?

texts:
I'm trying to find the centroid (center of mass) for several IUCN species range maps. Using the code below I've run into an issue in R.

library(rgdal)
library(rgeos)
shp <- readOGR("/Documents/data_1/data_0.shp")
class(shp)
plot(shp)
writeOGR(shp, ".", "shp", driver="ESRI Shapefile")
cents <- coordinates(shp)
cents <- SpatialPointsDataFrame(coords=cents, data=shp@data,
                                proj4string=CRS("+proj=longlat +datum=WGS84"))
points(cents, col = "Blue")
writeOGR(cents, ".", "cents", driver="ESRI Shapefile")
trueCentroids = gCentroid(shp,byid=TRUE)
points(coordinates(shp),pch=1)
points(trueCentroids,pch=19,col = "Red", cex = 1.8)

This is one of the resulting maps

The problem is that for some of the polygon maps I'm finding multiple centroids [as shown above in red] instead of just one (same result in QGIS). So, I'm wondering if there is a way to change my coding and find a single center of mass for the polygon?

I'd appreciate any help and/or explanation for my problem. Thanks.


texts:
I am trying to create a grid to use for zonal stats for a ag research plot. I have created some grids and rotated them and am getting decent fittings but I feel that there must be some other better way to fit the grid to the plot. Is there some way to define the boundaries of the plot with points and to fit the grid to those points or something along those lines (no pun intended)? What is the best way to create a grid for this type of analysis? Also and more importantly, I am thinking that I would like some empty space in-between the boundaries pf each grid cell, how would I create a such a grid?

I am working in qgis with imagery processed in pix4d. The plots are generally rectangular.


texts:
I am trying to classify 9 different classes of Land use deploying Resnet Architecture. Sentinel images come with 13 bands and in RS we normally combine different bands for different classes for example for agriculture color bands (11,8,3) are composites for agriculture, and so forth with RGB, True Color, NDVI.

I have tried combining different bands and so far accuracy is 0.83 with RGB. I tried SWIR and True Color combinations and both have lower scores.

In addition, I tried to have different and appropriate bands for each class, for example, Residential(Bands 4,3,2), River(8,11,4), and Crops(11,8,2). Well, this has improved the accuracy of the model up to 0.93.

My question is twofold:

Is this really an acceptable approach in Deep Learning practice? All the literature I have seen so far are using the same bands for each class.
If this approach is ok, then can I use different band combinations for the same class as data augmentation? this will increase data size and hence improve the accuracy.


texts: How to get drive distance polygon with ServiceAreaTask?

texts:
I'm using an ArcGIS Online ServiceAreaTask to calculate drive-distance polygons.
In the ServiceAreaParameters that I pass to ServiceAreaTask.solve(), there is an attribute named impedanceAttribute.
In the API reference, it's explained that this attribute can take "Time" or "Length" as value: https://developers.arcgis.com/javascript/jsapi/serviceareaparameters-amd.html#impedanceattribute

When I switch between Time and Length, I get exacly the same polygons. I guess that impedence attribute isn't taken into consideration by the task.

I'm using ArcGIS Javascript API and ArcGIS Online services. So, I have no access to the Network.

Am I doing something wrong?


texts: what version of jetty the geoserver runs?

texts:
In order to set CORS correctly in Geoserver, I have to download some jar files and then edit the web.xml file.

I have a really old version of Geoserver(2.1.3) that runs Jetty 6.8.1 and there are no jars for this version, at least according to this. I was thinking to upgrading to the latest Geoserver (2.8.0). What version of Jetty does it run? Will I find the right jars to implement the CORS setting?


texts: Attribute-spliter to parse a text file containing key=values into separate variables in FME

texts:
Context
Using FME, I want to read a text file containing several lines of variables I already defined, given in the following format:
var_name = "var value"

There are many lines, each one contains one single variable with one single value:
var_name0 = "var value0"
var_name1 = "https://path-to-API/v1"
var_name2 = "https://path-to-other-API/v1"
trigger1 = 1
trigger1 = 0
path_shapefiles = "/path/to/some/shapefiles_folder/"
path_geopackages = "/path/to/some/geopackages_folder/"
...

For the moment, I read the text file using an AttributeFileReader to which I connect an AttributeSplitter.
My results come as one single concatenated line in the output, with absolutely no separators which is obviously not what I want, such as:
var_name0 = var value0var_name1 = "https://path-to-API/v1"var_name2 = "https://path-to-other-API/v1"...
In the "Delimiter or Format String" field of the AttributeSplitter, I tried inserting \n \r or\r\n without any success:

I still got one unique row with all attributes concatenated together.
I also tried inserting a single newline by pressing one time the Enter key in the "Editor" as discussed here (2013):
https://knowledge.safe.com/questions/2361/use-attributesplitter-with-a-carriage-return.html without much success.
Question
Do you know what I am missing or what I'm doing wrong with that? Or any kind of workaround to end up with a table which fields are the variables in the value along with their value on the first row for example?
Environment
I'm using FME(R) 2020.0.0.1 (20200316 - Build 20202 - linux-x64) on Ubuntu 18.04.


texts: Defining geometry types (point, polygon, geometry) in PostGIS?

texts:
I am new to PostGIS. I have checked through a number of articles and posts on this, but still not finding a definite answer.

I am setting up a database for geo tracking purposes, which tracks assets (people, vehicles, equipment, etc).  In two of the source systems, geo objects are defined as points, circles, polygons and rectangles.

One of my tables (a generic object table) will have a mix of all types.  I can define the circles, rectangles and polygons as geometry(polygon, 4326) (for example), but what about the points?

Can I use geometry(geometry, 4326) for all geo types (polygons and points), or should I have 2 columns, one defined as geometry(polygon, 4326) and the other defined as geometry(point, 4326)?

To further complicate matters, future requirements could lead to linestrings being stored in the same table.

Two issues concern me: performance and future proofing the design.

Anyone with hard experience on this one?


texts: How to extract values from a raster using R?

texts:
I'm a beginner at using R software, and I'm attempting to perform a percent linear stretch on the following raster data:

> KBnew

#class: RasterLayer
#dimensions  : 27018, 26389, 712978002  (nrow, ncol, ncell)
#resolution  : 15, 15  (x, y)
#extent  : 541792.5, 937627.5, 5454772, 5860042 (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
#data source : /private/var/folders/sy/qh_q0tbs29l_j9wfndb6cw5r0000gp/T/RtmpDc7N7c/raster/r_tmp_2016-02-15_212425_11368_77493.grd
#names       : layer
#values      : 7.712379, 205.2743  (min, max)

... although (if I'm doing this correctly) from what I understand, I need to know the mean and standard deviation of the raster values (which would be contained within the range of 7.712379,205.2743 or min/max values). How do I extract these values, within the minimum/maximum range?


texts:
I'm working on a shapefile of Brussels (QGIS 3). This shapefile contains all the streets in Brussels, see below.

With the attribute table, I have an "X" and "Y" position for each line (debut and end),

I used this to create nodes at each debut and end of the line in Oracle Spatial.

Now I need to check if

all the nodes are not in double
They are all connected to another
There isn't one node so close to another that a travelling test would take the wrong way

EDIT :

For now I'm just Using the network analasys from Processing Tools, and I take one node on the left and one on the right (then N and S etc) and ask for the shortest path to see if it's OK. But There should be an easiest way that provides to do an analysis on all the Network !

==> SO I'm looking for
- a way to calculate a lot of shortest path easily (maybe with a script) or
- another way to see if each point is really connected to the rest.

Maybe I could take a point in the center and test it with all of the point on the external (like the boarder of my shp?)


texts:
I'm stack with a simple issue.

I have a raster (DEM) and what I want to do is to assign -9999 where the pixel value of the DEM is 500 and leave all the other values as the original ones.

I tried:

 (  ( "dem@1" = 500)  * (-9999)) OR  (  ( "dem@1" != 500)  * "dem@1" )

but I get a new raster with nan and 1.

Some idea?


texts: What is generally faster in OpenLayers: WMS or GeoJSON?

texts:
I work both with WMS (connection with postgres) and GeoJSON in OpenLayers. I have the impression that a WMS is faster when more feature need to be loaded. And I have the impression that GeoJSON is fast when there are only a small amount of feature. When more features need to load the application is much slower with GeoJSON

What is generally faster in OpenLayers?

And is there way to make a GeoJSON file load faster.


texts: Color different of two Sentinel 2A (L1C) images from the same sensing date

texts:
i'm new to this GIS Q&A websites and remote sensing as well, so sorry for this too basic question. It's about my Sentinel 2 (Level-1C) images. I downloaded 2 images with the same sensing date (09-09-2018). One is next to the other, meaning the first image is just one path away (like in Landsat) from the second one. However, some portion of the first image overlaps the second one.

I made a natural color composite out of it and thought that the color would be perfectly the same. But turns out, they're slightly different. I need to know what seems to be the cause of this color different. Is it just the visualization (histogram) or it goes back to the spectral value? Since i'm willing to do pixel based image classification to these image, this matter would be so significant to the result.

Additional information:

no correction (geometric and radiometric) has been applied to the images
it's not even a mosaic image, they're two different natural color raster i simply visualize in qgis


texts: Making radiometric calibration on Planet Image?

texts: How can make radiometric calibration on Planet Image? I tried to on ENVI, but I did not succeed. Cause software need gain and offset value for every band. But I did not find this value. Do you have any solution this matter?


texts: HRESULT 0x80010202 exception in ArcMap closest feature searching by ArcObjects?

texts:
My method must return closest feature attribute. But there's a problem, it finds closest feature but returns COM exception HRESULT 0x80010202

void selectStreet(IPoint point, IActiveView activeView)
{

..........

if (getLayerByName._GetLayerByName(settings.StreetsLayerName) != null) // Settings1 settings = new Settings1(); settings.StreetsLayerName = "Street"
            {
                var foundLayer = getLayerByName._GetLayerByName(settings.StreetsLayerName);
                IQueryFilterClass qFilter = new IQueryFilterClass () { WhereClause = string.Format("{0} IS NOT NULL", settings.StreetsName) };
                IFeatureCursor featureCursor = (foundLayer as IGeoFeatureLayer).Search(qFilter, false);
                while ((feature = featureCursor.NextFeature()) != null)
                {
                    double distance = ((IProximityOperator)spatialFilter.Geometry).ReturnDistance(feature.Shape);
                    IProximityOperator prox = spatialFilter.Geometry as IProximityOperator;
                    distance = prox.ReturnDistance(feature.Shape);

                    if (closesToId == -1 || distance < closestDistance)
                    {
                        closesToId = feature.OID;
                        closestDistance = distance;
                        resultFeature = feature;
                        var one = feature.get_Value(feature.Fields.FindField(settings.StreetsName)).ToString(); // StreetsName = "STR_NAME"
                        var two = feature.get_Value(feature.Fields.FindField(settings.StreetsStreetCode)).ToString(); //StreetsStreetCode = "STR_CODE"
                    }
                 }
                }
            }
            catch(Exception e)
            {
                MessageBox.Show(e.Message); // HRESULT 0x80010202 exception
            }
            return null;
        }
            Where could be a problem?


texts: Working with large datasets with no server

texts:
I am in a bit of a weird work environment. I'm a data-scientist type person who is now doing more and more GIS work on some very large datasets - spatial regression, some mapping, spatial correlation, etc.
My data-engineer team are quite old-school and so instead of querying our spatial data directly, they only give us access to pre-processed shapefiles (most of my organisation is on Esri) on a network drive, to avoid storing data locally.  Because of that, I tend to work in Python and have little SQL experience.
Managing large datasets with these is a pain - I tend to iterate over them in GeoPandas, and do most of my data manipulation in Python, which seems slow. I suspect ideally I'd do my initial querying in a database. I'm trying to find ways of improving that, but have some limitations!
I suspect the ideal solution would be a PostGIS SQL server, but we've been told that's impossible for now.  I've considered a few options below...

Spin up a local PostGIS server on my machine?
aggregate those files into a SQLite/Spatialite local database
aggregate them into geopackage db

Do these make sense?  Should I be looking at something else? Do any of those sound mad?


texts: How to add projected point coordinates to attribute table?

texts:
I want to get the projected coordinates of each point into my attribute table? The software which I want to use can only compute Cartesian/projected files, and not data in lat long or any other kind of degrees and minutes situation.

I have XY Coordinates as Lat Long in my attribute table. I wish to obtain the projected coordinates for all of my points in meters as shown in the status bar in the bottom right corner. They are also identified as 'location' when i use the info button to click onto each point.


texts: Does the OSM provide data for street names in different languages? E.g. Spanish and English. If yes, how can I download them?


texts: GIS analysis to find duplicate geometries

texts:
I have a big shapefile that contains all the buildings and houses of the town that I work in (approx. 90,000 features). The data of the buildings/houses are saved by the town's surveying engineers and due to bad practice and the access of different surveyors to that data, many buildings/houses have been saved twice and show up in the map as duplicates.

Some of them are exactly duplicated (they appear one over the other) while others are duplicated with a space between the two objects (like one object is inside the other - see the attached screen shot).

I want to clean that data so that I have only the correct buildings/houses in the town so my question is:

Is there any GIS analysis or SQL expression that I can run to find all the duplicated features (both the exact ones and the ones that are located inside others)? I have both ArcGIS and QGIS so I am open to all your suggestions.


texts: Viewsheds analysis with real life photos?

texts:
I'm doing viewsheds presentation for a non-software audience this week. One of the ways I will try to get closer to them is by putting photos from my personal observations from a location + the viewsheds analysis derived from personal computational methods in my presentation.

Is there another, alternative option that I could use?


texts:
I am trying to understand the difference between the two coordinates system. By so far I grasp that EPSG:3857 is used for a web-based map like Google Maps, OpenStreetMap.

I am getting a position (in lat,lon or EPSG:4326). Now my task is to figure out where it is next through its velocity(heading and speed) and time duration of the measurement.

  new_pos = old_pos + velocity * time

Realise that the equation needs to be in meter calculation.

I then proceed to

Converted lat,lon into EPSG:3857
Perform calculation
Converted them back into lat,lon
Plotted the points onto a map based on EPSG:4326

and the points were far off the next measured point

Then I repeated the calculation except using UTM coordinates. The results were much better (as shown above by blue point)

I reckon the distance measurement should be roughly similar(1unit to 1 m...am I right?

I am working near the equator).

So I am wondering what is wrong down here and which measurement should I take?


texts: Define Range when calling SAGA Kriging from QGIS

texts: Is there a possibility to define the range-value or display the Variogram when you call the SAGA-Kriging-Module from QGIS? Maybe by using the python console?


texts: UpdateCursor cannot see new field added to feature class

texts:
I needed to add a new field to a feature class to act as a unique key, but when I run:

fields = ['OID@', 'state', 'navkey']
ids = inspected_points #list of oid values
oidfield = arcpy.Describe(points).OIDFieldName #reference to ObjectID field name

where_clause_points = """{0} IN{1}""".format(arcpy.AddFieldDelimiters(points, oidfield),tuple(ids))

with arcpy.da.UpdateCursor(points, fields, where_clause_points) as cursor:

I get the error:

for row in cursor:
RuntimeError: Attribute column not found [42S22:[Microsoft][ODBC Driver 13 for SQL Server][SQL Server]Invalid column name 'navkey'.] [db.DBO.lines_split][STATE_ID = 0]

I have been able to add and write to other fields, so I'm not sure why this one is not working.  I know the field is there and available:

fields = arcpy.ListFields('lines_split')
    for field in fields:
    print(field.name)

This produces:

ProjectId
SurveyId
TechnicianUsername
TimestampUTC
GPSLatitude
GPSLongitude
GPSAccuracyInFeet
GPSAltitudeInFeet
Shape
OBJECTID
state
NEAR_FID
NEAR_DIST
NEAR_X
NEAR_Y
navkey

This field is just a text field.  I just want to populate it with values, but cannot for some reason access it.  The code above all works perfect if you remove the 'navkey' from the list of fields, but once you add it in it fails.

Ultimately I am looking to create and populate a field to contain GUIDS and act as a unique key, since you cannot rely on ObjectID or GlobalID as they can change and CopyFeatures does not copy the same keys.  I have the python code to generate GUIDS as a string, but I need access to this new field to populate them.


texts: GIS dataset with latitude, longitude, and elevation

texts:
Websites such as this one: https://gisgeography.com/best-free-gis-data-sources-raster-vector/
have many references to elevation data. They are all in *.jpg or other image format, however. What I need is simply a flat text file with lat/long/elevation as real numbers.
I suppose there is a way to find or extract such data. Any pointers?


texts:
I'm a beginner in using arcpy, and having trouble with the following code. I have 2 different input feature layers in different folder. I would like to create a new gdb in the output_folder (this works). I get the spatial reference from the Sentinel image, which I use to project and copy new files into the new gdb, but I get an error message in line: arcpy.env.workspace= ujgdb

How can I fix it?
infc1 = arcpy.GetParameterAsText(0)
infc2 = arcpy.GetParameterAsText(1)
Sentinel = arcpy.GetParameterAsText(2)
output_folder = arcpy.GetParameterAsText(3)
#spatial reference
arcpy.env.workspace= Sentinel
inraster = arcpy.ListRasters("*", "TIF")[0]
spatialRef = arcpy.Describe(inraster).SpatialReference
ujgdb= arcpy.CreateFileGDB_management(output_folder, "DEM.gdb")
arcpy.env.workspace= ujgdb
spatialkont= arcpy.Describe(infc1).SpatialReference
if spatialkont.Name == spatialRef.Name:
    arcpy.CopyFeatures_management(infc1, "kontur_utm")
elif spatialkont.Name =="GCS_WGS_1984":
    arcpy.Project_management(infc1, "kontur_utm", spatialRef)
spatialmag= arcpy.Describe(infc2).SpatialReference
if spatialmag.Name == spatialRef.Name:
    arcpy.CopyFeatures_management(infc2, "magassag_utm")
elif spatialmag.Name =="GCS_WGS_1984":
    arcpy.Project_management(infc2, "magassag_utm", spatialRef)


texts: How to find mean coordinates of a cluster of points from a point layer in QGIS?

texts: I have a set of coordinate points. The points are distributed in such a way that some points are clustered around a specific area. There is an option to find out mean coordinates(s) of the whole point layer. I would like to know if there is any way in QGIS to find out the mean coordinates of the points from a specific area from a point layer.


texts: ArcGIS Collector for Windows 10 Desktop Issue

texts:
I've got a few collector apps and they all seem to work fine on Android and iOS devices, but I can never get them to load into the desktop Collector for ArcGIS. Every time I get this error on the layer:

Invalid SpatialReference: The Wkid or VerticalWkid value is invalid

I can still for some reason add new features. As soon as I finish adding it, it disappears from the map though.

Notes:

These Collector feature layers are stored on my ArcServer, from an SDE that is versioned and tracking edits.
I've got points lines and polygons across various SDEs and apps, tried a lot of other workarounds, but nothing is working.
Layers are in WGS84 and CVGD2013


texts: Deciding on classification method and number of classes for mean temperature

texts:
I'm trying to decide on the classification method: natural breaks, equal interval, quantiles, or manual and the number of classes for a map showing the mean temperature of places in Canada.

The values range from -20 degrees to 10 degrees. The mean is -2 and the median is -1.7.

I'm not sure what to consider when deciding these things. The higher the number of classes, the better it seems to look, but it's hard to decide on the number of classes that is most useful to viewers of the map.


texts: gdal_rasterize: "Attempt to create 0x0 dataset is illegal"

texts:
I am new at GIS.

I want to burn all polygons from a particular shapefile (i.e. Gbg13004.shp) into a RGB TIFF file (i.e. temp.tiff). I have tried the following command:

gdal_rasterize -burn 255 -burn 0 -burn 0 Gbg13004.shp temp.tif

This results in the following error:

ERROR 1: Attempt to create 0x0 dataset is illegal,sizes must be larger than zero.
ERROR 1: Cannot create temp.tif

Running ogrinfo on the input file results in the following output:

INFO: Open of `Gbg13004.shp'
  using driver `ESRI Shapefile' successful.

Layer name: Gbg13004
Metadata:
  DBF_DATE_LAST_UPDATE=2018-01-25
Geometry: Polygon
Feature Count: 15401
Extent: (177898.268767, 218649.356570) - (186487.066994, 226178.919467)
Layer SRS WKT:
PROJCS["Belge_Lambert_1972",
    GEOGCS["GCS_Belge_1972",
        DATUM["Reseau_National_Belge_1972",
            SPHEROID["International_1924",6378388.0,297.0]],
        PRIMEM["Greenwich",0.0],
        UNIT["Degree",0.0174532925199433]],
    PROJECTION["Lambert_Conformal_Conic_2SP"],
    PARAMETER["False_Easting",150000.01256],
    PARAMETER["False_Northing",5400088.4378],
    PARAMETER["Central_Meridian",4.367486666666666],
    PARAMETER["Standard_Parallel_1",49.8333339],
    PARAMETER["Standard_Parallel_2",51.16666723333333],
    PARAMETER["Latitude_Of_Origin",90.0],
    UNIT["Meter",1.0]]
OIDN: Integer64 (15.0)
UIDN: Integer64 (15.0)
VERSIE: Integer (4.0)
BEGINDATUM: Date (10.0)
VERSDATUM: Date (10.0)
TYPE: Integer (4.0)
LBLTYPE: String (64.0)
OPNDATUM: Date (10.0)
BGNINV: Integer (4.0)
LBLBGNINV: String (32.0)
LENGTE: Real (16.2)
OPPERVL: Real (16.2)

My problem seems similar to the one described in gdal_rasterize error: "attempt to create 0x0 dataset is illegal..". Here, the problem was solved by reprojecting the source file into a different coordinate system and rasterizing the projected shapefile. In my case, the source file is in EPSG:103300. After projecting it (using ogr2ogr) into WGS84 and rasterizing the result, I still got stuck on the same error.

How do I solve this?


texts:
I would like to calculate the shortest path with dijkstra algorithm.

sample query
SELECT seq, id1 AS node, id2 AS edge, cost FROM pgr_dijkstra('
                SELECT gid AS id,
                         source::integer,
                         target::integer,
                         length::double precision AS cost
                        FROM ways',
                30, 60, false, false);

I imported an osm map, it has proper network. The name of the table is at_2po_4pgr, colums are
id, osm_id, osm_name, osm_meta, osm_source_id, osm_target_id, clazz, flags, source, target, km, kmh, cost, reverse_cost, x1, y1, x2, y2, geom_way.

in result I would like to see the shortest path and streets name in a row where route was planning.

Thanks


texts:
In ArcGIS for Desktop 10.1, is there a way to change the temp folder where data gets processed?

My temp folder is C:\Users\langerb\AppData\Local\Temp\arc5D6A. I want to change this to a different drive with more space.

I know there's an environmental setting for scratch workspace, and a setting for the display cache, but I can't find a setting to change the temp folder.


texts: Is it possible to import a PDF with vector layers into MapInfo as vector layers? If yes, how would I go about doing this?


texts:
I tried to find tools to delimitates a specific region area of a map
When the mouse is on a specific area, there is the region delimition which appears.

Like this website :
http://www.police.act.gov.au/crime-and-safety/crime-statistics.aspx
or even livelovely maps

I try use google api Fusion Tables, but it seems it can't do what i want...
http://code.google.com/apis/maps/documentation/javascript/examples/layer-fusiontables-styling.html

Maybe do you have some advices ?
Thanks


texts: Imagine 2014 displaying images with large no-data values

texts:
I'm starting to use Imagine 2014 after a recent upgrade.  The images I'm using have data values in the range of approximately -.02 to 350.0.  The no-data values bordering the data has a value of 2 billion. The data displays white, nan black.  It was easy to limit the values considered for generating a lookup table for display in previous versions but I'm not finding the tools now.  What is the best way to  get this data to display legibly in Imagine 2014?

Thanks.


texts: Determining surface that reflectance curve matches?

texts:
I was wondering to what surface that reflectance curve belongs?

It is definetely no water or snow since their reflectance drops at a much smaller wavelength.

I was thinking about a dense cloud but the decrease seems to be too constant.

Has anybody an idea what it could be?


texts: Getting 'list' object has no attribute using ArcPy for ArcGIS Pro?

texts:
I am using ArcGIS Pro with Python 3.6.

I'm posting this as a new question because I've made significant headway in narrowing down the problem. I believe the issue is how ArcPy is used with ArcMap versus ArcGIS Pro.

Here is the error:

 Query = FULL_STREET_NAME LIKE '%'
Traceback (most recent call last):
  File "\\coacd.org\dfs\SWS\Code Enforcement\ACD Admin Operations\GIS\Operations\BSC Maps\Scripts\Python3compatible\BSCMapGenerator_v7.py", line 41, in <module>
    for lyr in mp.ListLayers(): #set visibility and remove older layers
AttributeError: 'list' object has no attribute 'ListLayers'
[Finished in 11.2s]

Here is a sample of the script:

workspace = "G:\\Code Enforcement\\ACD Admin Operations\\GIS\\Operations\\BSC_Maps\\AutoLayerFiles.gdb"
arcpy.env.workspace = workspace
arcpy.env.overwritOutput = True #overwrite existing files
aprx = arcpy.mp.ArcGISProject("G:\\Code Enforcement\\ACD Admin Operations\\GIS\\Operations\\BSC Maps\\BSC_Template_Auto_Generate_PRO_v1.aprx")
mp = aprx.listMaps(aprx)
#df = arcpy.mp.ListDataFrames(aprx)[0] # address location data frame and turn on select layers
location = arcpy.GetParameterAsText(0)
location = location.upper() #convert to uppercase
caseNumber = arcpy.GetParameterAsText(2)
userName = getpass.getuser().upper()
date = time.strftime("%m/%d/%Y")  #added to map
expression = "FULL_STREET_NAME LIKE '" + location + "%'" #this is an SQL expression and can be modified as necessary
expression2 = "LAND_USE_ID LIKE '%'"
arcpy.AddMessage("Query = " + expression)
visibleStructuresLayers = ["TRANSPORTATION.street_segment","Structures","Subject Property","Parcels","Lakes"]
visibleLocationLayers = ["TRANSPORTATION.street_segment","Subject Tract","Parcels","Lakes"]
visibleRoadMapLayers = ["Subject Tract","Subject Property","Lakes","Basemap","OpenStreetMap"]

locationMapScale = arcpy.GetParameterAsText(3) #default 2000
roadMapScale = arcpy.GetParameterAsText(4) #default 5500
structureMultiplier = (arcpy.GetParameterAsText(5) # default 1.1 (110% of parcel size) #float(arcpy.GetParameterAsText(5))

#-----------REFRESH TEMPLATE------------------
#arcpy.AddMessage("Refreshing template...")
#arcpy.Delete_management("Subject Tract")  #removes instances from memory (not from map)
#arcpy.Delete_management("Subject Property")

for lyr in mp.ListLayers(aprx): #set visibility and remove older layers
    lyr.visible = False
    if lyr.name == "Subject Tract":
        arcpy.mp.RemoveLayer(lyr)
    if lyr.name == "Subject Property":
        arcpy.mp.RemoveLayer(lyr)
    if lyr.name in visibleLocationLayers:
        lyr.visible = True

For reference, it used to read like this in 2.7 were df was the arcpy.mapping.listdataframes(mxd):

arcpy.AddMessage("Refreshing template...")
arcpy.Delete_management("Subject Tract")  #removes instances from memory (not from map)
arcpy.Delete_management("Subject Property")

for lyr in arcpy.mapping.ListLayers(mxd, "", df): #set visibility and remove older layers
    lyr.visible = False
    if lyr.name == "Subject Tract":
            arcpy.mapping.RemoveLayer(df, lyr)
    if lyr.name == "Subject Property":
        arcpy.mapping.RemoveLayer(df, lyr)
    if lyr.name in visibleLocationLayers:
        lyr.visible = True


texts: Center a group of features all to the same point

texts: Given a group of features in a geojson object, is there a way to translate them, such that they all have the same center? I'd like achieve this in the actual coordinates of the files.