texts: Change length of string fields in featureclass using arcpy v10.1

texts:
I am trying to edit the length of all the text fields in a feature class.  The feature class has a lot of fields so I want to do this in python.  I only have access to v10.1 so cannot use the AlterField_management option.

The way I'm going about it is to use fc to fc conversion with an appropriate field mapping.

Below is my code but it does not alter the length of the text fields as intended.  The fc to fc conversion goes ahead but with no changes to any of the field properties.  I'm pretty sure the problem is with the field mapping but I can't figure it out.

import arcpy

infc =  r"C:\path to original fc"

out_space = r"C:\path to out gdb"
out_fc = "fc_newfieldlengths"

# set up field mappings object
fms = arcpy.FieldMappings()

# Create list of all non-string fields to skip when creating the fieldmaps

text_fields = arcpy.ListFields(infc, field_type = 'String')
all_fields = arcpy.ListFields(infc)

skip_fields = [x for x in all_fields if x not in text_fields]

# create field mapping with new length for the text fields
for field in all_fields:
    if field in skip_fields:
        pass
    else:
        fm = arcpy.FieldMap()

        fm.addInputField(infc, field.name)

        outfield = fm.outputField
        outfield.length = 150
        outfield.name = field.name

        fm.outputField = outfield

        fms.addFieldMap(fm)

# Copy featureclass with new field mapping
arcpy.FeatureClassToFeatureClass_conversion(infc, out_space, out_fc, field_mapping = fms)

# End of script

If I tried only looping through a list of string fields to generate the fieldmaps I found that all the fields were changed to string fields in the output, not just the fields I am interested in amending.

To avoid changing the data type for all the fields I created a skip list of fields not to alter - this manages to preserve the datatypes but has not amended the text field lengths (as was occurring previously but just also annoyingly changing all the datatypes as well).

I suspect I just don't understand field mappings well enough to know what I'm doing wrong.


texts: 16-day-MODIS NDVI values seem to be wrong compared to my soil moisture ground data!

texts:
I'm using MOD13Q1 product of MODIS for estimation soil moisture. As I know, this product is made up of the maximum value of NDVI of each pixel over 16-day calculated NDVI values. The problem is that my soil moisture in situ measurements don't correspond to the NDVI values. For example, where NDVI is 0.75, the calculated soil moisture is 0.2; in another case, where NDVI is 0.4, the soil moisture is 0.5, which is the highest value of soil moisture calculated. Of course, this problem doesn't occur for all the pixels.

For working with MODIS product, I use ENVI software. I first use MODIS Conversion Toolkit(MCTK) to transform its sinusoidal coordinate system into Geographic Lat/Long system. Then, I examine the NDVI values based on the location of my soil moisture ground measurements.

I know that my ground measurements are point-wise and the NDVI product covers around 250m * 250m on the ground, but soil water content over this area doesn't differ significantly. What's the problem in this case? where I'm making a mistake?


texts: Converting text in image to labels in ArcMap?

texts:
I have map like the below image and I want convert numbers to labels automatically.

Is there any method to do this?


texts:
I have converted a polygon file to raster (ArcMap 10.4), after setting extent and cell size environments to match my DEM (which is also set as snap raster). On the new raster, the extents are the same as the DEM EXCEPT for the right (top, bottom and left match). The new raster has one more column than DEM, which prevents me from using it in the model I want to use.

How do I fix this?


texts:
Is it possible to add a where condition to an update query in MapInfo?

At the moment I have to do a query to select the records I need and then do an update on that query.

I believe this is possible with SQL so I wanted to know whether it could be done in MapInfo.


texts: Clickable label (HTML, Link) in PDF export (QGIS)

texts:
I need to export a (GEO)PDF - overview map with bounding boxes - with clickable links in every bounding box.

I tried to simply export a PDF with html-links as text labels, but as you can see in the screenshot of my PDF below, the HTML syntax is not interpreted as link

I also tried simple links like: <a href>www.qgis.org</a>

This solution works on clickable text field links (rendered as HTML) but does not match my needs:
Exporting PDF from QGIS Print Layout including clickable hyperlinks

Even this post is not helping me: Clickable HTML link in QGIS print composer pdf export?


texts: How to store AutoCAD data in ArcSDE database?

texts:
Does anyone have an idea how i can create a database instance on my ArcSDE geodatabase to store my .dwg files?

I have over 1 million of these files and would be happy to automatically extract and migrate them into a specified schema in the Geodatabase.

I am using ArcSDE 10.1 on SQL Server2008 R2


texts: How to quickly match many high-scoring geocoding ties

texts:
I am geocoding Portland Metropolitan Area in ArcMap 10.2 based on Street Address, City, State and Zip. I get decent matches (about 93%). However, I get thousands of high-scoring ties. For example 2 or 3 matches scoring 100%. These tied values are all "separate" points (e.g. two different points within a property polygon, and on point on the street facing the polygon).

Is there a way to either:

non-manually match at random based on any one of the top scoring ties, or
take the average lat and average long of the top scoring ties?

I am a noob, so I will happily entertain ways to improve this question based on constructive criticism.


texts: How to uninstall QGIS 3.8 from MacOS in a clean way?

texts:
I have two installations of QGIS, 3.14 and 3.8, I want to uninstall 3.8. What is the proper way to do it? Do these instructions also apply to this newer version?
The following are the paths that the Python interpreter sees:
/Applications/QGIS3.8.app/Contents/MacOS/../Resources/python
/Users/my_user/Library/Application Support/QGIS/QGIS3/profiles/default/python
/Users/my_user/Library/Application Support/QGIS/QGIS3/profiles/default/python/plugins
/Applications/QGIS3.8.app/Contents/MacOS/../Resources/python/plugins
/Applications/QGIS3.8.app/Contents/Frameworks/Python.framework/Versions/Current/lib/python3.7/site-packages
/Applications/QGIS3.8.app/Contents/Frameworks/Python.framework/Versions/Current/lib/python3.7
/Applications/QGIS3.8.app/Contents/Frameworks/Python.framework/Versions/Current/lib/python3.7/site-packages/geos
/Applications/QGIS3.8.app/Contents/Frameworks/Python.framework/Versions/Current/lib/python3.7/lib-dynload
/Applications/QGIS3.8.app/Contents/Resources/python
/Applications/QGIS3.8.app/Contents/Frameworks/Python.framework/Versions/Current/lib/python37.zip
/Users/my_user/Library/Application Support/QGIS/QGIS3/profiles/default/python

Those that are not in /Applications/QGIS3.8.ap are common to both installations.


texts: Achieving Flat-End line buffers in QGIS?

texts:
I need to buffer a series of lines by 40m, however I only want the buffer to come off the sides and not off either end of the line.

I know in ArcGIS there is a "flat-ends" tickbox that achieves exactly this, however unfortunately I can't use ArcGIS at this current time.

Is there any way to achieve this in QGIS?

I've had a look at some of the work arounds, but due to the scale of my datasets (national road networks) none of the ones I have found so far are feasible.


texts: Select individual lines from feature class, convert to unique rasters: "raise e" error?

texts:
I have a feature class with 5 lines in it, and need select each one and turn them into its own raster -- for example, line 1 is raster 1, line 2 is raster 2, etc. This is the python script:

import arcpy

SightLines = "C:\\temp\\SightLines.shp"

index = 0
rows = arcpy.SearchCursor(SightLines)
for row in rows:
    RasterFileName = "C:\\temp\\fakeData.gdb\\Raster" + str(index)
    arcpy.AddMessage("Raster File " + RasterFileName)
    arcpy.FeatureToRaster_conversion(row, "FID", RasterFileName, "1")
    arcpy.AddMessage(" >> MADE SUCCESSFULLY")
    index = index + 1

Unfortunately, this is throwing the following error at the FeatureToRaster_conversion step:

    raise e
RuntimeError: Object: Error in executing tool

I can't figure out what raise e is trying to tell me. What exactly is going wrong?


texts: ArcGIS model and python script integration

texts:
I'd like to make a model with a last step, a python script, that will send an email, containing messages and status of the model. The email sending part is okay, nut the hard part for me is to query the status and the message. I found the "result" and the "GetMessages" things but could not succeeded.

I'd like to make it work on a really simple model than use it when its needed. The dummy model looks like this:

And the script:

import arcpy
arcpy.ImportToolbox(r"c:\arcGIS_2_email\toolbox\tulboksz.tbx")
arcpy.topmodel_tulboksz()
result = arcpy.topmodel_tulboksz()
print result.getMessages()


texts: Selection input for QGIS Processing models

texts:
QGIS 2.14.11 LTR. In a processing script we can define a selection as input parameter by ##auswahl_test=selection spam;egg;dummy resulting in

Background is, that I want to provide a selection of formulas (strings) as input for field calculator algorithm used in a model. By now, I'm typing these formulas in a string input field by hand, which is annoying, time-consuming and error-prone.

Now is there a way to obtain such input field in the auto-generated gui of a processing model? Or any other suggestions of how to provide a selection of strings as a model input?


texts: QGIS becoming unresponsive after I try and add (relatively small) layer

texts:
Each time I try and add a CSV as a delimited text file, I see this:

The file is not large. Only ~2k rows and 462 kb. I've certainly uploaded larger files to QGIS and never had any issues. I've force quit the app and tried to re-add the file, but that hasn't worked. I've started a new project to see if I can add the file by itself, and it still is unresponsive. I've waited for more than an hour and no dialog shows up after this white box. What can I do? Is it a RAM issue? I have more than 20 gb of storage on my machine. It's QGIS 3.10.


texts: How can I open, visualize and convert LandXML data?

texts:
I've received what I assumed was a DWG file but when I open it in AutoCAD Map 3D 2015 only data frames placeholders are shown. I was under the assumption that the person who sent the data forgot to include linked objects and drawings.

When I look as the files i received there is a single DWG file which i opened and long list of XML files.  I was told the project was created in Civil 3D and the XML files are some type of spatial 3D data.  Possibly LandXML?

I thought that a full features AutoCAD product such as Map 3D would be able to open any autodesk spatial data.

Ultimately I need to convert this data into an ESRI geodatabase feature class.

How can I open, visualize and convert LandXML in AutoCAD Map 3D 2015 or ArcGIS 10.3 (Standard, Spatial Analyst, 3D Analyst)

Note: I tried the LandXML to TIN tool but its giving me these errors:

  ERROR 000732: Input: Dataset
  T:\GIS\DATABASES\Rainy_River_Project_Mine_Site_Features\Received
  Data5-06-11 Water Management\Extracted\Clark Creek Dam.xml does not
  exist or is not supported ERROR 000732: Output TIN Folder: Dataset
  T:\GIS\DATABASES\Rainy_River_Project_Mine_Site_Features\Received
  Data5-06-11 Water Management\Extracted\CONVERTED does not exist or is
  not supported Failed to execute (LandXMLToTin).


texts: How to combine bands from different rasters using the global mapper raster calculator?

texts: I am wondering how to specify what raster you are getting the band from in the global mapper raster calculator? When you open it with a more then one raster open you can check the boxes to define which ones you use. But once you are in the actual calculator there is no way (that i can see) to specify which raster band 1 is coming from.


texts:
I am trying to polygonize one of my rasters with PyGIS. I do this with the below script.

sourceRaster =
gdal.Open('RASTER.tiff')
sr_proj=sourceRaster.GetProjection()
raster_proj = osr.SpatialReference()
raster_proj.ImportFromWkt(sr_proj)

