texts: Dealing with overlapping grid (and graticule) labels using ArcGIS Desktop?

texts:
I do have 3 grids made by Grids menu in data frame properties. I do need them to display in certain layout and also to behave accordingly to couple of rules (overlapping issues), so I can automate my production workflow. I figured that its not possible in the way I am trying to do now which is as per this Ask A Cartographer post entitled changing text characteristics of graticule labels.

So what are my alternatives? Is there a way to access grid labels in the existing mxd document and modify them according to my needs?

Probably its not possible with ArcPy, but maybe someone did it via ArcObjects (preferably in Python) and would like to share some experience?

I do not have any of the Production Mapping extensions, so I guess ArcObjects are my only option?


texts: Column in dbf datatype string, should be float?

texts:
This code takes a .csv file, organizes the data and outputs another csv file that is used to create a .dbf file for ArcMap 10.4. The code appends the data as a float into the list 'output' and that list is written to the csv. However, the column that is 'begMi' in the .dbf, ends up being of datatype string, not float. Why is that happening?

Here is the code:

env.workspace = 'C:\Users\U2970\Desktop\SEC_FILE_EXPORTS'
os.chdir('C:\Users\U2970\Desktop\SEC_FILE_EXPORTS')

def headerBody():
    output = []
    outfile = csv.writer(open('program_NORTH_data.csv','wb'))
    with open('NORTH_V2017.SEC') as line:
        for data in line:
            col = data.split('!')
            corrRB = col[4]+col[32]
            begMi = round(float(col[12]),3)
            endMi = round(float(col[13]),3)
            dir_ = col[16]
            lane = int(col[31])
            set_ = int(col[41])
            output.append([corrRB,begMi,endMi,dir_,lane,set_])

    columns = 'corrRB','begMi','endMi','dir_','lane','set_'
    outfile.writerow(columns)
    for line in output:
        outfile.writerows([line])

headerBody()

n = 0
inTable = r'C:\Users\U2970\Desktop\SEC_FILE_EXPORTS\program_NORTH_data.csv'
outLocation = 'Somewhere deep in a sharedrive'
while True:
    try:
        outTable = 'secfileDBF'+str(n)
        arcpy.TableToTable_conversion(inTable, outLocation, outTable)
        break
    except:
        n = n + 1
        outTable = 'secfileDBF'+str(n)
        continue

print outTable,'Has been created'

Here is a snippit of the .dbf file before I cleaned it up:

 OID    corrRB       length endMi  dir_lane set_
   0    C000001E    383.967 385.72  I   0   0
   1    C000001E    0.0 3.795   I   1   0
   2    C000001E    3.795   16.928  I   1   0
   3    C000001E    16.928  32.259  I   1   0
   4    C000001E    32.259  35.738  I   1   0
   5    C000001E    35.738  41.946  I   1   0
   6    C000001E    41.946  81.332  I   1   0


texts: GPS data cleaning and sampling of track data?

texts:
I have gotten a fairly large  polyline-dataset of GPS-recorded tracks. I'd like to have a sample along the tracks for different variables. The problem is that the data has overlapping tracks. This type of situation occures in example when same track has been walked (recorded) back and forth:

Points in the picture are 3 m apart from one another (measured along the polyline). As you can see, there are duplicate observations, which would corrupt my sample.

What is a convenient way of moving points (=observations) that are too close to each other and, most importantly, how do I keep one observation and delete the other?

I think I should clean the polyline data first and then take the sample, but I'm still figuring out how to do this.


texts: Selecting points that intersect only once?

texts:
I'm working on a river layer that requires some refinements. I have attached an image for reference. I have a geofabric layer (light blue) that I have used to extract rivers with a flow >5ML/day (purple). The problem I have is that this geofabric layer is not consistent with the river layer. At some points the geofabric has touched the edges of unwanted tributaries. I have created points from the intersections (light green) and my idea is to select only the tributaries that intersect just once with the geofabric layer and cross-reference them with dangle points i have created (red cross) in order to remove them. I just can't seem to isolate the points at which they cross only once. Does anyone have any suggestions as to how I can achieve this?

River system (purple) is split at intersections.


texts: Why does connected polyline (KML file) appear unconnect in arcmap as layer?

texts:
I created polygone in google earth and imported it to arcmap as layer. All things are OK but when I converted polygon to polyline, there some parts of polyline has not shown although its vertices appear?
How to fix this problem ?


texts: ogr2ogr Conversion Clockwise to Counterclockwise?

texts:
Looking for the right combination of ogr2ogr parameters to get a .shp bundle to GeoJSON and reverse the Coordinate pairs. I recognize that ESRI defines Polygons as clockwise (http://www.esri.com/library/whitepapers/pdfs/shapefile.pdf) but GeoJSON RFC defines Polygons counterclockwise.
Current command is
ogr2ogr -progress -overwrite -f GeoJSON -s_srs "EPSG: 3857" -t_srs crs:84 file0.json file0.shp

which makes GeoJSON that is invalid when consumed by other tools that follow the RFC at https://tools.ietf.org/html/rfc7946.
Any ideas? If this needs to be run twice to convert to an intermediate format that's fine, if there is another (command line or automation) tool, I'm open to that also.


texts: Delete overlapping features (select & delete hidden feature below)

texts:
I have a whole map, with a lot overlapping between features. The first image showed part of the map, with 100% opacity, so the most top features are shown.

And in the second image showed the same part, but with 70% opacity, and I see many parts of "hidden" features below them. I already tried the Check Geometries to search duplicates and features within features, some of them deleted, some couldn't (it's said the features are obsolete)

I don't know what to do, but I think the I have to delete the features that hidden below, and I don't know how to select those features.


texts:
Why are the options in Join Data missing for me?

The middle part of the table is missing compared to the one here:

https://community.esri.com/thread/188841-join-data-from-another-layer-based-on-spatial-location


texts: Intersect overlapping entities and their attributes

texts:
I have a shapefile containing 6 classes of polygons on QGis 3.14. Some of these polygons overlap each others, but they're not geometry errors or duplicates, just polygons with different attributes. However I need to have maximum one entity at any place. So, I'd like to obtain, for each overlapping area, a new entity containing the attributes of the polygons that overlap. I can't make an intersection on a single layer, and the union tool gives me two polygons with equal geometries.
What should I do ?


texts: I would like to have map widget in .NET application. I need basic functions (zoom, move map) and possibility to display WMS. Thats all. Is there any library that you recommend for that? I am JavaScript developer and OpenLayers does all that things for me, but I need to recommend something for .NET dev.


texts: Geofencing routes from point A to point B using Google Maps?

texts:
I'm working on a task where given the coordinates of a line from point A to point B, I need to create a geofence around this line (road route).

But based on what I read here https://developers.google.com/maps/documentation/tracks/limits, geofence has a limit of 100 vertices per polygon.

Is my understanding correct, that its not possible for me to create the polygon for geofence given all the vertices that will be generated to enclose the route above will be more than 100?


texts: Tips for minimizing shapefile size without compromising geographical detail?

texts:
I work for a state government's fisheries service, and lately we've been digitizing a lot of our regulatory boundaries within the state's waters that we'd like to service out on web maps.  Unfortunately since a lot of these boundaries are in a bay and follow the shoreline, the polygons contain a lot of vertices which greatly inflates the file size of the shapefiles.

I know I can use the generalize tool to simplify the polygons, but as these polygons represent legal boundaries, we need to maintain their accuracy as much as possible.  Are there any other tricks I could use to reduce the file size without risking spatial accuracy?


texts: Geotools throws "Cannot encode multilinestring (yet)" while Shape import on Linux

texts:
I've encountered a strange platform-dependent SHP import problem.

Everything works fine on Windows, shapefiles are converted fine to database.
There is although small little warning:

  [[ACTIVE] ExecuteThread: '9' for queue: 'weblogic.kernel.Default (self-tuning)'] WARN geotools.jdbc  - No mapping for com.vividsolutions.jts.geom.MultiLineString

Obviously SHP contains Multilinestring. So why the same code on Windows works, but on Linux we've got exception:

  UnsupportedOperationException("Cannot encode MultiLineString (yet)")

Would appreciate if someone could help.


texts: Points and polygons difference issue QGIS

texts:
I'm trying to remove points that are in polygons with the basic difference tool (Vector-> Difference).
What I see after the operation suits me because all points in polygons are deleted, but as I open the attributes table, there are the same count of point in initial and processed layers. Is there any issue with this process ?
Is there any other process I could use? I tried SAGA difference but it's not available for points.

EDIT: There seems to be no incorect features topology


texts: How to define the number of decimals when adding a new field as double to an attribute table with python in Qgis 2.1?

texts:
How to change the expression in order to define the number of decimals?

layer.dataProvider().addAttributes([QgsField("fieldname",  QVariant.Double)])

As I read in QGIS API Documentation:

Parameters
        name Field name
        type Field variant type, currently supported: String / Int / Double
    typeName Field type (eg. char, varchar, text, int, serial, double).

         len Field length
        prec Field precision

So I try:

layer.dataProvider().addAttributes([QgsField("fieldname",  QVariant.Double, 10, 3)])

But I get this error:

TypeError: arguments did not match any overloaded call:

QgsField(QString name=QString(), Type type=QVariant.Invalid, QString
typeName=QString(), int len=0, int prec=0, QString comment=QString()):
argument 3 has unexpected type 'int'

QgsField(QgsField): argument 1 has unexpected type 'str'


texts:
I have the same problem stated here:

https://geonet.esri.com/thread/98493

When I set the visibility of a field to HIDDEN, the field will not become hidden.

l = arcpy.mapping.Layer("myLayer")
desc = arcpy.Describe(l)
finfo = desc.fieldInfo
index = 0
finfo.setVisible(index,"HIDDEN")


texts: Creating popups in Folium (Leaflet) for GeoJSON polygon layer

texts:
Looking through the Folium documentation I can't find how to add a popup to a polygon GeoJSON layer in Folium. Looking through the Leaflet documentation it seems like its possible, but I don't know how to translate the JavaScript to Python. Below is a snippet of the code I'm using to render the map.
stand_boundary = folium.GeoJson(stands_gdf, name='Stands')

The stands_gdf variable is a GeoJSON read into a GeoDataFrame.


texts: I'm trying to create a labeled heat-map of incidents. Currently, I can label every town in the state -- which clutters up my map -- or I can have no labels at all. My boss would like me to instead label only the towns where we have an incident recorded. I've been trying to puzzle out a solution using conditional labeling, but I'm not very familiar with programming. I'm wondering if this is a) possible at all and b) if possible, how would I tell QGIS to only display the labels of certain towns?


texts: Gdal.warp outputs black image when cutting from json

texts:
I have the following task: to cut a polygon out of an .tif file.
my task is to: create a new image from the polygons I have
I perform the following steps:

get image and polygon data
reproject polygon to EPSG:4326 for json standard
make json
get polygon's reprojected bounds
crop image via gdal.warp

geom.Transform(transform)
json_ap = geom.ExportToJson()

gx_min, gx_max, gy_min, gy_max = geom.GetEnvelope()

