texts: Arcpy script for summing pixel values in loop can't handle large number of rasters

texts:
I want to sum the pixel values of large numbers of rasters of different extents (ArcMap v10.2).  16GB RAM.

I have the following code (with help from this site) which solves the different extent problem and works fine for up to about 1000 rasters (40-50 mb) and outputs a single raster of about 50mb correctly.  But when the number gets much larger ArcGIS closes.  It iterates through all the rasters but crashes at the save - I suspect it is accumulating data in memory throughout the loop and eventually runs out of memory.

At the moment the rasters are all in directory folders (5000 per folder - I need to find a better way to do this as well).

None of the file paths have spaces in them.

I have tried saving as tif rather than GRID but no difference.
I have tried to save within the loop but couldn't get that to work.
I've tried using a .gdb as the workspace and scratchWorkspace.
Ive tried some of the things suggested here https://geonet.esri.com/thread/101053 .
Next im going to try storing the rasters in a mosaic dataset and iterate through that but I expect that will take me awhile to figure out.

import arcpy
import os
from arcpy.sa import *

ws=arcpy.GetParametersAsText(0)
outputFolder=arcpy.GetParametersAsText(1)
outputName=arcpy.GetParametersAsText(2)
extentRaster=arcpy.GetParametersAsText(3)

arcpy.env.overwriteOutput=1
arcpy.env.scratchWorkspace = ws
arcpy.env.workspace = ws
arcpy.CheckOutExtension("Spatial")

eRaster=Raster(extentRaster)
arcpy.env.extent=eRaster.extent
arcpy.env.snapRaster=extentRaster

#create a list of rasters in the workspace
rasterList = arcpy.ListRasters("*", "GRID")

i = 0
#loop through rasters in list
for raster in rasterList:
    arcpy.AddMessage ("processing raster: %s" %i)

    #sum rasters together and convert no data to zero
    if i == 0:
        outSUM = Con(IsNull(raster),0,raster)
        i += 1
    else:
        outSUM = outSUM + Con(IsNull(raster),0,raster)
        i += 1

#save final output to the current workspace
arcpy.env.workspace = outputFolder
outSUM.save(outputName)


texts: Offsetting polygons by meters with clipper

texts:
I would like to create buffered polygons of locations (towns, villages etc) in order to use them for searching in radius.

This is what I would like to achieve (units for ilustration):

This is how I do it in pyclipper:

import pyclipper

coordinates = # Array of lat,lng tuples
clipper_offset = pyclipper.PyclipperOffset()
coordinates = pyclipper.scale_to_clipper(coordinates)
clipper_offset.AddPath(coordinates, pyclipper.JT_ROUND,
pyclipper.ET_CLOSEDPOLYGON)
scaled_coordinates = clipper_offset.Execute(1000.0)
scaled_coordinates = pyclipper.scale_from_clipper(scaled_coordinates)

Number 1000.0 is arbitrary and my question is - how do I calculate the right offset ratio for Execute method, so that the offset polygon will approximately be offset by 10,20 and 50km ?

Btw. is this the right approach to this problem ?


texts: Joining attribute of smaller polygon to attribute of biggest polygon when they are overlapping

texts:
I have two shapefile layers which include polygons. Each polygon on the first shapefile overlaps a polygon on the other shapefile.
I am trying to join the attribute of smaller polygons to the overlapping bigger polygons. I have used spatial join, but the attribute table of the bigger polygon is not showing the correct attributes of the smaller polygons inside it. I figured out there are some tiny polygons inside each polygon.
How can I join the attribute of the biggest overlapping polygon to the other shapefile?


texts: Missing processing toolbox in qgis 3.10.3 for mac

texts: I am trying to understand the basics of QGIS following their online training manual but I can't find the processing toolbox anywhere in qgis 3.10.3 for mac.


texts: Online/offline tool for collection of coordinates for simulation purposes

texts: Is there any good free online tool to collect coordinate positions given, say a trajectory; or a polygon. For example, retrieve all coordinates between two points, every 500 meters.


texts: Using for loop to execute query in ArcPy?

texts:
I'm trying to provide the breakdown of fire detection for an area by sorting of confidence levels. The dataset has 12 files (one for each month in a year) with fire records which all have a column called "confidence" with a value between 0 and 100. I have to use a for loop to query the files by 3 different confidence levels.
Group 1: >85% confidence
Group 2: 60-85% confidence
Group 3: <60% confidence

I HAVE TO USE a for loop to execute this. Any ideas? I started by appending all the files into one list. I'm told that I could use lists to set up lower and upper ranges and SelectLayerByAttribute_management, and use that in a loop, but don't know how to do that.

This is my final code that works. I decided to merge the files together into one big file.

arcpy.env.workspace = "***********/data/"

list1 = [];

for month in range(1,13):
    inf = "2007" + str(month).zfill(2)+ "_rfe"
    arcpy.DefineProjection_management(newpath + inf + ".shp", 4326)
    list1.append(newpath + inf + ".shp")

out1 = "big2007" + "_rfe" + ".shp"
layer1 = "layer1" + ".shp"
arcpy.Merge_management(list1,out1)
arcpy.MakeFeatureLayer_management(out1, layer1)

lower = [0, 60, 86]
upper = [60, 86, 101]

for i in range(len(lower)):
    arcpy.SelectLayerByAttribute_management(layer1, "NEW_SELECTION", """ "CONF" >= {} and "CONF" < {}""".format(lower[i],upper[i]))
    result = arcpy.GetCount_management(layer1).getOutput(0)
    result = int(round(result))


texts: Repositioning polygon Shapefile using ground control points and v.rectify?

texts:
I'm new to Grass.

I would like to reposition a polygon shapefile using ground control points and v.rectify seemed just the perfect tool to do that.

The raster from which the shapefile comes from was georeferenced using wrong control points for the KKJ1 projected coordinate system (EPSG:2391), now I obtained the correct ones and I am trying to find a way to reposition the vector without having to digitize everything again.

After running the command, the vector does move to the desired position, yet the original conformation is not retained and in particular some of the areas "lose" the centroid and become just empty boundaries.

Moreover I noticed that the areas whose centroids had not gone missing in the process, presented wrong attributes when inspected.

I had a look at the output and 3 errors show up (my output is in Italian so the translation might not be exactly the same):

WARNING: Number of incorrect boundaries: 1752
WARNING: Number of centroids out of the area: 139
WARNING: Number of duplicated centroids: 25


texts: Geomedia WebMap WFS to vector marker layer in Openlayers

texts:
Looking for minimal example/advice on getting a WFS service to show up as a vector layer with markers/icons in Openlayers map.

I have created the WFS service in Geomedia WebMap 6.01.11.19 and the Openlayers map.

a snippet from my map html code:

var wfs = new OpenLayers.Layer.Vector("myFeats", {
            minScale: 100000,
            maxScale: 99,
             strategies: [new OpenLayers.Strategy.BBOX()],
             protocol: new OpenLayers.Protocol.WFS({
                    version: "1.1.0",
                    srsName: "EPSG:900913",
                    url: "http://gisdev2/OL_WFS/Request.aspx",
                    featureType: ["myFeats"],
                    featurePrefix: "gmgml",
                    featureNS: "http://www.intergraph.com/geomedia/gml",
                    geometryName: "GEOMETRY_PT",
              }),
             renderers: OpenLayers.Layer.Vector.prototype.renderers
        }
    );

Opening the IE11 debug console it bombs near this line in the Openlayers.js code

{return this.xmldom?"string"==typeof a?this.xmldom.createNode(1,b,a):this.xmldom.createNode(1,b,""):document.createElementNS(a,
b)}

I suspect the problem is with the namespace in the WFS. I have read elsewhere to toggle the checkbox "on" for the 'Feature Collection in WFS namespace' setting in Geomedia WebMap. I have done this. I have also kept the 'Features namespace' value equal to 'http://www.intergraph.com/geomedia/gml' and the 'Features namespace prefix' equal  to 'gmgml'.

Still no luck getting the feature to show in the Openlayers map.


texts: How to find all ICommandItem instances of a custom button command in an ArcMap add-in?

texts:
I have an ArcMap 10.0 add-in in which I have a button that the user can place (using the Customize dialog in ArcMap) on any toolbar including user-created toolbars as well as standard ones or even 3rd party ones. They could even put multiple instances of this button on the same or many toolbars (what can I say, some people can be crazy).

I have some code I need to run on each ICommandItem instance (to set its FaceID property, see this question if interested) wrapping this button command.

How can I reliably find each of these instances?

Currently I only set the first (and only) one returned by ICommandBars.Find, leaving the rest alone. Not the worst thing in the world, but less than ideal.

I thought I had a lead with this snippet (List ArcMap Visible Toolbars) but this only lists standard ArcMap toolbars (those in the 'ESRI Mx Command Bars' component category in the registry). From that I could iterate over each command item in each toolbar.