band = sourceRaster.GetRasterBand(1)
print  band
bandArray = band.ReadAsArray()
outShapefile = "POLYGON"
driver = ogr.GetDriverByName("ESRI Shapefile")
outDatasource = driver.CreateDataSource(outShapefile+ ".shp")
outLayer = outDatasource.CreateLayer('polygonized', srs=raster_proj)
newField = ogr.FieldDefn(str(1), ogr.OFTInteger)
outLayer.CreateField(newField)

gdal.Polygonize( band, None, outLayer, -1, [], callback=None )
outDatasource.Destroy()
sourceRaster = None

print "finish"

When I run the above, a polygonized version of my raster file has been created. The polygon looks like expected, it at least seems to fit the raster.

The weird thing however is that the attribute table is empty and it looks like only the shadow is projected

also when I Print the band

    band = sourceRaster.GetRasterBand(1)
print  band

Below is printed

osgeo.gdal.Band; proxy of Swig Object of type 'GDALRasterBandShadow *' at 0x12a90a450
finish

When I use the QGIS GDAL Polygonize inbuilt function, the polygon looks like expected with a populated Attribute table.

Who can help?


texts: Event raised when selected layer changes in TOC

texts:
Is there an interface which exposes events for when a TOC selected item changes?

I can find snippets showing how to get the currently selected feature from the TOC.  However, I have not been able to find anything that allows you to track when a selection changes - only query the current state.

What would be the best approach (preferably other than polling) to track the user's current layer selected within ArcMap and ArcScene?


texts: Converting large CSV to line-delimited GeoJSON?

texts: What is a good way to convert a large (100MB) CSV file to GeoJSON on the command line, preferably without GDAL? I tried Mapbox's csv2geojson but it ran out of memory.


texts:
I'm trying desperately to write my first query in pgRouting using maps in osm format (as shown in http://www.carsonfarmer.com/?p=799).

Im using Ubuntu 10.10 with Postgresql 8.4 ang postgis 1.5.
After adding pgrouting functionality to my database and run osm2pgrouting command:

./osm2pgrouting -file /home/cfarmer/Downloads/ireland.osm -conf mapconfig.xml -dbname routing -user postgres -clean -passwd myPW

I get this message (that looks to be correct):

connection success
Trying to load config file mapconfig.xml
Trying to parse config

and I was expecting to find in my base tables such as ways, but I see only geometry_columns and saptial_ref_sys.
Moreover when I try to run the query:

select * from shortest_path
        ('select gid as id,
            source::int4,
            target::int4,
            length::double precision as cost
            from ways',
        52343, 39219, false, false);

I get this error: errore di sintassi vicino il simbolo non atteso "from".
If I run the same query on pgAdmin I get a more specific error i.e.:
relation "ways" does not exist.

Can anyone help me?
Thanks a lot
Laura

postgres@laura-VirtualBox:/home/laura$ cd /usr/bin
postgres@laura-VirtualBox:/usr/bin$ ./osm2pgrouting -file /home/laura/Scaricati/ireland.osm -conf mapconfig.xml -dbname routing -user postgres -passwd laura -clean
host=127.0.0.1 user=postgres dbname=routing port=5432 password=laura
connection success
Trying to load config file mapconfig.xml
Trying to parse config
postgres@laura-VirtualBox:/usr/bin$ select * from shortest_path ('select gid as id,source::int4,target::int4,length::double precision as cost from ways', 52343, 39219, false, false);
bash: errore di sintassi vicino il simbolo non atteso "from"
postgres@laura-VirtualBox:/usr/bin$ psql -U postgres routing
psql (8.4.8)
Digita "help" per aiuto.

routing=# select * from shortest_path ('select gid as id,source::int4,target::int4,length::double precision as cost from ways', 52343, 39219, false, false);
ERROR:  relation "ways" does not exist
RIGA 1: ...int4,target::int4,length::double precision as cost from ways
                                                               ^
QUERY: select gid as id,source::int4,target::int4,length::double precision as cost from ways
routing=#


texts:
This link suggests that you can saves edits (Currently Editor -> Save Edits) as the same time as a project.

Is it still possible to do something similar, change the option(s) so it saves edits AND the project using ArcMap?


texts: Error converting Json file to a feature class

texts:
I want to convert some json data being returned from a request to a feature class, but I keep getting a RuntimeError: Object: Error in executing tool. I made a test gdb in a folder 'C:/Workspace/Sandbox/ScratchTests/cslf.gdb' in which to populate my new feature class after conversion. To test that the request is correct and that I am returning Json data, I added a couple of print statements. Otherwise, everything is pretty straight forward. Does anyone see a problem with my code? I am following the arcpy [JSON to Features][1] directions from the documentation. The only thing I am doing differently is instead of using an actual file, I am just plugging in the variable 'cslfJson`.

import arcpy, sys, os, arcgis, requests
arcpy.env.workspace = "C:/Workspace/Sandbox/ScratchTests"

params = {'f': 'json', 'where': '1=1', 'geometryType': 'esriGeometryPolygon', 'spatialRel': 'esriSpatialRelIntersects','outFields': '*', 'returnGeometry': 'false'}
r = requests.get('https://hazards.fema.gov/gis/nfhl/rest/services/CSLF/Prelim_CSLF/MapServer/3/query', params)
print(r.url)

cslfJson = r.json()
print(cslfJson)

arcpy.JSONToFeatures_conversion(cslfJson, os.path.join("cslf.gdb", "cslf"))

Also, Here's the Traceback:

  https://hazards.fema.gov/gis/nfhl/rest/services/CSLF/Prelim_CSLF/MapServer/3/query?f=json&where=1%3D1&geometryType=esriGeometryPolygon&spatialRel=esriSpatialRelIntersects&outFields=%2A&returnGeometry=false
  Traceback (most recent call last):

  File "", line 1, in
      runfile('C:/Workspace/Sandbox/ScratchTests/CSLF.py', wdir='C:/Workspace/Sandbox/ScratchTests')

  File
  "C:\Users\jbridwell\AppData\Local\Continuum\anaconda3\envs\acrpro\lib\site-packages\spyder_kernels\customize\spydercustomize.py",
  line 668, in runfile
      execfile(filename, namespace)

  File
  "C:\Users\jbridwell\AppData\Local\Continuum\anaconda3\envs\acrpro\lib\site-packages\spyder_kernels\customize\spydercustomize.py",
  line 108, in execfile
      exec(compile(f.read(), filename, 'exec'), namespace)

  File "C:/Workspace/Sandbox/ScratchTests/CSLF.py", line 20, in

      arcpy.JSONToFeatures_conversion(cslfJson, os.path.join("cslf.gdb", "cslf"))

  File "C:\Program
  Files\ArcGIS\Pro\Resources\ArcPy\arcpy\conversion.py", line 403, in
  JSONToFeatures
      raise e

  File "C:\Program
  Files\ArcGIS\Pro\Resources\ArcPy\arcpy\conversion.py", line 400, in
  JSONToFeatures
      retval = convertArcObjectToPythonObject(gp.JSONToFeatures_conversion(*gp_fixargs((in_json_file,
  out_features), True)))

  File "C:\Program
  Files\ArcGIS\Pro\Resources\ArcPy\arcpy\geoprocessing_base.py", line
  506, in
      return lambda *args: val(*gp_fixargs(args, True))

  RuntimeError: Object: Error in executing tool

  [1]:
  http://desktop.arcgis.com/en/arcmap/10.3/tools/conversion-toolbox/json-to-features.htm


texts: Error with categorized render tool in model builder QGIS 3.4.3

texts:
In my model, I'm trying to color my polygon layer using an xml file through categorized render tool (under cartography toolset). However, when I run the model the output file has the attribute data but nothing is visible in the canvas.

In the same model, when I'm doing the same process for a line layer it gives me the color as specified in the xml file

Any idea how this can be solved? Please help me

PS: When I'm using the same xml file to color the polygon layer outside the processing modeler it is working perfectly fine


texts: Can I create a polyline with one long segment and three shorter underlying segments?

texts: We currently have a street centerline polyline feature class where the street segments correspond to our pavement management database.  Each segment has a unique ID that matches a record tracking pavement type/year.  Segment breaks occur at intersections, but can also occur mid-block if a subdivision was built in multiple additions. We would like to start adding Left From/To and Right From/To fields to the street centerline feature class for geocoding purposes.  I would like to be able to have a street segment start at an intersection and end at an intersection for address coding, but still keep the underlying multiple segments that correspond to the breaks in the pavement management database.  Is there an easy way to do this?  Do I need to have two seperate line files or is there way to have this in one feature class?


texts: QGIS doesn't display compound polygon (from a spatial oracle DB)

texts:
When I want to display the elements stored in a spatial oracle DB, the compound elements aren't display. There is no error message.

The type of these compound elements is "1005".

Are these elements managed by QGIS? Can you suggest me how to display this type of polygon?

I'm using Oracle database 11g (11.2.0.1.0) and QGIS 2.16.0


texts: How to update projection of layers in Geoserver using gsconfig

texts:
Im using following code to find projection of layer

from geoserver.catalog import Catalog
cat=Catalog("http:// localhost:8080/geoserver/rest")
myResource = cat.get_resources(store='storeName', workspace='workspaceName')
layers = cat.get_layers(myResource)
for layer in layers:
    layer.projection

I found that few layers are having projection "None". How can I update those layers projection with "EPSG:4326" using gs config.


texts:
I'm trying to clip a soil map layer and am getting invalid geometry that affects my clip result. The problem appears to be a self-intersecting polygon.

I've tried using the 'check geometries', and can identify the offending feature point, but this tool gives no solution.

Can anyone point me to a walkthrough or advise an easy-to-understand solution?


texts: Creating CON.shp in ArcPy gives ERROR 000732

texts:
I am creating  bunch of point shapefiles by using watershed boundaries in a loop

import arcpy
from collections import defaultdict

geo_db = 'C:\\CHIA_Routine\\CHIA_Watersheds.gdb'
smis_wq_resultsXY = geo_db + 'smis_wq_resultsXY'
huc12_370 = geo_db + '\\HUC12_370'

with arcpy.da.SearchCursor(huc12_370, ['WatershedD', 'CIAMaster', 'Code3', 'HUC_12']) as cursor:
    n = 1
    for row in cursor:
        watershed_report_dic[row[0]] = [row[1], row[2], row[3], n]
        watershed_shp_dic[row[1]] = [row[0], row[2], row[3], n]
        n = n+1

    for watershed in watershed_shp_dic.keys():
        code3 = watershed_shp_dic[watershed][1]
        code3_prefix = main_folder + '\\' + 'Watershed_DataSets' + '\\' + watershed_shp_dic[watershed][0] + '\\' + code3
        code3_path = main_folder + '\\' + 'Watershed_DataSets' + '\\' + watershed_shp_dic[watershed][0]
        arcpy.Clip_analysis(smis_wq_resultsXY, geo_db + '\\' + code3, code3_prefix + '.shp')

code3 is three letter unique code for each watershed. When code3 is 'CON' it gives me error:

ERROR 000732: Input Table: Dataset C:\CHIA_Routine\Watershed_DataSets\Contrary_Creek_Kentucky_River\CON does not exist or is not supported
Failed to execute (Clip)