output_options = gdal.WarpOptions(cutlineDSName = json_ap,
                                    outputBounds= [gx_min,
                                                    gy_min,
                                                    gx_max,
                                                    gy_max],
                                    srcSRS = g_ref,
                                    dstSRS='EPSG:4326',
                                    cropToCutline=True)

gdal.Warp('teste.tif', raster_path, options=output_options)

All images I export are blank!
I tried to create the crop with no reprojection, with a different reprojection, with a reprojected tif as source...
What is going on?


texts: Dissolve Shapefile and sum attribute in R

texts:
I have a shapefile with thousands of polygons and a separate csv of additional attribute data. I have joined these by

new_shp <- merge(old_shp, add_data, by.x="ID", by.y="ID")

However, now I want to summarise some attribute fields that have been added by any aggregate function. For instance I want new_shp@data to go from:

small_area  large_area    age
    1           A          50
    2           A          40
    5           B          40

to:

  large_area     avg(age)
    A              45
    B              40

with the boundaries now defined by large_area.

Is this possible using R and if not then by any other tool?


texts:
ogrmerge.py -single -o merged3.shp *.shp -src_layer_field_content {DS_BASENAME}

let's me painless merge several shapefiles into one and adds a new field with the basename of the shapefile. Is it possible to only write the first 9 characters of the basename to the shapefile?

Is:

2938_9282_myfile

should

2938_9282


texts:
I'm trying to visualize a 11-bit TIFF image using matplotlib. To do it, I load the image using tifffile to obtain a 3D numpy array. Next, I have to down-sample the image because the array contains values beyond the range [0, 255]. Here is my code:

import tifffile as ti
from scipy.misc import bytescale

x_3band = ti.imread(three_band_file_path)
# transpose to make the shape [M, N, 3]
X = np.transpose(x_3band, (1,2,0))
X = bytescale(X)
plt.imshow(X)
plt.show()

The problem is that the visualized picture is different with the picture visualized in QGIS (see the picture, the left picture is with matplotlib).

How can I show the picture like in QGIS?


texts: How to automatically calculate maximum cell values contained within polygons?

texts:
I've got a dataset made up of a large range of polygons (about 70 or so) which represent structures. I also have raster layers which detail a flood depth across an area. The floods depth raster has a 10x10m cell resolution, so for most polygons (structures) there are multiple flood depth cells contained within each polygon. I'm trying to create a model which will identify the maximum flood depth value in each polygon. Now most of the cells aren't wholly contained within the polygon, but these still need to be included in the assessment, as even part of the structures will be flooded to that maximum depth.

I don't want to do this manually as there are hundreds of cells and I have over 20 scenarios to repeat the process for, so this would take days to do manually.

I feel like this should be pretty simple but i've run up short so far.


texts: v.generalize not working QGIS 2.14 "The following layers were not correctly generated" polygons

texts:
I am trying to generalize a polygon layer using the v.generalize algorithm in QGIS 2.14. While it works for linear features, polygon features are refused by the program. Might anyone know, what the specific problem in this case is? The Data is referenced in DHDN GK 4.

I don't want to do the workaround using linear features, but rather would like to directly generalize (chaiken) the polygon features.


texts: ImportError: cannot import name QgsVectorFieldSymbolLayer

texts:
I am using QGIS 1.7.4., downloaded April 2012, but receive an import error loading the class QgsVectorFieldSymbolLayer. As i understand this class has been added in the Python bindings on 2011-11-15. Does anybody know whether this class has already been included and if so, what else can be wrong?

No, while executing a Python script.
First line:

from qgis.core import (QgsSymbolLayerV2Registry, QgsDataSourceURI, QgsVectorLayer, QgsVectorFieldSymbolLayer, QgsSymbolLayerV2AbstractMetadata, QgsCoordinateReferenceSystem, QgsMapLayer, QgsField, QgsComposition, QgsComposerLegend, QgsComposerLabel, QgsRectangle, QgsFeature, QgsMapLayerRegistry, QgsGraduatedSymbolRendererV2, QgsSymbolV2, QgsRendererRangeV2, QgsMarkerSymbolLayerV2)

All classes are imported okay, except QgsVectorFieldSymbolLayer.


texts:
Using Carto, I am trying to create shapefiles of regions using points(Link to Dataset in Carto ) with a shared column value to group by.   I've included 3 photographs. 1st photograph of the points outlining the shape. the second photo is my desired result which I have drawn out by hand. The third photo is an attempt that doesnt meet my desired result..

I am using the analysis tool provided by Carto, running a concave hull and have tried various target percent values and without the desired results. For this example I've included a target percent of 80 and the result is shown in a photograph at the bottom. The code is provided below

  SELECT row_number() over() as cartodb_id,
         territoryid as category,
         ST_ConcaveHull(
             ST_Collect(the_geom),
             0.8,
             false
          ) AS the_geom,
         count(1) as count_vals
    FROM (
         SELECT *
           FROM ( SELECT * FROM "st-braund".bd_regions_1 ) _camshaft_category_filter
          WHERE territoryid IN (3)
         ) _analysis_source
GROUP BY territoryid


texts: Error When Trying to Run Slope Analysis Script

texts:
I have written a script which is supposed to loop through all the rasters in a folder and perform a slope analysis on each raster and save the results in another folder.

I keep getting the following error:

  < type 'exceptions.RuntimeError'>: ERROR 000875: Output raster:
  C:\Slope_Outputs\burd_25m_dem_slope's workspace is an invalid output
  workspace.

Here's my code:

class LicenseError(Exception):
    pass

# Set desktop license used to ArcView
#
import arcview
import arcpy
from arcpy import env

try:
    if arcpy.CheckExtension("Spatial") == "Available":
        arcpy.CheckOutExtension("Spatial")
    else:
        # Raise a custom exception
        raise LicenseError

except LicenseError:
    arcpy.AddMessage("Spatial Analyst license is unavailable")
except:
    print arcpy.GetMessages(2)

from arcpy.sa import *

originLocation = arcpy.GetParameterAsText(0)
slopeMeasurement = arcpy.GetParameterAsText(1)
destinationLocation = arcpy.GetParameterAsText(2)

arcpy.env.workspace = originLocation

rasterList = arcpy.ListRasters("*")

for raster in rasterList:
    finalDestination = destinationLocation+"\\"+raster+"_slope"
    arcpy.AddMessage(finalDestination)
    outSlope=arcpy.sa.Slope(raster, slopeMeasurement)
    arcpy.AddMessage(outSlope)
    outSlope.save(finalDestination)

I'm not really sure why I'm getting this error.

Any ideas?


texts: K means clustering where centroid is a node

texts:
How would you approach a k mean clustering where the calculated centroids have to be on the same position as one of the nodes?

as you can see we have 2 centroids which are arguably the center of each cluster while being "ontop" of a node.

Is there already a name for this problem? I couldn't find anything related to it.


texts: How to extract surface from air photo avoiding trees and buildings?

texts:
I am working with Drone air images. These are geotiff images and clearly portray trees and buildings. I need to first classify the image according to land cover distribution, then need to remove the trees, grasses and buildings to get the original terrain model. I have another image of same area with elevation values and I need to reduce the elevation of objects over the surface.

I used ArcGIS 10.2 Spatial Analyst>Image Classification. I selected the Training Areas against my classification requirement and the generate signature file using automated toolset. I used this signature file for supervised classification and used Maximum likelihood for my classification.

The result of classification was awful. I found the big trees are classified as both grass land and trees as well as crops field also classified as big trees.

The image is RGB and it shows similar reflectance of big trees, crop fields and grass land. So, it is difficult for me to classify different classes separately to extract actual surface reducing tree height or crop height.

I am wondering if anybody has any suggestions regarding this issue.


texts:
I want to create and elevation hillshade but I found the script for QGIS 2.18. What should I use instead of processing.runandload? Here is my script on QGIS 2.99 python console:

import processing
rasterLyr = QgsRasterLayer("C:/Users/Mustafa Uçar/Desktop/Tutorial/qgis_data/dem/dem.asc", "hillshade")
rasterLyr.isValid()
processing.runandload("gdalogr:hillshade", rasterLyr, 1, False, False,1.0, 1.0, 315.0, 45.0, "C:/Users/Mustafa Uçar/Desktop/Tutorial/qgis_data/dem/hillshade.tif")
#to verify the output image
processing.runandload()


texts: Script Works in Console but not in QGIS Script Editor

texts:
I've recently started working with Python in QGIS and I've been trying to make a script that dissolves a polygon layer and saves the output in the same folder as the input, something I believe isn't possible to do with QGIS Graphical Modeller.

The code I've written is as follows:

processing.run("native:dissolve", {'INPUT':iface.activeLayer(),'FIELD':[],'OUTPUT':iface.activeLayer().dataProvider().dataSourceUri().split('2D', 2)[0] +"2D Simulation Elements dissolved.shp"})
iface.addVectorLayer((iface.activeLayer().dataProvider().dataSourceUri().split('2D', 2)[0] +"2D Simulation Elements dissolved.shp"),"","ogr")

It works fine in the console, however when I put it in the Script Toolbox it fails with the following error:

The code is confusing me cause, I obviously don't have 224 lines of code for it have an error there...
I've looked online and followed the suggestion to add the following on top of my code:

QgsApplication.initQgis()
app = QgsApplication([], True)

But it didn't make my script work as well.


texts: How to migrate a GeoServer setup to Tomcat?

texts:
I have some problems with Tomcat/GeoServer. I was able to create project in stand-alone GeoServer and all works fine.

This is my code.

<%@ page language="java" contentType="text/html; charset=UTF-8"
 pageEncoding="UTF-8"%>
  <!DOCTYPE html PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN"
   "http://www.w3.org/TR/html4/loose.dtd">

<meta http-equiv="Content-Type" content="text/html; charset=UTF-8">
<title> My first JSP </title>
<link href="http://localhost:8080/geoserver/openlayers/theme/default/style.css"type="text/css"/>

But when I try to put this in Tomcat nothing happens and I don't see a map. What do I need to do with Tomcat to run project which uses GeoServer?

I know my English is terrible.


texts: Combine multiple points to polygons by property

texts:
I'm working on a website with a Cesium map, and I'm using turf.js for geo calculations.
I need to find a way to merge a large set of points on the map (with a fixed distance between them) to polygons, by a certain property.
For example, here is a set of points, with two different properties (marked red and orange):

And the desired output should be something like this:

The way I tried to approach this is by transforming each point into a square, and then use turf's union on all squares that have the same property. This produced the second image above.
The problem with this solution is that I'm getting a large number of points (can be around 100K) and performing the union on them can sometimes take a very long time (up to 2 minutes), and I need to find a more performance efficient way of doing this.


texts: Find nearest line feature from point in QGIS

texts:
I have a line feature set and point feature set. I want to get the nearest line feature from every point. For this, the line features are converted to point using "Extract Nodes" tool. Then I use the "Distance to Nearest Hub (line to hub)" tool.
Here, many points got the nearest line but for some lines that don't have many nodes are not recognized by the nearest point. Here I marked one such case.
How to solve this?

