texts: Orthogonal Height Profile along line (sextante -> v.segment, r.profile?)

texts:
While trying to find a way to find "valleys/ridges" within a DEM, I'd
like to know if there are tools, plugins or "methods" known or
available to automatically generate height profiles along a given
path.

Or maybe a tool that (allows to) calculates the height profile
(UNIX preferred, no Windows at all) at a certain position, takin
coords and azimut as input?

--- / edit 19.12. --- /

Even though I'm yet not used to sextante and GRASS at all, I've searched and read a lot. Shouldn't it be possible to solve above discribed problem by using the sextante modeller?

If there is a single line within a line-shape plus a raster elevation map...

Shouldnt it be possible to...

Step: generate points using v.segment with direction (within line) as attributes for the resulting points P.
Step: generate 2 new points for every P within a fixed given distance, which mark the start and end of an orthogonal line (angle from "directon" at step-1), where P is the middle resp. point of intersection.
Step: calculate profile of this orthogonal using r.profile
Step: write/append elevation profile to a CSV with lat/lon/height x/y/elev, so that it can be further processed or simply beeing 3D illustrated via Speadsheet, gnuplot,..

The lan/lot is delete cause I didnt point out I would like to have it "straightned". Otherwise I could simply wrap the area from the DEM.

An please - as mentioned - I'm new to GRASS and even less when it comes to the modeller :-/   It's not too long ago I started using the sextante toolbox to improove the possibilities within QGIS.

Below show my two rasters examples. First one is straight ASTER elevation data overlayed by a rugh drawn line from a geological map, respresenting a tectonic fault.

Second on is the terrain convergence index (tci) of the DEM, where I get the impression that I can see this line.


texts: How can I download data from NYC open data into GIS?

texts:
I have a data from NYC open data website and I need to download them to GIS so it will be easier to analyze
https://data.cityofnewyork.us/Public-Safety/NYPD-Motor-Vehicle-Collisions/h9gi-nx95

How can I open and explore data from NYC open data into GIS when I only have longitude and latitude figures, and there is no coordination system given?


texts:
I'm new to ArcPy. I have coordinates of points

70°8'12.938"E  ,25°26'11.569"N

79°10'33.945"E,  25°43'58.898"N

80°25'46.695"E , 16°45'15.907"N

70°6'57.527"E , 17°10'51.245"N

I tried below  code:

>>> fc=r'D:\Arcpy Practise\arcpy.gdb\Points_UTM43N'
>>> file=open(r"D:\Arcpy Practise\xy.txt","r")
>>> read=file.readlines()
>>> with arcpy.da.InsertCursor(fc,["shape@"]) as cur:
...     for line in read:
...         x,y=line.split(",")
...         pt=arcpy.Point(x,y)
...         cur.insertRow(pt)

code shows RuntimeError: Point: Input value is not numeric


texts:
I am trying to map data by querying a table called xyz. I need to show the the_geom data on the map.

I have below query in my MapServer .map file:

DATA "the_geom from xyz using unique id using SRID=xxxx"

I need to generate a shapefile in order to get the data. I have

SHAPEPATH "../data" defined in my .map file.

But I am not sure command I should write for generating shapefile for above query, I tried doing

/usr/local/pgsql/bin/pgsql2shp -f /home/name/public_html/proj/data/xyz -h xxx.xxx.xx.xx -u username -P passwrd dbname "select the_geom from xyz"

This generated some files which I have it in by data folder. But I see blank page on my webpage.

Please correct me!


texts: Generating routes between every point in a set of points in QGIS

texts:
I'm looking for a way in which to generate a set of routes between every point in a layer.

To be clear, if I have four points [A,B,C,D] I want to get routes:

A - B
A - C
A - D
B - A
B - C
...

... and so on for a total of (4*3=) 12 routes in this example; in fact, I have 35 points in my dataset so it's a bit more than that!

I can't ascertain a way to do this programmatically in QGIS using the Road Graph function or the OSM Route plugin (although it does give good results for a single route), and web searches have turned up nothing (most answers relate to routes that follow a number of points).

Preferably, as this is for cycling purposes, I'd like to use a routing service that would take the method of travel in to account, but not strictly necessary as I'd be happy just to use my line vector layer that I have for the area containing my points.


texts: Can't retrieve geometry (Python, OGR, Rtree)

texts:
ISSUE FIXED SEE: '.Simplify' causes 'null geometry' errors (Python, OGR)

Im trying to create a spatial index on one of my layers for the purpose of performing an intersect analysis. However, I get the following error:

AttributeError: 'NoneType' object has no attribute 'GetEnvelope'

It seems that when I iterate over the features and try to use '.GetGeometryRef()' no geometry is returned. When I test '.GetGeometryRef()' outside of the loop it doesn't throw me an error. What could be the issue? Here is my code:

import sys, os, rtree
from osgeo import ogr

inFile = r'C:\Data\intermediates\class.shp'
driver = ogr.GetDriverByName('ESRI Shapefile')
inData = driver.Open(inFile, 0)
inLayer = inData.GetLayer()

feature = inLayer.GetFeature(0)
geom = feature.GetGeometryRef()
if geom is None:
    sys.exit(1)
feature = None
geom = None

index = rtree.index.Index(interleaved = False)
for fid1 in range(0, inLayer.GetFeatureCount()):
    feature = inLayer.GetFeature(fid1)
    geom = feature.GetGeometryRef()
    xmin, xmax, ymin, ymax = geom.GetEnvelope()
    index.insert(fid1, (xmin, xmax, ymin, ymax))

UPDATE

Apparently the 'class.shp' contains geometry errors according to the Check Geometry tool in ArcMap, (71 issues): 'null geometry'  and ' could not find spatial index'. I created class.shp from an original shapefile only by simplifying the features. I created the layer with this code, did I do something wrong here?:

if.os.path.exists('class.shp'):
    driver.DeleteDataSource('class.shp')
outFile = driver.CreateDataSource('class.shp')
outLayer = outFile.CreateLayer('class.shp', srs, geom_type = inLayer.GetLayerDefn().GetGeomType())

selectLayerDef = selectData.GetLayerDefn()
for i in range(0, selectLayerDef.GetFieldCount()):
outLayer.CreateField(selectLayerDef.GetFieldDefn(i))