If I changed both output shapefilename and watershed code then it works. I thought maybe 'CON' is a reserved word.

Is it possible to select a polygon from a polygon shapefile and export it as CON.shp without getting ERROR 000732?

.


texts: Geometry column naming convention - 'geom' or 'the_geom'?

texts:
I'm starting my first PostGIS project, and in various books and tutorials I've seen the geometry column labelled either 'geom' or 'the_geom'. Is one more conventional than the other?

Furthermore, is there a good reason to use 'geom'/'the_geom' instead of a more descriptive name for the geometry column? (E.g. 'centre_point' for the central point of a feature.) We don't label conventional db columns 'the_int' or 'the_string', so why label geometry columns this way?


texts: Calculating number of residents living within service areas using ArcGIS ModelBuilder?

texts:
I'm trying to build a Map of service areas of several facilities. In addition, I want to calculate the number of residents living within those service areas.

Therefore I have a Layer of Census-Blocks which I want to intersect with the new layer of my service areas. The problem is now that I get a new Layer of Blocks lying within my service areas but I don't get any attributes with information about shape length, shape area or amount of residents. Where is my mistake?

Picture 1: Model to build Service Area and intersect with census blocks to calculate residents within

Picture 2: Service Area

Picture 3: Intersection of Census-Blocks (RED) and Blocks within Service Area (GREEN)


texts: Calling variables within ArcPy functions using Calculate Field tool?

texts:
I have a defined function (ex. def my_function(x,y)).

def my_function(x,y):
  """algorithm"""
my_variable = my_function(filex, filey)

How can I call this variable within the code block section of the following ArcPy statement?

CalculateField_management (in_table, field, expression, {expression_type}, {code_block})

I want to include an if/elif/else statement based on the value of my_variable.  For example:

if field > my_variable:
  return my_variable / 2
elif field < my_variable:
  return my_variable ** 2
else:
  return my_variable

I think I can just include the my_variable = my_function() at the beginning of the code block but my function is huge.  Also, I don't know if CalculateField can see the function outside of the code block.


texts: Difference in destination location between pyproj and geopy

texts:
I am looking at ways to calculate (in Python) the destination location from a point given bearing and range.

Based on the results comparison from the 2 libraries in subject(geopy and pyproj), I noticed that there is an increasing difference in the final output. For instance, at 100 km is roughly of the order of decimeters. This is a minimal example of what I mean:

from __future__ import absolute_import, division, print_function

long_1 = -1.729722
lat_1 = 53.320556
bearing = 96.021667
distance = 124.8  #in km

# using geopy

import geopy
from geopy.distance import VincentyDistance

origin = geopy.Point(lat_1, long_1)
destination = VincentyDistance(kilometers=distance).destination(origin, bearing)
gp_lat_2 = destination.latitude
gp_long_2 = destination.longitude

# using pyproj

from pyproj import Geod
g = Geod(ellps='WGS84')
prj_long_2, prj_lat_2, prj_bearing_2 = g.fwd(long_1, lat_1, bearing, distance*1000)

# results comparison
print("        | pyproj        | geopy")
print("long_2   %.6f     %.6f" % (prj_long_2, gp_long_2))
print("lat_2    %.6f     %.6f" % (prj_lat_2, gp_lat_2))

print("> DELTA pyproj, geopy")
print("delta long_2   %.7f" % (prj_long_2 - gp_long_2))
print("delta lat_2    %.7f" % (prj_lat_2 - gp_lat_2))

I got these results:

        | pyproj        | geopy
long_2   0.127201         0.127199
lat_2    53.188432        53.188432

> DELTA pyproj, geopy
delta long_2   0.0000021
delta lat_2    -0.0000002

My main question is whether I am doing something wrong (both settings should be using WGS84).

If not, is the difference due to different formulas in use (Vincenty for geopy vs. Karney for pyproj)? E.g., the round-off error cited here.


texts:
I have a basic geojson source, which is

{
  "type": "FeatureCollection",
  "features": [
    {
      "type": "Feature",
      "id": "30cde2b6-767d-426d-8f41-c6b0fad78343",
      "geometry": {
        "type": "Point",
        "coordinates": [
          -75.9337031788254,
          38.8970092675465
        ]
      },
      "properties": {
        "Title": "Test Title",
        "DateCreated": "2018-07-10T13:38:59Z"
      }
    }
  ],
  "crs": {
    "properties": {
      "name": "EPSG:4326"
    },
    "type": "name"
  }
}

and want to show it on a simple map like so:

var vectorSource = new ol.source.Vector({
        features: (new ol.format.GeoJSON()).readFeatures(geojsonObject)
      });

      var vectorLayer = new ol.layer.Vector({
        source: vectorSource,
        style: styleFunction
      });

      var map = new ol.Map({
        layers: [
          new ol.layer.Tile({
            source: new ol.source.OSM()
          }),
          vectorLayer
        ],
        target: 'map',
        controls: ol.control.defaults({
          attributionOptions: {
            collapsible: false
          }
        }),
        view: new ol.View({
          center: [0, 0],
          zoom: 2
        })
      });

my actual code can be found here: https://codepen.io/vishnu4/pen/rqoYzp

Why is this displaying the point at 0,0?  I have the projection declared in the geojson, and putting it in the layer seems to have no effect.  what else am i missing?


texts: How to export only some feature classes to xml using ArcObjects?

texts:
I need to export some of the feature classes to an xml file ("export xml workspace") quite often so I'd like to automate the process. I've found an example here how to export a feature dataset to an xml file. I haven't been able to modify the code so that it would export only the feature classes I want (for example only the ones that start with a letter 'A'). I also need to export the tables and relationshipclasses (schema only).

Is this possible at all? I know there is an interface IFeatureClassName, but I don't know how to use it in this case.
Should this be done with arcpy instead?

I'm using ArcGIS 10.0 SP3.


texts:
I'm fairly new to QGIS and I'm trying to use Heatmap.  I should be doing every thing correct. But when I click OK for it to generate a heatmap, there is no Heatmap at all (nothing changes).

For the points I created a New Shapefile and then I used "Add Feature" button to add many points to the map.  The only thing I can think might be causing the problem is I'm adding points to map the wrong way perhaps.

I have the latest version of QGIS 2.8.1

When i zoom out to the world view I can see a giant 1/4 of a heat map circle at the bottom right of the map. So i was thinking the map coordinates (CRS) might be different, but they are really set to the same map coordinates (CRS). So I don't know what is going on.

EDIT:
It also was doing this in an older version I had, so last night I download the latest version (maybe should not have) hoping it would help and nothing changed.


texts:
I needed to load a shape in GRASS GIS to use the v.clean.snap function. The .shp was first viewed in QGIS, projected in 3857 or 4326, I don't remeber. After the adjustments in GRASS GIS I exported to shape.
Now when I load the corrected version of my shape in QGIS my features are not in europe but in afrika and whatever CRS I pick it is never right. Did I mess up my file in GRASS? How can I correct that?

The dataset is this one: http://www.unece.org/transport/areas-of-work/transport-statistics/statistics-and-data-online/e-roads/transmainwp6e-roads-census-2005.html


texts:
I have a shapefile containing polygons and a shapefile containing points. Now I would like to select the polygons that contain points.

How can I do this from the command line using ogr2ogr?


texts:
I need to create distance matrix between cities in a country in Python. First step I did was to download nodes with place=city|town|village from OpenStreetMap, then I was able to download street network with desired level (highways, first class ...), so I am able to route between any street-network nodes.

Problem is that place nodes are not part of street network - they are not connected - I can only search for closest street node.
But the bigger problem is how to determinate what are the city's neighbors. I don't need to know distance from every city to every city, but only direct connection between them.

Is there some way to filter the simplified network directly? For example through overpass API?


texts: What is the best ETL tool within ArcGIS standard to bring CSV to an existing schema with domains and other constraints?

texts:
We will be bringing some GPS related data from a business partner to our existing schema. We have provided them the sample data and data model from ArcGIS Diagrammer. There are bunch of data import options available in ArcMap and ArcCatalog. I am looking for an efficient approach to finish the ETL process with proper messaging about the exceptions of domains, datatype and mandatory/non-mandatory - so that any novice user of ArcMap can follow it. I don't want any third party tool or code development, but I can build a model, if there are 2-3 tools involved. I am new to data interoperability - has anyone used it in similar situation.
-Thanks


texts: Create individual (temporary) shapefiles for each object in source shapefile

texts:
Is there a way in QGIS or opensource GIS to split a shapefile containing point objects into multiple temporary shapefiles each containing 1 object (i.e. each point to a new scratch shapefile)?

Needed to insert in a routine where I will call each temporary individual point shapefile, as the routine isn't set up to cycle through multiple points within the original source file.


texts:
I'm trying to import OSM to PostgreSQL using the osm2pgsql tool.

osm2pgsql -cGs -d osm -U postgres -W postgres -H localhost -S /usr/local/share/osm2pgsql/default.style /Users/asmaa/SIG/Projet_SIG/campus.osm

I am getting this output with an error:

Setting up table: planet_osm_point
NOTICE:  table "planet_osm_point" does not exist, skipping
NOTICE:  table "planet_osm_point_tmp" does not exist, skipping
Setting up table: planet_osm_line
NOTICE:  table "planet_osm_line" does not exist, skipping
NOTICE:  table "planet_osm_line_tmp" does not exist, skipping
Setting up table: planet_osm_polygon
NOTICE:  table "planet_osm_polygon" does not exist, skipping
NOTICE:  table "planet_osm_polygon_tmp" does not exist, skipping
Setting up table: planet_osm_roads
NOTICE:  table "planet_osm_roads" does not exist, skipping
NOTICE:  table "planet_osm_roads_tmp" does not exist, skipping
Using built-in tag processing pipeline
Allocating memory for dense node cache
Allocating dense node cache in one big chunk
Allocating memory for sparse node cache
Sharing dense sparse
Node-cache: cache=800MB, maxblocks=102400*8192, allocation method=11
Mid: pgsql, scale=100 cache=800
Setting up table: planet_osm_nodes
NOTICE:  table "planet_osm_nodes" does not exist, skipping
Setting up table: planet_osm_ways
NOTICE:  table "planet_osm_ways" does not exist, skipping
Setting up table: planet_osm_rels
NOTICE:  table "planet_osm_rels" does not exist, skipping

***Reading in file: postgres
error while opening file postgres***

Any suggestions?


texts: Landsat Collection 1 Level-1 or Landsat Collection 1 Higher Level?

texts:
I want to use NDWI and NDVI index in my study for a time-series analysis using all products available (Landsat 5,7 and 8).

Can I use the product Landsat Collection 1 Level 1 for my analysis or
should I use Landsat Collection 1 Higher Level?

My doubt is because the calibration that use Landsat Collection 1 Level 1 is regardless of the sensor and the indexes that I want to use are dimensionless, however, always in this kind of studies we use surface reflectance (Landsat Collection 1 Higher Level).


texts: QGIS Map Composer Label Unit Multiplier - Feet and miles?

texts:
I have a map in the WGS84 (EPS:4326) coordinate system and I am putting a scale bar on the map with Scalebar units of Feet and label in miles.

When I select Feet from the Scalebar units dropdown menu, the automatically selected Label unit multiplier is 5419.95 and the Label for units is miles.

There are 5280 feet in a mile, not 5419.95.

Why does QGIS automatically use this value when feet is selected (and miles is automatically filled in for the label)?

Note that if I change to a different coordinate system (e.g. NAD83), the 5419.95 conversion remains.


texts: Expression to create and populate column based on filtered fields

texts:
I wanted to create a column and populate it with data based on data from other fields.

If column_1 contains fields which has similar string (eg. field_1 = "Country: UK", field_2 = "Country: US", field_3 = "City: London"), is there an expression to update fields which ONLY contain part of the string - in this case "Country"?


texts: Panning the map to certain extent javascript API

texts:
I want to limit map extent to the initial extent of the map and limit user from panning more than certain extent.

I tried following but nothing has changed:

map = new Map( "map" , {
        basemap: "gray",
        center: [-85.416, 49.000],
        zoom : 6,
        logo: false,
        sliderStyle: "small"
    });
    var imageParameters = new ImageParameters();
        imageParameters.format = "PNG";

    dojo.connect(map,"onLoad",function(){
        var initExtent = map.extent;

        dojo.create("div",{
          className: "esriSimpleSliderHomeButton",
          onclick: function(){
            map.setExtent(initExtent);
          }
        },dojo.query(".esriSimpleSliderIncrementButton")[0], "after");
      });

    dojo.connect(map, "onExtentChange", function (){
    var initExtent = map.extent;
    var extent = map.extent.getCenter();
    if(initExtent.contains(extent)){}
    else{map.setExtent(initExtent)}
    });


texts:
Trying to write the city name and point coords to a file. The following works with a polygon shapefile but not the point shapefile. I know I'm missing something basic.

import arcpy

# Set up the Environment
arcpy.env.overwriteOutput = True
arcpy.env.workspace = r"C:\temp"

fc = r"C:\Utemp\cities.shp"
fields = ['SHAPE@XY']
output = open(r"C:\temp\cities.txt", 'w')
COFilePath = r"C:\temp\cities.csv"
COFile = open (COFilePath, "w")

# For each row write citiy name & coords
cursor = arcpy.da.SearchCursor(fc, "SHAPE@")
for row in cursor:
    for point in row[0].getPart(0):
        COFile.write(str(point.X) + ", " + str(point.Y) + "\n")
COFile.close()

  Error: TypeError: 'Point' object is not iterable


texts: Using CASE expression in QGIS data defined properties

texts:
I try to use CASE expression in my data defined properties, like this:

CASE WHEN  "ASULKM1999" > 100000 THEN 5 END

ASUKASLKM is population information and "5" will be size of point symbol. However, when I click ok/apply, QGIS will turn this expression to funny, like

CONDITION:
- WHEN ASULKM1999 > 100000 THEN 5

And nothing happens. Sometimes result is "(no root)".

Is this a bug / feature or am I doing something wrong? Using QGIS 2.1 master on Win and Ubuntu.

Rgs

P

Add some snapshots:

Start from here:

Define expression for Size-property:

Expression seems to be allright:

Click ok and then re-opening Data Defined properties:

Expression is quite odd:


texts: How to format an attribute of type esriFieldTypeDate from a FeatureSet Graphic

texts:
Hoping someone can tell me if I'm going about this wrong, but I'm executing a QueryTask / Query on a map service and getting back a FeatureSet as expected. Within the attributes of each Graphic, I have two fields that are of type esriFieldTypeDate and their values are the long integer value representing the date. My problem is I'm not seeing it formatted in any way when it's displayed in my dynamically created DataGridColumn. I'd like to see it in MM/DD/YYYY format.

I was hoping to check the attribute's data type and if it was equal to 'esriFieldTypeDate' - run some formatting / conversion method on it. But, I can't even get the data type of the attribute. There is no spatial querying going on here, this is purely based on attribute data entered on a form in the application.

The last if statement of the sample below is an error - just showing my thought process.
Also, the "resultFields" variable is an array of field names. I use that to allow the field order to configurable.

Should I not be using a QueryTask for this???

function onResult(featureSet:FeatureSet, token:Object = null):void
{
...
    for each (var gr:Graphic in featureSet.features)
    {
        for (var i:int = 0; i < resultFields.length; i++)
        {
            var hg:HGroup = new HGroup();
            var lab:Label = new Label();
            lab.width = 200;
            var dgc:DataGridColumn = resultFields[i] as DataGridColumn;

            lab.text = dgc.headerText + ": ";
            hg.addElement(lab);

            var valLabel:Label = new Label();
            if (gr.attributes[dgc.dataField].type != "esriFieldTypeDate")
            {
                //Not a date - display text value
                valLabel.text = gr.attributes[dgc.dataField].type;
            }
            else
            {
                //Attribute is a date - format it to MM/DD/YYYY
                valLabel.text = gr.attributes[dgc.dataField];
            }

            // ... add value label to group, add group to column, add column to grid etc.


texts: Filter and highlighting feature using Atlas Generation in QGIS

texts:
I have a layer containing 200 features (constituencies). There are six regions and there's a variable called REGION in my shapefile too.

I want to filter on a region (by using query builder), and highlight each constituency within that region and save the map (and get 200 maps).

How can I automate generating these 200 maps? Maybe somehow using Atlas Generation? Or a Python script?

I'm using QGIS 2.16.1 - Nodebo.


texts: Temporal nearest neighbour analysis in QGIS

texts:
How do I conduct a temporal nearest neighbour analysis in QGIS?
Not sure if that's the right terminology, but basically, I have a dataset with ~100,000 house price transactions over a 5 year period. Over this period, there were new train stations being built. I know how to perform a basic nearest neighbour analysis, but I’d like to perform a temporal NNA to determine the nearest train station at the time of the transaction.


texts: Google Maps: find information about roads according to GPS data

texts:
I'd like to know if it is possible, using Google Maps, to get informations about roads, such as:

Name of the road;
If it's a highway or not;
Length of the road;
Altitude of the road, and even better: ranges of altitudes (example : from 15ft to 150ft)
Traffic data, like traffic lights, traffic usage, etc.

The input data I have are a bunch of GPS data (longitudes and latitude coordinates).

So basically, I'd like to find a way to do this :

I can imagine there isn't a software that does this, but maybe if you can give me some clues on how to do such a thing.


texts: Seeking API for Countries, States, Cities and Airports?

texts:
I'm looking for an API or data source (or multiple sources) that provide the following information:

Countries Information: ISO code, Region, Capital, population, languages, currencies and climate.
States Information: name and time zone
Cities Information: name and time zone
Airports Information: IATA, ICAO codes, longitude, latitude, altitude and time zone.
Hotels Information: how many stars, number of rooms, longitude, latitude and type of rooms.

I want this information to be highly accurate and include latitude, longitude and time zones.
also, I don't mind if it was a paid source.
I can't use any source that have an affiliated program

I have checked geonames.org restcountries.eu world bank APIs booking.com they either have missing information or inaccurate


texts: How to display a progress bar in an arcpy python toolbox?

texts:
I'm programming a python toolbox (*.pyt) and want to display a progressbar in arcgis 10.3. I'm using the same piece of code as this one, but I can't see the same dialog box as above.
I can only see the messages updating in the geoprocessing "results" window.
I also tried to add my Pyt to the toolbox window, but this didn't solve this issue.

Is this progressor class only working for python script, and not for python toolbox ?


texts: Access 2000 database via ArcCatalog displays decimals as integers?

texts:
I import an Access 2000 database through ArcCatalog. When I complete the process and check the attribute table, the values come not as decimals as they are in the database but as integer numbers. Why is that?

I have done this in the past successfully but now it is not working properly. Does it have to do with the decimal separator being a dot and not a comma? Does ArcGIS understand . as decimal point?

Process:

Excel file import to Access, save as Access 2000
ArCatalog connect database Microsoft Jet 4.0
open ArcMap and refresh ArcCatalog
take the database sheet and drop it on the table of contents
display xy data
data export data

And this is basically it. ArcGIS shows it as a double.

I am using ArcGIS 10.2 for Desktop.


texts: QGIS. Calculate length of lines in temporary draft layer

texts:
In QGIS, I need to calculate the distance between certain points and save this information as a line shape file, in which the attribute of each element (line) is its length (distance between points).
For this I have created a temporary draft layer of lines to draw the lines and I have obtained this:

Now, how do I calculate the lengths and save them as an attribute of the elements?


texts: Geocoding Multiple Addresses using ArcGis JS API

texts:
I am new to using the ArcGIS JavaScript API. I am currently trying to develop an application but I am having troubles with it.

The application I am developing uses a .csv file with a bunch of locations. The fields for each location are the address, zipcode, county, city, and state. It then geocodes each location and returns another .csv file with all the the address, zipcode, county, city, state, latitude, and longitude.

Using a .csv file to geocode the addresses isn't the important part, I could convert the .csv file into something else if necessary, what I am having trouble with is with the geocoding.

I looked at the samples online and the tutorials but did not find any that would geocode multiple addresses and that would return a file instead of just placing a mark on the map.

Could anyone please give me an idea of how this could be accomplished using the API? Tutorials/Samples would be extremely helpful.


texts: How to set conditional imports for a PyQGis script running inner or outer of QGIS?

texts:
Related to Import error for qgis.core when running OSGeo4w shell script

I try to detect if my script is run from a QGIS instance (through the plugin) or from a standard python console. This detection will enable conditional imports of the required librairies which are different from a case to the other.

For conditional imports, I found here

try:
    import module
except ImportError:
    import otherModule as module

and tried to adapt it for my use case but I can't raise a convenient error:

E.G.

try:
    import processing

except:
    # Load required libraries to run from python (!Unstable!)
    # See https://gis.stackexchange.com/questions/129915/cannot-run-standalone-qgis-script
    # for any improvements
    import os, sys, glob

    # Prepare the environment
    ...
    # See https://gis.stackexchange.com/questions/129915/cannot-run-standalone-qgis-script

Returns me following error :
QPixmap: Must construct a QApplication before a QPaintDevice

and does not instead enters the except clause what is sad.


texts: Extracting GPS coordinates from a large number of photos and then bring them into ArcMap as a shapefile

texts:
Long story short I'm doing surveying work and there is some debate on my team about whether or not pictures were taken in certain areas or if they were mislabeled and taken in the wrong sample area.

We use a GPS enabled camera and it is easy to check each photo individually by copying the lat&long from the exif data and loading it into GIS and seeing if it is within the Sample polygon that we were supposed to inspect.

Is there a way to do this process for hundreds of images and to make them into a shapefile with all the points?

I have tried searching for a program that can do this but have come up short. It seems like a pretty basic process, but I don't know much in the way of programming or anything like that. It would also be good to export them to an excel file with each row containing the filename of the image and the Lat&Long coordinates.


texts: What is the best way to create polygons on a map and download files associated with the polygon?

texts:
My client wants to build a tourist site and for that he wants a backoffice for the map and the ability to draw polygons on the map (that user can see).  When the user clicks on a polygon, he gets a list of files he can download (pdf , images etc...)

I'm using OpenLayers as a map client and haven't decided on the server side yet (probably GeoServer).

What is the best way to plan this? To where should I save the polygons and file list (shp file ? database ?)

thx


texts: Resources for creating a matrix of relationships between street segments

texts:
We are looking at examining traffic in a city using data coded to TMCs (which are long directional links, primarily major arterials and highways, where each link is the length of multiple blocks). In order to investigate impacts of incidents and bottlenecks on adjacent streets, I thought of constructing a table like the following:

{street1, street2, relationship, distance}

where the relationship would be one of [upstream, downstream, intersect, adjacent].

Any suggestions or resources for this endeavor? Preferably compatible with one of QGIS, PostgreSQL, Python. Does not have to be exclusive to TMCs, any form of street network.


texts: Equivalent method to arcpy.Raster(), but for a TIN

texts:
I am creating a Python script in a toolbox in ArcGIS 10.0. When the input dataset is a raster, I can use:

myRaster = arcpy.Raster('elevation')

What should I do when the input dataset is a TIN? Is there an equivalent for TIN datasets?

I only have found "CreateTin_3d", but I dont want to create another TIN.


texts: PostGIS St_snapToGrid doesn't work on Z-value

texts:
As the title says, it seems ST_SnapToGrid doesn't work on Z-values.

Example:

SELECT
id, ST_AsText(ST_Snaptogrid(ST_Transform(pos,2056),0.01)) as pos,
FROM tablename LIMIT 2;

Give me 2 WKT geometry as expected, X,Y are rounded, but Z-value comes with it's 11 or 12 decimals.

Here's my PostGIS version as given by SELECT PostGIS_full_version();:

"POSTGIS="2.2.1 r14555" GEOS="3.5.0-CAPI-1.9.0 r4084" SFCGAL="1.2.2" PROJ="Rel. 4.9.2, 08 September 2015" GDAL="GDAL 1.11.3, released 2015/09/16" LIBXML="2.9.3" LIBJSON="0.11.99" TOPOLOGY RASTER"


texts: Using IDW in ArcPy script gives ERROR 000824

texts:
Using ArcMap 10.4 and Python 2.

On ArcMap, I have the 3D Analyst and Spatial Analyst extensions checked out. I manually ran the IDW tool on a points feature class, and I was successfully able to create a raster.

After copying/pasting this python snippet into my python code, I am getting an error saying:

  ERROR 000824: The tool is not licensed. Failed to execute (Idw).

My code is as follows:

import arcpy
rock_points = r"C:\User\Path\Geodatabase.gdb\RockPoints"
idw_path = r"C:\User\Path\Geodatabase.gdb\IDW"
arcpy.CheckOutExtension("Spatial")
arcpy.gp.Idw_sa(rock_points, "Depth", idw_path, "0.32")

I ran this code in the python window of the ArcMap app, and I was successfully able to create the raster again, but when putting the exact same code into a python script, I'm encountered with the above error.


texts: get style item name of one marker symbol

texts:
i used next code to load a symbol with specific style..
now i need to do the opposite operation..
i mean, i have one symbol and i want to know it's style.

is it possible?

private ISymbol LoadStyleSymbol()
{
    try
    {
        IStyleGallery styleGallery = new ServerStyleGalleryClass();
        IStyleGalleryStorage styleStorage = styleGallery as IStyleGalleryStorage;
        styleStorage.TargetFile = @"C:\Program Files\ArcGIS\Engine10.0\Styles\3D Billboards.ServerStyle";
        styleStorage.AddFile(@"C:\Program Files\ArcGIS\Engine10.0\Styles\3D Billboards.ServerStyle");
        IEnumStyleGalleryItem enumStyleGalleryItem = styleGallery.Items["Marker Symbols", @"C:\Program Files\ArcGIS\Engine10.0\Styles\3D Billboards.ServerStyle", ""];
        enumStyleGalleryItem.Reset();
        IStyleGalleryItem styleItem = enumStyleGalleryItem.Next();

        while (styleItem != null)
        {
            if (styleItem.Name == "Red Pushpin 4") break;
            styleItem = enumStyleGalleryItem.Next();
        }

        ISymbol pSymbol = styleItem.Item as ISymbol;

        return pSymbol;
    }
    catch (Exception ex)
    {
        MessageBox.Show("error 8, " + ex.Message);
        return null;
    }
}


texts: Carto, Import API: Create/Delete datasets

texts:
I've been doing a lot of things lately with Carto, but I haven't been able to find the way to create or upload and delete datasets from my user.

In the API I just can see how to add a dataset to my user, but I cannot eliminate it. I need to do that instead of using the syncing option as i don't have available a CDN or similar and I'm using 8Gb of space in total.


texts: Sampling random cells from certain areas of a cropped raster?

texts:
I have a raster cropped to the extent of some glaciers and divided it into three masks to address inner glacial cells, border cells and non-glacier cells adjacent to the glacier. Data can be found here: https://drive.google.com/open?id=1ERFdsqDGLH1a_FbxwawE_gPm0Au0Q9vT

EDIT: Data preparation was as follows

glaciers <- readOGR("***", "Glacier_Clip_Polygon")
DEM <- raster("***Glacier_Clip1.tif")
glaciers_sf <- st_as_sf(glaciers)
glaciers_comb <- st_union(glaciers_sf)
glaciers_comb <- st_cast(glaciers_comb, "POLYGON")
glaciers_large <- glaciers_comb[st_area(glaciers_comb) > units::set_units(50000,"m^2")]
glaciers_large_sp <- as(glaciers_large, Class =  "Spatial")
glaciers_ras <- mask(DEM,glaciers_large_sp)

To divide the masked DEM into my three classes I used the function proposed by @obrl_soil:

outer_cells <- function(x = NULL, diag = TRUE) {
  focal(x, w = matrix(1, ncol = 3, nrow = 3),
        fun = function(win) {
          if(!is.na(win[5])) {
            return(NA_integer_)
          }

          if(diag == TRUE) {
            if(any(!is.na(win[c(-5)]))) {
              0L
            } else {
              NA_integer_
            }
          } else {
            if(any(!is.na(win[c(2,4,6,8)]))) {
              0L
            } else {
              NA_integer_
            }
          }

        }
  )
}

## presence/absence layer with no DEM-information
g_fp <- glaciers_ras
g_fp[!is.na(getValues(g_fp))] <- 1L

## glacier-adjacent cells
y_msk <- outer_cells(g_fp)
glacier_adj_cells <- y_msk + DEM
plot(glacier_adj_cells)

## glacier marginal cells
rcl <- data.frame('f' = c(1L, NA_integer_), 't' = c(NA_integer_, 1L))
g_fp_inv <- subs(g_fp, rcl)
g_msk <- outer_cells(g_fp_inv)
glacier_marg_cells <- g_msk + DEM
plot(glacier_marg_cells)

## inner glacier cells
lb_msk <- subs(g_msk, rcl)
lb_msk <- lb_msk + g_fp
glacier_inn_cells <- DEM + lb_msk
plot(glacier_inn_cells)

Now I have three raster layers with my three classes for various glaciers. From the inner cells mask of each glacier I now want to draw random cells (in order to lay a buffer and calculate the slope and do some more stuff). The ratio between the random cells and the number of glacier marginal cells should be approx. 0.3.

I want to draw the same ratio from each of my glaciers, so I would have to address them separately. Is there a possibility of doing this in a raster object, maybe by using the polygons as markers of some kind?

I hope, this edit makes my problem a bit clearer.


texts: ArcToolbox Create Replica tool advanced options

texts:
I am trying to create a two-way Replica of my SDE database to another SDE database in ArcGIS 10.0

If I use the Create Replica wizard from the 'Distributed Geodatabase' toolbar in ArcMap I am able to replicate for single tables both schema and data.
This is possible by setting option "Check out" option for this table to "All Records" value in Advanced Create Replica Options step.
BTW default for this option is "Schema only".

When I tried the Create Replica tool from ArcToolbox I haven't found such an option.

Is there any possibility to create replica for tables for both sctructure and data via ArcToolbox?


texts: How to add Opencarto to a Geotools project?

texts:
I've installed geotools through maven (and m2ecplise) but I'm not able to add opencarto library to my project. (As you can guess, I'm new to java)

I tried by adding opencarto's main folder as library through the Project panel, but still, maven is not proposing me opencarto as dependency when I try to find it.

Any idea of what I'm doing wrong?

(this is probably a question for you Julien. I sent you a mail through sourceforge. I'm not sure you got it ... anyway, it's a good thing to make noise around opencarto since it seems very interesting!)


texts: QGIS change raster tiles properties at once

texts: I have a set of several raster tiles and I need to change the properties of each one, by setting up a transparency for instance. I have found no other solution but changing properties of each one at a time... When I select all tiles by pressing ALT + click, the properties of the fisrt one only are modified. Is there a way to perform that?


texts: Splitting proportional/graduated symbols using QGIS

texts:
I am currently trying to make a split graduated/proportional symbol map on Classical music composers by country and musical era and I am having trouble with making the symbol in QGIS.
I had to make a .csv-table consisting of the totals by country and by era, and I am able to join my this table to a shape file of European countries, but the desired column (totals, representing the total composers of each country from the entire span of all musical eras) from the joined table I want to show represent with the graduated symbol does not show up in the drop down bar on either program.
I am able to get the split symbol using the pie graph feature on QGIS, but I still would like to be able to represent size changes of each pie graph showing the total composers from all eras.


texts: Fixing artifacts in rendered ocean polygons

texts:
I am trying to build a base map for North America with political boundaries, lakes&rivers, ecoregions, and oceans in QGIS.  I can retrieve shapefiles from opensource websites that provide the data that I want but they are in different projections and aren't displayed correctly (below is the image of the map when the files are imported directly as is to the map).

The North American data are in US National equal area (ESPG 2163) and the ocean shapefile is in WSG 84.  The current map projection is ESPG 2163 (which is better than WSG 84 at my map extent).

Attempts: I tried to reproject the ocean shapefile with "Save as... " and setting CRS to ESPG 2163 to the ocean shapefile but it still displays the oceans backwards in parts of the map at different scales (the polygon covers the landmass instead of the ocean).

Is there a reprojection that I can use to make the ocean shapefile "fit" properly?  Or is it a problem with the ocean shapefile itself that can't handle other projections in which case any suggestions for open ocean shapefiles?

Details on files:

Ocean shapefile:
WSG 84 (+proj=longlat +datum=WGS84 +no_defs)
http://www.naturalearthdata.com/downloads/10m-physical-vectors/10m-ocean/
(note: the ne_10m_lakes file from the same source does not have a wonky projection in the map)

North America Spatial data (shapefiles with political boundaries, lakes&rivers, ecoregions):
ESPG 2163 (+proj=laea +lat_0=45 +lon_0=-100 +x_0=0 +y_0=0 +a=6370997 +b=6370997 +units=m +no_defs)
https://www.sciencebase.gov/catalog/item/4fb545b6e4b04cb937751d3d


texts: How to set up graduated colors for a layer using ArcObjects?

texts:
I have to render a feature layer with a graduated colors renderer, using the ArcGIS API for Java (or VB, .NET).

As the feature layer data is created dynamically using webservices, I don't in advance the range of values, so I need to get the classification to be generated using the natural breaks algorithm.

Here is what I need to set:

In addition, I need to get the largest values overlaying the smallest:

I need to get n classes, using a AlgorithmicColorRamp from one given color to another color and symbolized by a simple line with a given width.

I tried to write this method, called after the layer data source is set:

    public static void applyGraduatedColorsLineRenderer(IGeoFeatureLayer geoFeatureLayer, String field, int count, IColor fromColor, IColor toColor, double width) throws UnknownHostException, IOException {
    IAlgorithmicColorRamp colorRamp = new AlgorithmicColorRamp();
    colorRamp.setAlgorithm(esriColorRampAlgorithm.esriCIELabAlgorithm);
    colorRamp.setFromColor(fromColor);
    colorRamp.setToColor(toColor);
    colorRamp.setSize(count);
    boolean[] ok = new boolean[1];
    colorRamp.createRamp(ok);

    ClassBreaksRenderer classBreaksRenderer = new ClassBreaksRenderer();
    classBreaksRenderer.setField(field);
    classBreaksRenderer.setBreakCount(count);
    classBreaksRenderer.setSortClassesAscending(false);
    classBreaksRenderer.setFlipSymbols(false);
    classBreaksRenderer.setColorRampByRef(colorRamp);

    for(int i = 0; i < count; i++) {
        SimpleLineSymbol simpleLineSymbol = new SimpleLineSymbol();
        simpleLineSymbol.setColor(fromColor);
        simpleLineSymbol.setStyle(esriSimpleLineStyle.esriSLSSolid);
        simpleLineSymbol.setWidth(width);
        classBreaksRenderer.setSymbol(i, simpleLineSymbol);
    }

    geoFeatureLayer.setRendererByRef(classBreaksRenderer);
}

The problem is that I don't get the right range for my data, and the classification method is set to "manual" instead of "Natural Breaks":

Moreover, I don't know what is the equivalent in ArcObjects of checking "Draw this layer using the symbol levels specified below".


texts: How do I add a marker symbol to an ArcGIS Engine MapControl at a specific coordinate?


texts: Grass 6.4.3 v.rast.stats: Calculate only one value

texts:
I am using the grass (6.4.3!) function v.rast.stats to calculate the mean value of rast values over polygons.

The tool is somewhat slow and I actually only need the mean value.
Is there any possibility/trick to reduce the calculation to the mean value (just like in grass 7.0)?


texts: Importing multiple stacked raster images in R?

texts:
I would like to import all multiple mosaic rasters (each mosaic is composed by 7 layers) from indicated folder into my R. Then access them as single multilayers rasters.

So I did:

# read all mosaics named "mos....img" in R
raster_data<-list.files(path=getwd(), pattern="mos.*.img$")
# read files as rasters
s <- stack(raster_data)
# check my imported rasters p.ex. raster n°8 from "s" raster stack
s[[8]]

and my raster s[[8]] contain only 1 layer, so not the whole mosaic was imported!

nlayers(s[[8]])
[[1]]

If I read each mosaic separately, it works:

# read 1 mosaic (composed by 7 bands)
mosaic1<-brick("mosaic1.img")
# extract one band
band4<-subset(mosaic1, 4)

Why "stack" tool doesn´t import whole mosaics but only one band of the mosaic and how it is possible to arrange it?


texts: Methodology for creating accurate drainage networks (and catchments) from high resolution LiDAR DEM?

texts:
This is not the first time I've stumbled over this issue; it seems that I am unable to generate a correct drainage network model and the resulting catchments from full resolution (1m cells) LiDAR data.

When I generalize the LiDAR dataset, convert it to an integer DEM and fill sinks, all is good and I can easily create what appears to be a very generalized model. However, I would like to produce a detailed site model for a large scale map and this is where I am having problems.

I should point out that most issues occur in flatter areas.

I would like the drainage network to accurately follow the terrain but when I use create the drainage network from an integer DEM input the resulting streams are very general and often "disconnected" in areas where it should not be. The streams are do not even closely follow the natural ridges in the terrain. There are also a lot of "orphan" or "go nowhere" segments. When I use a floating point DEM input, the resulting drainage network is detailed and accurate but very disconnected, clustered and "littered" with orphan streams.

I suspect my problem lies somewhere in the data preparation; integer vs floating point raster DEM input, filling sinks correctly, etc. Or could it be that I have to process the surface data somehow to first create a "hydrologically correct" input DEM?

Can someone describe the correct methodology for creating continuous drainage networks and catchments using high resolution LiDAR?

As it stands I have more success with creating the model from an integer DEM input. This however is not ideal for detailed large scale analysis:

The first attached image is a model produced from an integer DEM input. Several obvious problem areas are circled. Please note that there is actually a stream in what appears to be the main drainage channel. I added a very generalized version of the stream.

EDIT: As I already mentioned I have more success with creating the model from an integer DEM input. The following screen captures illustrate why that is. Even though the integer DEM input has many problems as can be seen above it still produces a drainage network that is less disconnected albeit not conforming to the terrain characteristics. As you can see on the image directly below using a floating point DEM input produces a very disconnected and clustered network full of small orphan segments.

Flow Accumulation raster produced from a floating point DEM

Flow Accumulation raster produced from an integer DEM

As far as i can deduct, both methods yield dramatically different results, both methods are unusable for a detailed model.

EDIT:  I apologize for making this post longer and longer (perhaps I am not expressing myself clearly in English)  To further illustrate the problem with using a floating point DEM for input I am attaching the resulting Stream Link output as well as the resulting watersheds. What I am expecting is a continuous Stream Network and a the entire area covered in basins that all flow into each other.

Stream Link produced from a floating point input DEM:

Watershed basins produced from a floating point input DEM:

Here is an example (nearby area, same data) of where the entire flow direction of a basin is changed due to the use of integer DEM input: Red arrow is the flow direction of the model and blue arrow indicates the direction of the actual flow. (blue lines - actual streams, red network is the LiDAR derived stream network Strahler order)

Link to data: https://www.yousendit.com/download/MEtSOGNVNXZvQnRFQlE9PQ
(Will expire May 13, 2011)


texts:
I recently installed QGIS 2.18.3. I loaded some WFS layers (1.1.0) I have previously loaded using 2.14.9 and noticed the record count had doubled (and one time they trippled but I can't replicate that). The duplication was consitent with a few layers I checked from the same URL.

I tried the same WFS URL using 2.18.0. I have one layer with just one point feature and 2.18.0 returns only one feature however, in other layers 2.18.0 seems to randomly duplicate features and drops other features. For example, 2.18.0 returns 1217 features in a layer with 1220 features - some missing and some duplicated.

I tried a 'select distinct * from table' query using 2.18.3 and get 1186 features returned (2372 without distinct). I can't see any duplicate values but there are obviously some missing.

Are there issues with WFS in 2.18.3?

Further information: I have confirmed that the WFS is being served with Geoserver 2.7.0. I have one layer that has an ID (it is not titled FID) that is unique. The IDs range from 1 to 48 (I'm assuming this meets the requirement of a WFS FID??). If I use ArcGIS 10.4 Interoperability extension, QGIS 2.6.0 or QGIS 2.18.0 then I get the correct number of features. If I use 2.18.3 then I get exactly twice the number of records. It seems the request is made twice and I get one record set appended to another.


texts: Adding a linestring by st_read in Shiny/Leaflet

texts:
I have a large shapefile. For a better perfomance I try to load data with st_read() from sf package. The official documentation states that it should be very easy to integrate the shapefile into Leaflet/Shiny. Nevertheless, it does not work.

Error message: missing value where TRUE/FALSE needed

It works with rgdal, but the performance is not that good. Link

I hope that someone can help me out.

small TestData

library(shiny)
library(leaflet)
library(sf)

shapefile <- st_read("pathToshape")

shinyApp(
  ui <-fluidRow(
    column(8,leafletOutput("map", height="600px"))
  ),

server <- function(input, output, session) {
  output$map <- renderLeaflet({
    leaflet() %>%
      addTiles() %>%
      setView(lng=16.357795000076294, lat=48.194883921677935, zoom = 15) %>%
      addPolylines(data=shapefile, layerId = shapefile$id, group = shapefile$zeitver, color="red", weight=3,opacity=1)
  })
})


texts:
EdI'm a completely new user to QGIS, using version 2.12.

I'm trying to built a map featuring a quite large river. The waterway layer displays all streams in the area, plus the edges of the river, which I'd like to fill in to make the map a bit clearer and have a neat border.

Anyone know how to do this?


texts: Issues calculating VARI from Sentinel L2A

texts:
I have a dataset of Sentinel-2 L2A imagery for which I am calculating the Visual Atmospheric Resistant Index (VARI).
The bands that I am using (B2 = blue B3 = green B4 = Red)
I am working in ArcMap. The bands have been divided by 10,000.0 to convert to float and decompress, I have also set the values to exclude value greater than 1.
The reflectances in each band seem fine, but the calculation result in images that are nearly all of one value, completely black or completely gray. I have done the exact same method with other dates and they come out fine. I have also calculated NDVI with the whole dataset with good results.
Any thoughts or suggestions for these errors?

The ranges of these bands are normal but they are rather dark. Not all of the dates I was having issues with were this dark.


texts: Render raster tiles from a Mapbox GL style

texts:
I'm struggling to put together the appropriate stack for rendering raster tiles (for display of slippy maps) from a given Mapbox GL style.

Style specifications: First there were the Mapnik XML stylesheets, which were then superseded by CartoCSS, which in turn are to be replaced with Mapbox GL styles. That's quite an evolution. It is costly to move from one to the other, so today the best option is to go with Mapbox GL styles, if possible.

Style conversions: To my knowledge, the only automated style conversion available is CartoCSS to Mapnik XML using TileMill.

Map editors: MapBox Studio Classic works with CartoCSS. It is replaced by Mapbox Studio, which instead works with Mapbox GL styles. Maputnik, a FOSS alternative to Mapbox Studio, also works with Mapbox GL styles. (Actually, prior to Mapbox Studio Classic there was TileMill, which also works with CartoCSS, and can export Mapnik XML.)

Obviously Mapbox GL styles are intended to be used with MapBox GL js, which uses WebGL. On new hardware this is very nice indeed. On older hardware, without graphics acceleration, MapBox GL js is not an option, however. Yet it desirable to use one style specification to render both with Mapbox GL js and with the old-fashioned raster tiles.

The first conclusion is that, today, we can safely discard MapBox XML.

This leaves the following options:

Use Mapbox Classic, work with CartoCSS: There is no path to Mapbox GL js rendered maps
Use MapBox Studio or Maputnik with Mapbox GL styles: There is no path to old-fashioned raster tiles

Am I missing something? Ideally it should be possible to render raster tiles with the MapBox GL toolkit, but I don't see how.


texts: Improving IDW and Sample Model using ModelBuilder?

texts:
I want to create a dataset of estimated values using inverse distance weighting. I am starting with a dataset of daily values of measurements for specific coordinates. I created a model (below) which iterates through daily values and creates raster layers of interpolated values (IDW tool). Then, values from a set of specific coordinates are extracted from each raster later and are output as a table (sample tool). Finally, tables are converted to .xls and I will use a different data management software to complete the final merge.

As far as I can tell, this model works. But, it takes some time and I have three decades worth of daily measurements to run this for.

How can I make this model more efficient?

Thoughts I had:

-delete the intermediate raster and info table files after they have been converted to .xls --> i tried selecting "delete intermediate data" but files weren't deleted.

-use collect values, can't figure how this works

-use the batch option to run this model in batches by 1..5..10 years.

-merge in arcmap


texts: Calculate distances between consecutive points in QGIS

texts:
I have a large dataset containing GPS points of various vehicles that are recorded every 10 seconds.

I have added a new id column to these GPS points in order to be able to connect them via the Points to Path tool by using the id field as order field.

car_id     LAT         LON         id
2891       4.899883    52.367300   1
2891       4.900583    52.366983   2
2891       4.901033    52.366900   3
2891       4.901100    52.366867   4
...        ...         ...
2891       4.901572    52.367967   82
5892       4.901794    52.369302   83
5892       4.901749    52.369442   84
5892       4.901829    52.361554   85

The resulting GPS trajectory dataset contains outlying observations (the straight lines that  jump across large distances not along with the road network) which I want to remove. These route points are impossible because it would require acceleration/speed that these vehicles can not achieve.

In order to detect these unrealistically travelled distances, I would like to calculate the distance between each point to the next point in sequence (id) for each unique car_id. This mean that the distance to the end point of each car_id should be zero, since there is no subsequent recorded GPS point (see below). Eventually I want to set a maximum distance depending on how far it is possible for a vehicle to get in 10 seconds (150 meters?) and filter out distances above this threshold.

car_id     LAT         LON         id      DISTANCE
2891       4.899883    52.367300   1       13
2891       4.900583    52.366983   2       2
2891       4.901033    52.366900   3       14
2891       4.901100    52.366867   4       10
...        ...         ...         ...
2891       4.901572    52.367967   82      0
5892       4.901794    52.369302   83      9
5892       4.901749    52.369442   84      12
5892       4.901829    52.361554   85      6

Is this the best way to get rid of these outliers, and achievable in QGIS?


texts: Identify the fraction of raster cells and listing all the cell values in every polygon

texts: In ArcGIS 10 I have a netcdf file that maps the fractions of cropland in the world from 1700 to 2007. I extracted the years that I needed as raster layer format (continuous data from 0 to 1) and now, using a shapefile of the world's countries, I would like to identify every raster's cell inside every country (possibly in a table format where I find id of the cell, country name, fraction cultivated). Is it even possible to do it? Moreover I would like to have also the fraction of cell belonging to a country and the fraction belonging to another one (i.e. in the boundary).


texts: How to access the icon of a feature retrieved from text in a vector-layer?

texts:
I`m using OpenLayers to show POIs on a map. These POIs are loaded from a textfile and I`m using a Vector Layer to get it onto the map - works great.
Now I want to use the cluster strategy to have points clustered to one when the user is zooming out and the points get to close (see this great tutorial). Works, too.  But only for a static image. I want to get the image that is defined in the textfile, but I`m not able to get this.
The text file is in the normal format for OpenLayers.Format.Text. The image is defined in the column named icon and as I told, it is shown when not using clustering.
According to the Source code of OpenLayers.Format.Text.js the value of the column icon should be stored in the externalGraphic-property. But this doesn`t seem to be the case. Here`s my code to try to get the Image:

var clusterStyle = new OpenLayers.Style({
    externalGraphic: "${icon}",
    graphicWidth: 24,
graphicHeight: 24,
    description: "${description}"
}, {
    context: {
        icon: function (feature) {
    if(feature.cluster)
    {
        feature.externalGraphic="/img/clusterimage.png";
    }
    return feature.externalGraphic;
        }...

I`ve also tried other constructs like

feature.cluster[0].attributes['icon']

or

feature.externalGraphic=feature.cluster[0].attributes['icon']

to get the icon, but neither worked.
So where is the icon stored, how to access it? Is there a qay to iterate the attributes-object?
I`m using OpenLayers 2.12


texts: How can I retrieve the features in a layer without adding it to a map?

texts:
I am using the ArcGIS API for JavaScript.

My layer of Point features is exposed through a REST service. I want to retrieve all of these features when the application starts up, then show different subsets of features at different times, depending on user input (such as text search) and client-side calculations (such as clustering calculations). If I just create a FeatureLayer normally, it won't actually retrieve any features before it's added to a map -- its graphics property is just [] when I look at it in Firebug. I tried calling FeatureLayer.refresh, but it just complains that map is undefined.

How do I retrieve the features without having to attach them to a map?


texts: QGIS length & scale dependent curved label visibility

texts:
I have my labels set on curved. Unfortunately quite often happens, that the curved label remains not visible because my line is too short on the map. I can't make my map larger, because I need some minimum area to cover. The example you can see below.

I have two lines, where the curved line is visible for this one, which is long enough. It's not visible for a shorter lines. In this event, I have to toggle between curved and parallel label appearance.
Is it possible to set the parallel label only for the item, which will be too short for the curved label to be visible?


texts: Sum of values based on specific field name or Sum of values based on specific field names from different layers using ArcGIS Desktop

texts:
I am looking for a way to sum the numerical values contained in a layer file based on a specific field name in ArcGIS 10.7

Let's say that I have a polygon shapefile that contains region names and the number of inhabitants. I would like to sum the number of inhabitants per region name to then illustrate population density on the map.

Is there any tool or code in ArcGIS that can generate a shapefile (not a dbase or a table) where starting from the first example (in the image) I can obtain the attribute table as in the final results?

Alternatively, is it possible to obtain the same results starting from 2 different layers?


texts: Cropping borders from images without invalidating georeferencing using ArcGIS Pro

texts:
I am using ArcGIS Pro2.4.
I have a large collection of georeferenced historical aerials.  They all have borders of varying dimensions and varying shades of grey.  These borders have to go somehow.  I'd like to crop out the borders using a raster editing program as this seems simplest and quickest, provided I can preserve the georeferencing.
Is there any way to do this crop that won't mess up the georeferencing?


texts: Making attribute features pop up in a static location

texts:
I am writing an ArcGIS application that has a 20% left pane div and an 80% map div. The left pane div has four tabs: one for layers, one for search, one for what I want to display information, and one for what will eventually become a geoenrichment display.

I have no problem making the popup happen. Everything displays correctly and it looks fine. However, I want the popup to only populate into the left pane.

I've seen the ESRI published tutorial and although it is what I need to do, I'm having trouble making it work.


texts: Rotate the map to GPS heading / direction in ArcGIS using ArcObjects SDK 10 for .NET

texts: Assuming I use ArcGIS or an ArcEngine application on a mobile PC. Is it possible to rotate the map in ArcGIS to the direction the device is travelling using a GPS receiver? If yes, how? I would like to use ArcObjects SDK 10 for .NET to accomplish this task.


texts: Check that a Variable is NOT Equal to Multiple Values

texts:
I am trying to use a check a field against several values (eg: [Name] != 'NONE', and [Name] != 'N/A'). This is for the label expression builder within ArcGIS and I am using it eliminate label for fields with NONE, @ or other invalid fields.

I have tried:

def FindLabel ([Name]):
      if [Name] != 'NONE' or [Name] != 'N/A':
         return None
      elif:
         return [Name]

I also tried separating these out:

def FindLabel ([Name]):
     if [Name] != 'NONE':
        return None
     elif [Name] != 'N/A':
        return None
     elif:
        return [Name]

The above results in the labels being printed for the invalid values still despite the if loops.

Is there a way to check a value against a list and eliminate labels, while printing those for valid values?


texts: Python Addin only works if Python Window open in ArcMap?

texts:
Using ArcGIS 10.2.1 and Python 2.7

My python addin works in ArcMap if I open ArcMap with the Python window open, but if I open ArcMap without the Python window open - I get a "MISSING" tool error.   See code below.

import arcpy
import pythonaddins

class IntersectionSearchCmb(object):
"""Implementation for IntersectionSearch_addin.combobox (ComboBox)"""
def __init__(self):
    arcpy.AddMessage("HERE I AM 1")
    print "Here i am 1"
    cur,row = None, None
    self.street_list = []
    self.street_dict = {}
    self.street_list.append("Battle St at 6th Ave")
    self.street_list.append("6th Ave at Battle St")
    self.street_dict["Battle St at 6th Ave"] = [688855.366,5616889]
    self.street_dict["6th Ave at Battle St"] = [688855.366,5616889]
    self.street_list.sort()
    #self.street = list(set(self.street))
    self.items = self.street_list

    self.editable = True
    self.enabled = True
    self.width = 'WWWWWWWWWWW'
    self.dropdownWidth = 'WWWWWWWWWWWWWWWWWWWWWWWWWW'
def onSelChange(self, selection):
    arcpy.AddMessage("HERE I AM 2")
    print "Here i am 2"
    self.mxd = arcpy.mapping.MapDocument('current')
    self.df = arcpy.mapping.ListDataFrames(self.mxd)[0]

    try:
        x = self.street_dict[selection][0]
        y = self.street_dict[selection][1]
        extent = arcpy.Extent(x,y,x,y)
        self.df.panToExtent(extent)
        self.df.scale = 1500
        arcpy.RefreshActiveView()

    except:

        pass

def onEditChange(self, text):
    pass
def onFocus(self, focused):
    arcpy.AddMessage("HERE I AM 3")
    print "Here i am 3"
    self.mxd = arcpy.mapping.MapDocument('current')
    self.df = arcpy.mapping.ListDataFrames(self.mxd)[0]
    pass
def onEnter(self):
    pass
def refresh(self):
    pass


texts: Making ArcGIS Geocoding prefer matches with street address than matches with zipcode?

texts:
I have around 700,000 address that I need to geocode. I want to use ESRI online service for the biggest portion (I have an organizational account with credits on it) and then put the bad matches in other online services.

Doing some tests, I realized that the ESRI service prefer matching the zipcode than matching the street and number. In my data, I'm certain that the street address is correct but I'm not certain that the zipcode is correct (in fact, a lot of the times is not).

Example (made up to avoid putting someones address in internet):

Original Address: 1250 THOMAS WILSON PKWY, HAMILTON, ONTARIO, J0K OP5, CANADA
Winner Match: J0K OP5, CANADA
Runner Up Match: 1250 THOMAS WILSON PKWY, HAMILTON, ONTARIO, P7W 3G9, CANADA

I really prefer the runner up (and Google Maps too!). Do you know a way to tell ArcMap to prefer matches with street address on it than matches with zipcode?


texts: MapServer very weak performance vs. GeoServer - Why?

texts:
I search for a good WMS server that support PostGIS at the end I found MapServer as the fastest and best choice.

Every performance test say it's faster than GeoServer because of its native compiled code.

Then when we developed our server side using MapScript  (the MapServer extension for PHP).

When we did some benchmark on server and compared it with GeoServer and

Wow! MapServer is 2 to 3 times slower than GeoServer.

We also repeat the experiment with MapServer CGI and results not changed.

Our experiment environments:

Windows 7 64 bit
MapServer version 6.0.3 (MS4W 3.0.6) (MS4W/PHP Mapscript)
GeoServer 2.7
PostGIS 2.1
PostgreSQL 9.4

The PHP code which used to render tiles and used as WMS server:

<?php
if(!isset($_REQUEST['layerId']))
    die('Layer Is Not Valid');

$layerId=(int)$_REQUEST['layerId'];

$map=new mapObj(null);
$map->setFontSet(realpath("lib\\server-side\\fonts.list"));
$map->setConfigOption("MS_ERRORFILE",__DIR__.("\\ms_error.log"));

$map->metadata->set('wms_enable_request','GetCapabilities GetMap GetFeatureInfo');
$map->metadata->set('wms_getmap_formatlist','image/png,png,png8,png24');
$map->setProjection("init=epsg:4326");

$layer = new LayerObj($map);

$layer->set('labelcache', MS_ON);
$layer->setProcessing('LABEL_NO_CLIP=FALSE');
$layer->setProcessing('FORCE2D=YES');

$layer->set('name', "MasterLayer");
$layer->set('status', MS_DEFAULT );

$layer->setConnectionType(MS_POSTGIS);
$layer->set('connection', "user=postgres password=psgpass dbname=map host=127.0.0.1 port=5432");
$layer->set('data', "geodata from geo_data as subquery using SRID=-1 using unique id ");

$layer->set('type', MS_LAYER_POLYGON);

$class = new ClassObj( $layer );

$style = new StyleObj( $class );
$style->color->setRGB(255,0,0);
$style->set('opacity',100);

$request = new Owsrequestobj();
$request->loadparams();
$request->setparameter('REQUEST', 'GetMap');
$request->setparameter('LAYERS', 'MasterLayer');
$request->setparameter('FORMAT', 'png8');

ms_ioinstallstdouttobuffer();
$map->owsdispatch($request);

$contenttype = ms_iostripstdoutbuffercontenttype();

if (!empty($contenttype))
{
    header('Content-type: $contenttype');
    ms_iogetStdoutBufferBytes();
}
else
    echo "Fail to render!";
ms_ioresethandlers();
?>

DataBase Data Projection: EPSG:4326

Does anybody know why GeoServer is much faster when it isn't compiled natively? Am I missed something?


texts: Lines on map when I change projection on QGIS

texts: I changed the projection of my base to WGS 84 PDC but I get these lines on top. This also happens to a few other projections. I am trying to make the Pacific the center of the map.


texts: Why does Python if/else statement only return else condition when trying to reclassify field using ArcMap field calculator?

texts:
Using ArcGIS 10.4. I'm trying to reclassify values from a field (ECOREGION_DIVISION), which is in string format, into new values using the following Python syntax:
def Reclass(ECOREGION_DIVISION):
  if (ECOREGION_DIVISION == "Desert"):
    return "Not forest"
  else:
    return "Forest"

    Reclass(!ECOREGION_DIVISION!)

The result is that all the values in the new_class field are returned as "Forest" (the if statement does not return "Not forest" when value = "Desert"). The new_class field is also string format. Does this syntax seem correct?


texts: Different results using ST_Length and ST_LengthSpheroid in PostGIS

texts:
I have a 3D MULTILINESTRING geometry in my database contains lat/long and EPSG:4326:
'MULTILINESTRING Z ((50.9872438056729 35.0857686001374 906.97,...)

and I calculated the length through two methods:
conn = psycopg2.connect(...)
cur = conn.cursor()

_spheroid = 'SPHEROID["WGS 84",6378137,298.257223563]'
cur.execute("SELECT ST_LengthSpheroid(geom,'{spheroid}')/1000 from ...;".format(spheroid=_spheroid))
#68.030 KM

cur.execute("SELECT ST_Length(geom, false), id from ...;")
#68.001 KM

Where does this difference come from?


texts: Querying over multiple bands at one point

texts:
I have a multiband raster in PostGIS. I can obtain the value for a specific band at a specific point with
SELECT ST_Value(rast, band, ST_SetSRID(ST_MakePoint(x, y), srid))
FROM table;

I have 12 bands in my raster. I would like to efficiently query the values for all 12 bands at one specific point, returned in, say, an ordered array. Is this possible? Or do I have to do the ST_Value query 12 times?


texts: I am struggling to find or import a base map without street labels in QGIS


texts: Using ArcGIS Pro Graticule Ticks and Labels?

texts:
I want to display my Latitude ticks and labels along the left and right of my data frame and the longitudinal ticks and labels along the top and bottom of the data frame. In ArcGIS Desktop this has never been an issue. In ArcGIS Pro 2.1 it seems to want to display the longitudinal ticks and labels on every side.
The image below shows how the longitude labels are replacing some of the latitude ones when I tick 'South' as visible.

Without 'South' visible the latitude labels come back on each side but all the longitude ones are gone along the bottom.

Is there a way to indicate which ones I want visible on each side?


texts: Looking for open source component to create PDF from WMS in a web application

texts:
I want to create a PDF from WMS in a web map application.

Please can anyone recommend a robust component? I'd like free open-source under an relatively business-friendly licence (no GPL, BSD/MIT etc. would be ideal, LGPL would be OK).

Our WMS layers are produced by ncWMS and UNM MapServer. We will also want to overlay some layers from external WMS providers.
EDIT: We are developing our own JavaScript client based on OpenLayers.
EDIT: We have some backend code in Java. We would also consider PHP or other platforms to create the PDF.

MapFish would be excellent but it is GPL.


texts: Sum up same objects with the QGIS composer atlas filter

texts:
I have a line layer.
Atlas shall filter by a number column.
Many of the lines have the same number like this: 500, 500, 500, 520, 530, 530.

Now Atlas makes an own image for any 500, when I try to save my work.
How can I sum up all 500 to get only one image instead of three?
Can I use somehow count_distinct?


texts:
Yes, I know the math behind this is simple, but I'm looking for a way to automate the process.

In ArcGIS 10.2.2. is it possible with the Create Random Points data management tool to generate a set of random points based on a percentage of the area of the constraining polygon? More specifically, I have a 1 meter square sampling plot I want to apply so as to cover 10% of the total area of any constraining polygon.

I can think of multiple ways to accomplish this, but they all involved manual intervention. I'm looking for a generic solution that is automated.


texts: Protect our WMS map tiles to only accessible from our site

texts:
We are using geoserver (2.3) to serve WMS map tiles generated from mapping data we licenced. We use leaflet on the client side to show the map (part of a ruby on rails site running on heroku). Because of the licence of the mapping data, we should not allow anybody to use these map tiles, it should be only accessible by our users when they are using the map on our website.

Is there any way to avoid somebody just taking the WMS URL from our javascript code and using it on another site?
Or bulk downloading all the map tiles?

We can configure the geoserver to require username/password or certificate, but then we would need to expose it to the client side.

Thanks!


texts: Cannot correctly iterate IFeatureCursor (calling ArcObjects from Python)

texts:
I have started recently using some of the ArcObjects inside my Python modules. Having all the useful posts and insights shared by @matt wilkie et al, I was able to get started pretty quickly (installing the comtypes with pip and downloading the 10.2 snippet from Pierssen and changing "10.2" to "10.3" everywhere).

I am trying to iterate IFeatureCursor and get all the features inside a feature class. However, I am getting back only the latest feature (with the highest ObjectID value).

There are 6 features in the feature class hence xrange(6) to keep it simple.

from comtypes.client import GetModule, CreateObject
from snippets102 import GetStandaloneModules, InitStandalone

# First time through, need to import the “StandaloneModules”. Can comment out later.
#GetStandaloneModules()
InitStandalone()

def iterate_features():

    # Get the GDB module
    esriGeodatabase = GetModule(r"C:\Program Files (x86)\ArcGIS\Desktop10.3\com\esriGeoDatabase.olb")
    esriDataSourcesGDB = GetModule(r"C:\Program Files (x86)\ArcGIS\Desktop10.3\com\esriDataSourcesGDB.olb")

    # Create a file geodatabase pointer
    file_gdb_pointer = CreateObject(progid=esriDataSourcesGDB.FileGDBWorkspaceFactory,
                                    interface=esriGeodatabase.IWorkspaceFactory)

    file_gdb = file_gdb_pointer.OpenFromFile(r"C:\GIS\arcobjects\MyData.gdb",hWnd=0)

    #access contents inside gdb
    feature_workspace = file_gdb.QueryInterface(esriGeodatabase.IFeatureWorkspace)

    in_fc = feature_workspace.OpenFeatureClass("Warehouses")

    def enum_features(in_fc):
        """returns pointers to IFeature objects inside Feature Class"""
        cur = in_fc.Search(None,True)
        for i in xrange(6):
            feature_obj = yield cur.NextFeature()

    feats = [feat for feat in enum_features(in_fc)]
    print [f.OID for f in feats]

iterate_features()

The line print [f.OID for f in feats] returns [6, 6, 6, 6, 6, 6].

What am I doing wrong? The same logic with generator/yield (def enum_features()) works fine when iterating feature classes inside the feature dataset.

the feats_OIDs = [feat.OID for feat in enum_features(in_fc)] will give correct results, [1, 2, 3, 4, 5, 6], without me making any modifications to the code. The problem seems to be in that when I create a list of features [feat for feat in enum_features(in_fc)], they all refer to the same feature (because when I explore each of them later, each of them have the same OID).


texts:
I developed a plugin in QGIS and I have attached a function on button clicked but it is calling the function two times and I don't know the reason why.

I am registering click function in plugin run event:

    def run(self):

        """Run method that performs all the real work"""
        AttachBtnClicked(self)

        # show the dialog
        self.dlg.show()
        # Run the dialog event loop
        result = self.dlg.exec_()

def AttachBtnClicked(self):
    self.dlg.btnDataSourceEncodingChangeAllLayers.clicked.connect(lambda: btnDataSourceEncodingChangeAllLayersClicked(self))

while the method is:

   def btnDataSourceEncodingChangeAllLayersClicked(self):
        print 'DataSourceEncodingChanageAllLayers : '# + encodingName
        LayerDataSourceChange()

it runs two times.


texts: Blank error while using load objects in ArcMap 10.2

texts: I'm trying to load data from an excel file to a feature class using load objects. After matching all the fields I get the following blank error and no data is added to the feature class. Does anyone know what this error means?


texts:
I'm new to Carto, and using any of the file types Carto relies on, so treat me as a complete beginner. I'm trying to create a heat map for the number of road casualties in 2015 in England.

I used a shapefile from the ONS which includes the geometrical data of all the different authorities in England, and uploaded that to Carto. This seemed to work, because the map is displaying the borders between each authority

However it wouldn't work with my CSV file, which included the data for the casualties in each authority. So, I inserted the geometrical polygon data from the shapefile for each authority on my CSV file manually, and deleted the shapefile. The above photo is displaying from just using the CSV file with the added polygon data.

Below is a photo of the spreadsheet in Carto.

While the polygon geometry data is working, it's not properly reading the Casualties column, and isn't allowing me to create any analysis from this to make a heatmap. I'm not sure if I'm just not understanding the interface properly, because it is reading the casualties to some extent, shown by the numbers in the right locations below.

What I'm asking, is how do I get these numbers into something more like a heatmap, as shown in the first photo, but without the values between each colour making the majority of the map purple. I can't seem to find a way to change how Carto quantifies the numbers.


texts: Calculate moving average (low pass filter) using fft in GRASS

texts:
I am interested in calculating a moving average for a raster in GRASS/QGIS. In other words, each cell needs to receive as value the average of its current value and the neighboring cells in a square (let's say 3 cells x 3 cells square).

I understand I can achieve that with r.mfilter, but I found this solution is pretty slow for large rasters or bigger filters (larger neighborhoods).

I notice GRASS has i.fft for fast fourier transform, which is super fast. I heard in theory one can calculate moving averages using fft, but I could not find an example for a 2d image, and I cannot wrap my head around the math to figure out on my own.

Does anyone can conceive way to use to calculate moving average with a 3x3 filter, using i.fft?


texts: Split line by polygon where polygons overlap ArcGIS

texts:
I would like to split line features to according to boundaries of individual polygons within a layer.

Where a line crosses two polygons, I would like the line split into two lines

I would normally use the intersect tool to do this, but because my polygons overlap the intersect tool is splitting the line crossing two polygons into 4 lines, rather than 2.

I could re-merge the lines, but not sure if there is a common field I can use to do this?