Red circle marked showing the point which falsely recognizes the nearest line due to less node of nearest line.


texts:
I am a beginner at mapping and I am wondering how best to structure a dataset I am working on. Ideally, I would like one data point for each city on my map, but I would also like to include alternative names for these cities. I would also like these alternative names to be searchable on the map.

Any advice on structuring my data?

I will either be using Carto or QGIS for this project.


texts: I have a point layer and a polyline layer. I want to select all points in the point layer that are not snapped to a vertex in the polyline layer. How can I accomplish this? I'm using MapInfo and hoping that some SQL select statement can do the trick, but I'm happy with solutions using other GIS packages/program languages as well.


texts: Syntax error during calculate field between two long integer fields

texts:
I keep getting an error when trying to run a calculate field section of my code. I'm trying to simply add values from two long integer fields (Speed and Lanes) into a newly added field (Stress_Calc). I know that it has to be with the syntax of my expression. I have tried several different formatting combinations to try and get the code to work.

Below I will list the different combinations of formats of the expressions I have tried and the results that come out of them :

"!Speed! + !Lanes!" = ERORR 999999
'!Speed! + !Lanes!' = ERORR 999999
"'!Speed!' + '!Lanes!'" = ERORR 999999
'"!Speed!" + "!Lanes!"' = No error but the calculation does not occur

Code listed:

arcpy.AddField_management(MergeNetwork, "Stress_Calc", "LONG")
StressCalcExp = '"!Speed!" + "!Lanes!"'
StressCalc = arcpy.CalculateField_management(MergeNetwork, "Stress_Calc", StressCalcExp)


texts: Assertion Error with arcpy.mapping.UpdateLayer

texts:
I"m working with arcpy, where I'd like to write a script that examines the selection and if only 1 element is selected creates a buffer. I want to get the symbology from another layer, as the buffer needs to be hollow, becouse at the further part of the script I export the layout for printing. But I get the following error message:

Runtime error
Traceback (most recent call last):
  File "<string>", line 9, in <module>
  File "c:\program files (x86)\arcgis\desktop10.2\arcpy\arcpy\utils.py", line 181, in fn_
    return fn(*args, **kw)
  File "c:\program files (x86)\arcgis\desktop10.2\arcpy\arcpy\mapping.py", line 1869, in UpdateLayer
    assert isinstance(update_layer, Layer)
AssertionError

My original code looks like:

if "something.mxd" in path or "otherthing.mxd" in path:
     if len(arcpy.Describe("mylayer").FIDSet.split('; ')) == 1:
         mxd =arcpy.mapping.MapDocument("CURRENT")
         df = arcpy.mapping.ListDataFrames(mxd)[0]
         arcpy.Buffer_analysis("Kutak", "Buffers.shp", "1000 Meters")
         updateLyr = arcpy.mapping.ListLayers(mxd, "Buffers.shp")
         sourceLayer = arcpy.mapping.Layer(r"style_ref.lyr")
         arcpy.mapping.UpdateLayer(df, updateLyr, sourceLayer, True)

Any ideas?


texts: OSM administrative polygons : what is the current worldwide coverage?

texts:
Does OSM have a summary view of its administrative polygons coverage worldwide for each administrative level ?

NaturalEarth have Countries (L0), and subdivisions (L-1) with complete (100%) coverage.

GADM.org have a very clean summary of its coverage for L-1 to L-5, it's super easily readable thanks to a nice visualization (see previous link and following snapshot) :

Does OSM have a similar coverage summary somewhere ?


texts:
I've got an FME Workbench (version 2015.0) which searches for shapefiles in a folder tree.  The shapefiles are contours which are split into tiles for a particular area.  The workbench is fairly simple, it filters based on certain attributes then combines the output into a merged shapefile containing all the contours.

However, the workbench fails because the resulting shapefile is too big.  I get the following error:

  Error - Cannot write out a file larger than 4.295GB. Split your data
  into two or more shape files

Of course, I didn't necessarily know beforehand that my output would be so large.  Is there any way for FME to detect that the file will be too big and somehow automatically write the output to two or more shapefiles?  Perhaps with the suffix "_1.shp", "_2.shp"?


texts:
I have a layer file that has a few polygons in it, each of these polygons is required to be a different color based a column in the DBF file.

I normally color layers like this:

ISimpleRenderer simplerenderer = geolayer.Renderer as ISimpleRenderer;
SimpleFillSymbol simpleFillSymbol = simplerenderer.Symbol as SimpleFillSymbol;
simpleFillSymbol.Color = selectedColor;

But this doesn't have a mechanism to use multiple colors on the same layer.  How does one accomplish this?

ArcEngine 10, C#, VS2010


texts: Why would "import arcpy" hang Python script run from DOS/cmd prompt?

texts:
I have a very small test script (test.py) that I use to see if Python is installed and ready to use ArcPy. All it contains is:

print "Starting"
import os
print "OS imported"
import arcpy
print "ArcPy imported"

I’m running it under Vista SP1 with ArcGIS Desktop 10 SP1 and the only version of Python present is 2.6.5 which was installed by ArcGIS Desktop 10. ArcMap 10 SP1 is starting fine and the "import arcpy" and other arcpy functions work fine from the Python Window.

However, when I run “python test.py” from the DOS/cmd prompt, the first two print statements execute as expected, but the last one does not, so it seems to be failing/hanging on “import arcpy”. For a brief time on two separate days the script actually ran fine, but then it seems to stop working (which to me is bizarre). To try and figure what the problem is I have been reading Importing ArcPy and I think I have the system variables Path and PYTHONPATH set to appropriate values:

Path=C:\Program Files\weoapp;C:\Program Files\Common Files\Microsoft Shared\Windows Live;C:\Windows\system32;C:\Windows;C:\Windows\System32\Wbem;C:\Program Files\Intel\WiFi\bin\;C:\Program Files\ATI Technologies\ATI.ACE\Core-Static;c:\Program Files\Common Files\Lenovo;C:\Program Files\Common Files\Roxio Shared\10.0\DLLShared\;C:\Program Files\Common Files\Roxio Shared\DLLShared\;C:\Program Files\Common Files\Roxio Shared\DLLShared\;C:\Program Files\Common Files\Roxio Shared\10.0\DLLShared\;C:\Program Files\ThinkPad\ConnectUtilities\;C:\Program Files\Lenovo\Client Security Solution;C:\Windows\System32\WindowsPowerShell\v1. 0\;C:\Program Files\QuickTime\QTSystem;C:\Program Files\Windows Live\Shared;C:\Program Files\QuickTime\QTSystem\;C:\Python26\ArcGIS10.0
PYTHONPATH=C:\Python26\ArcGIS10.0

but I’ve also tried:

PYTHONPATH=C:\Python26

When I first noticed the problem I had both Python 2.5.x (C:\Python25) and 2.6.5 (C:\Python26) present, so I went to considerable trouble to remove both, and uninstall ArcGIS Desktop 10 SP1 before reinstalling ArcGIS Desktop 10 (which installed Python 2.6.5) and then its SP1 again.

Does anyone have an idea as to what may be causing “import arcpy” to hang for me most of the time, or have any ideas as to what debugging steps may be available to me?


texts: Does building a raster pyramid have benefits even if rasters only displayed zoomed in?

texts:
Standard advice about handling large raster datasets on QGIS (and elsewhere?) usually suggests building a raster pyramid. If I understand correctly this process creates reduced resolution (i.e. resampled rasters with smaller dimensions) copies of the raster images so that when zoomed out the GIS software can display these instead of the full resolution rasters.

When using professionally produced tiled raster mapping datasets, for example rasters from the UK Ordnance Survey for use as basemapping, it is normal to set up a QGIS project to use multiple different raster layers depending on zoom level. The consequence is that when zoomed out a different raster layer is used altogether.

When working in this way is there any gain from producing a raster pyramid, or will this actually be complicating things?