outLayerDef = outLayer.GetLayerDefn()
for i in range(0, selectData.GetFeatureCount()):
    inFeature = selectData.GetNextFeature()
    outFeature = ogr.Feature(outLayerDef)
    for i in range(0, outLayerDef.GetFieldCount()):
        outFeature.SetField(outLayerDef.GetFieldDefn(i).GetNameRef(),
                            inFeature.GetField(i)
    geom = inFeature.GetGeometryRef()
    geomSimple = geom.Simplify(10)
    outFeature.SetGeometry(geomSimple.Clone())
    outFeature.SetField('FEAT_TYPE', 'LAV')
    outLayer.CreateFeature(outFeature)
    inFeature.Destroy()
    outFeature.Destroy()


texts: ArcMap Calculate Field..if "text" in !Field!

texts:
Is it possible to use Field Calculator in ArcGIS to add text to a new field based on an if in statement?

Example: In an attribute table I have a some fields called !X!, !Y!, !Z! where "A (xx%)", "B (xx%)", "C (xx%)", "D (xx%)" are the data.

The data can be found in any field but does not have to.

I would like to use python in some form so that once I have created a New Field, I can calculate if A in !X! copy A in !X! to New Field then continue this through !Y! and !Z!.

I believe outside of Field Calculator this will involve If in statements followed by arcpy.da.insertcursor() but I'm pretty stumped at the moment.

I hope this is clear...thanks for any help and comments!


texts: How to set date as default value in Arcmap?

texts:
I need to set default value of the field "ExpireDate" to 31/12/9999. When I'm trying to set it it changes to the number 2958465.

What is the problem here?


texts: Join and export .csv file as .shp in QGIS?

texts:
Here's what I'm trying to achieve:

Load the required GIS layer (*.shp file)
Load the .csv file in QGIS
Perform the join based on common attributes in the two files
Export the mapped .csv file as a .shp file

This is where I'm having trouble. The error that I'm getting is

"Invalid Data Source: C:/GIS Projects/Services.shp is not a valid or recognized data source"

I am so used to doing this easily in ArcGIS and MapInfo I don't understand why I cannot get this to work on QGIS. I have write permissions and everything and I can see 3/4 files being created in the specified folder just the .shp file is missing in the output folder.


texts:
I am counting the number of features in a geodabase. I  am triying to insert the feature name and the count inside a dbf table which has only two fields: "FEAT" (Text field which containts the feature name) and "COUNT" (short field which contains the feature count).

This is my code:

  listlayers= arcpy.ListFeatureClasses()
  for feat in listlayers:
    count = 0
    count = arcpy.GetCount_management(feat)
    feat = str(feat)
    if count > 0:
      print feat + " has " + str(count) + " entries"
      with arcpy.da.InsertCursor(destinationgeneraltable, ("FEAT","COUNT")) as icursor:
         icursor.insertRow((feat,count))
    else:
      print feat + " is void"

I use an insert cursor to insert these records inside this table but it fails without any error message.


texts: CSV to Geodataframe : How to have valid geometry objects?

texts:
I'm writing a script with Geopandas. I try to use a csv of blocks to make a spatial join. So I convert it as a Geodataframe. But when I want to set geometry column it returns me Input geometry column must contain valid geometry objects.

Here is my code to import csv file :

csv_df = pandas.read_csv(csv_file)
csv_gdf = gpd.GeoDataFrame(csv_df)
csv_gdf = csv_gdf.set_geometry('geometry')

Here is csv_gdf.head() before I try to set geometry column :

    id      name shortName  accountId  isMonitored  varietalId
ranchId  \
0  14633.0    HC4bas      HC4b      346.0        False         4.0
855.0
1  14634.0   HC3haut      HC3h      346.0        False         4.0
855.0
2  14637.0      HC12      HC12      346.0        False         2.0
855.0
3  14638.0  HC11haut      HC11      346.0        False        72.0
855.0
4  14641.0    HC9bas      HC9b      346.0        False         4.0
855.0

inRowDistance  betweenRowDistance  \
0            1.2                 1.5
1            1.2                 1.5
2            1.2                 1.5
3            0.9                 1.5
4            1.2                 1.5

                                        geometry       ...         \
0  POLYGON ((-0.1642995066034836 44.9397295596186...       ...
1  POLYGON ((-0.1634854066129132 44.9405302549332...       ...
2  POLYGON ((-0.1624824342183362 44.9398350833047...       ...
3  POLYGON ((-0.1592356652491378 44.9399712591478...       ...
4  POLYGON ((-0.1610166332996532 44.9391145465108...       ...

slopeInclination  slopeOrientation  soilType  rowOrientation
dashboardId  \
0               NaN               NaN       NaN             NaN
NaN
1               NaN               NaN       NaN             NaN
NaN
2               NaN               NaN       NaN             NaN
NaN
3               NaN               NaN       NaN             NaN
NaN
4               NaN               NaN       NaN             NaN
NaN

canopySystemId  canopyWidth  topWireHeight  clusterWireHeight  \
0             NaN         0.45            NaN                NaN
1             NaN         0.45            NaN                NaN
2             NaN         0.45            NaN                NaN
3             NaN         0.45            NaN                NaN
4             NaN         0.45            NaN                NaN

pruningSystemId
0              NaN
1              NaN
2              NaN
3              NaN
4              NaN

[5 rows x 27 columns]


texts: What happened to buildings on USGS topo maps?

texts:
I have been using the National Map website to download some topographic maps and I can't help noticing a big difference between the old and new. Historical USGS topographic maps have building outlines on them, but current ones do not. Current maps have things like schools and cemeteries, but no outlines depicting general built structures.

Is there a reason for this omission?

Is there a data source where building information is available?

My apologies if this is not the correct forum for these questions.


texts: Error keeps showing while I try to merge shapefiles

texts:
I am trying to merge shapefiles (I have multiple of them in one folder) but it keeps returning this error. I'm not sure what that means. I didn't change the path or the data.
import arcpy, os
arcpy.env.overwriteOutput=1
arcpy.env.workspace="C:/Users/data"

fc = arcpy.ListFeatureClasses()
for shapefile in fc:
    if arcpy.ListFeatureClasses("","Polygon"):
       arcpy.Merge_management(fc,"merge_poly.shp")
    elif arcpy.ListFeatureClasses("","Point"):
       arcpy.Merge_management(fc,"merge_point.shp")
    elif arcpy.ListFeatureClasses("","Line"):
       arcpy.Merge_management(fc,"merge_line.shp")
    else:
        print("cannot merge")

ERROR 002948: The shape type of C:/Users/data\obce.shp is not the same as that of previously entered datasets.
Failed to execute (Merge).


texts: Pass attribute to python script with QGIS command line launcher

texts:
I didn't find the answer to this question : is there a way for passing attributes to python script when you try to run it from QGIS command line launch :

os.system('qgis --nologo --project project.qgs --code code.py [arg1; arg2; ...]')

Because on official documentation there is no word about this possibility. Otherwise do you know a way for passing a layer path as attribute to this code.py ?


texts:
I am downloading polygons from one of our feature layers on ESRI Online. I am creating sqlgeometry polygons in a table in MSSQL server.  Then I am downloading some points that come from another feature layer an storing those points in another table in MSSQL server.  However when I try to find the polygon for a point, I keep coming back with a null value.  From what I have read, it usually seems to be a problem with the order of points used to make the polygon.  Here is the order I am using.

Can anyone confirm if I am going in the correct direction or not?

I have tried to create the polygons the two different directions, but to no avail.  Maybe if I know I am going counter clockwise, I can focus on something else.

my polygon is constructed this way

-9963016.632 5551661.8461
 -9963013.4085 5551537.4852
 -9962845.196 5551391.701
 -9962564.1822 5551377.1886
 -9962477.7671 5551462.9439
 -9962380.3512 5551669.254
 -9963016.632 5551661.8461


texts: Call FWTools (ogr2ogr) from ArcGIS python script

texts:
In ArcGIS, I want to use a Python script (and a standard ArcGIS script tool UI) to call functionality from FWTools.

How would I make the call to FWTools?  I tried using os.system(C:\WINDOWS\system32\cmd.exe /K "C:\Program Files\FWTools2.4.7\setfw.bat") which is what the windows shortcut calls, but that is not it.  It seems I am probably missing something rather simple, but I can't figure out how to make the right call through python to get FWTools & ogr2ogr operating on the string I pass.

My UI looks like this:

and the pseudo-code behind it is something like:

'ogr2ogr -f "' + outformat + '" "' + outfile + '" "' + infile + '" ' + ogrvars


texts:
Digitizing a utility network is feasible in QGis, but I face great difficulties to digitize the huge number of annotations and dimensioning attached to the network. I try to achieve this operation in the QGis by creating blank layer with labels but it is not convenient for dimensioning. Also, working with the composer is very tedious. Note that this operation is easy in CAD software (Autocad map for example).
Is there a plugin that could achieve easily this operation in QGis ?


texts: Joining tables results in null values within joined part using ArcGIS Pro

texts:
I have tried various geoprocessing tools, read through the stack exchange, googled questions, and I still have not figured this out.

I am trying to join a table based on GEOID field. I made sure both tables have the same text data types, same field name, and same data type length. However, every time I try to join the table it results in null values. I am not sure what is causing this.

Table 1:

Table 2:

What the join looks like:


texts: Raster Cell Size differences between ArcMap and R

texts:
I'm having issues when projecting habitat suitability rasters (WGS84) to a projected coordinate system (Africa Albers Equal Area Conic). The reason I reproject the rasters is so that I can calculate the area of suitable habitat which is intended to be used for conservation-related purposes.
Considering that I have about 900 individual rasters that need to be reprojected for area calculation purposes I decided to automate the process in R. However, I realised that there is a difference in 'cell size' values for the same rasters between R and ArcMap 10.7.
For example, when I use 'Raster A' (cell size: X: 0.00833333376786619; Y:0.00833333376786619) and project it to Africa Albers Equal Area Conic the cell sizes are different between R and ArcMap and therefore also the area of suitable habitat...
After reprojecting in ArcMap the cell size of the rasters are X: 879.789038737608
Y: 879.789038737608
But in R the 'cell size' is different: X: 853
Y: 902
The resampling method used was bilinear when reprojecting the raster.
This turned out to be a massive headache to me because I do not know which area calculated is the correct one as both R and ArcMap provide area estimates vastly different from one another. Area calculated from R = 18000km2 and in ArcMap 14000km2.


texts: How can I create a polygon from points and assign it a color in QGIS using python?

texts: I want to create a new vector layer, containing a polygon from a set of given points, and giving it a color.


texts: No module named "SSUtilities" within Optimized Hot Spot Analysis tool

texts:
Background

I'm trying to execute Optimized Hot Spot Analysis tool in ArcGIS 10.3 but receive the following error:

Additional Information

I've had a look at some other tools within the same toolbox and have noticed I'm also missing the "WeightsUtilities" module within the Cluster and Outlier Analysis tool. In this particular question, I am only interested in finding the "SSUtilities" module, but it would be grand if the method used to find it can be applied to other modules which may be missing in the future.

Question

Where and how can I find the missing "SSUtilities" module needed within the Optimized Hot Spot Analysis tool?


texts: How to get the contributors on OpenStreetMap

texts:
I have put my osm file in base thanks to osm2pgsql and I am looking for a way to get the contributor list. I have searched in my database and haven't found so much... and the osm wiki about the osm2pgsql schema does not contain a lot of details...

Is it possible?


texts: Copy/Paste Text box ends up way off screen in ArcMap

texts:
I'm trying to copy/paste some drawing text box elements in data view to preserve the formatting and for some reason it's pasting the element way out of my view area and way outside the scale of my drawing.

Is this a projection issue?


texts: I want to make velocity map in this format. I have points with vector attributes (u, v and velocity), but my arrows have only heads!


texts:
To get rid of "Mixed content" warnings I'm trying to access stamen map tiles via https. For http the used url is:

http://{s}.sm.mapstack.stamen.com/(toner-lite,$fff[difference],$fff[@23],$fff[hsl-saturation@20])/{z}/{x}/{y}.png

It's working fine with leaflet. After reading the notes on maps.stamen.com the protcol independet endpoint is http(s)://stamen-tiles-{s}.a.ssl.fastly.net/

I replaced the url by https://stamen-tiles-{s}.a.ssl.fastly.net/(toner-lite,$fff[difference],$fff[@23],$fff[hsl-saturation@20])/{z}/{x}/{y}.png
But this solution won't work. I got only file not found errors
After poking around https://stamen-tiles-{s}.a.ssl.fastly.net/toner-lite/{z}/{x}/{y}.png succeded but isn't the wanted sytle.

Any ideas how to solve this problem?


texts: Making button clickable to perform function in Leaflet?

texts:
I am trying to create buttons that will initiate a flyTo. I can click the buttons on the page but they do nothing

action functions:

<script>
      function Action(){
      map.flyto([43.372000,-80.987697], 13);
      }
</script>
<script>
      function Action1(){
      map.flyto([43.370627,-80.998735], 13);
      }
</script>
<script>
      function Action2(){
      map.flyto([43.377328,-80.963505], 13);
      }
</script>
<script>
      function Action3(){
      map.flyto([43.388931,-80.972062], 13);
      }
</script>

Buttons

<div class = "leaflet-bottom leaflet-left" style="pointer-events: auto">
 <button onClick="Action()" block>Douglas Flood Risk Area</button>
 <button onClick="Action1()" block>Avondale Flood Risk Area</button>
 <button onClick="Action2()" block>Millenium Park Flood Risk Area</button>
 <button onClick="Action3()" block>Mornington Flood Risk Area</button>
</div>

I have added 4 buttons into my leaflet map and none of them can actually be clicked.

Is there a way to make them clickable?


texts: ogr2ogr merge multiple shapefiles: What is the purpose of -nln tag?

texts:
The basic script in order to iterate recursively over sub-folders and merge all shapefiles into single one is:

#!/bin/bash
consolidated_file="./consolidated.shp"
for i in $(find . -name '*.shp'); do
    if [ ! -f "$consolidated_file" ]; then
        # first file - create the consolidated output file
        ogr2ogr -f "ESRI Shapefile" $consolidated_file $i
    else
        # update the output file with new file content
        ogr2ogr -f "ESRI Shapefile" -update -append $consolidated_file $i
    fi
done

Hoverer in vertaully all examples around the web I noticed that for the case where I update the output file, -nln tag is added, for example:

ogr2ogr -f "ESRI Shapefile" -update -append $consolidated_file $i -nln merged

According to the documentation it says:

  Assign an alternate name to the new layer

And I noticed it creates a temporary shapefile called "merged", and in the end of the loop the file is identical to the last shapefile I merged.

I don't understand why I need this? Because I succeeded to merge successfully without this tag.


texts: Extract rows and columns from a raster in Google Earth Engine

texts:
I am trying to extract column and row information from a raster in Google Earth Engine, but am having some trouble doing so. Here is sample code.

var ROI = ee.Geometry.Rectangle(-82.56277, 35.58935,-82.53436, 35.59996); // define area

var collection = ee.ImageCollection('LANDSAT/LC08/C01/T1')
    .filterDate('2014-05-01', '2015-10-01');

var composite = ee.Algorithms.Landsat.simpleComposite(collection);
var image = composite.clip(ROI);

Map.centerObject(ROI);
Map.addLayer(image,{},"image");
print(image.arrayLength(image,0));

However, image.arrayLength(image,0) returns the error message

Too many (3) arguments to function: Image.arrayLength(input, axis)

Returns the length of each pixel's array along the given axis.

Args:

  input (Image): Input image.
  axis (Integer): The axis along which to take the length.

I have tried to enter image.arrayLength(0) as well, but that returns the whole image. According to the arrayLength help page, it seems like it should return an integer with the number of pixels in the specified axis.

The information is contained if you drill down into the image object (Image->bands->0->dimensions->0 or Image->bands->0->dimensions->1), and the object appears to be a json object, but I can't figure out how to access this information.


texts: Creating a new integer field based on text in a different field - ArcGIS 10.1

texts:
I need to create a new integer field based on another text field. All instances of a the same text string in the existing field should be assigned the same number (starting at number 1). Can I do this in field calculator? If so, how?

I'm a beginner code-wise.

EDIT: Here's a sample

Each place name should be assigned a number. The third column shows what I'd like as a result. This is so I can use Hawths tool BATCH FIXED KERNEL DENSITY ESTIMATOR which only works with number categories:

http://www.spatialecology.com/htools/bkde.php

NOTE: Although I'm using 10.1, I've tried the non-batch kernel density estimator and it works well, no hitches so expect / hope the batch kernel estimator will work too.


texts: What is the raster merge criteria algorithm on rasterio.merge?

texts:
I'm creating a mosaic from many sentinel images with a bunch of overlapping areas, although I'm trying to understand what is criteria. For example:

This is a small piece totally clear of clouds

And this is a larger piece with many clouds

This is the result of the merge using rasterio and gdal
I used rasterio.merge and gdal_merge. Both brings me the same result, but don't understand the algorithm method yet. How can I tell the tool to preserve the image with fewer clouds as possible? Which method are these algorithms using to preserve the larger raster and not the cleanest (no clouds)?
It would be absolutely awesome if I could select the raster with fewer clouds to be merged for this mosaic. Is there a way to choose the raster I want to keep?


texts: ArcGIS historical traffic data with Impedance attribute of miles

texts:
I'm trying to solve a route using an impedance of miles (to get the shortest route) through ArcGIS Sever 10.3.1 rest API.

However, I would also like to get the travel time of this route based on historical traffic data for a specific date/time. It seems like ArcGIS only supports historical traffic data if the impedance is TravelTime.

The message returned by the server is:

type: 50,
description: "The start time was ignored because the impedance attribute is not time-based."

Is it even possible to use traffic data on non time-based routes?

EDIT:

I did manage to find this in some documentation:

  To use traffic in the analysis, set impedanceAttributeName to TravelTime, and assign a value to startTime...
  ...Traffic is supported only with the driving time impedance or travel mode.

Which seems to indicate this isn't possible, but that surprises me. I would have thought this was not an uncommon request. Does anyone know of a workaround or alternative method?


texts: Don't allow to save edits in the table untill the field will be populated in ArcGIS with pythonaddins

texts:
I have a table and some fields there(OBJECTID, UN_ID, for example).Table is not empty. how to prevent an edit session from stopping or saving until  the field will be filled.
I wrote an extension onStopEditing(). It displays a message that the field is empty.

And how to make the module does not complete the editing?

Any suggestions as I am new to this.

import arcpy
import pythonaddins

class Uniq_ID(object):
    def __init__(self):
        self.enabled=True
    def onStopEditing(self, save_changes):
        mxd = arcpy.mapping.MapDocument('current')
        with arcpy.da.SearchCursor(table.name,['OBJECTID','UN_ID'],"UN_ID IS NULL") as cursor:
             for row in cursor:
                 pythonaddins.MessageBox("UN_ID cannot be NULL","Error Message")

        with arcpy.da.SearchCursor(table.name,['OBJECTID','UN_ID'],"len(UN_ID) <>9") as cursor:
             for row in cursor:
                 pythonaddins.MessageBox("UN_ID must have 9 symbols","Error Message")
    del mxd
    pass


texts: When exporting, mask behind text becomes garbled (Data Driven Pages)

texts:
I am using Data Driven Pages with different rotations for a project. When I export to PDF or PNG, the mask behind the text moves and rotates. It does not seem to rotate in the direction of the data frame. We have tried MXD Dr, Repair Geometry tool, and re-creating the MXD. I am able to export fine until I insert the data driven page index. We also created a new data driven page index.shp and had the same problem.

What it should look like:

9/25: Tried exporting again; masks w label engine (in red) did not show up, Fischer Rd (annotation) is still rotated/ wonky.


texts: How to customize logic for defining default speed limits for roads in osm2po

texts:
I need to customize osm2po logic in defining possible speed limits of ways. osm2po takes into account only class of road while I need to consider one more parameter: is road inside or outside some city, i.e. for example: speed limit for higway=trunk inside city = 60km/h while outside it should be equal to 90km/h. What is a proper way to implement such logic? Is it feasible at all to do withing osm2po? Or may be it better/easier to load graph data into Postgres/PostGIS using osm2po generated sql-script and do speed/cost corrections there (assuming that we already have all city borders loaded into Postgres/PostGIS db) and use pgRouting as routing engine?
Any advices highly appreciated.


texts: Are there any tools / plugins to develop mapbooks using QGIS?


texts: Seeking urban area shapefiles for France, Germany and Belgium?

texts:
I'm trying to make a map like this for the Low Countries/Rhineland, where I basically show how spread out all the urban areas would be there if the urban areas had similar densities to American urban areas. I have all the data I need for the UK and for the Netherlands, but I can't find all the data I need for France or Germany or, in particular, Belgium. For each of those countries, I'm looking for:

a shapefile of the boundaries of all the built-up/urban areas in that country (as opposed to metropolitan areas or administrative city boundaries), which includes data of the population and area (or density) of each of those built-up areas
a shapefile of the boundaries of the smallest subdivisions possible for that country (akin to Local Super Output Areas for the UK or Buurten for the Netherlands), which includes data of the population and area (or density) of each of those subdivisions, and preferably also the population centroids of each of the subdivisions. I've found a 1km^2 grid for Germany, but not for the other two countries.

Does anyone know where I could get these shapefiles?

I would prefer shapefiles based on 2011 census data. I've looked on Eurostat, but they don't have built-up area shapefiles, and they only do detailed shapefiles for the whole of the EU put together, which my computer isn't powerful enough to cope with.


texts: Adjusting sampling points to closest land surface QGIS

texts:
I have many sampling points in mangroves of Southeast Asia and am trying to extract climatic data for each using the WorldClim dataset. However, many of the points fall on areas classified as "water" for which the WorldClim dataset values do not exist.

How do I adjust each of the point long/lat values to the closest land value in QGIS? Alternatively, if there's a way of doing so in Program R, I'd also be interested in that.


texts: Fix an export for loop that recently stopped working?

texts:
I have created several for loops in my GEE script.  These loops worked fine up until a week ago (about 4 months old).  I receive the following error in my console when the code tries to run the for loop to export the rasters as GeoTIFFs:

For loop starts on Line 219
Array.prototype.indexOf called on null or undefined
https://code.earthengine.google.com/75c7c2b933b631e8a61ab4a47d39fabe

I am very new to coding so I don't know too enough to really understand and implement the information in the GEE guide about Server side vs. Client side looping.

Can someone look at my export code and help me figure out how to write new code to export these rasters (1988 to 2017 time series)?


texts: Safe default value for missing street address number when reverse geocoding?

texts:
Context

We are building an app that reverse geocodes a user location to retrieve a street address. The street address is then used to request transportation pick-up at the person's location.

Problem

The transportation company needs a street address including street address number (the '123' in '123 Main St.') in order to dispatch a ride. However, the street number is not available in 100% of reverse geocoding results.

Question

What are some options that might be safe to consider for default street number value when no data is available?


texts: How to create a sequence for a table which was created by selecting data from another table?

texts:
I have created a new table by selecting data from an existing table:

CREATE TABLE polygon2 AS
SELECT gid, geom
FROM polygon1;

CREATE SEQUENCE polygon2_gid_seq;
SELECT SETVAL('polygon2_gid_seq', (SELECT MAX(gid) FROM polygon2));

ALTER TABLE polygon2
ALTER COLUMN gid SET NOT NULL,
ALTER COLUMN gid SET DEFAULT nextval('polygon2_gid_seq'::regclass),
ADD CONSTRAINT polygon2_pkey PRIMARY KEY (gid);

Unfortunately neither QGIS nor pgAdmin are dropping the corresponding sequence when dropping the table. This issue doesn't appear when dropping tables which were created using the QGIS 'Create table' dialog. So I suppose there is something missing in the way I've created the sequence. Can anybody please help me?


texts: Populate Alternate ID Fields in Turn Feature Class in ArcGIS 10.5

texts:
I have a Network Dataset with turns and need to edit some of the edges (I found out some were not planarized yet). Using the Populate Alternate ID Fields tool I keep getting the Error 030030.

The Editor is off and no feature class is in use. I get the error both in ArcMap and ArcCatalog.


texts:
I attempting to make flow graph in QGIS 3.4, some thing like

I have found several 'ways to do it' in QGIS 2.xx eg QGIS2.0 method

I known its possible to get a line with arrows by just change the line type.

But I am not sure how to get a curved line

The Arrow and Flow mapper plugin only seems to be a QGIS 2.xxx feature.

Any suggestions?


texts: Labeling lines that are not within polygons using field calculator

texts:
I have two features. One line feature and another polygon feature. I want the line feature label not to appear when it is within a polygon similar to the one in the figure. That is, I want the label to only appear when the line is outside the polygon. Is it possible to do this in QGIS 3?

I tried the field calculator with the within function but it didn't work very well.

Below is the code I tried to do in the QGIS field calculator:

CASE
   WHEN within('Lines_4e1a1089_ba02_495e_8c79_d32b79cf3b54','polygon_406a7fae_06b4_45a1_a20e_d64b8682e501') THEN " "
   ELSE "Label"
END


texts: How to search vector layer features in Mapbox GL JS

texts:
I need to add a search bar that will search vector layer attribute features which are hosted by mapbox. I'm using the Mapbox GL JS API. I know this can be done fairly easily with Mapbox js/Leaflet. I'm totally stuck as to how to do the same in Mapbox GL JS. I want to search for a zip code, which corresponds to a feature attribute in my vector layer (please note, I don't wish to geocode an address or zip code. I am aware that there is a geocoder plugin for Mapbox GL. This is not what I'm looking for). I'm then hoping to update a html info panel with various feature values from the vector layer based on the matching search entry.

I need to know how to marry search functionality with the vector layer features served by Mapbox (This question specifically relates to Mapbox GL, not the older Mapbox JS library).


texts: How to write a point-in-polygon query for Postgis?

texts:
I'm trying to write a query in Postgis which will return all rows where a dynamically changing longitude/latitude point resides within the geometric polygon of that row (in the geom column). This is the query I'm trying to use:

SELECT * FROM "geospatial-data"
WHERE ST_Contains(
    "geospatial-data".geom,
    ST_Transform(
        ST_GeomFromText('POINT(-90.56210 38.60369)', 4326), 4269
    )
)=true

I'm basing this query off of an older post here: http://codemagician.wordpress.com/2009/09/12/point-in-polygon-in-postgis/

I'm not getting any errors but I can't understand what I'm doing wrong. How should I write/modify the query in order to query for a point-in-polygon?


texts: Can I reproduce this classic cartographic hillshade relief effect (hachures) electronically using modern DEM data?

texts:
I own a number of old cetral-European maps and I always liked the classic hand-drawn cartographic hillshade relief presentation that these maps feature. I am not wanting to reproduce the shading (which is similar to swiss shading) but the accompanying lines. These are not contours but rather down-slope "gradient" lines that (often with a degree of imagination) create a hillshade effect. I don't have my maps on me right now so I can't scan a good detailed example of this but I found some examples online.  As far as I know, and at least in Czechoslovakia, these types of maps were being produced up until the second world war. I would like to produce some modern maps with an old feel and featuring this type of effect.  Is there any way that this could be done in a modern GIS (Preferably ArcGIS) using an existing DEM?  Illustrator or Photoshop solution would also be acceptable.

Not a part of the question but I would also be interested to find out:
What is this particular hillshading called?
Where could I get some more information about how it was done and what rules it adheres to?


texts: Remove decimal points from string field using field calculator

texts:
I have a string field populated with shape lengths in my attribute table but they have 12 decimal points.
How do you remove decimal points from a string field populated with numbers using the field calculator?


texts:
I am trying to order one my  features based on some of the attributes from ascending to descending. It works for some of my attribute columns but not others. Are there any requirements for addOrderBy to work?

for feature in layer.getFeatures(QgsFeatureRequest().addOrderBy('Length ID',True)):

The above doesn't seem to work, but

for feature in layer.getFeatures(QgsFeatureRequest().addOrderBy('Circuit',True)):

the above works. I have tried to change my Length ID as float, int, and string, but nothing has worked.


texts: Transform SPOT5 images to natural color images

texts: We're currently using PixelSense to convert raw SPOT5 images to natural color ones for display and printing. Do you know a good algorithm to perform the same task as we want to ingest raw images in an image server and provide an image service with on the fly natural color conversion.


texts: Is it possible to add an atlas variable to the SQL query when creating a virtual layer?

texts:
I'm trying to do something like this:

select *
from table_a
inner join table_b
  on table_b.fid = table_a.id
  and table_b.code =  @atlas_featureid

I have tried this using QGIS 2.18.16 and QGIS 3.0 but it doesn't seem to work. Am I missing something? Is there another approach to achieve this?

The alternative I could figure is to add a virtual layer without the last condition and symbolize just the ones that match the @atlas_featureid but this layer has about 100 times more features (≈ 300k) and it takes ages to work with it, it's not a viable solution.


texts:
I have a set of functions I am using to read shapefiles and produce some geojson files for parsing:

def read_shapefile(path, buffer = []):
    # read the shapefile
    reader = shapefile.Reader(path)
    fields = reader.fields[1:]
    field_names = [field[0] for field in fields]
    #buffer = []
    for sr in reader.shapeRecords():
        for i in range(len(sr.record)):
            if type(sr.record[i]) == bytes:
                sr.record[i] = str(sr.record[i])
        atr = dict(zip(field_names, sr.record))
        geom = sr.shape.__geo_interface__
        buffer.append(dict(type="Feature", \
        geometry=geom, properties=atr))
    return buffer

def write_geojson(buffer):
    # write the GeoJSON file
    #print(buffer)
    print(type(buffer))
    geojson = open("pyshp-demo.json", "w")
    geojson.write(dumps({"type": "FeatureCollection",\
    "features": buffer}, indent=2) + "\n")
    geojson.close()

def read_shapes():
    """ read all of the shapes in our shape directory and write them to geojson"""
    files = get_file_list(SHAPE_DIRECTORY, [".shp"])
    print(files)
    buffer = []
    for file in files:
        buffer = read_shapefile(file, buffer)
    write_geojson(buffer)

which works fine but sometimes I am getting files in X\Y instead of lat long. I see How to convert projected coordinates to lat/lon using Python? which is useful for converting single points but is there any way inside read_shapefile I could easily apply the conversion to the whole geometry object instead of parsing through my end result line by line?


texts:
I would like to transform KNMI MSGCPP geoview netCDF file to WGS84 netCDF file with gdal_translation and gdalwarp. I got error msg.
I shared the geoview netCDF file (140 MB):  https://z001.s3.amazonaws.com/test.nc
szabot@XSRV:/mnt/c/gdal/bin/gdal/apps$ gdal_translate -a_srs WGS84 NETCDF:"test.nc":cldmask test001.nc
ERROR 5: OSRCalcInvFlattening(): Wrong input values
Input file size is 3712, 3712
Warning 1: Several drivers matching nc extension. Using NETCDF
0...10ERROR 5: OSRCalcInvFlattening(): Wrong input values
ERROR 5: OSRCalcInvFlattening(): Wrong input values
ERROR 5: OSRCalcInvFlattening(): Wrong input values
ERROR 5: OSRCalcInvFlattening(): Wrong input values
...20...30...40...50...60...70...80...90...100 - done.


texts: IBufferConstruction create buffer with different geometry types (C# ArcObjects)

texts:
I am using the following method in C#

public void ConstructBuffers (
    IEnumGeometry pInputGeometries,
    double distance,
    IGeometry CollectionpOutputBuffers);

to create buffers.

In the method above if I have a IEnumGeometry with same types of multiple geometries I can dissolve the overlaps to get a single geometry. But when I have a IEnumGeometry with different geometry types (points, lines and polygons) and try to dissolve the overlaps I don't get a single geometry but my output (CollectionpOutputBuffers) is the same number of geometries as input. Anyone has any idea what may be wrong.


texts: How to call gdal_retile from Python code?

texts:
I am trying to call gdal_retile from my Python code like this:

import gdal_retile

and in main:

gdal_retile.main("-v -r bilinear -levels 4 -ps 2048 2048 -co \"tiled=YES\" -targetDir pyramid --optfile files.txt")

but I get this error in my Debug Probe window:

Unrecognised command option: -
Usage: gdal_retile.py
    [-v] [-co NAME=VALUE]* [-of out_format]
    [-ps pixelWidth pixelHeight]
    [-ot  {Byte/Int16/UInt16/UInt32/Int32/Float32/Float64/
           CInt16/CInt32/CFloat32/CFloat64}]
    [ -tileIndex tileIndexName [-tileIndexField fieldName]]
    [ -csv fileName [-csvDelim delimiter]]
    [-s_srs srs_def]  [-pyramidOnly] -levels numberoflevels
    [-r {near/bilinear/cubic/cubicspline/lanczos}]
    [-useDirForEachRow]
    -targetDir TileDirectory input_files

Any idea what I am doing wrong?


texts: Python crashes when I try to import GDAL after re-installing it

texts:
A while back I went through the installation process for GDAL by installing the core components and then the Python/GDAL bindings.  Everything was fine when I would type in from osgeo import gdal.  However, I had to uninstall and reinstall everything due to a problem I was having.  Now when I type in from osgeo import gdal, Python crashes and I get the message that it has stopped working.  When I click on "Close the Program" Python restarts.

I've looked to see if there was anything left over from the previous install, but couldn't find anything.

Any ideas on something I could try?


texts: Selecting all polygons where neighbor vertex points are below distance threshold in ArcPy?

texts:
I have a shapefile with over 100k polygons (some of them have holes in them but apart from that they are simple polygons - no multipart).  I want to loop through all of these polygons and select only the ones where at least two neighbor vertices (within the same polygon) are within a distance threshold.  Another way of putting it is:  if two neighbor vertices (ex:  id #3 and id #4) within a polygon have a distance < threshold_distance, then select them for further processing.  Any idea how to do this efficiently?

I am a beginner with arcpy (python). Conceptually, the brute force approach I see is:

Use a search cursor to loop through all polygons
For each polygon, use "getPart" to create an array of points
Use a for loop to cycle through all points and calculate distance between
    point i and point i+1 (except for last point). If the distance is less than distance_threshold, then stop the loop, copy the polygon and insert it into a new shp file. (or maybe just copy the polygon's unique identifier and save that into a txt file)

Since I need to cycle through + 100k polygons, is there a more efficient way of doing things?


texts: getDestinationGeographicPoint() in geotools returns wrong coordinates in South/West

texts:
I am trying to calculate a coordinate given an azimuth, distance and start point.
If my code uses an start coordinate in North/East quadrant, all goes fine.
But, by using 0,0 coordinates, an azimuth -90 or 180 give me an inverted result movement.

Here is my sample code:

CoordinateReferenceSystem crs = CRS.decode("EPSG:4326");
GeodeticCalculator calc = new GeodeticCalculator( crs );

calc.setStartingGeographicPoint(0, 0); // Starting from Lon/Lat 0,0
calc.setDirection( -90, 2000 );

Point2D dest = calc.getDestinationGeographicPoint();

System.out.println(calc);
System.out.println( dest.getX() + " " + dest.getY() );

This is the output:

Coordinate Reference System:  WGS 84
Ellipsoid:                    WGS 84
Source point:                 00°00,0'E 00°00,0'N
Target point:                 00°01,0'W 00°00,0'N
Azimuth:                      -90°00,0'
Orthodromic distance:         2000 m

-0.017966305682390427 1.107533167043945E-18

I'm using this site http://earthexplorer.usgs.gov/ to show the points.
As you can see, a movement in azimuth -90 ( to west ) results in West coordinate of target be grow (Target point: 00°01,0'W 00°00,0'N) but the numeric representation (degrees) is not correct ( the second argument is the latitude and MUST be ZERO and not 1.107533167043945E-18, since my movement was to WEST).

What I'm doing wrong ?


texts: Integrating OpenLayers.Layer.ArcGIS93Rest layer on OpenLayers 2 map where base map uses different projection

texts:
OpenLayers 2 has a layer subclass called OpenLayers.Layer.ArcGIS93Rest which is great but it is passing BBOX parameter in Google Projection (EPSG:900913) which is what my base layer uses, instead of transforming it from the map projection to the layer's projection (in this case 102656).  Even if I change the layer's projection, the URL generated tells the server of its correct projection, but the BBOX is NOT in this projection :-(

I have no control over the projection of the server maps.  So I cannot change the server.

Here is an example - it is using the CT State Plane projection...

    var layerB = new OpenLayers.Layer.ArcGIS93Rest("My new layer",
        myMapserverExportURL, //Note this ends with 'MapServer/export'
        {
            layers: "show:6,7,9",
            srs: 'EPSG:102656',
            TRANSPARENT: "true"
        },
        {
            isBaseLayer: false,
            visibility: false,
            resolutions: this.getScales(6, 11)
        });
    layerB.projection = new OpenLayers.Projection("EPSG:102656");
    map.addLayer(layerB);


texts: Adding layers with subtypes and feature templates using ArcPy?

texts:
I am trying to clip a feature from one geodatabase into another and retain the editing feature templates I have assigned using subtypes for the original feature class in the original geodatabase. After clipping, the subtypes and editing feature templates are retained, but only when I drag and drop the layers into the map document. When I use arcpy.mapping.AddLayer(), the editing templates are lost.

How do I use python to add layers such that the editing feature templates are retained?

The clipped feature is named "Feature_Clip" in the code below.

workspace = arcpy.env.workspace
clippedFeature = os.path.join(workspace, "Feature_Clip")

mxd = arcpy.mapping.MapDocument("CURRENT")
df = mxd.activeDataFrame
layer = arcpy.mapping.Layer(clippedFeature)
arcpy.mapping.AddLayer (df, layer)
arcpy.ApplySymbologyFromLayer_management(layer.name, layerFile)
del mxd, df, layer

note layerFile points to a .lyr file I have created from the original feature class in hopes of using it's symbology, but it doesn't seem to have any effect.


texts: Releasing features from a recycling feature cursor

texts:
I have a recycling feature cursor that is used to fetch a single row. Is it my responsibility to release that returned feature, or will it be cleaned up by the cursor as it would be if I consumed the entire cursor?

public object GetFirstValue(IFeatureClass featureClass, string field)
{
    IFeatureCursor cursor;
    IFeature feature;
    try
    {
        cursor = featureClass.Search(null, true);
        feature = cursor.NextFeature();
        return feature.get_Value(feature.Fields.FindField(field));
    }
    finally
    {
        if(cursor != null) Marshal.ReleaseComObject(cursor);
        //do I need to release the feature here?
    }
}


texts:
I'm new in Geographic Information Systems and programming.

I have some points referencing Barcelona, Spain. These points are comming from City Hall's opendata site. They say they are in ETRS89 reference system. I don't really know what it means. I need them in Latitude Longitude system so I can pass them to CARTO, OpenStreetMaps, GoogleMaps...

Here some points example:

430344431;4582166733
430344431;4582261851
430344431;4582060238
430344431;4582301912
430344431;4582201841
430344431;4581863161
430344431;4580885282
430344431;4582098182
430344431;4582313751
430344431;4581900735
430344431;4581281417


texts: Calculating centroids without adding field to table using ArcPy?

texts:
I haven't been able to find any bit of code that does this successfully.
I need to delete some fields in a table that I'm modifying using a script I'm writing. I don't want the fields to be left on the table, but I'm pretty sure that I need the fields to calculate the centroids. I haven't found a way of calculating parcel centroids that works without adding a field to a table. I want the xCentroid and yCentroid fields to be calculated, and I want to delete those fields at the end of the script. Here the code I'm using.

 try:
                # Set environment settings
                arcpy.env.workspace = "U:/JackBuildingFootprints.gdb"

                # Set local variables
                inFeatures = "parcelsCopy"
                xCentroid = "xCentroid"
                yCentroid = "yCentroid"
                fieldPrecision = 18
                fieldScale = 11

                # Add fields Remember to delete fields later. Get help with that
                arcpy.AddField_management(inFeatures, xCentroid, "DOUBLE",
                                          fieldPrecision, fieldScale)
                arcpy.AddField_management(inFeatures, yCentroid, "DOUBLE",
                                          fieldPrecision, fieldScale)

                # Calculate centroid
                arcpy.CalculateField_management(inFeatures, xCentroid,
                                                "!SHAPE.CENTROID.X!",
                                                "PYTHON_9.3")
                arcpy.CalculateField_management(inFeatures, yCentroid,
                                                "!SHAPE.CENTROID.Y!",
                                                "PYTHON_9.3")

except Exception:
                e = sys.exc_info()[1]
                print(e.args[0])


texts: How can I launch a QWidget after pressing a submenu (QAction) in QGis3.4?

texts:
I need your help
I am developing a dialog box (which I will call toolraster) in QGis to interact with raster layers. I want it to display when I click on a menu I insert in the QGis interface.
I created Toolraster as a class that extends QtWidgets, a QWidget.
Create the QAction and connect it to the toolraster instance function. Create the submenu with iface.pluginMenu().addSeparator() and iface.pluginMenu().addAction(activar). Example below.
The problem: when I run the menu the dialog box opens and closes instantly. It doesn't report any errors
When I perform the same procedure with a QInputDialog or a QMessageBox it doesn't present any problems.
But with a class like this it is a problem:
from qgis import PyQt
from qgis.PyQt import QtWidgets

class test(QWidget):
    def __init__(self):
        super().__init__()
        self.setWindowTitle("clase de prueba")
        layout = QVBoxLayout()
        widgets = [QCheckBox,
            QComboBox,
            QDateEdit,
            QDateTimeEdit,
            QDial,
            QDoubleSpinBox]
        for w in widgets:
            layout.addWidget(w())

        self.setLayout(layout)

def activar():
    h=test()
    h.show()

test_action=QAction("Menu de prueba")
test_action.triggered.connect(activar)
iface.pluginMenu().addSeparator()
iface.pluginMenu().addAction(test_action)

Note: when I execute the function activar() in the console the same thing happens, if I insert the show() inside the class too


texts: Consuming Forecast.io Radar Data in an API or ArcMap

texts:
I want to consume forecast.io's radar data as an ArcGIS Server REST or JSON to use in an API or ESRI's ArcMap GIS Software.  MapBox is doing it here. I don't think forecast.io is currently making the radar data service available as REST or JSON, it looks like they are using TileMill to serve out 16-bit grayscale GeoTIFF images of the radar data and applying some MapBox magic.  I want to do the same as a service for ArcGIS Server 10.2. I have looked into Cutting Tiles for ArcGIS Server Using TileMill but have not tested.
It's rumored that ESRI will build functionality into server to be able to consume MBtiles and that may be a potential solution, but for now I am looking for a workaround.
Forecast for Developers


texts: Comparing Google Earth Enterprise to ArcGIS for Server and GeoServer?

texts:
I'm working on a web app that includes mapping capabilities (currently just markers and KML overlays on an embedded Google Map). We're starting to outgrow the limitations of Google's free stuff, and Google Earth Enterprise initially seemed to be the natural upgrade path. But at least judging from the tags on this site, it doesn't seem like it gets a lot of use in the GIS community.

Has anyone used GEE who could compare it to the more widely used server platforms?


texts: Overwriting features without first removing their layer using ArcPy

texts:
I'm currently using arcpy.env.overwriteOutput = True because I want to be able to play around with variables and not have to delete the file from my gdb every time. But whenever I rerun the script with a changed variable it runs fully but then my changes were not implemented even if I drag the feature back onto my map.
I found the only way to get it to overwrite is to delete the layer from the dataframe first (see image) then run script and then it appears with my changes.

Is there a way around this?


texts: Can ArcGIS Operations Dashboard manage routine data collection at same geographic points?

texts:
I have started working an organisation that uses ArcGIS to manage its GIS data generally.  They use a wide range of other software and methods to also manage collection of routine data, such as:

Water Meter readings - the meter locations are stored in GIS, but the usage data is stored elsewhere
Water quality sampling - there are a range of locations where the samples are collected weekly/monthly/quarterly for laboratory analysis and stored in an excel spreadsheet.

The current system is fragmented and I thought that ArcGIS could possibly manage this better by consolidating the data.  I am a novice in ArcGIS, but thought that possibly each "round" could be stored as a new feature, with geographic location linking results and the attribute table recording the specific round data (eg date, time, analysis results, etc).  In addition, I thought GIS system may enable capturing the spatial variation of results - for example areas of high water usage.

Is this a common use of ArcGIS?  I understand that ArcGIS can create "dashboards" and would this be a suitable delivery method?


texts:
this is a pretty noob question but how do you save a project without losing the layers? I was doing a hillshade analysis and saved the file as .qgz but only loads the .tif file. There is a dialog box that pops up and says "bad layers."

Using:

Macbook Air Sierra
Python 3.6
QGIS MacOS 3.2

How can I save the project when I'm not done yet?


texts:
Is there a way to select features in a layer that have an attachment?

I want to identify all features on a map that have an attachment associated with it


texts: Join a table to a pre-joined field in QGIS

texts:
I'm trying to join a .csv (X) to a .shp layer (Y), but X and Y don't share an ID field.

Another .csv, Z, has TWO ID fields--one shared with X and one shared with Y.

In any other language, I'd first join Z to Y, then join X to Z/\Y--but after I join Z, the fields of Z are unavailable as "target layers" in the join dialog.

Is there any way around this without editing the underlying .shp files?


texts: Adding feature to remote geodatabase takes a long time

texts:
I have a program that I must add some features to a feature class that is in a sql server remote database. But when I want to open connection and add them, it takes a long time, about one or two minutes and it is not good. I tried to add them manually using SQL statements and I cannot do that if I enable geodatabase on my database. So if I want to add features manually, I should not enable geodatabase. On the other hand, I want to edit shapes using arcobject so I must enable geodatabase.

Here is the question:

If I dont enable geodatabase, I can load and add features very fast and it is great but I cannot edit them using arcobject.
If I enable geodatabase I can edit features but it is too slow.

What are your suggestions? Is there an alternative way? Are there any libraries that I can use to edit shapes directly from none geodatabase?


texts:
I'm not a Linux person; I'm finding it very difficult to install CartoDB 2.0 on my virtual machine.
I believe there was a V1 AWS VM that was created at some point.
Are there any VMs that are available for experimentation using version 2.0 of CartoDB?


texts: How can i add coordinates of Polygon nodes on a Polygon data table?

texts: I'm new in this domain and I want to know if I can add to data table of a Polygon.shp, coordinates of node. I have this table pkuid/ name_poly/ area/ lenght/ (here nr_point ex:1,2,3,4 it' about order/ x or long/ y or lat/ and other things about poly. I need in this format because I want a simple format for daily use by people who doesn't know how to manipulate layers of QGIS. At this moment I have already 2 layers one with points (coordinates) and one with polygons.


texts: Selecting specific counties using selectlayerbyattribute in arcpy

texts:
I have a shapefile for counties containing many fields, among which Pop_Descrp is one.

I need to select the counties with Highly Increased attribute under the field and export it.

I am not able to make correct expression for query. Can anyone help me with this?

import arcpy
from arcpy import env
env.workspace=r"Z:\Ash Tree Project\Shapefiles_Arkansas"
env.OverwriteOutput=True

arcpy.MakeFeatureLayer_management("County_AR.shp","County_layer")
arcpy.SelectLayerByAttribute_management("County_layer", "NEW_SELECTION", "[Pop_Descrp]='Highly Increased'" )
arcpy.CopyFeatures_management("County_layer", "HighPopR_counties.shp")


texts: Difference between planet_osm_polygon and planet_osm_point tables

texts: What's the difference between the planet_osm_polygon and planet_osm_point tables? I find similar columns in both tables.


texts: Why is file geodatabase so large when feature classes inside do not add up to anywhere near as much?

texts:
I recently noticed when publishing maps for use on ArcReader from ArcMap that when I create the map package my file size is now 18GB when it used to be around 800MB. I have not added any new feature classes to my file geodatabase that I am using to create feature classes for the map packages. I have no idea why there is this huge spike in data. When I create a map package it creates a pmf folder where the actual map resides and a data folder where all the geodatabase info is being stored. There are three cw folders, cw00, cw10, cw20, and under cw00 there is a file name a00000075 and it says it is a GDBTABLE File that has 16.6 GB in it but I have no idea how to tell what feature class it is coming from.

Here attached is what the file looks like


texts: I'm using QGIS 2.6 and regularly update a shapefile layer to indicate felling progress. However, if I edit the nodes of a polygon, save the layer out then look at the Attributes table the area shown doesn't change.


texts: Calculate path which follows a given grade or grade/slope relation

texts:
How is it possible to draw a path, which follows a given grade or a trail grade / slope relation?

I know least cost path where I was able to calculate least cost path using the slope and avoid steep areas (or steep fall line trails), as far as I understood it. That's what I did for now.

What I need is a trail, which is either not steeper than 10 %, independent from the slope, or whose trail grade is not more than half of the slope.

I try to set up a method for a bike trail suitability analysis based on multiple indicators. These two criteria are conditions for sustainable trails regarding erosion and water runoff. The trail would have a start and end point or several possible start and end points.

If I was not clear: Trail grade is the steepness of the trail, slope the steepness of the hill.

I'd also be glad with a corridor instead of one specific narrow path.

Further, I'd prefer a QGIS solution, but could also work with ArcGIS.


texts: I have setup a Haicom bluetooth gps on a machine running Windows8 and have used GPSGate to create some virtual COM ports. I have valid data arriving at the COM port, but then I run QGIS and try to connect to the specified COM port and each time QGIS crashes. Can anyone suggest what is going on please?


texts: Adding labels to raster with levelplot and spplot

texts:
I am trying to add district labels for a country, which I am plotting using levelplot (to use worldclim precipitation data) and adding boundaries using spplot.

I am having trouble adding text labels of the districts onto the map. I have the following code:

tan <- getData("GADM", country="Tanzania", level=1)
alt_world <-getData('worldclim', var='alt', res=10)

tanalt.sub <- crop(alt_world, extent(tan))
tan.sub <- mask(tanalt.sub, tan)

p1 <- levelplot(tan.sub, col.regions=colorRampPalette(brewer.pal(9, 'Blues')),margin=FALSE,main ='Tanzania rainfall (mm/year)')
p2 <- spplot(tan)

p <- p1+p2

I have tried the following to add labels:

labels <- layer(sp.text(coordinates(tan), txt = tan$NAME_1, pos = 1))
p2 <- p2 + labels

But get an error regarding 'inherited function for coordinates'.


texts:
I'm trying to generate shapefiles from png image and isolate polygons by colors. I'm using gdal_translate to georeferencing image to GTiff and gdal_translate.py to convert GTiff to shapefile. Converting the raw image,  the output shapefile is correct, but isolating one color, using imagemagick, the output shp have a colored background.

Is there any tool that would help me in this task ? And why converting to transparent is breaking my shapefile?

PS: Sorry, but I could not post the shapefile image without isolating the colors because I do not have enough reputation to post more than 2 images.

Raw file:

Shapefile with one color isolate:

Scripts i'm using:

      #Remove o yellow E red (d1)
    /usr/bin/convert $RAWDIR/$IMG -fuzz 2%  -transparent "#6A6A00" -transparent "#6A0000"  $COLORDIR/$D1

    #Remove green e red  (d2)
    /usr/bin/convert $RAWDIR/$IMG -fuzz 2%  -transparent "#008200" -transparent "#6A0000" $COLORDIR/$D2

    #Remove green E yellow (d3)
    /usr/bin/convert $RAWDIR/$IMG -fuzz 2%  -transparent "#008200" -transparent "$#6A6A00" $COLORDIR/$D3

    gdal_translate -of GTiff -a_ullr -82 12.04 -31.96 -45 -a_srs EPSG:4326 $COLORDIR/$D1  $GEODIR/$D1
    gdal_translate -of GTiff -a_ullr -82 12.04 -31.96 -45 -a_srs EPSG:4326 $COLORDIR/$D2  $GEODIR/$D2
    gdal_translate -of GTiff -a_ullr -82 12.04 -31.96 -45 -a_srs EPSG:4326 $COLORDIR/$D3  $GEODIR/$D3

    #Converte para shpefile
    gdal_polygonize.py $GEODIR/$D1 -f "ESRI Shapefile" $SHPDIR "${IMGNAME}_d1"
    gdal_polygonize.py $GEODIR/$D2 -f "ESRI Shapefile" $SHPDIR "${IMGNAME}_d2"
    gdal_polygonize.py $GEODIR/$D3 -f "ESRI Shapefile" $SHPDIR "${IMGNAME}_d3"


texts:
I need to calculate how much (in percent) each of thousands of spatial polygons shrink after some events.

I use R and the very useful sp-library.

Example:

Projection: "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84
+towgs84=0,0,0"
poly.1 is a Spatial Polygon. It has an area-slot with a value
area.1
poly.2 is the poly.1 after shrinking. It has an area-slot with a
value area.2.
I need to calculate area.2/area.1 = percent after shrink.

I know that to get the square-meters etc. of the polygons, they must first be transformed to the appropriate projection. But I hope to avoid to do the transformation, to save calculation-time.

To calculate the percentage-difference, can I use the projection above? I’m not sure what unit the areas are in, square-degrees? Each polygon is less than 200x200 meters.


texts: permissions problem with creating a new text file from the python interpreter in ArcMap?

texts:
I'm attempting to create a .txt file as the output from a larger python script I'm running within ArcMap. The code I'm having trouble goes like this: x = open('adfd.txt', 'w'), and my error message says:

  Runtime error
  Traceback (most recent call last):
    File "<string>", line 1, in <module>
  IOError: [Errno 13] Permission denied: 'adfd.txt'

I've tried this on ArcMap installations on other computers, as well as on non-Esri python interpreters on my own computer (Anaconda's Spyder), and it works fine there, so I'm pretty sure it's something going on with how my installation of ArcMap handles permissions. Does anyone know what might be going on here? I'd like to avoid reinstalling if possible, because that's a big pain with my student license.

Edit prompted by @AlexGIS: When I try to specify a different existing folder, I can change folders fine using arcpy.env.workspave, but I still get the same error message from the open command. When I try to make a new directory within ArcMap using os.mkdir I get a new error message:
os.mkdir('asdfadfa') Runtime error Traceback (most recent call last): File "<string>", line 1, in <module> WindowsError: [Error 5] Access is denied: 'asdfadfa'

Maybe a new clue?


texts:
While using/trying to make a local relief in GRASS with the plugin r.local.relief I stumble across a error.

ERROR: G_realloc: unable to allocate 15000000 bytes of memory at lib/vector/diglib/cindex.c:113

Does this mean that I don't have enough memory(got 16 GB) or that GRASS can't handle this size of data(the tiff is 477 MB) Or is it something totally different. I'm using Windows.

I discovered that the map lib/vector/diglib/cindex.c:113 doesnt exist.

my region is:

g.region -p
projection: 99 (sterea)
zone:       0
datum:      rd18
ellipsoid:  bessel
north:      506250
south:      500000
west:       215000
east:       220000
nsres:      0.5
ewres:      0.5
rows:       12500
cols:       10000
cells:      125000000

Does someone know what the problem is?

This is the whole traceback:

Smoothing using r.neighbors...
Subtracting smoothed from original elevation...
Finding zero contours in elevation difference map...
Reading data...
Displacing data...
Writing vector contour (one level)...
WARNING: 890973 crossings found
Writing attributes...
Building topology for vector map <r_local_relief_7784_vector_contours_tmp@berend>...
Registering primitives...
Current region rows: 12500, cols: 10000
ERROR: G_realloc: unable to allocate 15000000 bytes of memory at lib/vector/diglib/cindex.c:113
Cleaning temporary maps...
Traceback (most recent call last):
  File "C:\Users\berend\AppData\Roaming\GRASS7\addons/script
s/r.local.relief.py", line 260, in <module>
    main()
  File "C:\Users\berend\AppData\Roaming\GRASS7\addons/script
s/r.local.relief.py", line 184, in main
    overwrite=gcore.overwrite())
  File "C:\Program Files (x86)\GRASS GIS
7.0.0\etc\python\grass\script\core.py", line 375, in
run_command
    return handle_errors(returncode, returncode, args,
kwargs)
  File "C:\Program Files (x86)\GRASS GIS
7.0.0\etc\python\grass\script\core.py", line 310, in
handle_errors
    returncode=returncode)
grass.exceptions.CalledModuleError: Module run None
['r.contour',
'input=r_local_relief_7784_subtracted_smooth_elevation_tmp',
'levels=0',
'output=r_local_relief_7784_vector_contours_tmp'] ended with
error
Process ended with non-zero return code 1. See errors in the
(error) output.
(Fri May 22 11:46:33 2015) Command finished (10 min 1 sec)

Does someone knows what the problem is?


texts: Different SVG markers for different attributes

texts:
I am using QGIS 3.0.

I want to use different angle arc SVG markers depending on the data in the attribute table. For example, some points require a 60-degree marker, other a 35-degree marker etc. to give an idea of a field of view from a given point.

Just to clarify, I have 3 different fields of view, depending on the angle view of the photograph. So I want to assign a different angle arc marker depending on the angle specified in the attribute table. View_angle 62, 30, 20.

For each I want an arc marker to match the angle.


texts: Reading a specific overview layer from GeoTIFF file using GDAL / Python

texts:
I know how to create the overviews to a specific geoTIFF file using GDAL ( i have tested it in QGIS), but I cannot find any good information on how to read these overviews individually in Python.
gdal.SetConfigOption('HFA_USE_RRD', 'YES')
outData = gdal.Open("test2.tiff",gdal.GA_Update)
outData.BuildOverviews(overviewlist=[2,4,8])

i would expect the answer to be something like:
array = outData.GetOverviews(overviewlist = 2)


texts: How to create map with graticule with geoserver

texts: I wonder if it is possible to render maps with graticule like this one using geoserver, if not what kind of FOSS4G could do the task in server side.


texts: CartoDB: Center and zoom to stacked layer items/clusters with setView?

texts:
Is there any way to setView on a location with several layer items on top of each other?

Mapbox example:
http://www.zeit.de/zeit-magazin/2015/16/regional-einkaufen-deutschlandkarte

Google Maps example:
http://www.kaffeemaschine-gastronomie.com/kaffeemaschine-gastronomie--Hersteller--s-rancilio.aspx


texts: Why is QGIS 1.8 crashing when exporting to PDF?

texts:
I want to export for example this:

but qgis crashes. Sometimes with an error message std::bad_alloc. But: My Mac has 8 GByte RAM and "just" 5,... something (of 8) GByte are used, moreover there is > 11 GByte free space on the disc.

Sometimes it crashes and you can have a look on an error report, like this one:

Now I noticed that everything works fine if I disclaim labels like "Unserer Lieben Frau". Is this a bug? It seems to be independent on the font used. Or... what can I do to export the results to PDF? Thanks a lot!!

Harald


texts: How to re-project graphics returned from a geoprocessing service on the fly?

texts:
I have a geoprocessing service returning a feature-set in WGS_84. I want to show that graphics on the map using Arcgis JS api.
As I can see my geometries are not re-projected on the fly, since my base-map is in a projected coordinate system and my geometries keep their geographic coordinate system (wgs_84). How can I tell the graphic layer to be re-projected in the fly so as to follow my base-map's spatial reference ?
In other words I want to find a way to reproject a graphics layer in a desired projected system, from inside the arcgis JS api using javascript, and after I have received my results (graphics )from the GP service.


texts: Subset polygon data of SpatialPolygonsDataFrame by Distance to a Set of Coordinates in R?

texts:
I am learning and I made a very similar question when I hadn't learnt what a Spatial Polygons Data Frame was about 10 days ago: https://stackoverflow.com/questions/49193867/select-raster-in-ggplot-near-coastline.

Now, I have discovered the magic of SPDF and choropleth maps and have in essence the same question but with different types of files. I am still wrapping my head around S4 objects and I can't figure out how to subset certain mini-polygons MUNICIPI from my data set.

To the point!

Context:
I have a SPDF that contains these data:

Aim:

I would like to subset all MUNICIPI that are within a certain distance from the coast, clear data for the rest of the MUNICIPI and just work with the my new subset yet plotting the entire geometry.

A choropleth map of what I want looks like this.

Let's say 20km as @urbandatascientist set in his answer to my first question and create a choropleth map of, for example, upper_trees with the subsetted MUNICIPI.

From the https://stackoverflow.com/questions/49193867/select-raster-in-ggplot-near-coastline post I also have the list.of.polygon.boundaries that we'll substract MUNICIPI coordinates from.

Data is here.

Once I subset the coastal MUNICIPI, I'm hoping the map will look something like the green shaded region here.  I have also tried to make sure that coordinates are in the same between list.of.polygon.boundaries.

So far, here's my chloropleth map for the entire region using upper_trees as an example:

tm_shape(catasense2)+
  tm_fill(col="upper_trees",n=8,style="quantile")+
    tm_layout(legend.outside =TRUE)


texts: How to transform WRF Spherical Earth to WGS 84?

texts:
Weather Research Forcasting Model assumes a spherical earth model with a radius of 6370km. I'm looking to display the data onto the WGS-84 Ellipsoid model as that is the model my system uses for its simulations.  We are looking to incorporate real-world changing weather so we are hoping to use WRF's output as we don't want to write our own weather simulator.

While I'm primarily concerned with the continental US I will eventually want to do this for the entire globe.

I've found information, most notably here: http://computing.io/wp/2016/01/wrf-lati ... nd-datums/
That suggests that WPS and WRF may just be treating the WGS-84 information in the land-use data as spherical earth information and this suggest using the XLAT and XLONG output directly from WRF and putting it onto the WGS-84 Ellipsoid. This doesn't sound 100% correct however, and certainly wind information would need to be transformed slightly differently as well well some other information, in particular I would expect the geopotential height may need to be adjusted for the change in surface shape.

I know that to do transforms between Ellipsoids I'd need to use one of the Helmert, Molodensky-Badekas, or Molodensky methods. But it seems like I need reference points on both which I don't actually have for WRF.  Has anyone done such a conversion before or could you provide some insight into the steps to take?


texts:
I'm not sure we can get a straight answer from ESRI on this question at his point because there is so much uncertainty with Microsoft in this area right now but "What is the future of ArcGIS Mobile for Windows?".

The concern is the diminishing future of Windows Mobile 6.x, which is probably the main platform ArcGIS mobile is designed for.  When Microsoft recently announced they are changing the name of Windows Mobile to Windows Embedded Handheld, they have said a new version of this will be released in second half of 2011 and will utilize the latest development tools such as Visual Studio 2010 and Silverlight/XNA.  They did say there is a clear migration path for 6.x apps but that sounds pretty vague to me.

That said, along with the fact you currently need Visual Studio 2008 and WinForms for ArcGIS Mobile development, I don't see it sticking around too much longer.  I see ESRI efforts going more toward the iPhone, Android, Blackberry, and WP7 API's.  Maybe the WP7 API will be the future for the new Windows Embedded Handheld OS?  Any thoughts, predictions?


texts: Can I select a datasource for all layers in QGIS?

texts: I received a zip file with a QGZ map and all the layers that go into the map. When I open the map, it is not recognizing the sources for any of the layers, so I have to go through and individually select the datasource for each layers. It's A LOT of layers, so it would save me a lot of time if I could just select the folder where all the layers are and have it figure out the rest. I'm pretty sure you can do this is ArcMap. Is there a way to do this QGIS?


texts: Error transforming EPSG:4326 to EPSG:90013

texts:
When I'm trying to run a map on geomajas an error has occurred:

  o.g.internal.service.GeoServiceImpl - Problem during transformation
  EPSG:4326->EPSG:900913 of MULTIPOLYGON

The multipolygon coordinates are x = -170 and y = -85.0511... .


texts: Is there a way to disable usage of cached feature types in GeoServer?

texts: Because of thread locking issue connected with code which  locks threads when working with cached feature types I would like to disable this feature. Root issue GeoServer 2.9.2 database connection pooling/thread locking issue


texts: Creating time series chart in Google Earth Engine?

texts:
I'm facing an error in my code when I try to generate time series chart.
Anytime I run the code I get this error :

  Error generating chart: No features contain non-null values of
  "system:time_start"

I have seen post here which suggests to add this code to the function:

  return collection.filterDate(startDate, endDate)
      .reduce(ee.Reducer.mean().combine({
        reducer2: ee.Reducer.minMax(),
        sharedInputs: true
      }))
      .set('system:time_start', startDate.millis());

but it didn't function (I changed the collection into the name of the image collection I have but it didn't work).

here is the code I have :

/**
 * Function to mask clouds using the Sentinel-2 QA band
 * @param {ee.Image} image Sentinel-2 image
 * @return {ee.Image} cloud masked Sentinel-2 image
 */
function maskS2clouds(image) {
  var qa = image.select('QA60');

  // Bits 10 and 11 are clouds and cirrus, respectively.
  var cloudBitMask = 1 << 10;
  var cirrusBitMask = 1 << 11;

  // Both flags should be set to zero, indicating clear conditions.
  var mask = qa.bitwiseAnd(cloudBitMask).eq(0)
      .and(qa.bitwiseAnd(cirrusBitMask).eq(0));

  return image.updateMask(mask).divide(10000);
}

// Map the function over one year of data and take the median.
// Load Sentinel-2 TOA reflectance data.
var dataset = ee.ImageCollection('COPERNICUS/S2')
                  .filterDate('2015-06-23', '2019-07-02')
                  // Pre-filter to get less cloudy granules.
                  .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 20))
                  .select('B2','B3','B4','B8','QA60')
                  .filterBounds(geometry)
                  .map(maskS2clouds);

var rgbVis = {
  min: 0.0,
  max: 0.3,
  bands: ['B4', 'B3', 'B2'],
};

var clippedCol=dataset.map(function(im){
   return im.clip(geometry);
});

//test if clipping the image collection worked
Map.centerObject(geometry);
Map.addLayer(clippedCol.median(), rgbVis, 'RGB');

// Get the number of images.
var count = dataset.size();
print('Count: ',count);
// print(clippedCol);//here I get the error messege "collection query aborted after accumulation over 5000 elements
// print(dataset,'dataset');//the same error here

//function to calculate NDVI
var addNDVI = function(image) {
  var ndvi = image.normalizedDifference(['B8', 'B4']).rename('NDVI');
  return image.addBands(ndvi);

};

//NDVI to the clipped image collection
var withNDVI = clippedCol.map(addNDVI);

// // Test the addNDVI function on a single image.
// var ndvi1 = withNDVI.select('NDVI').mean();

var colorizedVis = {
  min: 0.0,
  max: 1.0,
  palette: [
    'FFFFFF', 'CE7E45', 'DF923D', 'F1B555', 'FCD163', '99B718', '74A901',
    '66A000', '529400', '3E8601', '207401', '056201', '004C00', '023B01',
    '012E01', '011D01', '011301'
  ],
};

// Test the addNDVI function on a single image.

// var ndvi1 = withNDVI.select('first').mean();
// var ndvi2 = withNDVI.select('second').mean();

Map.centerObject(geometry);
print(withNDVI,'NDVI Images');
// Map.addLayer(ndvi1,colorizedVis,'test1');
// Map.addLayer(ndvi2,colorizedVis,'test2');

//Long-Term Time Series
var NDVIalltimes = Chart.image.seriesByRegion(withNDVI,geometry, ee.Reducer.mean(),'NDVI',
30, 'system:time_start').setOptions({
title: 'NDVI for Royal Exchange',
vAxis: {title: 'NDVI'},
});
print(NDVIalltimes);

The end goal is to get chart that shows the mean NDVI value for each Image in the image collection.


texts: Are shapefile coordinates unaffected when layer exported from Data Frame with different CRS?

texts:
I would need to deliver a result shapefile in WGS 84 (GCS: Geographic –lat/lon- coordinate system), the shapefile is originally in WGS84; but in the process of creating this result shapefile I need to work with planar data.  What I am considering is to set the Data Frame coordinate system to a proper projected coordinate system for that zone, work with the shape layer  (in WGS 84 (GCS) whereas the Data frame is in a Projected coordinate system) and once finished working and obtained the result shape layer, Export the shapefile by doing Data > Export Data and, for the option “Use the same coordinate system as:”, choose “this layer’s source data”.

Is the resulting .shp going to be fully operative and to keep its geographic coordinates unmodified even if the Data frame coordinate reference system is a projected one?


texts:
I'm looking into possibilities to align track data from a GPX file with known roads, perhaps coming from OSM.

So far I'm able to show the tack as such on a map powered by leaflet. Since the app should be as easy to use as anyhow possible I look for web based and if possible fully automated solutions.

The GPX files are stored in a standard CouchDB, so is the App itself.

I think about an approach in which users can upload their GPX file. Once they do the map goes into an edit mode that shows the track and may already suggest a "road only" version of the track. In a best case scenario all the user has to do is to accept the suggestion.

How could something like that be done? Are there any libraries or plugins doing that?


texts: More than one result for coordinates from PostGIS?

texts:
I'm using a modified version of pyris to get IRIS codes from France.  I have a PostgreSQL/PostGIS database based on a shapefile.

My problem is, when I query the result for a coordinate, I get more than one result. That should not. An IRIS code should be unique for one result.

For example:

SELECT
iris ,nom_com ,depcom ,nom_iris ,dcomiris ,typ_iris
FROM geoiris
WHERE geom && st_makepoint(4.808744, 45.980649);

iris |     nom_com     | depcom |    nom_iris     | dcomiris  | typ_iris
------+-----------------+--------+-----------------+-----------+----------
 0000 | ARS-SUR-FORMANS | 01021  | ARS-SUR-FORMANS | 010210000 | Z
 0000 | MISERIEUX       | 01250  | MISERIEUX       | 012500000 | Z
 0000 | SAINTE-EUPHEMIE | 01353  | SAINTE-EUPHEMIE | 013530000 | Z
(3 rows)

The second result is right, but the two others not. The others are areas next to the right one.

Is there something wrong with my query?

I'm using PostgreSQL 9.3 and shp2pgsql for the shapefile.


texts:
What are all the folders I have to share with others who would like to use my data ?
I gave workspaces and styles folder. But it is not working in Geoserver 2.4.4. Since layergroup is inside the workspace I ignored the layergroup folder


texts: Exporting to PDF from print composer in QGIS with two maps on layout?

texts:
I'm using QGIS Desktop 2.10.1, and I created a print composer with two different maps.  The maps are based on the same shapefile, but each has a different variable and thus a different color scheme (one blue, one green).  To create them, I duplicated my shapefile in the main window and applied a different variable to each copy.  When I attempt to export my map to pdf, the resulting pdf looks fine except that instead of the two different maps, it has simply duplicated the first one where the second should be.  However, the print composer still looks correct and has two distinct maps.  I have locked both map layers.

Does anyone have an idea of what I'm doing wrong?


texts:
Using osm2po to extract traffic signal in the vertices used to provide the value 1 in the clazz column of the vertex table:

Add vertices costs to postgresql routable network

For some reason I keep getting the value -128 now. Is that supposed to be like this?


texts: Image sizes for VHR remote sensing image

texts:
I need to compare three or four different segmentation algorithms on very high spatial resolution remote sensing images, with a range between [0.31,1.00] meters/pixel.

What image size could you recommend to me?
Maybe 100x100, 500x500, 1000x1000 or other size?


texts:
I'm trying to compile QGIS 3.0 on Windows. I want to port it to Windows RT (Windows on Arm)

I keep getting this error:

Looking for gaiaStatisticsInvalidate in C:/OSGeo4W/lib/spatialite_i.lib
Looking for gaiaStatisticsInvalidate in C:/OSGeo4W/lib/spatialite_i.lib - not found
CMake Error at cmake/FindSpatiaLite.cmake:75 (MESSAGE):
  Found SpatiaLite, but version is too old.  Requires at least version 4.2.0
Call Stack (most recent call first):
  CMakeLists.txt:251 (FIND_PACKAGE)

How do I fix this? I have tried setting up OSGeo4W with spatialite 4.3, but it still thinks that my version is too old.

I have been following the install guide in the source code here

EDIT 1: According to CMAKE I get two errors: https://hastebin.com/woherijadi.tex


texts: How to generate DTM from a DSM (having 12.5 meter accuracy)?

texts:
The data that I have is an open-source DEM having a spatial resolution of 12.5 meters.
Now I want to generate the DTM out of it.
Is there any algorithm or model available to extract DTM out of DEM (using open-source or ArcGIS software)?


texts:
I have a scanned map in which roads are shown in red color, the RGB value being R>250, G=0 and B=0.

Is there a way to auto digitize the roads by specifying RGB criterion or snap criterion?

My scanned map only has the option of RGB composite or stretched under display options


texts: How to remove layers of detail from an osm file?

texts: I want to remove some layers from the osm map of an area. It has all the details like POI, offline routing and all. But I just want the map data. How can I remove this extra layer from the map?


texts: Configurable Map Viewer (CMV) - Print Widget - Proxy Issue

texts:
I'm not sure this is necessarily an issue with the CMV as much as it is with how I may have the proxy configured.  I'm not sure exactly what to ask so I'll start with describing the symptom.  When I am zoomed to a full extent and have few features being rendered the print widget generates a pdf just fine.  When I'm zoomed in and many features are rendered the print widget generates an error.

Is this a result of the url length and the proxy coming in to play?

The error is:

  POST
  http://myurl/myfolder/proxy/proxy.ashx?http://myurl2/arcgis/rest/services/Utilities/PrintingTools/GPServer/Export%20Web%20Map%20Task/execute
  500 (Internal Server Error)

When I copy and paste the url in my browser the export web map task opens fine, through the proxy.

I have <serverUrl url="http://myurl2" matchAll="false"/> in proxy.config.  The error message in the stack says ""Error: Unable to load proxy/proxy.ashx?http://myurl2/arcgis/rest/services/Utilities/PrintingTools/GPServer/Export%20Web%20Map%20Task/execute status: 500."  Should this show the full url to the proxy?

I have 2 operational layers, one of which has a minimum scale defined.  I noticed that the print works fine until the the minimum scale is reached on that layer.  Once the minimum scale is reached, the error above is thrown.  If you zoom out beyond that scale, the error is thrown until you refresh the page.  I replaced my operational layer that had the minimum scale with one that that didn't (the same feature class) and the print works fine.  I'm curious if this is expected (can others duplicate this) and what I might do to resolve it?


texts: How to manually input decimal value into z-coord in QGIS

texts: If I want to manually input a decimal value to Z-coord in QGIS, there's no effect. Only integer values are saved in nodes editor. Operating system Ubuntu 18.04, QGIS 3.6, locale: poland, data and money format: polish.


texts: Display points from CSV table (feature collection) in Google Earth Engine

texts:
Unfortunately, I cannot display the point data I uploaded as csv file. I specified latitude and longitude, the rest is unknown so I kept it at default. Can somebody help?
Link to asset


texts:
I set up GeoServer cluster with ebmedded JMS-brokers (basically as described here -  https://geoserver.geo-solutions.it/edu/en/clustering/clustering/active/index.html) in Multiple Masters - Multiple Slaves mode. Until yesterday everithing worked fine (~1 month).
But then we get just stuck in any changes made to catalog. Log was hanging on

Creating an instance of: class org.geoserver.cluster.impl.handlers.catalog.JMSCatalogModifyEventHandlerSPI

In my setup all configuration changes sending to just one (primary) GeoServer instance, based on haproxy load balancing.
I spent some time and found that GeoServer embedded broker database (kahadb) grew to 1Gb on this problem instance, and that default setting for disk usage in GeoServer config for it is exactly 1Gb -
activemq.broker.systemUsage.storeUsage=1 gb. Deleting kanadb and restarting problematic instance did the trick - it works again.

And now the question - should each GeoServer instance have own copy of kahadb? Or should it be the same kahadb for each instance?

If each GeoServer instance should have own kahadb, why log-files is not truncated? ActiveMQ documentation says, by default they will be cleared every 30s, but seems it is not working properly. More strange - except this kahadb disk space issue, cluster works just fine, replicating all the changes to each node.


texts: How can I assign attribute IDs to points in an 85000 point spread according to a grid scheme?

texts:
I am working on a pretty complicated grid system and cannot figure out a way to name each of the points in my shapefile.

The X of the grid originally started with A and ended with AX.  We have since had to broaden it both east and west to include YE through ZZ, then A through CT (so 146 columns in the grid).  The Y of the grid is 1 through 147.

To further complicate things, there is a great deal of overlap with the points that meet at the vertices of the grid.  This is because we have to name each corner of the square according to the grid AND the cardinal direction (e.g. A1_SE, which overlaps with B1_SW, A2_NE, and B2_NW).  Each of these will be associated with the projected coordinates (e.g. A1_SE at 700000, 300000).

I have tried playing around with some excel formulas, autofill lists, etc, but cannot figure out how to assign these IDs.  I am using ArcGIS 10.3 and don't yet have experience with python and extremely limited experience with R.

A simpler question:
How can I populate the attributes of the grid squares (in the polygon shapefile) with the grid name info?  e.g. A1, ZX146, CD56


texts: Undo option to 'move edits to base' ArcSDE

texts:
I have an ArcSDE Geodatabase where a number of layers are registered as versioned. Having spent some time troubleshooting an issue I'm having with syncing data on Collector, I have noticed that the layers which are not loading are those which are registered as versioned with the option to move edits to base.  I suspect this is my problem and I am looking for a way to remove the option easily.

Presumably, if I reconcile and post all changes to the DEFAULT version, then unregister as versioned, then register as versioned again (being careful not to tick the box with the move edits to base) this will solve my issue without losing any data? If I do this, will it be applied to all feature classes or can I specify the ones I want to unregister?

Is this a viable option or are there better ones? What are the risks?


texts:
Is there a way to make GeoPandas work with geographies?
If GeoPandas were used with pygeos, is there a way to force all calculations (ST_Distance, ST_Length, etc.) to be geodesic?  How do I tell it not to use Shapely, but PyGeos (or whatever) and that it should do geodesic math?
I like using GeoPandas, but in order to do work with spherical geographies it seems necessary to put the data into PostGIS and have it do spherical calculations.  Can GeoPandas be forced into geoid operations?  Is it just a matter of setting the proper CRS on a dataframe?


texts: Are there Bookmarks with selection highlighted in QGIS?

texts:
I have a search query on the attribute table which highlights the features and then zooms to that particular layer. I take a bookmark for this particular view with the selection being highlighted.
But later when I see the same bookmark it does not highlight those selected features.Any settings or any other plugin for this.


texts: How could I extract TMS and WMS link from ArcGIS Server?

texts:
I need link of TMS and WMS from our custom ArcGIS maps in order to use in qlikview.

How could I extract this link?


texts: Wrong position of tiles when using Proj4Leaflet for a tiled raster layer

texts:
Using Leaflet I am trying to display a tiled raster layer described at the following url : https://mpt.svp.sk/server/rest/services/podkladove_mapy/orto_2019_12v/MapServer

I understood that the layer was not available in the web mercator projection EPSG:3857 from that server. So I have tried to display it in the original projection EPSG:102067 .
My Leaflet map currently displays the layer tiles but at the wrong positions. Here is a Jsfiddle that demonstrates where I am standing:
https://jsfiddle.net/uhcbpz73/

Here is the JS code I am using:

var crs = new L.Proj.CRS('EPSG:102067',
  '+proj=krovak +lat_0=49.5 +lon_0=24.83333333333333 +alpha=30.28813975277778 +k=0.9999 +x_0=0 +y_0=0 +ellps=bessel +units=m +no_defs',
  {
    resolutions: [
      529.1677250021168, 264.5838625010584, 132.2919312505292, 66.1459656252646, 26.458386250105836, 13.229193125052918, 6.614596562526459, 2.6458386250105836, 1.3229193125052918, 0.5291677250021167, 0.26458386250105836, 0.13229193125052918, 0.06614596562526459
    ],
    origin: [-600000.0,-1120000.00]
  });

var map = new L.Map('map', {
    crs: crs
  }).setView([48.5, 19], 8);

L.tileLayer('https://mpt.svp.sk/server/rest/services/podkladove_mapy/orto_2019_12v/MapServer/tile/{z}/{x}/{y}', {
  maxZoom: crs.options.resolutions.length-1,
  minZoom: 0,
  continuousWorld: true,
  tms:true
}).addTo(map);


texts: What are the highest and lowest precision Coordinate Reference Systems for US East Coast?

texts: And does it matter what system your data layer begins with if you choose the right system for your project?  Do some applications use more precise projection systems than others?


texts: I am looking for California state parcel data such as owner name, no. of stories, year built etc., for each parcel. I can get to the level that I can browse them in property appraiser site but don't see data download options.


texts:
I have layer of water points and a layer of houses. I need to get the number of houses that are within 500 m from a waterpoint. After this I do a spatial join (join one to one) and then I get number of houses that are within 500 m of a waterpoint.

However many of the buffers are overlaying and I am wandering if the result I get in the attribute table includes all the double counts, and if so, how do I get rid of double counting?


texts: ST_DistanceSphere can't deal with the prime meridian?

texts:
I have a table of weather stations, and I wanted to find stations near London. I ran this query:

select id, name, location, ST_DistanceSphere(location, ST_GeomFromText('POINT(51.51 -0.12)')) as dist from stations order by dist;

The query seems to work fine, up until the 13th row:

id          name        location                    dist
USC00500252 AMCHITKA    POINT (51.3833 179.2833)    94097.52204081

This is a weather station in the Aleutians. Needless to say, it is quite a lot more than 94km from London (and I expect the weather there is not especially similar to London's). How can I get ST_DistanceSphere to account for longitude correctly?

NB: I've tried adding the SRID: ST_GeomFromText('POINT(51.51 -0.12)', 4326). This doesn't make any difference.


texts: Stacking multiple selected points in ArcGIS

texts: I would like to find a way to stack multiple selected points in ArcGIS, either using the field calculator or python script. For example, I want to snap all points representing houses on a street into one point. Is there a way to reference one of the existing coordinate sets to use as the values for all selected points?


texts: Contour lines and grid intersection points in QGIS

texts:
I'm looking for a way to get a list of altitude for de point in a grid and I think that is possible using intersection points between contour lines and the grid, in a similar way that the tool MapEditor of Wasp do (https://www.wasp.dk/dataandtools#map-editor, see that image)
The desired output is something like:
[East Coord] [North Coord] [elevation]
x1 y1 z1
x2 y2 z2
....
xn yn zn

Do you know some tool to do that?


texts: Splitting shapefile with counties into single boundary lines

texts:
I'm a beginning QGIS user and have a shapefile with polygons of all the county borders in Holland.

When I style this shapefile with a non-fill and a dashed line, all the adjacent (overlapping) lines are almost solid instead of dashed. This is because adjacent county borders are duplicates.

How do you get rid of those lines? tried to set the polygon shapes to lines and tried several geometry tools. Is this even possible? Or do you have to remove, clip all the double lines by hand?