This approach could work if I could get it to list all toolbars including user-defined ones (which are not in any category as far as I can tell), and, admittedly less-importantly, 3rd party ones (which may or may not always be in the same category used by the snippet, I'm not sure).

I found this crusty old unanswered post from 2003 on the ESRI forums asking how to find user toolbars: "Detect user defined toolbars programmatically" -- don't you just hate dead ends?

Lastly I found an answer in this post, "looping through all commandbars", (from our own GIS.se member Kirk Kuykendall, no less) which looks promising but I can't see how to translate the CreateObject call to .NET, and crucially, it doesn't look like you can actually get an ICommandBar reference from its items.

Any pointers would be greatly appreciated.


texts: How to process geo-referenced PDF files?

texts:
I need to process a batch of 3rd party PDF files to

extract the bounding lat/long box
convert them to a .mbtiles format

It looks like the gdal is the go to tool for this but non of the gdal command recognise my pdf as geo referenced. For instance:

gdalinfo 00989I11.PDF
ERROR 4: `00989I11.PDF' not recognized as a supported file format.
gdalinfo failed - unable to open '00989I11.PDF'.

I installed gdal 2.4.1 on macOS using brew install gdal
When I open the same PDF in Acrobat Reader, I can see the lat/long when moving my cursor:

You can download the sample PDF file here


texts: Change point symbol in mapbox.js/leaflet/esri-leaflet.js

texts:
I am completely new to these libraries.
Until now, everything shows up very fine.
I got my mapbox as a basemap and added some point data provided by an ArcGIS Server.
To get the AGS data onto the map I had to use the ESRI Leaflet API.
My first question is:
1.) How can I change the symbol to a different point marker?
and for later:
2.) How could I add some more fancy pointer style (like an CartoCSS)?
There are some examples on the ESRI site but I don't get everything.
Here's the code:
    <!DOCTYPE html>

<meta charset=utf-8 />
<title>A simple map</title>
<meta name='viewport' content='initial-scale=1,maximum-scale=1,user-scalable=no' />

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

It also seems like the script is loading the "red map" also in the background even though "pencil map" is activated.
Is that true?


texts: How can I make PROC GMAP (SAS) run is less time?

texts:
I am trying to create a map of Australia, by Local Government Area (LGA), with shaded regions displaying the population for each LGA. The shape file can be downloaded from the Australian Bureau of Statistics (ABS) Website. However, the program will not run in a reasonable time period. If I use a subset of the map (say, 100 observations), the code does run, although I get something that is obviously incomplete. If I run the the code with all observations for the map, SAS crashes. There are around 3 million observations in the map, which is too much for SAS to handle.

How can I make Proc GMAP run is less time?

Here is my code:

proc mapimport datafile="[insert path]\LGA11aAust.shp" out=mymap;
run;

proc sql;
create table population_by_lga as
    select distinct lga_code_2011, sum(total) as total
    from caseload2
    group by lga_code_2011
    order by lga_code_2011;
quit;

proc sql;
create table mymap as
    select *, lga_code11 as lga_code_2011
    from mymap
    order by lga_code_2011;
quit;

proc gmap data = population_by_lga map = mymap;
  id lga_code_2011;
  choro total;
run;
quit;


texts: Filling only the deepest sinks/depressions in a DEM

texts:
I am trying to come up with a workflow that will extract the location of surface mines from a DEM. My original thought was to fill the depressions in the DEM and then subtract the original DEM from the filled DEM to produce a grid of the sinks in the landscape. This worked well, and the surface mines stick out like a sore thumb because they are very deep sinks, but it also includes shallower natural depressions that were filled in the landscape.

Is there a way I can only fill sinks that exceed a certain depth such that when I subtract the original DEM from the filled DEM the resulting grid only includes the deepest of sinks (which will, hopefully, be all of the mines in the landscape).

I'm using ArcGIS 10 but I recently discovered SAGA GIS so any solutions using SAGA would be welcome as well.


texts: How to speed up custom processing tool on large dataset

texts:
I have developed a processing tool in QGIS to quality check the data in 26 layers. The tool goes to each layer and checks so many conditions. when I run the tool for small data (less than 500 features in each layer) the tool is quite fast, only takes 3-4 minutes to run. But when I test this to large data (around 10000 Features on each layer) the tool takes forever to run (last test run 27 hours)
The tool currently accepts both GeoJSON and TAB files as input layers and the layers need to be opened in QGIS.
Can someone guide me on how can I increase the speed of this tool?
I cannot attach the code as it is 9000 line of code but the code of simple if conditions.
    def initAlgorithm(self, config):

    parameter=QgsProcessingParameterMultipleLayers(
            'qgs_layers',
            self.tr('Input layers')
            )

    parameter.setMinimumNumberInputs(25)
    self.addParameter(parameter)

    parameter.setFlags(parameter.flags() | QgsProcessingParameterDefinition.FlagAdvanced)
    self.addParameter(parameter)

    self.addParameter(
        QgsProcessingParameterFolderDestination(
            'Output_folder',
            self.tr('Output Folder'),

        ))

def processAlgorithm(self, parameters, context, feedback):

    for layer in parameters['qgs_layers']:
        layer_name=layer.name()

        try:
           exec('self.' + str(layer_name).replace(' ','_')+'_layer_obj=layer')
        except:
            pass

    self.QA_layer_1(feedback)

    self.QA_layer_2(feedback)

def QA_layer_1(self,feedback):
    features=self.structurepoint_layer_obj.getFeatures(QgsFeatureRequest().setFilterExpression (filter))
    for f in features:
        if feedback.isCanceled():
            sys.exit()

        if len(f['id'])!=18 and f['structure_type']!='Pole':
            self.error_list.append([table_name,f['id'],'Invalid ID',f['id'],'The StructurePoint ID should be 18 Character','point',f.geometry()])
            new_error=CriticalError()
            new_error.update({
             'table_name':'structurepoint','Feature ID':f['id'],'Error Type':'Invalid ID','value':f['id'],
             'description':'The StructurePoint ID should be 18 Character','geometry_type':'point','geometry':f.geometry()
            })

            self.error_dict_list.append(new_error)
            feedback.pushInfo('\t{} current value is ({})'.format(str(new_error['description']),str(new_error['value'])))

def QA_layer_2(self,feedback):
    pass


texts:
One of our guys is using the "Functional Surface->Add Surface Information" feature in 3D-analyst tool in ArcMap 10.1 to add elevation data to an XY list.  It's the only feature we're using from the 3D-analyst.  Our trial period has ended and we're trying to figure out if we need to buy the whole 3D-analyst tool for $2500 just for this feature, or if there's a cheaper way to do this task since we don't use anything else in 3D-analyst.  Since it's a pretty simple task, I might just develop a tool myself to do this using the USGS web service given the cost of 3D-analyst.

Is there a cheaper / easier way to do this?


texts: How to access accuracy of binary probability in google earth engine?

texts: I want to use SVM model to do species distribution model in google earth engine. Therefore the output would be a binary probability layer, but the tutorial of earth engine didn't talk about how to access it. Is it possible to do ROC curve in earth engine? I really don't want to export large data and import it to R, then go to the traditional workflow.


texts: How to open and add SQLite connection to browser with PyQgis

texts:
I have a code to open the layers of a SQLite database. When executing this code, the layers are loaded correctly (1 & 2). However i also need to connect to the database and show the connection in the browser toolbar.

How doe load the connection and let it show in the database browser (3)?

uri = QgsDataSourceUri()
uri.setDatabase(r'C:\\path_to_db\\andorra.sqlite')
schema = ''
tables = ['points', 'lines', 'multilinestrings', 'multipolygons', 'other_relations']
display_names = ['points', 'lines', 'multilinestrings', 'multipolygons', 'other relations']
geom_column = 'geometry'

for i in range(len(tables)):
     uri.setDataSource(schema, tables[i], geom_column)
     iface.addVectorLayer(uri.uri(), display_names[i], 'spatialite')


texts: ArcMap ArcPy - Place signs with no lat/lng on map

texts:
I'm new to ArcPy and ArcMap but I do have programming experience. What I am trying to accomplish is to take over 22,000 signs from a database and place them on a map in reference to a street layer that is there. The problem is that there is no lat/lng but the database specifies the road it is on, and the two intersections it is inbetween and the distance it is from the first intersection. I found an example of a trace tool and I am attempting to use that with little luck. I have a mainline layer of the road that I can trace and I have intersection points to reference but again, I am new to ArcMap and ArcPy so any help in the right direction would be great.

I am thinking that I will put in a loop to select the intersection that corresponds with the sign and then trace the distance from that intersection towards the other intersection point. The loop is not a problem but the trace function is rather foreign to me so I was hoping someone could provide me with a place to find some good examples or just a basic overview of how to implement a trace function. The ArcGIS site wasn't much help. Thanks a lot!


texts: Determining if species point data is upstream or downstream of artificial fish passes using ArcGIS Pro

texts:
I am working with ArcGIS Pro.
I am working on a project to determine how to prioritize restoration of artificial fish passes. I have been supplied with species point data, point data of the artificial fish passes and river polyline data. Now I need to determine whether the species points are located upstreams or downstream from the fish passes. Do anyone have any suggestions?
I've tried to create a utility network, but have not been able to get any results.


texts: How to convert 1d numpy array to 2d numpy array

texts:
I used FeatureClassToNumPyArray like this:

arcpy.MakeFeatureLayer_management(datapath, "lyr_fixpoints", "area_in_m2 > 8000")
arr = arcpy.da.FeatureClassToNumPyArray("lyr_fixpoints", ('land_type', 'asset_cat', 'name', 'maintenance'), null_value=-999)
print arr

and get the following 1d array:

[(5, 0, 380, 3) (1, 4, 369, 3) (5, 0, 421, 2) (0, 7, 425, 1)]

Now I want to convert this array into a 2d array of this form, because the 2d array is very suitable for my further operations.

[[5 0 380 3]
 [1 4 369 3]
 [5 0 421 2]
 [0 7 425 1]]

Can someone help me?


texts:
I'm fairly new and entirely self-taught in using QGIS 3.4. I'm looking at turtle nest sites and trying to run a nearest neighbor analysis to analyze their distribution and establish whether or not they are clumped, dispersed or random. I'm doing this for a number of different years to determine whether nest sites are becoming more clumped as a result of environmental changes.

However, i'm not sure i understand the results that are being output. the z scores are just astronomical and the expected values are effectively 0s. How is QGIS establishing expected distances?

Maybe i'm just doing something fundamentally wrong, but i'm selecting the vector analysis -> Nearest Neighbor Analysis -> and selecting my vector layer containing the nest site points. These are some of the outputs I'm getting for different years:

Observed mean distance: 5.451281160546463
Expected mean distance: 0.00026372940677887315
Nearest neighbour index: 20669.97847197658
Number of points: 125
Z-Score: 442084.1070632401

Observed mean distance: 5.40372033661974
Expected mean distance: 0.00019992349282925835
Nearest neighbour index: 27028.941222203965
Number of points: 134
Z-Score: 598544.2290450602

Maybe i'm just doing something fundamentally wrong? Or maybe theres just a better way of achieving what i'm trying to do. I'm not sure.


texts: Placing data frame as inset on only one page of Data Driven Pages map book?

texts:
I am using data driven pages to make a 5-page mapbook of an area. There is also one far-away site that I wanted to include an inset of on the first page of my mapbook (near to where it is). However, I don't want that inset map to be in the following four pages.

How do I make a data frame appear in only one page of my data-driven map book and not in all of the pages?


texts: Calculating angle between polyline segments and assign this value to point?

texts:
I have a polyline feature with many segments (polyline split at vertices).

I plan to record angle values to attribute table per each polyline segment and assign it to point features.

I need to calculate angle value between segments of a polyline (Get value calculated angle less 180 ). .

Here almost is what I need:
Calculating incidence angle of polyline?


texts:
Related to this question asked here on gis.se and this thread on postgis-user, has anyone worked out a good solution for in-database rasterization of a vector layer using PostGIS?  It sounds like the necessary functions are available and it is possible to get a valid output, but not one that is readable by GDAL et al. because it has irregular tiles.

To be more specific, I would like to rasterize two distinct vector maps of the US, counties and NBCD mapping zones, to a particular resolution and both aligned with an existing raster.  It is possible to do this with gdal_rasterize, but with some trial and error as I recall -- mainly because the extents of the three inputs do not match, I expect.

Any thoughts?


texts: Identifying polygons that intersect with another layer using QGIS

texts:
Using QGIS, I have two polygon layers, a vegetation layer and a treatment layer. I want to add a field to the attribute layer of the vegetation layer to indicate if any part of the polygon has been treated (overlaps with a treatment layer).

For the above example - I want an attribute table that ends up looking like

I do not want to intersect the two layers - I need to keep the full polygons in the vegetation layer, I just want to know if they've been treated.  There are too many polygons to do this manually.
I have looked through the vector tools to try identify one that does this, this was not successful.  I've tried to find a workflow via rasters instead of polygons, but cannot figure out what would be needed for this to work.  I have searched online but not been able to find anything that answers this question, likely because I do not know what this is called so do not know what terms to search for.
The only solutions I have found are manual, but there are several hundred polygons so I would like to avoid that approach unless there is no other option.


texts:
As shown in the image is there any possibilities to avoid unwanted dangle errors in arcgis ?!

1:

Update: As you can see in the image above, the four red ones which are close to each other have been determined as dangle errors in ArcGIS, but they are not. They are just dead-end alleys, for example. How to force ArcGIS not to determine that as errors?


texts:
I am trying to use QGIS3.4 by using Emacs (elpy) on Mac OS Mojava. I checked following sources:

https://docs.qgis.org/3.4/pdf/en/QGIS-3.4-PyQGISDeveloperCookbook-en.pdf
http://osgeo-org.1560.x6.nabble.com/Error-importing-qgis-in-python-td4141685.html
Calling QGIS from MacOS python

Then, I use the following code:

import os
import sys
sys.path.append('/Applications/QGIS3.4.app/Contents/Resources/python/')
sys.path.append('/Applications/QGIS3.4.app/Contents/Resources/python/plugins')
sys.path.append('/Applications/QGIS3.4.app/Contents/MacOS/lib')
import qgis.core

But, I got the following result.

ImportError: dlopen(/Applications/QGIS3.4.app/Contents/Resources/python/qgis/_core.so, 2): Library not loaded: @executable_path/../Frameworks/qgis_core.framework/qgis_core
  Referenced from: /Applications/QGIS3.4.app/Contents/Resources/python/qgis/_core.so
  Reason: image not found

Could you help me to resolve this issue?


texts: Enable multiple GDAL threads for qgis2web?

texts:
I'm using qgis2web plugin for a project that has a raster base map and many layers. It works well but exporting is slow, taking 16 minutes. gdalwarp.exe seems to use only one thread (the CPU has eight). I have set GDAL_NUM_THREADS to 7 in the general settings and also project settings but it doesn't seem to do a difference.

Starting from a command prompt and setting GDAL_NUM_THREADS doesn't seem to help either. Using OSGEO4WIN shell doesn't help either.

Is it possible to enable multiple GDAL threads for QGIS plugins?

Windows 7 Enterprise

QGIS 2.18.7 (not newest, but that's what is installed)

qgis2web 2.33.1


texts: Does the number of decimal places in coordinates affect the speed of rendering of the features on a Leaflet map?

texts: I am trying to improve the performance of a webpage that renders hundreds of points on a Leaflet map. I noticed that the coordinates of these points are saved and retrieved with unnecessary accuracy (14 decimal places). I wonder if limiting the decimal places to 5, which is enough for my purpose, should (and if, why?) impact the speed of rendering of the map with features.


texts:
I have created a web based GIS using Geoserver-Postgis-Openlayers...

Every thing was working fine till I changed the server path from 'Testing to LIVE.'

WMS POP up Function has stopped working in LIVE Server Environment & showing the following error message..

  "Message: Access is denied.

  Line: 177 ; Char: 9 ; Code: 0 ; URI:
  ....../geoserver/OpenLayers-2.8/lib/OpenLayers/Request/XMLHttpRequest.js
  "

WMS Pop Up code is as follows (Working perfectly at my Test Server)

info = new OpenLayers.Control.WMSGetFeatureInfo({
            url: '..../geoserver/LT_WEB/wms',
            title: 'Identify features by clicking',
            queryVisible: true,
            eventListeners: {
                getfeatureinfo: function(event) {
                    map.addPopup(new OpenLayers.Popup.FramedCloud(
                        "",
                        map.getLonLatFromPixel(event.xy),
                        null,
                        event.text,
                        null,
                        true
                    ));
                }
            }
        });
        map.addControl(info);
        info.activate();

        map.addControl(new OpenLayers.Control.LayerSwitcher());
        map.zoomToMaxExtent();

SERVER - Microsoft Windows 2003 R2.
Web Browser - Internet Explorer 8.0.

In both Servers (Testing & LIVE) GeoServer, PostGIS and OpenLayers 2.8 installed and all paths verified.

Awaiting for your valued response.


texts: I've got a txt/csv file with the utm coordinates and the hights. Now I want to upload this into qgis. This is not a problem so far. The problem is to create a shapefile out of the data with an attribute table for adding other values to work with. Has anyone an idea how to do this?


texts: Loading and rendering a .ecw raster with PyQGIS?

texts:
I am trying to load an .ecw file with pyqgis on a pc. Usin version 2.14 of QGIS. I tried the code i could find in the cookbook about loading and rendering raster from it I interpred this code:
from qgis.core import QgsRasterLayer
import os
os.chdir(r'C:\some\path')
file = 'FooBar.ecw'
fileInfo = QFileInfo(file)
baseName = fileInfo.baseName()
rlayer = QgsRasterLayer(file, baseName)
if not rlayer.isValid():
    print "Layer failed to load again!"
else:
    print rlayer.renderer()
    print rlayer.renderer().type()
    if hasattr(rlayer, "setCacheImage"):
        rlayer.setCacheImage(None)
    rlayer.triggerRepaint()

I will get a dialog window where I get to define the projection for the raster. But th file will not draw.
The file is fully functionall to load normally. How can I make the raster draw with pyqgis?


texts: How can I superimpose a geopandas dataframe on a raster plot?

texts:
I have a geopandas dataframe countries with country polygons, and a raster dataset raster (read with rasterio). I've already aligned the two as follows:

countries = countries.to_crs(raster.crs.data)

I can see they're aligned when I plot each separately. However, I have been unable to display the country borders overlaid on top of the raster. The following just shows the raster image, no country borders:

import matplotlib.pyplot as plt
from rasterio import plot as rioplot
plt.figure(figsize=(15,15))
base = rioplot.show(raster)
countries.plot(ax=base, color='red');


texts: Problem installing OTB for QGIS 2.8.15 in ubuntu

texts:
I'm experiencing a problem in installing OTB for QGIS 2.8.15 in my ubuntu 16.04.

I'm following the instructions in
https://wiki.orfeo-toolbox.org/index.php/QGIS_access_to_OTB_applications

but when I do "make" I obtain this error

/usr/bin/ld: /usr/local/lib/libpython2.7.a(abstract.o): relocation R_X86_64_32S against `_Py_NotImplementedStruct' can not be used when making a shared object; recompile with -fPIC


texts:
I am trying to select the latest date by expression in QGIS.

I tried to use this code "Date"= maximum("Date"). But nothing was selected.

I wonder what is the right expression code to select the latest date in general?

The date is formatted as follows for example: 28-4-2018

And could I select the latest date for a specific "ID"?


texts:
I'm trying to work with a layer that uses a PostGIS geometry column of 3D polygons, but hit detection seems to be broken when trying to use the select tools. The geometry renders fine, but when I click on a feature using the highlight or identification tool the status bar in the bottom left reads "No features at this position found."

I can use ST_Force2D on the geometry, but I'd like to avoid having an extra calculated column on the original table. A view would work, but I need to be able to edit attributes of selected features, and I'm not sure how to do that or even if it can be done if my geometry column is from a view.

What's a good solution that will allow me to select features?


texts: How to intersect more than 2 polygons using Turf.js?

texts:
We use Mapbox GL JS 1.12.0 in our Vue.js project. I need to create intersection objects if I have 2 or more polygons. But, I'm able to do it with 2 polygons only. The code looks as following:
import * as turf from '@turf/turf';

export function createIntersection(features) {
  // features = 3 polygons on the screen bellow
  const intersection = turf.intersect(...features);

  return intersection;
}

The screenshots:

Unselected polygons:

Selected polygons

After code execution

As you can see there is only 1 intersection object created.
How to do it with more than 2 polygons?


texts: Conversion of CSV file in MyGeoData Converter - XY columns not recognized

texts:
Recently I have tried to import various CSV files with XY coordinates to MyGeoData Converter. In some cases I got a message:

  Selected data does not have a spatial part!

For example, this file works fine:

id,name,coordx,coordy
1,Start,18.63,49.5
2,End,18.65,49.42

Within this file I am able to see correct extent in the map, assign coordinate system and transform to desired coordinate system and format, but I got a wrong result with this file (no spatial part detected):

id,name,koordinatenx,koordinateny
1,Start,18.25,41.16
2,End,18.26,41.12

Seems like columns with XY coordinates were not recognized.

Any idea?


texts: Clipping a shp returns another with an area bigger than the mold used on the clip

texts:
The problem at hand includes 3 shapefiles. 2 Land Use classified satellite images and 1 regular polygon shapefile that covers the area of a National Park in my country. The software used in this analysis is QGIS 2.14.7. All shapefiles are using the SIRGAS 2000 CRS (EPSG: 4674)

My assignment was rather simple, I had to calculate the area for each class on the Land Use inside the park. The mold's area (the shapefile with the boundaries of the national park) is of 870.000 hectares.

I made sure both satellite images didn't overlap and merged them together. After that, I clipped the merged shapefile using the boundaries shapefiles as the mold. The output shapefile was created with the exact shape of the mold. Apparently everything was ok.

At last, I use the statistics for categories tool (located on the toolbox) to get a small table with the sums of each class. When I sum all the sums, the area goes to approximately 938.000 ha.

Any clues on why this might have happened?

I tried the topology checker plugin. I found 20 errors close to the limits of the satellite images, but they are small and wouldn't make up for almost 60.000 ha. Is it possible that just by the presence of errors, QGIS goes cuckoo?

I know I could dissolve the satellite shapes, but I'm not on a very powerful computer and the processing takes too long. Would it make a big difference anyway?

They're all on the same CRS, including the project.


texts: Using Python dictionary as parameter and variable in Python Toolbox tool of ArcPy

texts:
I'm building a feature class merging and comparing toolbox in arcpy, ArcGIS 10.7.
I have a working command line utility with fixed featureclass values, and now I'm building this toolbox, I just copy the relevant sections:

        param2 = arcpy.Parameter(
        displayName="My filter list:",
        name="fc_dict",
        datatype="GPString",
        parameterType="Required",
        direction="Input",
        multiValue = True)
    params = [param0, param1, param2]
    return params
def updateParameters(self, parameters):
    fdict = {'key1':'value1','key2':'value2','key3':'value3'}
    exportList = []
    for key in fdict:
        exportList.append(key)
    parameters[2].filter.list = exportList
    return
def execute(self, parameters, messages):
    """The source code of the tool."""

    import os
    exportList = parameters[2].valueAsText
    for search in exportList.split(";"):

        arcpy.AddMessage("Exporting " + search)

With this I can populate my GPString with the keys and use them, what's nice. But I have to use my values also as variable. But I don't understand the accessing approach.


texts: Grid convergence in ArcGIS Desktop gives conflicting values?

texts:
I calculated the grid convergence of a feature in ArcGIS (using geographic rotation and none coordinate system). The returned value is negative, which is conflicting to my understanding.

The easting of the calculated location is 570000m, east of central meridian -117, northern hemisphere. My understanding about grid convergence is any position east of central meridian should be positive.

How ArcGIS gave a negative grid convergence value?


texts: Updating PostGIS table AND setting area_id

texts:
I have a table with lots of points which I'd like to update and set a consecutive area_id to every point "group". I define a point "group" as points that are within a 500 meters radius.

The raw table is something like:

point_id    desc_string       lat            lon          gid        geom
A           Glion           25.057939     55.187833       1
B           Glion           25.061244     55.185937       2
C           Springs         46.433853     6.926215        3
D           Springs         46.431851     6.923763        4
E           Hobby Center    24.221813     55.780285       5
F           Ponte           -12.643624    14.864234       6
G           Ponte           -12.644966    14.866923       7
H           Ponte           -9.093045     17.658906       8

You can now see that there are different pairs or groups of points.
With the following SQL I can identify the points which by definition belong to these groups:

SELECT distinct on (point_id) t1.*
FROM ciss_poi AS t1 JOIN ciss_poi AS t2 ON St_dwithin(t1.geom, t2.geom, ( 0.5 / 111.111 ))
WHERE t1.desc_string = t2.desc_string AND t1.gid != t2.gid

What I want is to group these pairs and add an area_ or group_id:

point_id    desc_string       lat            lon          gid        geom     group_id
A           Glion           25.057939     55.187833       1                   1
B           Glion           25.061244     55.185937       2                   1
C           Springs         46.433853     6.926215        3                   2
D           Springs         46.431851     6.923763        4                   2
F           Ponte           -12.643624    14.864234       6                   3
G           Ponte           -12.644966    14.866923       7                   3

I have no idea how to write the query. I can identify the groups but I don't know how to add the area_id.

Can someone help me?


texts: Splitting polygons and dividing attribute value proportionally using ArcGIS Pro?

texts:
I am using ArcGIS Pro 2.3.0.

I have two feature classes (pictured below), RED polygons and BLUE polygons. The RED polygons have a numeric attribute, Value1.

How can I split the RED polygons with BLUE polygon features and have Value1 be divided proportionally?

i.e. RED feature A has Value1 of 100. If a BLUE polygon split a RED polygon in half, each new RED polygon B & C would have a Value1 (or new field) of 50.


texts: How can I have my map open in the center of my visitors IP address?

texts: I have a map of all restaurants in Florida. I would like to have the center of the map be the address of the visitors IP. How can I do this?


texts: URL " http:// /ArcGIS/rest/service/GPTask/MapServer/jobs/" not found in the ArcGIS Server?

texts:
I am Using ArcGIS Server 10.03 SP3. I created a geoprocessing service in the server which returns a result map service. I successfully created the gp service and map service with same name. I need to display the result in my application using ArcGIS Javascript API v3.3.
The logic works like the sample application  http://developers.arcgis.com/en/javascript/samples/gp_resultmapservice/

I tested the gp service in the REST Web interface and successfully executed.The result map is obtaining after completion of the task.

But when I use it in the JS API front end , the fire bug tool shows success execution of the gp task , but when it executes

map.addLayer(gpLayer);

the response shows

dojo.io.script.jsonp_dojoIoScript9._jsonpCallback({"currentVersion":10.03,"error":{"code":400,"message":"Unable to complete  operation.","details":["Invalid URL"]}});

When I tried the url  http:// /ArcGIS/rest/service/GPTask/MapServer/jobs/

The page not found shows.

When I changed the code with

gp.getResultImageLayer(jobinfo.jobId, "output", imageParameters, function (gpLayer) {
 map.addLayer(gpLayer)
});

I got  the result .

Can anybody help to find the exact reason? I need to include the service via ArcGISDynamicMapServiceLayer

Here is the My sample code that Did not work

function applyGPTask() {
var gpServiceUrl = "http://localhost/ArcGIS/rest/services/GPTasks/GPServer/GPTask";
gp = new esri.tasks.Geoprocessor(gpServiceUrl);
var params = {
    param1: "SUBDIVID=202",
    param2:"value2",
    param3: "value3"

 };
  gp.setOutputSpatialReference({ wkid: 4326 });
  gp.submitJob(params, gpJobComplete, gpJobStatus, gpJobFailed);
 }

function gpJobComplete(jobinfo) {
var mapurl = "http://localhost/ArcGIS/rest/services/GPTasks/MapServer/jobs/" + jobinfo.jobId;
var gpLayer = new esri.layers.ArcGISDynamicMapServiceLayer(mapurl, {
    "id": "gpoutput",
    "opacity": 0.7
});
    map.addLayer(gpLayer);
 }


texts: Heatmap is rendering different radii (radius) for the same data in different layers

texts:
Background: I am plotting car crashes along a road and creating a heatmap of crash frequency. I'm using multiple layers that were pulled from the same source data. The layers just sort things by whether they were northbound or southbound on this particular stretch of road.

Issue: When I create a heatmap with the heatmap plugin, different layers return different sized circles dispite the radius setting in the heatmap setup being identical in each of these layers. Oddly, the cell size changes for these different layers as well even though they come from the same source.


texts:
I am new to using QGIS and I have been trying to build a Map which will eventually be put on my firm's website with overlaying layers for different indicators in the local area.

I have made a draft map with OS open background and 16 interchangeable layers using QGIS2Web. However, it is taking around 6-8 minutes to load and is very slow.

I think this is due to the size of the underlaying OS map? Does anyone know how to crop the extent of the map in QGIS2Web?


texts: "Horizontal" and vertical datums terminology

texts:
I have heard people referring to as a horizontal datum to models such as WGS84 and as a vertical datum to geoids. The former is primarily used for determining latitude and longitude but also the height (and thus X, Y Z can be calculated) above the ellipsoid while the latter is used for determining height above an equal-potential surface only.

I wonder why one would call "horizontal" a full 3D datum. The fact that you have a pure vertical datum (i.e. geoid) does not mean you have to call "horizontal" the remaining one.

Here are two articles where you can find the term:

http://geography.about.com/od/geographyintern/a/datums.htm
http://www.esri.com/news/arcuser/0703/geoid1of3.html


texts: QGIS Installation on Mac - Mac doesn't accept QGIS as a trusted developer

texts: I am trying to install QGIS on a MacBook, using MacOS High Sierra Version 10.13.6. For some reason, having installed Python 3.6, my Mac refuses to use the installer saying that QGIS is not a trusted developer. Is there a solution to this?


texts:
I have a GRD-File where all pixelvalues which are less than 0.05 should be converted to nodata.

I know how to do this with ArcMap, but the process should be done automatically, for example with cronjobs.

Is there a possibility to do that with gdal (maybe gdal_warp or gdal_translate)?


texts: Selecting MultiPart Polygons if there are buildings in at least two of the parts

texts:
my data situation:

shape of polygons (some are MultiPart Poylgons)
shape of building locations (Points)

What I need:
Select all polygons that have two ore more parts and contain one building or more in at least two of the parts.

What I did so far:

I added and calculated a field that tells me if a poylgon has multiple parts.
Field calculator Python: !shape.isMultipart!
I calculated of how many parts each MultiPart Polygon consists.
Field calculator Python: !Shape!.partCount

Now I have all the MultiPart Polygons with two or more parts selected.
From this point, how can I select only the polygons that contain a building in at least two of their parts?

PS: I'm working with ArcGIS


texts: Replacing attributes of raster cells with other raster values?

texts:
I have two raster layers.  The first is a terrestrial habitat raster that includes an open water value.  The second is a lentic habitat raster.  I would like to "burn in" the lentic habitat raster values to the terrestrial habitat raster anywhere that it overlaps a cell with an open water value.  I am close to completing this with the conditional tool, but all open water cells in the terrestrial raster that are not overlapped by cells in the lentic raster are being replaced by no_data.  I understand why this is happening, but I am having trouble figuring out how to modify my workflow to get this done. Please let me know if you have any suggestions for a modified workflow or a different method.  Here's what I'm currently doing:

1. Expand tool on the lentic raster to ensure lentic habitats are at least the width extent of the open water classes
2. Conditional tool with parameters shown in first figure

And below is a snapshot of the output with black cells representing cells that originally had open water values but now have no_value.


texts: Fixing non-noded intersection problem using PostGIS

texts:
I'm using a PL/R function and PostGIS to generate voronoi polygons around a set of points. The function that I am using is defined here. When I use this function on a particular dataset I get the following error message:

Error : ERROR:  R interpreter expression evaluation error
DETAIL:  Error in pg.spi.exec(sprintf("SELECT %3$s AS id,
st_intersection('SRID='||st_srid(%2$s)||';%4$s'::text,'%5$s')
AS polygon FROM %1$s WHERE st_intersects(%2$s::text,'SRID='||st_srid(%2$s)||';%4$s');",
:error in SQL statement : Error performing intersection: TopologyException: found non-noded
intersection between LINESTRING (571304 310990, 568465 264611) and LINESTRING (568465
264611, 594406 286813) at 568465.05533706467 264610.82749605528
CONTEXT:  In R support function pg.spi.exec In PL/R function r_voronoi

From examining this part of the error message:

Error performing intersection: TopologyException: found non-noded intersection between
LINESTRING (571304 310990, 568465 264611) and LINESTRING (568465 264611, 594406 286813)
at 568465.05533706467 264610.82749605528

This is what the problem listed above looks like:

I initially thought that this message might be caused by the existence of identical points, and tried to solve this using the st_translate() function, used in the following way:

ST_Translate(geom, random()*20, random()*20) as geom

This does fix the problem, but my concern is that I am now translating all points up to ~20m in the x/y direction. I also can't tell what an appropriate translation amount is need.

For example, in this dataset through trial and error a 20m * random number is ok, but how can I tell if this needs to be bigger?

Based on the image above I think the problem is that the point is intersecting with the line while the algorithm is trying to intersect the point with a polygon. I'm not sure what I should be doing to ensure that the point is within a polygon, rather than intersecting with a line. The error is occurring on this line:

"SELECT
  %3$s AS id,
  st_intersection(''SRID=''||st_srid(%2$s)||'';%4$s''::text,''%5$s'') AS polygon
FROM
  %1$s
WHERE
  st_intersects(%2$s::text,''SRID=''||st_srid(%2$s)||'';%4$s'');"

I have read through What is a "non-noded intersection"? to try to better understand this problem.


texts: Missing buildings on London Mapzen OSM Metro Extract

texts:
Please could someone explain why about 25% of the buildings are not showing on this map of Central London (eg between Hyde Park and the Thames). I have just downloaded the data, as a shapefile, and styled the buildings purple.

I am fairly new to QGIS (and to Stacexchange). The map was produced by:

downloading the data as a Shapefile from mapzen.com
using Layer>Add Vector Layer to add the shapefile to the project
selecting the Buildings layer and changing the style to purple.

I would like to know how to show all the buildings.

[Note: I have tried several ways of making a map with OSM data (Geofabrik, OSM download etc) and am finding that each method has problems I can't resolve. My previous GIS experience was with ArcView, quite a long time ago, and I do not remember problems of this type].


texts: Dynamically assign feature properties on L.marker popup interactions

texts:
I have a Leaflet based map on which a user can click to create a marker. Then, if one clicks on the marker a popup is opening, where one can change text in a textarea, activate the camera of the mobile device and save all markers with feature properties into a geoJSON.

I initialize an empty L.featureGroup(), generate dynamic HTML popupContent inside an addMarker(e){} function and add all this data to an object, that is being created via .toGeoJSON(), where I can access feature.properties. Because it is overwritten each time, I have to write all properties into another array infoTextArray

Is there a efficent way to directly write a feature's properties on creation?

edit: to further specify my question:

I want to assign my geoJSON features properties when they're created. I tried something like:

let layerGroup = L.geoJSON([], {
pointToLayer: function (feature) {
    feature.properties.infoText='';
    feature.properties.imageURL='';
},
onEachFeature: function (feature) {
    feature.properties.infoText='';
    feature.properties.imageURL='';
}
}).addTo(map);

layerGroup.toGEOJSON() still returns an empty feature properties object. I'm guessing my functions aren't triggered, when a new features is added to the layer.


texts:
Has anyone some examples of  info-windows showing highcharts.js in it ??

something like this but using carto data . I also found  this example but I need an interactive chart .


texts:
Please help a beginner in R

I am trying to do some regression analysis in R  between pairs of precipitation and NDVI rasters. I got the regression to work, and now I would like to loop this over all the rasters, pair-wise. All my data is labelled NDVI0101, pre0101, NDVI0102, pre0102.....

Is there a way to loop so that the regression is preformed between two rasters with the same date, and which also saves the output residual-file with the date?

Here is my code:

library(raster)
library(sp)
library(rgdal)

#set working directory
setwd("C:/users/neija/documents/Rworkingfolder")

#load raster
ndvi <- raster("rdata/ndvi/ndvi0101.tif")
pre <- raster("rdata/precip/pre0101.tif")

#extract values from raster into a vector
y <-getValues(ndvi)
prere <- resample(pre, ndvi, method="bilinear")
x <-getValues(prere)
y[which(is.na(x))]<-NA

#combine in data frame
f<-data.frame(x,y)

model<-lm(x~y,data = f)
summary(model)

#extract residuals
r <- ndvi
r [] <- residuals.lm(model)
writeRaster (r, file = "res0101.asc")

Any suggestions?


texts: How could I avoid extracting incorrectly closed contours from a raster with angled sides?

texts:
I'm updating a Python script which currently calls gdal_contour as a subprocess to extract contours. I'd like to replace gdal_contour with skimage.find_contours.

However the output of the GDAL script differs slightly from scikit-image's output when used on a raster with angled sides. Here are a few images to illustrate:

Before extracting contours, my source raster is first masked (via rasterio.mask.mask) with a rectangular polygon, producing something like this:

When I use skimage.find_contours, it produces open contours (which is what I want) when the desired height disappears out of bounds, as in the bottom of the right edge below:

But when the mask is rotated, like this:

The contours produced by skimage.find_contours now follow the edge of the raster when they pass out of bounds and are closed, as in the bottom-right corner.

(Contours produced by GDAL from the same raster are open in the same corner.)

I guess I could add a margin, then clip the contours to the correct bounds, but is there another way to deal with this situation?

(I'm open to alternatives to scikit-image, see my previous question here.)


texts: Impossible to save a layer as a shapefile (shp)

texts:
I'm trying several spatial statistics analysis in my QGIS 2.18.3 project (convex hull, random points inside polygons, point sampling tool, basic statistics for numeric fields, distance matrix, nearest neighbour analysis, mean coordinates, polygon centroids).

Since I used the nearest neighbour analysis, which is saved under a .html in my data folder, I cannot save any other layer as a shapefile (shp). The proposed type of files are just .dbf, .xlsx, .csv and .ods.  Encoding is set as System. Pre-selected type of file is .dbf.

Result is that mean coordinates and polygon centroids are not visualised in the project.

I tried to create a fix buffer that I already did this morning and even that layer was not able to be saved under a shapefile, but in the morning it was possible.


texts: Dissolve overlapping polygons according to percent overlap using PostGIS?

texts:
I have a shapefile with convex hull polygons, all in the same layer/shapefile. Some of the polygons overlap (not duplicate polygons). I would like to dissolve/union the polygons which overlap with more than a specified percentage (65%).

Any idea how to go about solving the issue using PostGIS?


texts: Mapbox GL Expression to filter array in property

texts:
I have GeoJSON Feature properties that look like this:

{
  "type": "FeatureCollection",
  "features": [
    {
      "type": "Feature",
            "properties": {
        "id": "someID",
        "lines": [
          {
            "name": "Foo",
            "anotherProperty": 1
          },
          {
            "name": "Bar",
            "anotherProperty": 2
          }
        ]
      },
      "geometry": {
        "type": "LineString",
        "coordinates": [
          [
            5.9765625,
            48.45835188280866
          ],
          [
            40.42968749999999,
            47.517200697839414
          ]
        ]
      }
    }
  ]
}

And I'd like to select features that match 'name': 'Foo' in their lines array. But I'm struggling to write a match expression that will check each entry in 'lines'.


texts: Creating polygon outside buffer according to its area using R

texts:
I have a lot of polygons (shapefile) and need to creat their own buffer according to the area of each polygons,despite I have found some infomation in Creating OUTSIDE_ONLY buffer around polygon using R?.
# create basic function
outerBuffer<-function(x, dist){
  buff<-buffer(x, width =dist - 1, dissolve = F)
  e<-erase(buff,x)
  return(e)
}

library(maptools)
library(rgdal)
# create data
p1 = shapefile('42gp.shp')

# apply function, create only outside buffer
a2<-outerBuffer(p1,500)

# highlight isobath with the buffer and add outline
plot(buf, outline=FALSE, n = 10, col = 2, lwd=.4)
plot(buf, lwd = 0.7, fg = 2)
# }

# tradaaa ! :)
plot(a2, col = "white", main= "Outer buffer only")

What can I write the R script for creating their own buffer according to the area of each polygons?


texts:
If there are a way to do

call("gdal_contour -i 10 -f GeoJSON".split(' ') + [file1, file2])

as Python call method of some library? Not via command line call and intermediate file storing?


texts:
I have been trying to grid a point cloud, version1.2, LiDAR datasets which are intended to be used for hydraulic simulation of flood propagation. Already, I have used Arc* to work on the traditional (.las) format, to create lasdataset. Then, a raster dataset was created from the lasdataset and finally interpolated by means of a nearest neighbour algorithm.
The interpolated raster data were exported as txt files.

Having gone through this routine, I expected a result with header file as well as regular values in X, Y and Z orthogonality. Instead, the result I got was distorted sets of (real) numbers I could not make sense of.  I tried using Surfer* to grid txt files created from (.las) point clouds, by means of lastools*, but the result was a binary file and I couldn’t read it. Lastools* gave a reasonable result when used to grid the point clouds, only that the cell size varies inversely with the number of cells that result in the output grid.

My challenge is how to grid the point cloud LiDAR and specify the number of rows and columns to be contained in the resulting gridded data, giving a fixed cell size.


texts: Adding Not Operator in Definition Query of ArcGIS Pro?

texts:
I am using ArcGIS Pro 2.2

I want to apply a definition query saying

NOT (name contains 'road%')

apparently there is no option to add Not Operator in query neither there seem to be any option saying not contains.

Is there a way to add Not Operator in definition query using ArcGIS Pro?


texts: implement aspect in least cost path analysis to account for downhill and neutral ground movement

texts:
I am creating a least cost path model for large species that migrate over large distances across a rugged landscape  (elk, deer, etc.)

To calculate movement cost I am using a slope raster but I realize that I should also take into consideration aspect to figure out if an animal is moving downhill, uphill or just along a contour. The idea is to have no penalty for moving downhill.

I know ArcGIS has a least cost path analysis but I am not sure how to implement direction of travel and aspect. This paper http://www.sciencedirect.com/science/article/pii/S0168159100001015 explains it well but does not give any details on how it was implemented ArcGIS. Is there something already built-in ArcMap or R that I can use?


texts: Using Python script to control GRASS GIS from outside?

texts:
I am working with GRASS 6.4.3 and python 2.7.6 under WIN 8.1.

I am a novice with computer science and coding, and I have seen many posts including those:

Is it possible to use GRASS GIS in Python stand alone scripts?
GRASS and Python, in order to control GRASS from outside by using python.

Here's my python code, which wants to use the GRASS module g.list to list my raster data in GRASS. I type those code in Python IDLE.

import os
import sys

gisbase = os.environ['GISBASE'] = 'C:\OSGeo4W64\apps\grass\grass-6.4.3'  #GISBASE needs to point the root of the GRASS installation directory
gisrc = 'C:\Users\Heinz\Documents\grassdata'
gisdbase = 'C:\Users\Heinz\Documents\grassdata'
location = 'newLocation'
mapset = 'TC'
LD_LIBRARY_PATH = 'C:\OSGeo4W64\apps\grass\grass-6.4.3\lib'
PATH = 'C:\OSGeo4W64\apps\grass\grass-6.4.3\etc';'C:\OSGeo4W64\apps\grass\grass-6.4.3\etc\python';'C:\OSGeo4W64\apps\grass\grass-6.4.3\lib';'C:\OSGeo4W64\apps\grass\grass-6.4.3\bin';'C:\Python27';'C:\OSGeo4W64\apps\Python27';'C:\OSGeo4W64\apps\msys'
PYTHONLIB = 'C:\Python27'
PYTHONPATH = 'C:\OSGeo4W64\apps\grass\grass-6.4.3\etc\python'
GRASS_SH = 'C:\OSGeo4W64\apps\msys\bin\sh.exe'

sys.path.append(os.path.join(os.environ['GISBASE'], 'etc', 'python'))

import grass.script as grass

And I have set every path in my code as environment variable in the windows control panel.

As I run module, I got the error:

Which line of my code is wrong or not necessary and what code should I add to my python script, or can't I use Python IDLE to do this?

@ustroetz I have found no file in my computer named .grassrc6 , so I create a new file named .grassrc6(which is empty), adding a line in my code:

gisrc = 'C:\Users\Heinz\.grassrc6'

But still got the same error.

Did I misunderstand the instruction in the post？

I create a text file named grass.pth in this folder C:\OSGeo4W\apps\python27\lib\site-packages, and type those two lines in it:

Here is the list of my environment variables I have set in windows penal:

GISBASE = C:\OSGeo4W64\apps\grass\grass-6.4.3
GISRC = C:\Users\Heinz\Documents\grassdata
Path =
C:\Python27
C:\OSGeo4W64\apps\grass\grass-6.4.3
C:\Users\Heinz\Documents\grassdata
C:\OSGeo4W64\apps\grass\grass-6.4.3\lib
C:\OSGeo4W64\apps\grass\grass-6.4.3\etc\python
C:\OSGeo4W64\apps\grass\grass-6.4.3\etc
C:\OSGeo4W64\apps\grass\grass-6.4.3\bin
C:\OSGeo4W64\apps\Python27
C:\OSGeo4W64\apps\msys
C:\OSGeo4W64\apps\grass\grass-6.4.3\scripts
C:\OSGeo4W64\bin
C:\OSGeo4W64\apps\Python27\Lib\site-packages

then I ran my python script, and got the same error:

You can see that I have added this line:sys.path.append('C:\OSGeo4W\bin') in my script(it's because the script didn't work before I added it), but it still got the same error.

@Martin, I have found that I have three places where python installed:

C:\OSGeo4W64\apps\Python27
C:\Python27
C:\Program Files (x86)\Quantum GIS Lisboa\apps\Python27

and I also found site-packages files in those places, so I copied grass.pth to those places, but my script still can't work.

I also found this python-related file:

C:\Users\Heinz\Downloads\http%3a%2f%2fdownload.osgeo.org%2fosgeo4w%2f\x86_64\release\python

I don't know if this keeps my script useless, should I delete this file, or it doesn't matter?

I have tried this in cmd:

I will keep trying to make my script work!


texts: Seeking shapefile of 2km DINTY Tetrad squares for UK?

texts:
does anyone know where I can download a shapefile of 2km DINTY Tetrad squares for the UK please?  A google search wasn't very helpful so I thought I'd try here.

A Tetrad is derived from the British National Grid, but is not the same as the British National Grid, so this is not a duplicate question.


texts: Deleting many fields from many feature classes using ArcPy?

texts:
I'm currently formatting our GIS data for use in a mobile app. The purpose of altering the data is to reduce the output file size that's imported into the app. I need to delete many fields from many feature classes. I will need to run this script about every two months to provide updated data.

The first part of my script copies the feature classes from our SDE geodatabase to a file geodatabase using arcpy.FeatureClassToFeatureClass_conversion(). This function has a field_mapping parameter which I'm not very familiar with.

What would be the quickest way to delete many fields from many feature classes?

Using the field_mapping parameter, or just using arcpy.DeleteField_management()? I would think field_mapping.

The real problem is that I want to remove around 30 fields from each feature class, and only keep around 4-5.

With field_mapping, (I haven't researched this yet) is there a way to simply state only the fields I want to keep without listing all the fields I want to remove too?


texts: Creating overview map in ArcGIS Pro that refers to extents of map frames in other layouts

texts:
I'm putting together a map book and would like an "overview" map which will show the extents of the maps that will follow it. I've looked at the Extent Indicators help page here:
https://pro.arcgis.com/en/pro-app/latest/help/layouts/extent-indicators.htm
The help page only gives instructions on how to set up Extent Indicators that are referring to another map frame on the same layout.
Is there a way to refer to a map frame on a separate layout in the same project or any workarounds to achieve the same outcome?


texts: ArcGIS JS API: Convert ESRI Date Type into String

texts:
I am following this ESRI tutorial to display information from a related table for the selected features. The example in itself for my project works fine, however, my data include two Date fields and when I display them on the right side grid, I see something like "1237507200000" as starting date and "1237161600000" as ending date.

I understand I have to format this esriFieldTypeDate so I can display it as a text field like "YYYY-MM-DD", but I don't know how. It is important to highlight that this related records are stored in a dojo.data.ItemFileReadStore as in the example.

Please, does anyone know how to do this?

Thanks!


texts: Location Allocation (lines) to Polygon map?

texts:
I have a solved "Location-Allocation" Analysis, providing a map result like this.

What methods do I have, to help transform this into a segmented polygon areas?  Current doing it by hand, but have the entire nation to do!


texts: Adding perpendicular lines to points along a polyline

texts:
I'm looking for a method to add perpendicular lines at points along a polyline. Essentially I'm searching for a solution in QGIS or Python (I'm aware of an ArcGIS solution available here: http://www.arcgis.com/home/item.html?id=3f05ef15203e4a3ca88fc1d92347512b.

A similar question (How to add a perpendicular line to the start & end nodes of every polyline?) addresses what I'm looking to achieve, but again, I need a QGIS/Python solution.

Does anybody have any advice? Currently I'm using the QChainage plugin to generate equally spaced points across lines (see below), but can't find a method to generate perpendicular lines from the points.


texts: gdal_rasterize failure with QGIS version 3.6.3

texts:
I have both QGIS 3.6.3 and 3.0.1 running on Windows 10 64bit. In 3.6.3 I am unable to rasterise a point layer, while the same operation seems to work in 3.0.1. I'm having to switch between versions just for this one operation....

Logfile from 3.6.3:-
QGIS version: 3.6.3-Noosa
QGIS code revision: 0c5774c068
Qt version: 5.11.2
GDAL version: 2.4.1
GEOS version: 3.7.2-CAPI-1.11.0 b55d2125
PROJ version: 520
Processing algorithm…
Algorithm 'Rasterize (vector to raster)' starting…
Input parameters:
{ 'BURN' : 0, 'DATA_TYPE' : 5, 'EXTENT' : '163952.5,185652.5,9144881.5,9171981.5 [EPSG:XXXXX]', 'FIELD' : 'field_3', 'HEIGHT' : 25, 'INIT' : None, 'INPUT' : 'file:///Y:/Projects/2019/XXXX/Analysis/ExportData.txt?type=csv&delimiter=%5Ct&useHeader=No&trimFields=Yes&skipEmptyFields=Yes&xField=field_1&yField=field_2&crs=EPSG:XXXXX&spatialIndex=yes&subsetIndex=no&watchFile=no', 'INVERT' : False, 'NODATA' : 0, 'OPTIONS' : '', 'OUTPUT' : 'Y:/Projects/2019/XXXXX/QGIS/CoverageV6A.tif', 'UNITS' : 1, 'WIDTH' : 25 }

GDAL command:

gdal_rasterize -l None -a field_3 -tr 25.0 25.0 -a_nodata 0.0 -te 163952.5 9144881.5 185652.5 9171981.5 -ot Float32 -of GTiff file:///Y:/Projects/2019/XXXXX/Analysis/ExportData.txt?type=csv&delimiter=%5Ct&useHeader=No&trimFields=Yes&skipEmptyFields=Yes&xField=field_1&yField=field_2&crs=EPSG:XXXXX&spatialIndex=yes&subsetIndex=no&watchFile=no Y:/Projects/2019/XXXXX/QGIS/CoverageV6A.tif

GDAL command output:

FAILURE: No output file specified.

Usage: gdal_rasterize [-b band]* [-i] [-at]

{[-burn value]* | [-a attribute_name] | [-3d]} [-add]

[-l layername]* [-where expression] [-sql select_statement]

[-dialect dialect] [-of format] [-a_srs srs_def] [-to "NAME=VALUE"]*

[-co "NAME=VALUE"]* [-a_nodata value] [-init value]*

[-te xmin ymin xmax ymax] [-tr xres yres] [-tap] [-ts width height]

[-ot {Byte/Int16/UInt16/UInt32/Int32/Float32/Float64/

CInt16/CInt32/CFloat32/CFloat64}] [-optim {[AUTO]/VECTOR/RASTER}] [-q]

<src_datasource> <dst_filename>

'delimiter' is not recognized as an internal or external command,

operable program or batch file.

'useHeader' is not recognized as an internal or external command,

operable program or batch file.

'trimFields' is not recognized as an internal or external command,

operable program or batch file.

'skipEmptyFields' is not recognized as an internal or external command,

operable program or batch file.

'xField' is not recognized as an internal or external command,

operable program or batch file.

'yField' is not recognized as an internal or external command,

operable program or batch file.

'crs' is not recognized as an internal or external command,

operable program or batch file.

'spatialIndex' is not recognized as an internal or external command,

operable program or batch file.

'subsetIndex' is not recognized as an internal or external command,

operable program or batch file.

'watchFile' is not recognized as an internal or external command,

operable program or batch file.

Execution completed in 0.18 seconds
Results:
{'OUTPUT': 'Y:/Projects/2019/XXXXX/QGIS/CoverageV6A.tif'}

Loading resulting layers
The following layers were not correctly generated.

Y:/Projects/2019/XXXXX/QGIS/CoverageV6A.tif

You can check the 'Log Messages Panel' in QGIS main window to find more information about the execution of the algorithm.


texts: In ModelBuilder, can I reclassify the first set of the result of previous step?

texts:
I am trying to create this model and run it on multiple data sets. It's about interpolation of certain signals (electromagnetic data). So what I want to do is to get the first set of the interpolation results: the lowest numerical value set, my targeted set. Normally Arcmap will make it into  sets, so the first set can be accurate enough for my interpolation.

Here is a glance of my model.

As a matter of fact, after I got the first set of the reclassify result, which is basically isolation, I want to get the centroid of that area.

But now the biggest problem seems to be: the reclassify option will always retain the value of last time. As shown below


texts: Calculate band using data from existing bands and look-up table in Google Earth Engine

texts:
I have a 4 band composite image and a csv table consisting of regression coefficients. Each row in the table represents a regression specific to a day of year (doy). The 1st column in the table is doy and Band 2 of the composite is doy

What I would like to do is calculate a new band using the coefficients specific to the doy of each pixel.

The process should be something like this:

Find doy value in pixel.
Look-up the appropriate row in the coefficient table with matching pixel doy.
Extract appropriate values from columns in coefficient table.
Extract appropriate values from other bands in pixel.
Plug values into equation to calculate new band.

The Equation can be something like:

Col2 + Col3*Band3 + Col6*Band3*Band1

How can I implement this in Google Earth Engine?


texts: What are the most detailed and reliable address databases for the United States?

texts:
Here in Britain people can license the Ordnance Survey's AddressBase product, which provides a list of almost all the country's addresses, along with individually surveyed building locations. Like all geodata products it's got some errors and omissions, but it's generally good - even in remote areas. This is useful for address suggestion and geocoding (some of the major parts of what I want to achieve).

I'd like to broaden my knowledge beyond just the UK. What's the best data available for the United States, and how good is it?

I'm aware that the USPS has a "zip+4" file that contains 30 million addresses (presumably not a full list of addresses in a country of 300 million people); that Navteq's map data contains street number ranges (in the UK their accuracy is OK but not the best); that Google have a geocoding API (I would prefer a database rather than an API); and that TIGER data contains some address information (I haven't fully investigated this - is it any good?). Should I be looking at one of these options in more detail, or are there other sources I've missed?


texts: Refreshing map canvas in for loop in QGIS 3

texts:
I have a console script where I'm trying to iterate through saved bookmarks from the open project to have the map pan/zoom to the bookmark extent.  I would like to eventually save each map canvas as pdf file.  I can manually set the right parameters to cycle through the bookmarks and the canvas refreshes each time I set the extent.
However, if I cycle through the saved bookmarks and set the canvas extent in a for loop (as below), the canvas does not refresh with each extent change, it only refreshes once the script has run it's course.  I've added sleep commands in the for loop to allow for redraw time, but that didn't work.  Am I missing something simple?
bmarks = QgsApplication.bookmarkManager()

bmarks.bookmarks()

for b in bmarks.bookmarks():
    canv = iface.mapCanvas()
    canv.setExtent(b.extent())
    iface.mapCanvas().refresh()


texts: Storing metadata in Geopackage for QGIS?

texts:
Is there a way to store metadata in a Geopackage along with the layer data?

I'm using QGIS 2.18.24. When you add or edit metadata via the "Metadata" pane of the layer's Properties, this is saved in the project file, not the geopackage.

There is an open issue that states the "metadata data panel is still a work in progress" which then links to an enhancement proposal on github. Perhaps this issue is already addressed in QGIS 3. Regardless, is there some kind of workaround? Maybe there's another program that can open the gpkg and add metadata in a format that QGIS will still read?


texts: Why QGIS make request to WMS in wrong projection?

texts:
On GeoServer I have published a layer in custom projection I called it EPSG:60900. When I open this layer in Leaflet or OpenLayers it's works fine. In GeoServer log I see that client sends requests in my projection:

Request: getMap
    Time = []
    Buffer = 0
    Format = image/png
    Filter = null
    Filters = null
    Height = 752
    Width = 768
    RemoteOwsType = null
    RemoteOwsURL = null
    FormatOptions = {}
    CQLFilter = null
    Elevation = []
    FeatureId = null
    StartIndex = null
    ViewParams = null
    Palette = null
    Tiled = false
    SRS = EPSG:60900
    FeatureVersion = null
    Styles = [StyleImpl[ name=line]]
    Layers = [org.geoserver.wms.MapLayerInfo@2d3b64a0]
    MaxFeatures = null

Now I want to open this WMS layer in QGIS. At first I add SRS definition in QGIS. And open my layer. And here a problem. In moment where I add WMS layer QGIS does not allow me set custom SRS to layer, so I choose Mercator.

After in layer properties I set layer and project SRS as my custom coordinate system.

But QGIS sends request in another SRS:

 Request: getMap
    Time = []
    Buffer = 0
    Format = image/png
    Filter = null
    Filters = null
    Height = 1195
    Width = 1236
    RemoteOwsType = null
    RemoteOwsURL = null
    FormatOptions = {DPI=96}
    CQLFilter = null
    Elevation = []
    FeatureId = null
    StartIndex = null
    ViewParams = null
    Palette = null
    Tiled = false
    SRS = urn:x-ogc:def:crs:EPSG:3857
    FeatureVersion = null
    Styles = [StyleImpl[ name=line]]
    Layers = [org.geoserver.wms.MapLayerInfo@2d3b64a0]
    MaxFeatures = null

What can be wrong?


texts: Simple PostGIS ST_Distance query for two points from the same geometry table?

texts:
Note: This question was published at Stackoverflow at first.

I have a geometry table points_table with SRID 4326 and a few points in it.

I need to find a distance between a pair of points (features) from this table.

I would like to select one point with id = X and another point with id = Y and then find a distance between them with ST_Distance function.

I know how to query these two points from the database:

select * from points_table where id = X or id = Y;

But this query returns a result, which I do not know how to pass to ST_Distance, because ST_Distance takes 2 parameters.

I have tried Select ST_Distance(select the_geom from points_table where id = X or id = Y);

But it is obviously wrong.

How can I find the distance between two geometry points (features) from a table?

P.S. All the examples I have found so far show the ST_Distance usage with points hardcoded into the arguments of the ST_Distance function. I have not found an example, where ST_Distance is used to calculate a distance between two points from the same table.


texts: GRASS GIS: From multiple lines to single line

texts:
This seems like a simple process but somehow I am not able to solve this.
The situation is as follows:
I have three separate lines which form visually 1 line. As I want to further process these lines they have to be "merged" into 1 line from a topology side of view.

The lines are build from a thinned raster. Nodes numbers are counted anew for each line.

The attribute table is as follows:

Layer 1
cat: 1-3
value: length

I have tried the following tools:

v.edit to try to merge
v.build.polylines
v.dissolve

I might have used to wrong parameters or did not select certain options in the right way. I tried to search for similar topics and tried several solutions but as of yet to no avail.


texts: OpenLayers Circle Style Update During zoom

texts:
I am using a stylefunction to dynamically draw an 'image' based on the feature geometry. The style function is dependent on the resolution (zoom level) so that the drawn items are a certain consistent real-world size. Here is an example of what I am doing. Is there a way to have the orange circles scale with the map during zoom (or at a minimum hide during zoom while other features remain).

Here is a code sample:

var styleFunction = (feature, resolution) => {
  return [
    new Style({
      stroke: new Stroke({
        color: "blue",
        width: 3
      }),
      fill: new Fill({
        color: "rgba(0, 0, 255, 0.1)"
      })
    }),
    new Style({
      image: new CircleStyle({
        radius: 500000 / resolution,
        fill: new Fill({
          color: "orange"
        })
      }),
      geometry: function(feature) {
        // return the coordinates of the first ring of the polygon
        var coordinates = feature.getGeometry().getCoordinates()[0];
        return new MultiPoint(coordinates);
      }
    })
  ];
};


texts: Non-random sampling of point vector layer in QGIS

texts:
I have a point vector layer with 47.000.000 points representing elevation. It is a regular grid with constant spacing. I aim for rasterization but it is very slow.
I know how to do random sampling but have not found a grid sampling. (not on this stack either)

How to create a more spaced grid, based on the initial layer's values?


texts: How to map and color localities (provinces) with an external CSV using Quantum GIS

texts:
I use a Swiss shape file map (called PLZO_PLZ) available on http://www.cadastre.ch/internet/cadastre/fr/home/products/plz/data.html (ZIP file Shape LV03)

PLZO_PLZ contains a column called PLZ which correspond to a swiss locality. For instance number PLZ 1202 is in Geneva. By the way, PLZ can also be called NPA (same meaning).

My aim is to get a colored map according to the NPA frequency coming from an external CSV file: http://pastebin.com/raw.php?i=iAjV88KQ

For instance NPA not appearing in the CSV should be not colored (white) and NPA coming few time should appear in light blue and the more frequent the NPA is the more dark blue should be the map. According to my CSV NPA 1203 should be the darkest.

Is it possible to obtain what I am looking for with Quantum GIS? (I tried to use "Joins" function in proprieties but it's not working)

Edit:
I try to join like this :


texts:
I have this original image

But when i open it in TileMill, it converts image like this. White areas are transparent in the image

This is more like, it removes all no data values i.e. 255 in my original image. How can i set custom no data value with TileMill and also set RGB bands null value to null. I don't want RB


texts: Looping through a table to allocate kids to nearest school

texts:
The general problem is to allocate town's kids to their nearest school and get the distance and the geometry.

Tools in use is Postgresql 10.6, PostGIS 2.4.3 and PGRouting 2.5.2

There are three tables:
school_locations (n=5)
kids_locations and (n=1500)
roads

So far the schools and the kids have be allocated to their nearest network nodes and the nearest school and the distance have been calculated.

The pairing results are stored in a table that can be seen below.
To get this result pgr_dijkstraCost has been used. The drawback is that it does not return the geometry. So another algorithm has to be used to create a new table or update a geom field in the result table mentioned above.

Here my SQL so far:

CREATE TABLE shortest_path AS (
SELECT
    kids.id AS kid_id
    , subquery.start_vid AS source
    , school.id AS school_id
    , subquery.end_vid AS target
    , subquery.agg_cost
FROM (
    SELECT
        row_number() OVER (PARTITION BY
                start_vid
            ORDER BY
                agg_cost ASC) row_num
            , *
        FROM
            pgr_dijkstraCost ('SELECT id, source, target, length as cost FROM roads WHERE length >= 0'
                , (
                    SELECT
                        array_agg(node_id)
                    FROM
                        kids_location)
                    , (
                        SELECT
                            array_agg(node_id)
                        FROM
                            school_location)
                        , FALSE)
                ORDER BY
                    start_vid
                    , agg_cost) AS subquery
                INNER JOIN kids_location AS kids ON (kids.node_id = start_vid)
                INNER JOIN school_location AS school ON (school.node_id = end_vid)
            WHERE
                row_num = 1);

The resulting table looks like this:

kid_id | source | school_id | target | agg_cost
   523 | 164933 |         1 | 156849 |     5202
   522 |  43687 |         1 | 156849 |     4470
   521 | 115352 |         1 | 156849 |     4088
   520 |  64630 |         2 | 130234 |     4748
   519 | 114828 |         2 | 130234 |     4138
   518 |  82010 |         2 | 130234 |     5443
   517 |  87980 |         2 | 130234 |     6237

The problem is that I know how to calculate the geometry for one pair but I do not know how to do it for the entire table. I guess looping should solve the problem, but I do not have the knowledge to do so.

SQL for one pair of kid<->school:

SELECT
    seq
    , edge
    , b.geom
FROM
    pgr_dijkstra ('
               SELECT id, source, target,
                       length as cost FROM roads'
    , 43687
    , 156849
    , FALSE) AS a
INNER JOIN roads AS b ON (a.edge = b.id)
ORDER BY seq;

This would be the hardcoded version. I could construct n version of this SQL but looping through the table and setting source and targed would be the better solution. I imagine this is common problem but I could not find a solution so far. If anyone could help me with code please.


texts: Code to print in the Qgis console the values of the attribute table that are selected in the shapefile

texts:
I found a code that prints all data that a shapefile has in the attribute table in the qgis console. But i wanted to have a code that just printed the data that i select in the shapefile. Can someone help me?

I'll leave here the code that prints all the data.

lyr = iface.activeLayer()

features = lyr.getFeatures()

featList = []

selected_features = layer.selectedFeatures()

for feat in features:

featList.append(feat.attributes())

featList.sort()

for f in featList:

print f


texts:
I have been trying to get to some county data from geocomm.com but the website has been down for several days.

Does anyone have information regarding their status?


texts: Where is it possible to download a DEM for the Romanian east carpathians?

texts: I'm quite new to GIS and have only done a few small tasks for Projects at University and the required data has always been provided. I am now struggling to find the specific DEM I need. Also other terrain data such as presence of glaciers (or former glaciers) and river networks would be useful.


texts: Getting data driven pages to apply gamma stretch for rasters from current display extent?

texts:
I'm trying to export a single, multipage PDF using Data Driven Pages in ArcGIS 10.3. The background images for my maps are 15 large, high res overlapping rasters that are sometimes poor quality due to being underexposed or overexposed, but only in parts of each raster. I've dealt with this problem by going into layer properties, Symbology, and choosing a "Percent Clip" stretch and applying a gamma stretch from the current display extent. That way, my raster colors appear quite good on each page. When I click through the data driven pages manually, the rasters appear nicely rendered, with bright whites and dark blacks for each image.

However, when I export the map using data driven pages using the GUI (I'm a python coding noob), those preferences are not being respected in the final pdf document after exporting via data driven pages.

I thought the solution on the bottom of this page might work, because it forces each page to refresh before printing it to the final pdf, but my scripting knowledge leaves a lot to be desired and I'm having trouble applying the code to my own document.

It sounds like the solution might lie in using python scripting to force data driven pages to draw each page separately (using mxd.dataDrivenPages.refresh() ?), then append all the pages together in a final document.

I've messed around with script I found on these questions,

How to Prevent Data Driven Pages From Hanging on Subsequent Export?
https://support.esri.com/en/technical-article/000014465
https://community.esri.com/thread/19604

but I just can't get it to work.

Has anyone has been able to get the stretched renderer to work in data driven pages?


texts: Reconstructing MODIS time-series applying Savitzky-Golay Filter with Python/Numpy

texts:
I want to apply the 'Savitzky-Golay' (savgol) filter to my time series, MODIS dataset, to remove noise (i.e., cloud pixels, etc.) in my data. MODIS have quality flags that indicates the reliability of each pixel values or if the pixel is possibly affected by clouds. So I would like to incorporate these quality flags in my filter by putting less weight or ignoring those pixel values and let savgol filter predict the optimal pixel value. I am testing np.NaN/np.nan/isnull but it seems that it removes the element in the array, and consequently savgol filter also skip those values. I would like my resulting data to be like in the attached figure.
 (https://matinbrandt.wordpress.com/2014/12/02/smoothingfiltering-a-ndvi-time-series-using-a-savitzky-golay-filter-and-r/)


texts: Writing gridded NetCDF format CF Convention from a 4 dimensional Array in Java

texts:
I am using Unidata NetCDF Java library to write NetCDF, its being written in the classical format, but I want to write it as a gridded NetCDF file (as THREDDS / WMS are unable to read a classic NetCDF file for serving WMS requests)

I have a 4 dimensional Array of dimensions like this

Temperature[TIME][DEPTH][LAT][LONG], how can I write it using NetCdfCFWriter and using GridDataSets?


texts: Does IDW interpolation in ArcGIS Geostatistical analyst work as exact interpolation method?

texts:
As we know, Inverse distance weighted (IDW) is an exact interpolation method, which means the interpolated values of the known samples should be the same with its original observed values. But when I used the IDW in Geostatistical analyst (ArcGIS 10) to interpolate, I found that the interpolated minimum value and maximum value are not the same with the observed minimum value and maximum value.
Actually I did get a Predication Map using the interpolation method of Geostatistical analyst, its values seemed to be the same with the observed ones, but when I changed it into a raster layer, the values all changed, why?
In ArcGIS software, is the IDW interpolation in Geostatistical analyst the same as that in Spatial analyst? Are their equations the same with each other？

Forgive my poor English.


texts: Adding annotation programatically in QGIS 3.10

texts:
I'm able to add an annotation in QGIS 3.10 it works with this code

from qgis.core import QgsTextAnnotation, QgsPointXY, QgsCoordinateReferenceSystem, QgsAnnotationManager
from qgis.gui import QgsMapCanvasAnnotationItem
from PyQt5.QtGui import QTextDocument
from PyQt5.QtCore import QSizeF,QPoint

annot = QgsTextAnnotation()
annot.setFrameSize(QSizeF(100, 50))
annot.setMapLayer(layer)
annot.setFrameOffsetFromReferencePoint(QPoint(10, 10))
annot.setDocument(QTextDocument('TEST'))

# X and Y are defined previously
annot.setMapPositionCrs(QgsCoordinateReferenceSystem(layer.crs()))
annot.setMapPosition(QgsPointXY(x, y))

QgsMapCanvasAnnotationItem(annot, qgis.utils.iface.mapCanvas())

The annotation is added but I can't modify it in the GUI or delete it. How can I do it? There is an other thing to add in my code?


texts: Bounding Box to list of included Tile Quadkeys

texts:
In all places where it's possible to find formulas related to Tiles & Quadkeys, it seems that the function translating bounding box (rectangular on exact Tile borders - cannot divide a Tile, however not necessary square) into contained Quadkeys is missing.

Is there a closed formula for this?


texts: Shortcut to remove current layer in QGIS

texts: I wonder if there is a shortcut / hotkey to remove current layer in QGIS? And if no - is it possible to assign a custom one?


texts: Is there limit of 5,000 rasters stored in folder by ArcGIS for Desktop?

texts:
I have a python script that creates batches of ~6,000 cost distance rasters. These rasters are intermediate data that will be summarized in another process. Output is stored in plain old folders. However, at output raster 5000, my python script stops with the error msg something like: 'forcible close by a remote host.'

If I redirect the output for files 5000+ to a new folder, all works as before. I'm not sure if it's a folder size limit, since I can add shapefiles or other folders to the output folder (already containing 5000 rasters) without issue. The error is repeatable with different output and different folder locations.

Has anyone had this issue before and does anyone know where to change the folder limit?


texts: I'm using ArcGIS 9.3 to find centroid of some polygon features. I'd like to do this in Python, but so far I can only find solutions with VBA. Can anyone provide some leads here?


texts: Counting Water Crossings using ArcGIS Pro?

texts:
I've got a polyline network and i'm trying to calculate the number of water (polygon) crossings for each polyline.

I'd like to end up with a count associated with each unique polyline.

Is there an easy way to do this using ArcGIS Pro?


texts: Does transforming postGIS geographic data to EPSG:4326 , makes them geographic?

texts:
I am reading this about "Clustering on Indices" and in 25.2 says that

The geohash algorithm only works on data in geographic
(longitude/latitude) coordinates, so we need to transform the
geometries (to EPSG:4326, which is longitude/latitude) at the same
time as we hash them.
CREATE INDEX nyc_census_blocks_geohash ON nyc_census_blocks (ST_GeoHash(ST_Transform(geom,4326)));

I was under the impression that EPSG:4326 is available to both geometry and geography data. Also I thought that you can have geographic data by setting the column type as "geography" when creating a table or by geometry casting using the ::geography
Does the transformation to EPSG:4326 automatically turns the geometry to geography data?


texts: Missing output from ArcGIS Split Raster tool?

texts:
I have a number of raster datasets from different sources sitting in one GDB. In the same GDB I have a polygon feature layer I am using to split each raster layer using the Split Raster tool, the outputs of which are being placed in a folder in TIFF format.

For some reason, with some of these raster layers (and there seems to be no common theme) the process is running, giving me a completed message, GeoProcessing results indicate everything has worked correctly, but they are nowhere to be seen in the destination folder. For more than half of the raster layers, following exactly the same steps, the process has worked without fault. All rasters have the same projection & resolution.

I have tried exporting the problem datasets from the GDB to TIFF and then splitting, with the same result. There must be a property common to these datasets which is causing the error, but from comparison with those datasets that have worked correctly I can't identify what it might be.

Has anyone experienced a similar error?

Screenshots for reference:


texts: Get resolution and whether right/left looking from sentinel metadata

texts:
Is it possible in Sentinel hub products to know the resolution of the product and whether the SAR was right or left looking at the scene?

I know that it's possible to know if the SAR was operating in a descending or ascending mode from the /annotation/*.xml files (from the `...' value), but I couldn't find the aforementioned parameters inside those files!


texts: If a point (lat long) lies within any polygon of a shapefile

texts:
i checked different places to finally come up with a code which used ray casting method to check if the given lat, long was inside one of the many polygons of my shapefile. this code worked fine when i tested with a lat long inside the polygon but when when the point was not in any of polygon  it still printed TRUE

from shapely.geometry import Polygon, Point, MultiPolygon

import shapefile as shp

import matplotlib.pyplot as plt

polygon = shp.Reader("D:/New folder1/21_8_2016.shp")

polygon = polygon.shapes()

shpfilePoints = [ shape.points for shape in polygon ]

print shpfilePoints

polygons = shpfilePoints

for polygon in polygons:
    poly = Polygon(polygon)
    print poly

lat=15.139078
lon=120.763779

def point_in_poly(x,y,poly):

    n = len(poly)
    inside = True

    p1x,p1y = poly[0]
    for i in range(n+1):
        p2x,p2y = poly[i % n]
        if y > min(p1y,p2y):
            if y <= max(p1y,p2y):
                if x <= max(p1x,p2x):
                    if p1y != p2y:
                        xints = (y-p1y)*(p2x-p1x)/(p2y-p1y)+p1x
                    if p1x == p2x or x <= xints:
                        inside = not inside
        p1x,p1y = p2x,p2y

    return inside

print point_in_poly(lon,lat,polygon)

the lat lon given in the code dont lie in any of the polygon where the following lie inside of the polygon
lat: 15.094997
lon: 120.836563


texts: Can't close custom Web AppBuilder widget programmatically

texts:
I am creating a custom widget in ArcGIS Web AppBuilder.

This page shows how you can use the WidgetManager to close a widget programmatically, saving the user from having to click the X button manually.

The difference between the example code given on the above page and the code in my widget is I have a declare portion of my class like so:

define(['dojo/_base/declare', 'jimu/BaseWidget', 'jimu/WidgetManager', 'dojo/dom-style'],
function(declare, BaseWidget, WidgetManager, style) {
    return declare([BaseWidget], {
        ....
    });
});

Within an onClick function in the above, I attempt to call WidgetManager.getInstance.closeWidget(this.id), as recommended by the code sample on the ESRI page. This does not work, and the failure is silent.

I'm a Dojo/Web AppBuilder newbie, but I vaguely  understand that declare creates a new instance of a class, whereas the ESRI website seems to be using WidgetManager in a more abstract way. Does this have anything to do with it?


texts: Loaded SQL Point Data not displaying correctly in QGIS map

texts:
Having some issues getting my SQL data to display correctly in my QGIS map layout.  I am using Postgresql/Postgis, QGIS 2.4, Windows 7 (64bit OS).  I created a database with some point data that is supposed to outline the values collected at a specific geographical location in Colorado, but when I load the data, it is displayed as 1 single point in the middle of the Atlantic off the coast of Africa.  I will attach a screenshot of what I am seeing.

enter image description here


texts:
I have tried to create a fishnet(10 by 10) using the below code. But the output is as shown below. What is that I'm missing in the below code?
import geopandas as gpd
from shapely.geometry import Polygon
import numpy as np
ras= rasterio.open(r"My Raster file")

xmin,ymin,xmax,ymax = ras.bounds

length = 1000
wide = 1000

cols = list(range(int(np.floor(xmin)), int(np.ceil(xmax)), wide))
rows = list(range(int(np.floor(ymin)), int(np.ceil(ymax)), length))
rows.reverse()

polygons = []
for x in cols:
    for y in rows:
        polygons.append( Polygon([(x,y), (x+wide, y), (x+wide, y-length), (x, y-length)]) )

grid = gpd.GeoDataFrame({'geometry':polygons})
grid.to_file("grid.shp")
df= gpd.read_file('grid.shp')
df.plot()


texts: I am searching for a data source on past hurricanes that includes the coordinates of the path at different times including wind speed pressure and category. Any other added data would also be helpful.


texts:
I'm working with data containing multiple lat/long entries within a single row. I am using PostGreSQL to process the data. Is there a way to take the WKT lat/long field and convert it into a line geom field to then bring into ArcMap/QGIS?

Field example:


texts:
I have a table with geometry column called point, saved in SRID 4326. I want to get the geometry by its coordinates, but something didn't work. I suspect database configuration or something.

Here is what I do

select point from geom_data limit 1

And I get "0101000060E610000000000020558A5EC0000000C0F3B8424000B0CFCFA4717342"

This as text is

select AsText(point) from geom_data where point = '0101000060E610000000000020558A5EC0000000C0F3B8424000B0CFCFA4717342'

And the result is "POINT(-122.161445617676 37.444938659668)"

Trying this

select * from geom_data where ST_Contains(point, ST_GeomFromEWKT('SRID=4326;POINT(-122.161445617676 37.444938659668)'))

Didn't give a result. Tried this

select ST_PointFromText('POINT(-122.161445617676 37.444938659668)', 4326)

and get "0101000020E61000000F000020558A5EC0040000C0F3B84240" which is different.

Selecting within bounding box didn't worked either

select * from geom_data where
st_x(point) >= '-122.161445617670' AND st_x(point) <= '-122.161445617679'
AND st_y(point) >= '37.444938659660' AND st_y(point) <= '37.444938659669'

Any ideas?