(I'm investigating how to improve performance in the situation described - where we use several VRT files with QGIS, each accessing a tiled raster dataset, some with a large number of tiles. I'm also be struggling with one dataset with so many tiles referenced by the VRT that my PC/QGIS can't cope, even when zoomed in to show only one tile.)


texts: Using ModelBuilder to tabulate area of overlapping buffer zones in ArcMap

texts:
I'm looking into how much open and close canopy falls within given buffer zones. My buffer zones are a polygon layer created from point sites, but canopy data is a raster with 2 types of data, open or closed canopy. It was pulled from PASDA and has a fairly fine resolution, which is why some other solutions (Converting to polygon) did not work.

This is my first time using ModelBuilder, but I'm attempting to use the tabulate area tool in ModelBuilder by iterating through each row of my buffer zones, tabulating area, and then collecting the new tables. Finally, I use a second model (With the iterating model as a sub-model) to merge these tables together into one output (Added images of models).

However, I'm getting error 001156 (Attached screen shot of error) at the end.

Does anyone have a suggestion to get around this?

It suggests that the string length of columns is limited (This doesn't seem to be the case), or that there is data of one type being written to a column that contains a separate data type. All of the columns from Tabulate area are the same, however.

Unfortunately, I cannot include a screenshot of my tables as they contain locality information for sensitive species. The tables from Tabulate Area have 4 columns (Object ID, Name, Closed Canopy, Open Canopy)

https://imgur.com/a/lpJVe8x


texts:
I'm not in GIS, and I'm working on a nutty redball express project that involves matching some zip+4 data to congressional districts for the upcoming 113th congress. We have found very few sources of this data as it hasn't been officially published by the census bureau. One of the ones we did find was ZipInfo.com. The problem is I'm not sure if we can trust their data -- I'm not horribly familiar with the field, and my google fu leads to no smoking guns either for or against.

So has anyone had experience with their data? Is it worth the price of admission?


texts: Is there a lower bound on the std deviation for v.kernel in GRASS?

texts:
I'm using the v.kernel function in GRASS to create heat maps. I am able to generate an acceptable heat map using a std deviation of .2 and a several thousand rows and columns. After the map is generated, I look at it in QGIS.

v.kernel input=data output=heatmap stddeviation=.2

However, if I go to a std deviation below .2, all I get is a black rectangle. A std dev of .19 or .1 or .01 results in no variation being shown on the heat map, just all black. Does anyone know why this might be? I want to be able to increase the granularity and precision of the maps I'm creating, and lowering the stddeviation would seem to be the solution.


texts: Parsing CSV data and updating a dictionary and creating a polyline

texts:
I am a little stuck on my code.  I am new to Python and am learning as I go, but have usually been able to figure it out till now.  The code was working up till I put the dictionary part in and now this is where its failing.  When I run the program, it only goes through once and doesn't appear to be storing the data.  I have found people in the discussions have had the same problem, but they never said what helped them resolve the issue.  Below is a breakdown of the objective of the program:

The rhinos in the spreadsheet appear in no guaranteed order, and not all the rhinos appear at the beginning of the spreadsheet. As I parse each line, I must determine which rhino the reading belongs to and update that rhino's polyline track accordingly. I am not allowed to sort the Rhino column in Excel before I export to the CSV file. My script must be "smart" enough to work with an unsorted spreadsheet in the order that the records appear.

I do not immediately know how many rhinos are in the file or even what their names are. Although I could visually comb the spreadsheet for this information and hard-code each rhino's name, your script is required to handle all the rhino names programmatically. The idea is that I should be able to run this script on a different file, possibly containing more rhinos, without having to make many manual adjustments.

Sample Data:

Observer,X,Y,Rhino,Comments
Ben,26.99391,-19.10447,Bo,
Ben,27.00071,-19.1089,Tulip,Bathing
Ben,26.9919,-19.10511,Bo,
Ben,27.00071,-19.1059,Tulip,
Ben,26.96809,-19.09578,Patches,
Ben,26.97808,-19.11016,Dinky,
Ben,26.99213,-19.10395,Bo,
Ben,27.00083,-19.10326,Tulip,

I am lost and confused....

#Import a csv file into ArcMap and plot the points for each Rhino
import arcpy
import csv

arcpy.env.workspace = "C:\\GEOG485\\Lesson4\\Project4"
arcpy.env.overwriteOutput = True

targetFolder = "C:\\GEOG485\\Lesson4\\Project4"
rhinoShape = "C:\\GEOG485\\Lesson4\\Project4\\Rhino.shp"
rhinoFile = "C:\\GEOG485\\Lesson4\\Project4\\RhinoObservations.csv"
spatialRef = arcpy.Describe(rhinoShape).spatialReference

try:
    #Open CSV file and read the headers

    trackPoints = open(rhinoFile, "r")
    textLine = trackPoints.readline()
    headerLine = textLine.strip("\n")
    pairList = headerLine.split(",")
    print pairList
    polylineArray = arcpy.Array()
    #Get index for header fields
    xIndex = pairList.index("X")
    yIndex = pairList.index("Y")
    rhinoIndex = pairList.index("Rhino")
    print rhinoIndex

    #Created an empty dictionary
    dictionaryRhinoTracks = {}

    for line in trackPoints.readlines():
        line = line.rstrip("\n")
        splitLine = line.split(",")
        print splitLine
        #Created variable for needed field items based on index location
        pointX = splitLine[xIndex]
        pointY = splitLine[yIndex]
        rhinoName = splitLine[rhinoIndex]

        #Create a point object

        currentPoint.x = pointX
        currentPoint.y = pointY
        currentPoint = arcpy.Point(currentPoint.x, currentPoint.y)
        print currentPoint

        #Check if the dictionary does not contain name
        if rhinoName not in dictionaryRhinoTracks:
            rhinoTrackArray = polylineArray.Add(currentPoint)
            #Add the point object to the Array
            dictionaryRhinoTracks[rhinoName].Add(currentPoint)

    cursor = arcpy.da.InsertCursor(rhinoShape, ['SHAPE@'])
    cursor.insertRow([polyline])

    for rhinoName in dictionaryRhinoTracks:
        row.Shape = dictionaryRhinoTracks[rhinoName]
        row.Rhino_Name = rhinoName
        cursor.InsertRow(row)

except:
    print ("error again")


texts: Group by variable and create layers from that in QGIS

texts:
I have a large shapefile with one layer with several variables. The shapefile has US block groups that I'll use for analysis. However, as I want to use this for leaflet, I would like to create different layers/files for each state.

State Value1 Value2 Value3 ...
  NY    x      x      x
  CO    x      x      x
  NY    x      x      x

Is there a way of grouping per state and generate a layer from that, that I can later save as a shapefile?

I know that I can do it one by one, by using the SQL query tool, but as I have 50 something states, plus different files, that would take me quite of a long time.


texts: looping through a folder with QGIS script

texts:
I'm attempting to write my first Python script for QGIS. I have no programming experience and am having some difficulty.

The aim is to loop through a folder and check to see if the file name contains the string "rsa2". If it does, copy the file path to an existing attribute called "rsa2_document" for each feature in a shapefile.

I've been using this as an example of how to change attributes with a script but I'm getting the following error: name 'folder_containing_files' is not defined See log for more details. Since I'm completely new to this, I can't tell if the error message is not being helpful and it's a problem with my code or if there is some sort of issue with using a folder type input. I looked here for an idea of how to use the folder type input, but doesn't say if the folder type input gives a string or an object.

My code so far is below. I'm using QGIS 2.8.2 on Windows 7.

##inlayer=vector
##folder_containing_files = folder

from qgis.core import *
from PyQt4.QtCore import *
import os

for file in os.walk(folder_containing_files):
    if string(file).contains("rsa 2"):
        path = file

layer = processing.getObject(inLayer)
for feat in layer.getFeatures():
    outlayer.changeAttributeValue(feat.id(), rsa_2_document, path)

outlayer.commitChanges()


texts: Determining the number of vertices of a selected object in QGIS 3

texts: How can I determine how many vertices/nodes a selected feature has in QGIS 3?


texts: How to convert a GeoTIFF from WGS84 to EPSG:3857 with gdalwarp?

texts:
I'm trying to convert images with WGS84 projection to EPSG:3857 but warp is not giving me a output I would expect.
This is my command:

gdalwarp -s_srs "+proj=merc +datum=WGS84" -t_srs "EPSG:3857" SRTM_1km.tif SRTM_1km_warped.tif

If I overlay the warped image on top they looks exactly the same.

This is the gdalinfo of my source image:

Driver: GTiff/GeoTIFF
Files: SRTM_1km.tif
Size is 43202, 14402
Coordinate System is:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4326"]]
Origin = (-180.008333333767950,60.008333328098033)
Pixel Size = (0.008333333333000,-0.008333333333000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (-180.0083333,  60.0083333) (180d 0'30.00"W, 60d 0'30.00"N)
Lower Left  (-180.0083333, -60.0083333) (180d 0'30.00"W, 60d 0'30.00"S)
Upper Right ( 180.0083333,  60.0083333) (180d 0'30.00"E, 60d 0'30.00"N)
Lower Right ( 180.0083333, -60.0083333) (180d 0'30.00"E, 60d 0'30.00"S)
Center      (  -0.0000000,  -0.0000000) (  0d 0' 0.00"W,  0d 0' 0.00"S)
Band 1 Block=43202x1 Type=Int16, ColorInterp=Gray
  NoData Value=-32768
  Metadata:
    LAYER_TYPE=athematic


texts: I am trying to make my scanned and georeferenced TIFF file translucent so I can see the layers underneath, however when I do so using the global opacity tool, it becomes pixelated, as shown in the image. How do I fix this?


texts:
I've recently moved to a new city (and country) and I'm really eager to explore it. I want to gamify the process, so I thought about a scratch map for city's streets, districts and even country's regions but actually didn't find any good free apps. What I thought maybe is a self-hosted OSM where I could edit data myself? Thus I have a few questions:

Can I host an OSM (database + some UI) for a specific country only? How much space will it need for a county like Germany?
Is it possible to somehow "scratch" places in OSM and display a web page with that?
Any good guides for complete noobs on these topics?


texts: Executing QGIS tools in Python from Linux Shell

texts:
I would like to execute a Python script which uses QGIS Processing tools.
On Ubuntu 20.x I installed qgis, qgis-common, python3-qgis and python3-qgis-common.
My objective is to execute the tool Exploe HStore Field from shell within a python script.
The following lines are working well within the QGIS-GUI Python Console
from qgis.core import processing
from qgis.core import *
processing.run("qgis:explodehstorefield", {
    'INPUT' : '/home/myacc/geodata.gpkg|layername=points',
    'FIELD' : 'other_tags',
    'EXPECTED_FIELDS' : '',
    'OUTPUT' : 'ogr:dbname=\'/home/myacc/geodata_explode.gpkg\' table=\"points\" (geom)' })

My problem is, that processing.run is not recognized outsite the QGIS-GUI Python console.
I tried to create a pre-build script from the modeller, but there is the same problem.


texts: Get Intersecting roads using ArcGIS Server with Silverlight API?

texts:
I am developing a tool where when I click on a road in the map I have to list the roads intersecting the selected road
using ArcGIS Server10 with silverlight API

to develop this tool
I got the roadID of the selected road using IdentifyQuery
Then I got the envelope of that road using QueryTask
I passed that road envelope to IdentifyTask again and queried the intersecting roads

but its lists me all the roads inside the envelope but I need only intersecting roads
also its very very slow :(
I am stuck here.


texts: How to compute NDSI in ArcGIS for Desktop?

texts:
How do I compute a Normalized Difference Snow Index (NDSI) in ArcGIS's Raster Calculator? I have a Landsat 8 scene (hdf4 format) that contains derived surface reflectance band-by-band. Each band is a grid 16-bit signed format.

In Raster Calculator, I tried to compute NDSI using the following code: (float(b3)-float(b6))/(float(b3)+float(b6)), where b3 is Band 3 (Green) and b6 is Band 6 (SWIR). This results in an output of 32-bit and floating point pixel type. However, the resulting NDSI is not scaled from -1 to 1. Instead, the values range from, in this case, -367 to 251. The original surface reflectance for each band ranges from -1000 to 12018.

Does this have to do with rescaling? If I rescale the NDSI output (dividing the output results by 1000), the values still seem low for an ice and snow-covered surface (now, the output values are -0.367 to 0.251).


texts:
Shapefile A was joined with a new variable x from table C. Then I saved Shapefile A under a new name Shapefile B. The joined field (variable x) appears in Shapefile B as a "native" variable (which was my aim).

I tried to edit Shapefile B. Editing is possible, however, I cannot edit variable x in Shapefile B. I do not really understand why and how I can make the respective variable editable again.

All of this takes place in QGIS.


texts:
I would like to keep my legend on as I export a multipage map set to PDF, but the legend appears to lose any layer items and just keep the title ("Legend"), which is what appears in my PDF document.

My map set in Composer has 12 pages, 8 of which contain similar content: a map window with different views of the same extent for each page, a legend that is set to the correct map window on each page, a scale bar for each map window, and title and descriptive text. Each legend has a sequential ID value from 1 to 8 ("Legend 1" ... "Legend 8") and is filtered to viewable layers.

I am running QGIS 2.8 on Windows 7. All my layers are shapefiles.


texts: Closing AddIn in ArcGIS with Visual Studio 2013?

texts:
I was working in ArcGIS 10.2 with VBA but now I work with Visual Studio 2013 to create addins.

I need to close all processes, I did this with the end statement in VBA but I don't find a command to do the same process in Visual Studio.


texts: Inconsistent results using Zonal Statistics by Table for Min in ArcGIS Desktop?

texts:
Using ArcGIS Desktop 10.5 I am looking to find the minimum elevation in x number of polygons using a DEM.
Using Zonal Statistics method:

import arcpy
from arcpy import env
from arcpy.sa import *
env.workspace = "Path/Default.gdb"
outZSaT = ZonalStatisticsAsTable("Tests", "FID", "dem.img", "output_table","DATA","MIN")

In another approach, I convert the polygon to a raster, convert the raster to points, and then extract values to the points from the DEM. When I do this I see elevation values that are lower than the absolute minimum that Zonal Statistics as Table returns. (RASTER VALUE is the elevation). These are the pt values in my polygon with FID sorted by RASTERVALU in ascending order. Note the min values found by Zonal Statistics are somewhat far down this list.

I would expect the MIN for FID 0 here to be 5.380977

What may account for this discrepancy and how can I fix it?


texts: Programmatically fetching current Add-In version using ArcGIS Pro SDK

texts:
The version of a custom ArcGIS Pro SDK Add-In is held in the config.daml of the Add-In project:

<AddInInfo id="{some guid}" version="1.0" desktopVersion="1.1.3068">

How do I programmatically access the current version of the Add-In  (without reading the xml file directly at runtime using an xmlreader)?

In the Add-In framework for ArcMap using ArcObjects SDK, the Add-In version could be accessed like so:

My.ThisAddIn.Version

Has anyone discovered the ArcGIS Pro SDK way of doing this? I would imagine it's possible but can't find the method or property to retrieve it.

I guess the version number is also held in AssemblyInfo.cs, but i'd rather use the value from config.daml.


texts: ArcGIS Geodatabase raster mosaic "could not build dataset item" error 80042019

texts:
For ArcGIS 10.3.1 with new 10.3.1 file geodatabase raster mosaics we randomly produce an 80042019 error when building overviews (which fail to build). In some circumstances the participating rasters in the mosaics are removed after the build overview fails.  We approached Esri Inc technical support in USA but they were unable to replicate the issues.

Trolling technical forums has furnished me with a very small check list. All of which we have accounted for in testing, including:

Path file names length do not exceed 256 characters
File names do not break windows naming conventions
Raster's are visible in the mosaic and exist (double checked via repair)
There's no permission issues

The issue has been observed on a raster mosaics with only two rasters (one built, the other did not), and on raster mosaics with many rasters (15 odd raster's 36 did not build).

The rasters are added to mosaics on an ad hoc basis by different users. I suspect this the underlying reason, and it's probably to do with the mosaics properties conflicting with loaded rasters. My next tests will concentrate on rebuilding raster statistics and then mosaic histogram before the build overviews. But I'm flying blind. Has anyone see this issue before?


texts: What is the best web editor and server for PostGIS?

texts:
I'm looking for a web editor (similar to OpenStreetMap) that works with a map server (PostGIS). This should permit editing and creation of new structures on a map.

What is the best option?


texts: How to return plpgsql function as table?

texts:
I have three input tables and would like to return the result in a new table. Below is the code

CREATE OR REPLACE FUNCTION vrg.VRGTest(_tbl1_VRG regclass,_tbl2_urban regclass,_tbl3_postcodes regclass)
RETURNS table (
    PostcodeName text,
    vrgid double precision,
    urban_gid integer,
    vrg_area double precision,
    urban_area double precision,
    weights double precision
    )AS
      $BODY$
      BEGIN
      EXECUTE format
         ('SELECT
             postcodes.name as PostcodeName,
             fishnet.fid as vrgid,
             urban.gid as urban_gid,
             st_area(fishnet.geom) as vrg_area,
             ST_AREA(ST_INTERSECTION(urban.geom, fishnet.geom)) as urban_area,
             ST_AREA(ST_INTERSECTION(urban.geom, fishnet.geom))/st_area(fishnet.geom) as weights
         FROM '   || _tbl1_VRG       ||  'AS fishnet,'
              || _tbl2_urban     || 'as urban,'
              || _tbl3_postcodes || 'as postcodes
         WHERE
             ST_INTERSECTS(urban.geom, fishnet.geom) and ST_INTERSECTS(fishnet.geom, postcodes.geom)')

       END;
       $BODY$ LANGUAGE PLPGSQL;

The query executes successfully but when I call the function

select vrg.VRGTest( 'bb_sumatrapart_grid_1000m','ghs_pop_indonesia_v2_reclass_poly_union_v2_projected','gfk_2016_id_5_digit_pcd_areas2013_projected')

I get the error "relation does not exist" I even tried giving the schema name as prefix but no luck.


texts: Defining variable for output location based on feature class attribute?

texts:
I'm pretty new to python. I have a feature class of points (Points_Layer) that covers the state and a feature class of municipal boundaries (Municipalities).

My script is going to loop through each municipality in the Municipalities layer. First the script will use the Select Layer by Location tool (overlap_type = intersect) to select points from "Points_Layer" within the municipality then Feature Class to Feature Class to create a layer of the selected features. It will loop through every municipality in the state. I want to save the output to different geodatabases/feature datasets based on attributes from the Municipalities layer (county name and municipal code).

My files are set up like this:

C:\Counties\
    CountyA.gdb
        Town1_0001
        Town2_0002
    CountyB.gdb
        Town3_0101
        Town4_0102
        Town5_0103

The Municipalities layer has a "County" field that has the name of the county the municipality is in, and a "Muni_Code" field that has the 4 digit municipal code. To define the output location my plan is to do something like:

outLocation = os.path.join(r'C:\Counties', gdbName, FDName)

I just don't know how to define gdbName and FDName. I'm thinking that arcpy.SearchCursor, the List function (eg arcpy.ListGeodatabases), or even arcpy.da.Walk would be useful but I'm not sure how to use them in this case. I'm especially struggling with how to get FDName - I need to match the Muni_Code field to the last 4 characters of the FDs. The gdbName should be easier because the name of the gdb exactly matches the "County" field.


texts: Will there be a new ArcGIS API for Silverlight that will be updated for ArcGIS 10.4?

texts:
Will there be a new ArcGIS API for Silverlight that will be 100% compatible with ArcGIS Server 10.4 (I know a lot of 3.3 API works, but I also know that older APIS never completely work and have to be adjusted).

I also know that Esri is retiring 3.x on 2016-06-01, I'm thinking there's a good chance that there will be no new release. Can anyone confirm this?


texts:
Pgr_drivingdistance function works as :

SELECT seq, id1 AS node, cost
    FROM pgr_drivingDistance(
            'SELECT id, source, target, cost FROM edge_table',
            id, distance, false, false
    );

How could I get this to work substituting x, y for id?

SELECT seq, id1 AS node, cost
    FROM pgr_drivingDistance(
            'SELECT id, source, target, cost FROM edge_table',
            x,y, distance, false, false
    );

Checking the pgrouting functions in Postgis I have found SQL code to change but not for pgr_drivingdistance.


texts: Plotting a certain range in azimuthal projection using ggplot2

texts:
I would like to limit a map to a particular range of values (i.e. northern hemisphere) under an azimuthal or orthographic projection in ggplot2. When I set ylim, I receive an error message that relates to geom_path() and certain portions of the map are missing (note the triangle cut out of the United States). Any ideas on how to plot a particular region using a non-cartesian projection without losing portions of the shapefile?

world<-map_data('world')

p<-ggplot(data=world,aes(x=long,y=lat,group=group)) +
ylim(30,90) +
geom_polygon(fill='skyblue') +
geom_path(color='gray60',aes(x=long,y=lat),size=0.1) +
coord_map("orthographic")
p

Warning message:
Removed 276 rows containing missing values (geom_path).


texts: MySQL GeoJSON query for polygons format error

texts:
Trying to take MySQL polygons to GeoJSON. My SQL query is Okay and works in PhpMyAdmin and I'm trying to modify some code that works for points. I think I'm close but stuck.

In PhpMyAdmin the records are returned and the geom field returns:
{"type": "Polygon", "coordinates": [[[-111.4754251.....

I'm having trouble inserting this in my formatted geometry array.

 $sql = "SELECT `FID`, `STATE_NAME`, `SUB_REGION`, `STATE_ABBR`, `POP2010`, `POP10_SQMI`, `MALES`, `FEMALES`, `SQMI`, ST_AsGeoJSON(`geom`) FROM `usa` WHERE 1";

    $result = mysqli_query($connection, $sql) or die("Error in Selecting " . mysqli_error($connection));

    while ($res = mysqli_fetch_array($result,MYSQLI_ASSOC)) {

        $msg [] = array(
                'type' => 'Feature',
                'geometry' => array(
                    $res[ST_AsGeoJSON('geom')],
                ),
                'properties' => array(


texts: Advantage(s) of using underscores in ArcGIS file naming conventions?

texts:
In an introductory GIS course, it was insisted that all spaces within file and folder names be replaced with underscores.

Why is this taught?

Is there a clear benefit of doing so?


texts: Spatial Reference not working with If-Else

texts:
I'm attempting to use a script to check a folder of data sets to see if they match a projection of another data set in a target folder. I need to be able to re project any that don't match the same projection as the file in the target folder and then copy the datasets to a projected folder. If the projections do match I need to copy this to the projected folder without re projecting.

I'm coming unstuck with my if/else statement at line 28 where:

if srcDesc.SpatialReference.Name == targetDesc.SpatialReference.Name:

When i run the script via a custom toolbox in Arcmap 10.4 i get the following error:

***line 27, in <module>

 if srcDesc.SpatialReference.Name==targetDesc.SpatialReference.Name:
AttributeError: DescribeData: Method SpatialReference does not exist***

My code is not perfect in other areas and I'm hoping someone may be able to help me work out where I've gone wrong?
Please excuse the over excessive commenting also :)

I'm thinking my error is something to do with describing the two datasets and with the spatialReference usage.

# Import arcpy module
import arcpy
import os

# Get input, target projection and output locations
datasetFolder = arcpy.GetParameterAsText(0)
targetProjectionDataset = arcpy.GetParameterAsText(1)
projectedDataset = arcpy.GetParameterAsText(2)

# set workspace to the dataset folder
arcpy.env.workspace = datasetFolder

# returns a list of files in the dataset folder
for inputFC in arcpy.ListFeatureClasses():

    # Use AddMessage to print feature classes in dataset folder
    arcpy.AddMessage("Datasets in the dataset folder are: " + str(inputFC))

    # describe the dataset in the dataset folder
    srcDesc = arcpy.Describe(inputFC)

    # describe the target dataset
    targetDesc = arcpy.Describe(targetProjectionDataset)

    outputFC = projectedDataset

    # test for matching coord systems
    if srcDesc.SpatialReference.Name == targetDesc.SpatialReference.Name:
        arcpy.CopyFeatures_management(inputFC, outputFC)
    else:
        # set the output name to match the input name
        outputNewFC = os.path.join(projectedDataset, inputFC)

        # Set projection to be used
        outCS = arcpy.SpatialReference(targetDesc)

        # run projection tool
        arcpy.Project_management(datasetFolder, outputNewFC, outCS)

        # check messages
        print(arcpy.GetMessages())


texts: How to obtain an extent of a whole shapefile?

texts: I have got two shapefiles. One of them is a point shapefile and the second one is a polygon shapefile. I need to find out which of them has a larger extent and then set the environment extent to include both shapefiles.


texts:
I am currently working on a project where I am responsible for calculating the mean temperature and precipitation data across the USA.

My calculations are made for metro census data and county census data and it would be greate for me to add some information about population in those areas. The problem is that I can't find any data for metro areas that are older than 1970 and I need to have a data from 1950.

I am wondering if there is some dataset with demographic data on metro or county level that will have an information about population starting from 1950 or even older?


texts: Per-pixel (statistical) calculations on a raster stack using GDAL

texts:
In R, it's relatively trivial to perform per pixel calculations based on a raster stack (e.g., get std.dev for each pixel on a 12 layer GeoTIFF). Unfortunately, the speed is less than desirable when working with large rasters.

Recently, I've been trying to use GDAL for this purpose, but I can't seem to get the expression down. I've tried the following, but it doesn't seem to work as advertised (or how I think it should work):

gdal_calc.py -A tmean_monthly_1980.tif --calc="std(A)" --outfile=tstdev_1980.tif --allBands=A --calc="std(A)" --NoDataValue=-9999

In the above example tmean_monthly_1980.tif is a raster composed of 12 layers (Jan-Dec) of mean monthly temperatures in 1980 (stacked using gdal_merge). Although the GDAL command produces an output, the results are definitely not the standard deviation per pixel.

Another option might be to use single band files, but gdal_calc.py appears limited to "only" 26 inputs (A-Z; i.e., "-A file -B file -C file etc."), and the resulting command line/script would be fairly cumbersome when working with large, multi-layered rasters (e.g., time series).

Is there a specific syntax for performing per-pixel raster calculations on a multi-band image using GDAL? If only individual files/bands can be specified using gdal_calc.py, please provide a scripted example using Bash, Python, etc.

Specifically, my goal is to calculate the standard deviation and cv for each pixel in a multi-layer raster stack, or using individual files that make up the same raster stack.

Note: I'm holding out for a GDAL specific answer, but other FOSS solutions are welcomed and will be upvoted, especially if they are relatively compact and adaptable to large number of inputs.


texts: Forming continuous network out of discontinuous lines in ArcMap

texts:
I have river data, which consists of lines (narrow sections of streams) and polygons (wider sections of streams). I would like to have all data in line form.

First, I tried Collapse Dual Lines to Centerline tool. For that I converted river polygons to lines and deleted parts from both ends to form dual lines.

Outcomes were varied a little, depending on settings, but results were not really satisfying.

For second option, I tried Vectorization to Centerline in ArcScan toolbar. For that, I turned polygons to rasters.

This time results were somewhat better.

Next, I would like to connect newly formed lines to existing line network. Data includes 3842 polygons, so manual editing is not practical.
How can I handle this and more unusually shaped sections?


texts: Adding lines connecting labels to points in QGIS?

texts:
I'm working on a map that has a lot of points in a small area, so I'm adding lines from my labels to their points using a Geometry generator in the symbol designer.  The problem I'm having is that the lines always come from the bottom left corner of the label.

How do I get them to come from the part of the label that is closest to the point?

I'm using the following Expression and really don't have a clue where to go from there.

    make_line( make_point( $x,$y), make_point( "X","Y"))


texts: Concave Hull with points that have to always intersect

texts:
I use Concave Hull(Alpha shapes) tool in my set of points in QGIS, it works fine, but I want add some extra points in the beginning, which line of Concave Hull always intersect them.

In this image, we see a brown polygon in the background, which shows what it should look like. Orange polygon is output of Concave Hull and through the green points I want it to intersect.
Did anyone have this problem too? I think it is necessary to redesign the algorithm in QGIS in some way, but I don't know how.

EDIT: Separate images for clear view. It shows just clipping, the entire data are large:
Brown polygon, which shows what it should look like:

Output of Concave Hull with points, which I want intersect:


texts: Proper basemap/projection for displaying Alaska state boundary without distortion in ArcGIS Pro

texts:
This is a simple question - when working with basemaps in ArcGIS Pro, the default projection centers around the lower 48 states of the US. I am trying to do a basic visualization of a watershed within the state of Alaska, but I notice when I use shapefiles downloaded from the National Map database for this purpose, they appear distorted, looking like the top part of Alaska is curving away.
I tried reprojecting the shapefile to UTM Zone 5N, but it appeared the same. What is the best strategy to achieve this goal?? Will it always appear distorted as long as the basemap is centered around a different region?
This curved projection makes sense but doesn't serve my purpose - I would like a flat rendering of the state boundary. Here is an image of what I mean:

I would like something that looks more like this:


texts: QGIS and ArcGIS Pro give different results on spatial join

texts:
I have two point layers that one has 2885 and the other has same 2885 plus 1748 more features. The second one is a result of "v.net.centrality" analysis and now I want to join the result spatially into the first point layer in QGIS 3.4 LTR.

2885 of the features have exact same coordinates and reference system so I expect all of intersecting points to join spatially to give analysis' result to first point layer but in QGIS only 335 of them are intersecting. Then I tried to do it by using PostGIS with

SELECT first.id, first.geom, second.result
FROM first
LEFT JOIN second ON ST_Intersects(first.geom, second.geom)

but the result is same with QGIS "Join attributes by location" tool, only 335 of them intersecting. I also tried "Equals, touches" rules on QGIS and "ST_Equals, first.geom = second.geom conditions on PostGIS but the result didn't change.

Then I pass the data into ArcGIS Pro and do spatial join on it. It gives exactly what I want, all of the attributes (analysis result) from second layer are added to the first layer with "Intersect" rule.

I'm wondering why this happening? What is the reason behind it and how can I resolve in QGIS?


texts: Thread safe method to create print layout from within a processing script in QGIS

texts:
I have a PyQGIS script, which dynamically creates a print layout and opens it in the layout designer, when finished (example below). It works without problems when executed from the console. Now I want to integrate it into my scripts in the processing toolbox. Executed from the processing script, which was previously also running well, the same procedure causes problems and ultimately crashes QGIS. I have opened an issue on Github for this and have received the answer that this is not a QGIS error but that my script is not thread safe.

Does anybody know a way to accomplish my objective in processing scripts?

project = context.project()
manager = project.layoutManager()
layout = QgsPrintLayout(project)
layout.initializeDefaults()

layout.setName('Layout Name')

pages = layout.pageCollection()
pages.beginPageSizeChange()
page = pages.page(0)
page.setPageSize('A4')
pages.endPageSizeChange()
page_center = page.pageSize().width() / 2

title = QgsLayoutItemLabel(layout)
title.setText('Layout Title')
title.setFont(QFont('arial', 36))
title.adjustSizeToText()
title.setReferencePoint(1)
title.attemptMove(QgsLayoutPoint(page_center, 20))
layout.addItem(title)

map = QgsLayoutItemMap(layout)
map.setRect(20, 20, 20, 20)
ms = QgsMapSettings()
ms.setLayers([catchment])
rect = QgsRectangle(ms.fullExtent())
rect.scale(1.0)
ms.setExtent(rect)
map.setExtent(rect)
map.attemptResize(QgsLayoutSize(180, 90))
map.setReferencePoint(1)
map.attemptMove(QgsLayoutPoint(page_center, 50))
layout.addLayoutItem(map)

table = QgsLayoutItemTextTable(layout)
layout.addMultiFrame(table)
cols = [QgsLayoutTableColumn(), QgsLayoutTableColumn(), QgsLayoutTableColumn()]
cols[0].setHeading('Heading 0')
cols[1].setHeading('Heading 1')
table.setColumns(cols)
table.setContents([['Hello', 'world!']])
frame = QgsLayoutFrame(layout, table)
frame.attemptResize(QgsLayoutSize(30, 30), True)
frame.setReferencePoint(1)
frame.attemptMove(QgsLayoutPoint(page_center, 150))
table.addFrame(frame) # crash point 1

manager.addLayout(layout)

iface.openLayoutDesigner(layout) # crash point 2


texts: gdal_grid performance issue -- alternatives?

texts:
I use gdal_grid within a Python script, but for my large files the calculation time is enormous. Apparently this is a gdal_grid issue http://trac.osgeo.org/gdal/ticket/2411

Does anyone have an advice / some code which can serve as a simple alternative to gdal_grid?


texts: QGIS dissolve by attribute and retain primary polygon attributes

texts:
I have a parcel shapefile. Many of the parcels have been combined over time, but the shapefile still represents the original parcel boundaries. The attribute table has a field for the original parcel number as well as, for those that have been joined, the current 'parent parcel'.
I created a new attribute that replaces the historic parcel id with the current parent. When I dissolve by this I get a shapefile of that represents the current geographic conditions, but my attribute data is not accurate. In some cases QGIS dissolves and retains the parcel record that was merged into (and not currently recognized) with no associated attribute data.

I could just join up the data again after dissolving; but there must be a way to tell the dissolve function to retain a specific 'primary' polygon. I have a dummy attribute marking the 'primary' as 1, but I don't know how dissolve chooses which record to retain; sorting doesn't seem to make a difference. Is there a way to change how dissolve decides which record to retain, and feed it my dummy? Is there an entirely different approach I should try? Or should I just shut up and rejoin all of the data?

I'm hesitant to rejoin, because there are half a million records, and most operations are pretty taxing on my poor computer.


texts: Performing viewshed analysis with Dutch AHN files in QGIS?

texts:
i downloaded the new (dutch hight data) AHN files, and was trying to perform a viewshed analysis at some point near UTRECHT.

The plugin (grass plugin) works ok, i tried it before with a dem map.

But when i use the AHN file it shows not the correct visibility area.

so my question is

Does someone know how to setup the AHN files, for using in a viewshed analysis
Does someone know how to make the real height information visible on a map


texts: QGIS 3 Graphical Modeler pre calculated values using Advanced Field Calculator

texts:
I have a new model that tries to create a new column using the Advanced Field Calculator and pre-calculated values. I just seem not able to get the syntax right on how to populate the attribute using the values stored in the Variable. Here is my example model.

The input data is depicted below

Basically what I am trying to do is fix up some bad geometries and then reclassify the layer based on attribute values.
I have tried the following formula:
 Case when attribute(get_feature(@layers,@InputFields,'WC')) then 1 ELSE NULL END

Or
Case When layer_property(@layer,@InputFields) = 'WC' then 1 else 2 end

All of the above formula do not give a clear answer. I have searched the browser for different answers but I haven't got any hints why I am getting an error when I execute the model.
Any hints or links to documentation.


texts: Dissolve lines putting segments out of order QGIS

texts:
I'm using a polyline river layer in QGIS to calculate river miles, but am ending up with points out of order.  Starting with a river layer that is broken into reaches (short sections), I'm dissolving into a new layer based on the stream name.  I am then calculating the distance along the line (I've tried two separate plugins, "Locate points along lines", and "Qchainage").
Unfortunately, some segments seem to be out of order, so the distance along the line jumps around, and sometimes go in opposite directions.
I end up with river miles that look like :

The segment with miles 9-11 is between the segments with 0-5 and 6-8
I suspect this crops up in the dissolve stage because of some feature of the original data.  Is there a way to dissolve the lines into sections that are continuous lines without misaligned segments?
Here is a link to the shapefile I'm working with.  It's in a zip on dropbox. This version has been converted to single parts, with all the lines going in the same direction.
zipped shapefile on dropbox


texts:
I'm interested in displaying data in Leaflet from PostGIS, and would particularly like to have my data as GeoJSON when passing to Leaflet. I can currently do this via:

# set connection
conn <- dbConnect("PostgreSQL",user="postgres",password="pass",port=5432,dbname="postgis_user")

# set query and reproject to EPSG:4326 for GeoJSON creation
# (the real query is more involved with a bounding box etc)
qry <- "SELECT ST_Transform(table.geom, 4326) as geom FROM table"

# i then read as an sf object
pols <- st_read_db(conn, query=qry, geom="geom")

# and convert to GeoJSON
pols.js <- geojsonio::geojson_json(pols)

This is definitely reasonable for smaller datasets, but the bulk of processing time is in the conversion of simple features to geoJSON -
 geojsonio::geojson_json(pols)

How then do i extract geoJSON directly from PostGIS? And i presume this has to be done as a FeatureCollection? my attempts and reading have not got me to a solution:

attempt 1: simplistic ST_AsGeoJSON

qry1 <- "SELECT ST_AsGeoJSON(ST_Transform(geom,4326)) FROM table"
pols1 <- dbGetQuery(conn,qry1)
## no error messages but no display in Leaflet (incorrectly constructed geoJSON i assume, returned as a data.frame)

attempt 2: building featureCollection as per this answer

qry2 <- "SELECT json_build_object(
'type', 'FeatureCollection',

'features', json_agg(
    json_build_object(
        'type',       'Feature',
        'id',         gid,
        'geometry',   ST_AsGeoJSON(ST_ForceRHR(st_transform(geom,4326)))::json,
    'properties', jsonb_set(row_to_json(table)::jsonb,'{geom}','0',false)
)
)
)
FROM table"
    pols2 <- dbGetQuery(conn,qry2)
    ## hangs and errors: "DETAIL:  Cannot enlarge string buffer containing 1073740688 bytes by 3167 more bytes."

attempt 3: another try as per this answer:

qry3 <- "SELECT row_to_json(f) As feature
  FROM (SELECT 'Feature' As type,
                ST_AsGeoJSON(ST_Transform(l.geom, 4326))::json As geometry,
                row_to_json((SELECT l FROM (SELECT id, 'l' AS type) As l))
As properties
  FROM table As l) As f"
pols3 <- dbGetQuery(conn,qry3)
## this gives a similar output to attempt 1 above but gives a warning "RS-DBI driver
## warning: (unrecognized PostgreSQL field type json (id:114) in column 0)", which
## i'm guessing is bcos i've got the syntax wrong.

Can i return geoJSON directly from PostGIS in R without having to convert a spatial feature?


texts:
In a project in which I am using a specific shapefile with data from a specific date multiple times (around 100 times) with different legends each time, every few months, I save the shapefile with new data, and hence have my maps updated.

I now want to keep the old data still available and am looking for a way, in which I can easily exchange the one shapefile I use multiple times in one project, in one batch, so as not having to change it 100 times.

I am using QGIS 2.14.14 Essen.


texts:
In my web app I have a search side bar which can return results as either a bounding box, or an X and Y coordinate with a zoom level. The sidebar pops out the side of the screen and partially overlays the map, so I want the results centred on this slightly offset map. I can do this no problem when a bounding box is supplied to the search results, as I can simply use the padding FitOption like so
//leftPadding is a pixel value equal to the width of the sidebar
//result.bbox is an Openlayers Extent
olMap.getView().fit(result.bbox, { padding: [100, 100, 100, leftPadding] });

However, I cannot (to my knowledge) use .fit() with a coordinate and a zoom level.
This is how I zoom to a point at the moment (uses animation but doesn't have to)
olMap.getView().animate({
        duration: 1000,
        center: convertedCoord,
        zoom:result.zoom
    })

result.zoom is an integer indicating zoom level to zoom to
convertedCoord is an Openlayers Coordinate
I could use setView() like this
let view = new View({
    center: convertedCoord,
    padding: [100, 100, 100, leftPadding],
    zoom: result.zoom
});
olMap.setView(view);

However, this shifts the actual view of the map, so rotating now rotates around that point, and when I close the sidebar I would need to reset the padding every time.
What I ideally want is some way is to setCenter() with padding, but I can't seem to find a way to do that. Ideally I would want the animation to still work, but if it can't so be it.
I feel there is a possibility of using coordinate.add() to move the center point, but I can't quite figure out how to determine the delta required to move it to get the same padding.


texts: Generating global view Image of where Lat-Long location is?

texts:
I would like to automatically generate an image which highlights a location from a global view.

Here is an example of what I would like to achieve:

I am debating using python utilities such as mpl_toolkits.basemap,which examples can be seen here. But I was hoping to find an easier way from you GIS professionals.

This is for a small webpage, and it would be nice to have it automatically generated.


texts: What is the purpose of the KML XMLNS in a KML's source?

texts:
I am looking over some KML documents source code and have noticed a difference in some of them.
Some have:

kml xmlns="http://earth.google.com/kml/2.1">

and others:

kml xmlns="http://www.opengis.net/kml/2.2">

I was would like to know the purpose of this string (as a complete novice), does it have any implication to how a KML is generated / operates?


texts: Do I need IMU data for image stitching and orthomosaic generation?

texts:
Background: I am working on low-cost raspberry pi based camera system to integrate it on drones as an external module. So, I am using Pi camera and UBLOX Low cost GPS receiver which logs data in CSV file.
So, is this data sufficient to generate orthomosaic for using in GIS applications?
Also, Does the system need IMU data for orthomosaic generation? If no, What is the use of orientation x/y/z or roll/pitch/yaw  in photogrammetry?


texts: Pasting features into a feature class doesn't honour subtypes

texts:
I'm encountering an issue where features that I'm copying from a feature class in a file geodatabase are all going to a single subtype when I paste them in the destination feature class in an enterprise geodatabase (Versioned). Subtypes are identical in both feature classes. The paste operation is not honouring the subtype of the original data.

I'm using ArcGIS 10.3.1.

My workflow:

Start an edit session on a version of an enterprise geodatabase.
Select 15 features in a file geodatabase feature class called SewerMaintenanceHoles.
Copy these features and paste them into the corresponding feature class in the enterprise geodatabase. It is also called SewerMaintenanceHoles.
Because the destination feature class contains a subtype (which is identical to the source feature class subtype), I get the following option when I paste:

You can see in the attribute table the source values are shown and there are 3 different possible values (Maintenance Hole, Maintenance Shaft and TEP). I don't want to paste all the features into a specific subtype. I want the paste to honour the original values in the source feature class. So, I chose the top level (Sewer_Maintenance_Holes). I thought that would paste the features into their corresponding subtype. What is happening instead is that all features are pasting to the default subtype which is Maintenance Hole.

Note, USE is an alias for the Subtype field. But, they're all coming in as Maintenance Hole.

How can I get the original subtype values to come through to the destination feature class?


texts: Changing resolution of rasterbrick using R

texts:
I have a rasterbrick « env1 » with following informations :
class      : RasterBrick
dimensions : 3008, 2913, 8762304, 8  (nrow, ncol, ncell, nlayers)
resolution : 0.00331, 0.00207  (x, y)
extent     : -73.38788, -63.74585, 48.13865, 54.36521  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs

I wanted to change the resolution to have the same x and y with :
res(env1)<- 0.00331

But went I wanted to plot it, it shows the following error :

Error in .plotraster2(obj, col = col, maxpixels = maxpixels, xaxt =
xa,  :  no values associated with this RasterLayer

Is there any way to change the resolution in order to have the same x and y?


texts: Cut linestings at nearest point junctions

texts:
Using the "function to cut linestings and multilinestrings at nearest point junctions" found in chapter 11 of the book PostGIS in Action by Regina Obe and Leo Hsu (Second Edition)

It returns this error when running in PgAdmin 4:

ERROR:  syntax error at or near "END"
LINE 40:    END IF;
            ^
********** Error **********

ERROR: syntax error at or near "END"
SQL state: 42601
Character: 1064

Code lines around line 40:

        ELSE
         var_lset[j] := kutt(var_lset[j],var_pset[i])
        END IF;

What is wrong?


texts:
I want to add a context menu on right click for various elements on my GeoJSON layer (I'm doing a road map so on a right click on the road at any part I want to show the context menu).

I've managed to get the left click working fine by using the onEachFeature and doing the following

function onEachFeature(feature, layer) {
        layer.on({
            click: showAssetInfo,
            contextmenu: contextreg
        });

    }

 function showAssetInfo(e) {
        AssetMouseClick(e.target.feature.properties.objectid, e.latlng);
    }

For the context menu I have followed the example here  . The context menu library is found here

I have the following that gets called on the document ready (jQuery)

$.contextMenu({
    selector: 'path.leaflet-clickable',
    zIndex: 99999,
    callback: function (key, options) {
        var m = "clicked: " + key;
        window.console && console.log(m) || alert(m);
    },
    items: {
        "edit": { name: "Edit", icon: "edit" },
        "cut": { name: "Cut", icon: "cut" },
        "copy": { name: "Copy", icon: "copy" },
        "paste": { name: "Paste", icon: "paste" },
        "delete": { name: "Delete", icon: "delete" },
        "sep1": "---------",
        "quit": { name: "Quit", icon: "quit" }
    }
});

I've tested it and the selector does return the GeoJSON features, also if it attach the same menu to something else it works correctly.

Is there something I am missing here?

Also is there a good way to pass in the objectid to the menu when it starts up so I can use it when calling the different options of the menu

I have created this fiddle to demonstrate http://jsfiddle.net/Q3L4c/22/


texts: I am doing a research project on the north south divide, and have asked people to draw on a simple outline of England where they feel the north south divide should be in England. However, I am unsure how to now accurately plot this into Arcmap.


texts: Convex hull that follows irregular points distribution

texts: Based on an irregular distribution of points (longitude,latitude), I obtained the convex hull below (I use saga_cmd with polygon convexity in 1). Would you have a way to create a convex hull that closely follows this irregular shape using one-liner commands or similar? I call commands from a bash terminal in Linux.


texts: Adding new features with attributes to an existing layer in PyQGIS?

texts:
I have a vector layer, populated with points, and I want to programmatically add to this layer.

The part related to the geometry works just out of the box, but I was having trouble with adding attribute values to the new features.
The problem is, when I add a new feature to the layer, I do not manage to specify fields. The feature would complain on exceptions.KeyError.

        featureList = []
        ## now add the computed points to the layer
        for p in [p for p in points.values() if p['computed']]:
            x, y = p['coordinates']
            feature = QgsFeature()
            layerPoint = transf.transform(QgsPoint(x, y))
            feature.setGeometry(QgsGeometry.fromPoint(layerPoint))
            feature['id'] = p['id']
            featureList.append(feature)

        layer.dataProvider().addFeatures(featureList)
        layer.commitChanges()

And before you ask, obviously I'm sure it's the feature['id'] which is giving trouble, not the p['id'].

So the short version of the question was "how do I specify fields of a feature before it was ever part of a layer?"


texts: OpenStreetMap's Nominatim indexing restarted after rank 29 to rank 0

texts:
I'm trying a new installation of OpenStreetMap's Nominatim under Linux Debian on a small server that has 4 threads and sufficient space for full planet import.
I used the following command to resume indexing after rank 26 and when it came to end of rank 29,
it suddenly restarted to rank 0...
./utils/setup.php -v --index --create-functions --index-noanalyse --create-search-indices --osm2pgsql-cache 10000

Is it normal or did I make something wrong? Should I restart the import?


texts: Multiple Animated Layers on GIS Map User Movements

texts: I'm relatively knew to GIS, so I figured this might be a smart way to go about this. I have a data set with multiple users in it I was going to use Torque from cartodb to show in realtime the movements of users simultaneously, until I found out that I can only have one torque layer. So I'm wondering is there a way to do this is cartodb?


texts:
How can I find depression contour lines?
I already have find depressions via Fill-method, but now I try to detect them by using the contour lines, anyone have an idea ?


texts: Combine Raster and Feature Class selection

texts:
I've been working on a quite complex analysis with ArcGIS 10.2. But now I'm stuck at a problem...
I have a raster file and a polygon feature.
The Raster file stores information about Zmax, Zmin and the Range.
The Polygon feature stores a height value

What I need: I need to select Polygons with a specific value (above or below Zmax) that are within a certain distance to the raster cell with the specific Zmax value.

I know I can work around and use raster to point and then make an SQL selection (by location and by attribute); select all the points = specific Zmax + select polygons within the certain distance + select out of this the ones with the specific value

But I was wondering whether there is a shorter way/possibility to automate the whole process? I've been working with model builder before but I'm not sure if model builder can be used to automate this task


texts: Customizing a layer from public map server in ArcGIS Online

texts:
I added a layer from the web via publicly available map server in my ArcGIS Online Web Map. Now that I have this layer in my map, I'd like to limit the extent in which it appears. This layer covers the entire contiguous 48 states of the US and I really just need it for one county.
Is there a way to limit the extent, either in my Web Map, or via Web App Builder?


texts:
I have precipitation data from NOAA in a NetCDF format. I imported that into ArcGIS as a raster layer. Now I see that I can only get the precipitation value for one year at a time.

I need to calculate the annual average precipitation for the decade (1980-1990).

Is there any way to do this?

I am very new to GIS and don't know how to use Python. I have dabbled with the ModelBuilder a bit.


texts: Calculating Centroid of all lands of Earth using SQL Server Spatial?

texts:
I want to find the center point of all lands of the earth.

What I have tried up to now is downloading shapefile of all countries and polar parts of the earth. Then I imported that shapefile into a SQL Server 2008 database and then tried to run some queries to find that point but all of them were in vain.

How can I calculate the Centroid of all lands of the earth using SQL Server?


texts: What Open Source (command line/API) method is there to orthorectify and georeference an aerial image given telemetry data?

texts:
How can I take an image and orthorectify it automatically with open source software?
The information associated with each image is:

camera specs
latitude, longitude, elevation of the camera
orientation of the camera: facing down off a UAV, roll pitch and yaw are given

I want to output a geotiff file and hopefully end up stitching the images together to create a map.


texts: QGIS 3.2: Where are layout files stored?

texts:
I am changing computers at work and I want to be able to bring all of my old layouts (the ones accessed through the layout manager) over to my new computer. Each layout is set up on a particular location in the exact place I want it. If I have to recreate these in the future, I don't think I could make them 100% the same and I would have to redo hours of work.

I searched through the QGIS files and I didn't see anything that was definitely the save point for each layout.

Does anyone know where these are saved?


texts: Producing multiple PDF plots of mapsheets specified by shapefile key from ArcGIS Desktop?

texts: I'm looking for a product that will produce multiple pdf plots of mapsheets specified by a shapefile key (say, 11 x 17) from ArcGIS 9.3 using a single .mxd template, include and orient a common keymap, etc.  The DS Mapbook extension from the ESRI developer samples is one of these, but it's an old VB6 product and am wondering if there's anything else similar, newer and maybe commercially available.


texts: overwrite feature class in arcgis online using python

texts:
I have access to Arcgis for Server and am trying to overwrite a couple feature service layers I have hosted on arcgis online. I usually overwrite these manually by uploading zipped shapefiles from a folder on my desktop. I'd like to simply overwrite these feature service layers using python possibly with arcrest but don't know how to do this. I've found a lot of mxd to sd file conversion workflows with python but I want to use zipped shapefiles to upload to arcgis online not mxd's. I also am willing to use arcgis for Server if this is more efficient or secure. So far, my upload attempt has failed using the code in the answer below. with a
    KeyError: 'id'
Hoping to eventually upload/publish/overwrite several zipped shapefiles with for loops but for now I'm just looking to figure out how to do one at a time.


texts: Calculating surface area above a certain point on a mountain or calculating surface area above a certain elevation

texts:
I have a station point somewhere (hydrostation) on a mountain, I need to calculate the complete surface area/Landmass of the mountain above this station point.

There are approximately 50-70 station points.


texts: create dynamic geometry point from WKTElement geoalchemy2

texts:
I have a function

def get_centres(origin_lat, origin_long):
     pt = WKTElement('POINT({0} {1})'.format(origin_long, origin_lat), srid=4326)
     centres = Centres.query.with_entities(Centres,functions.ST_Distance_Sphere(Centres.location,pt).label('distance')).all()

The query it creates when I run the above function is

 `[SQL: 'SELECT kare.centres.id AS kare_centres_id, kare.centres.name AS kare_centres_name, kare.centres.address1 AS kare_centres_address1, kare.centres.address2 AS kare_centres_address2, kare.centres.city AS kare_centres_city, kare.centres.postcode AS kare_centres_postcode, ST_AsEWKB(kare.centres.location) AS kare_centres_location, ST_Distance_Sphere(kare.centres.location, ST_GeomFromText(%(ST_GeomFromText_1)s, %(ST_GeomFromText_2)s)) AS distance \nFROM kare.centres'] [parameters: {'ST_GeomFromText_1': 'POINT((-1.3686658,) (50.9965899,))', 'ST_GeomFromText_2': 4326}]`

This query throws an error

sqlalchemy.exc.InternalError: (psycopg2.InternalError) parse error - invalid geometry
HINT:  "POINT((-" <-- parse error at position 8 within geometry

The problem is WKTElement is converting geometry point to

'POINT((-1.3686658,) (50.9965899,))'

When postgress sql is expecting below format.

`POINT(-1.3686658 50.9965899)`

I copy the query from above and run it on a Postgres Sql client and it throws the same "invalid geometry" error.

When I format the point and run the query I get the results.

How do I pass the right geometry point


texts: ArcGIS JS API 4.4 - mouse over (hover) event

texts:
I'm trying to write a handler for a mouse-over event on a FeatureLayer in ArcGIS JS API 4.4. It seems in the previous API this was built in, but I'm not finding it anywhere in the new API. Am I missing something or was this just lost in the overhaul? If so, is there a workaround?

This thread mentions an option but the link doesn't take you there, I'm guessing this was for the old API.


texts: I found a way to extract data from Copernicus SciHub, but I'd love to find a way to automate to pre-processing in a Python script for both S2 and S1 data without having to use SNAP Desktop. Is there any way to use the algorithms like sen2cor so I can have a resulting geoTIFF?


texts: Splitting regular polygons into same-height segments

texts:
I have vector layer which looks like this:

It is result of raster classification, vectorization and finally Minimum Bounding Box Geometry. What I need is, that maximum "height" of one polygon is for example 20 cm. It means the small rectangles stay the same as they are smaller than this threshold, but the "long ones" divide into several smaller polygons. After I will create centroids and count number of these polygons in one row, so first I need to have the small polygon segments.


texts:
Why ST_AsText(location(ST_IsValidDetail(the_geom))) from LINESTRING return blank? From POLYGON, all ok... "POINT(...)"

For example:

SELECT ST_Astext(location(ST_IsValidDetail(ST_GeomFromText('POLYGON((0 0, 1 1, 0 1, 1 0, 0 0))')))) as location

--> return one row with "POINT(0.5 0.5)"

location
|"POINT(0.5 0.5)"|

but

SELECT ST_Astext(location(ST_IsValidDetail(ST_GeomFromText('LINESTRING(0 0, 1 1, 0 1, 1 0, 0 0)')))) as location

--> return one row with nothing

location
||

Documentation:
http://postgis.net/docs/ST_IsValidDetail.html

 --simple example
 SELECT * FROM ST_IsValidDetail('LINESTRING(220227 150406,2220227 150407,222020 150410)');

 valid | reason | location
-------+--------+----------
 t     |        |


texts: Change openlayers map div options dynamically

texts:
I have two HTML div on my page,

I want to move my map from "map" to "map2" but not working. I am using map.setOptions method.

var map = new OpenLayers.Map({
    div: "map",
    center: new OpenLayers.LonLat(0, 0),
    minResolution: "auto",
    maxResolution: "auto",
    layers: [
        new OpenLayers.Layer.OSM()
    ]
});

map.setOptions({
    div: "map2",
    center: new OpenLayers.LonLat(0, 0),
    minResolution: "auto",
    maxResolution: "auto",
    layers: [
        new OpenLayers.Layer.OSM()
    ]
});

DEMO is here.


texts: Alternative to Copy Features in ModelBuilder to display a result

texts:
I have created a model that takes in a raster, converts it to a point shapefile, then extracts a set of points based on location, then further selects a set of points based on a SQL expression. I tried to run it once but it was taking absolutely forever. I would have thought this step to be akin to 'Export Data' in the regular window. Is this the right method? Current model creation below.

The goal of this procedure is to find a region of high density crops in a selected country. I will eventually need to determine the coordinate footprint of the area so I can export that information to fine tune a forecast model to focus on that particular area.


texts: Converting Worldview 2 imagery to TOA reflectance for Vegetation Change detection in ArcGIS?

texts:
I am working on a change detection study to look at effects of herbicides on vegetation growth. Therefore I have Worldview 2 (4 band +pan band) imagery for two subsequent years that I would like to use for a change detection analysis, but need some advice on appropriate pre-processing steps.

My proposed techniques to evaluate the change is first is a simple comparison of NDVI between the two years and secondly to run a Random Forest Classification in R. I have been trying to follow the guidance outlined in the Digital Globe Radiometric Use of Worldview 2 Imagery document, but am having a bit of trouble deciphering exactly what formulas I need to use to make my imagery comparable between years.

Initially, my workflow was to:

1) Convert to DN to Top-of-Atmosphere Spectral Radiance using the second formula from section 5.1 of the document (see equation 1 below) using ArcGIS raster calculator.

2) Convert results from step 1 to Top of Atmosphere band-average reflectance using the last equation from section 7 of the technical document (see equation 2 below) using ArcGIS raster calculator.

3) Run the NDVI and R forest classification using the appropriate Top of Atmosphere band-average reflectance layers.

However, after re-reviewing the Digital Globe technical document, I am wondering if I also need to complete the radiometric balancing for multiple scenes outlined in Section 6 (See equations 3 or 4 below) or is step already accounted for in the formula for converting to Top of Atmosphere band-average reflectance (Equation 2)?

Equations From Technical Document Proposed: