Friday 14 February 2020

arcpy - Changing output name when exporting data driven pages to JPG?


Is there a way to save the output JPG, changing the output file name to the page name, instead of page number?


I mean changing the script found in this post?


import_path = r"..."   # Path of .mxd
export_path = r"..." # Path of output file
field_name = "Name" # Name of field used to sort DDP


mxd = arcpy.mapping.MapDocument(import_path)
for i in range(1, mxd.dataDrivenPages.pageCount + 1):
mxd.dataDrivenPages.currentPageID = i
row = mxd.dataDrivenPages.pageRow
print row.getValue(field_name)
arcpy.mapping.ExportToJPEG(mxd, export_path + "." + row.getValue(field_name) + ".jpg")
del mxd

Answer



If you define a new variable: pageName = mxd.dataDrivenPages.pageRow.Name You can then include the variable in the export path: str(pageName)


import_path = r"..."   # Path of .mxd

export_path = r"..." # Path of output file
field_name = "Name" # Name of field used to sort DDP

mxd = arcpy.mapping.MapDocument(import_path)
for i in range(1, mxd.dataDrivenPages.pageCount + 1):
mxd.dataDrivenPages.currentPageID = i
row = mxd.dataDrivenPages.pageRow

pageName = mxd.dataDrivenPages.pageRow.Name


print row.getValue(field_name)
arcpy.mapping.ExportToJPEG(mxd, export_path + "." + row.getValue(field_name) + str(pageName) +".jpg")
del mxd

gdalwarp - jpg to geotiff transformation with gdal_translate and gcps fails


I am trying to create a geotiff from a jpeg by inputting gcps with gdal_translate:


gdal_translate -of GTiff -a_srs EPSG:4326 -gcp [1616.0 0 -121.459869532 38.5822533549] -gcp [0 0 -121.460081357 38.5831052365] -gcp [0 1080.0 -121.460807872 38.5829948432] -gcp [1616.0 1080.0 -121.460596039 38.582142963] pict20140910_131103_0.JPG tmp.tif


followed by gdalwarp to set the geotiff projection:


gdalwarp -s_srs EPSG:4326 -t_srs EPSG:4326 tmp.tif geo_pict20140910_131103_0.tif

however the gdalwarp process fails with:



ERROR 1: Failed to compute GCP transform: Transform is not solvable.



The gcp coordinates are the image footprint vertices and were calculated using the image location, altitude, orientation and camera attributes. The gcps are in lat/lon (wgs84 epsg:4326), and the image size is 1616x1080.


Here is the gdalinfo for the file after gcps have been added, but before gdalwarp:


Driver: GTiff/GeoTIFF

Files: tmp.tif
Size is 1616, 1080
Coordinate System is `'
GCP Projection =
GEOGCS["WGS 84",
DATUM["WGS_1984",
SPHEROID["WGS 84",6378137,298.257223563,
AUTHORITY["EPSG","7030"]],
AUTHORITY["EPSG","6326"]],
PRIMEM["Greenwich",0],

UNIT["degree",0.0174532925199433],
AUTHORITY["EPSG","4326"]]
GCP[ 0]: Id=1, Info=
(0,0) -> (-121.459869532,38.5822533549,0)
GCP[ 1]: Id=2, Info=
(0,0) -> (-121.460081357,38.5831052365,0)
GCP[ 2]: Id=3, Info=
(0,1079) -> (-121.460807872,38.5829948432,0)
GCP[ 3]: Id=4, Info=
(0,1079) -> (-121.460596039,38.582142963,0)

Metadata:
AREA_OR_POINT=Area
EXIF_BrightnessValue=(4.275)
EXIF_ColorSpace=1
EXIF_ComponentsConfiguration=0x1 0x2 0x3 00
EXIF_CompressedBitsPerPixel=(3)
EXIF_Contrast=0
EXIF_CustomRendered=0
EXIF_DateTime=2014:09:10 13:11:02
EXIF_DateTimeDigitized=2014:09:10 13:11:02

EXIF_DateTimeOriginal=2014:09:10 13:11:02
EXIF_DigitalZoomRatio=(1)
EXIF_ExifVersion=0230
EXIF_ExposureBiasValue=(0)
EXIF_ExposureMode=0
EXIF_ExposureProgram=4
EXIF_ExposureTime=(0.001)
EXIF_FileSource=0x3
EXIF_Flash=16
EXIF_FlashpixVersion=0100

EXIF_FNumber=(4)
EXIF_FocalLength=(16)
EXIF_FocalLengthIn35mmFilm=24
EXIF_GPSAltitude=(66)
EXIF_GPSAltitudeRef=00
EXIF_GPSImgDirection=(79)
EXIF_GPSImgDirectionRef=T
EXIF_GPSLatitude=(38) (34) (57.45)
EXIF_GPSLatitudeRef=N
EXIF_GPSLongitude=(121) (27) (37.22)

EXIF_GPSLongitudeRef=W
EXIF_GPSVersionID=0x2 0x3 00 00
EXIF_Interoperability_Index=R98
EXIF_Interoperability_Version=0x30 0x31 0x30 0x30
EXIF_ISOSpeedRatings=2000
EXIF_LightSource=0
EXIF_Make=SONY
EXIF_MaxApertureValue=(2.96875)
EXIF_MeteringMode=5
EXIF_Model=NEX-5T

EXIF_Orientation=1
EXIF_PixelXDimension=1616
EXIF_PixelYDimension=1080
EXIF_ResolutionUnit=2
EXIF_Saturation=0
EXIF_SceneCaptureType=0
EXIF_SceneType=0x1
EXIF_Sharpness=0
EXIF_Software=NEX-5T v1.01
EXIF_WhiteBalance=0

EXIF_XResolution=(350)
EXIF_YCbCrPositioning=2
EXIF_YResolution=(350)
Image Structure Metadata:
INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left ( 0.0, 0.0)
Lower Left ( 0.0, 1080.0)
Upper Right ( 1616.0, 0.0)
Lower Right ( 1616.0, 1080.0)

Center ( 808.0, 540.0)
Band 1 Block=1616x1 Type=Byte, ColorInterp=Red
Band 2 Block=1616x1 Type=Byte, ColorInterp=Green
Band 3 Block=1616x1 Type=Byte, ColorInterp=Blue

Answer



The gdal_translate document page http://www.gdal.org/gdal_translate.html may indeed give an impression that the values for Ground Control Points should be closed between brackets



[-gcp pixel line easting northing [elevation]]*



However, that is not the case as close reading reveals and correct syntax for your case is



gdal_translate -of GTiff -a_srs EPSG:4326 -gcp 1616.0 0 -121.459869532 38.5822533549 -gcp 0 0 -121.460081357 38.5831052365 -gcp 0 1080.0 -121.460807872 38.5829948432 -gcp 1616.0 1080.0 -121.460596039 38.582142963 pict20140910_131103_0.JPG tmp.tif

That gdal_translate accepts the input with brackets without sending an error message and makes on invalid output as a result may be considered as a bug. Write a mail to gdal-dev mailing list and ask an opinion.


A hint for the future: The best format to use for the interim result is GDAL Virtual raster .VRT http://www.gdal.org/gdal_vrttut.html. It is a small text file that has only a reference to the physical image file. Because image data is not copied you will get the result faster and save disk space. In the next step you just warp the VRT file.


gdalwarp -s_srs EPSG:4326 -t_srs EPSG:4326 tmp.vrt geo_pict20140910_131103_0.tif

data - Seeking shapefile for postal codes/zip codes of India?


I have been searching for free/paid (preferably free) sources for postal code shapefiles for India.


Does anyone know where I can find one?





How to prevent small polygon WMS features to be hidden when zoomed out (GeoServer + OpenLayers)


Currently I'm developing a Web-GIS app where I use GeoServer 2.61 and OpenLayers 2. The map data is served from PostgreSQL+PostGIS via GeoServer (WMS).


I'm having a problem when I'm viewing polygon map. Some maps have a very little polygon feature, which do not show at certain zoom level when I zoom out (because the feature will be too small). Is there any workaround to prevent this to happen? I need to see ALL the polygon features at the lowest zoom level.



Answer



Polygon layers by default do not stay the same size when zooming further out because they represent an area. The only aspect that you will see for smaller polygons when zoomed out is just the outline. If you still want to represent the feature when zoomed far out one option would be to add a within the sld that will turn on at x scale (defining MIN and MAX scales for the polygon and point symbolizers). Refer to this Q/A for how to set this up:


Use SLD to conditionally display points or polygons based on zoom level?



arcgis desktop - Setting single output location for multiple files in ModelBuilder?


I have created a model in ModelBuilder. The first step is for the user to specify the location of the File Geodatabase that all the outputs will be saved in. What I want is for all the outputs (from the different tools within the model) to be saved in that FGDB; however, it could be called anything and located anywhere.


So how do I cause the output to be saved in the location that is specified in the first step?



Answer



This page on Esri's site should give you all the information you need to do this within ModelBuilder. Essentially you create a variable for the output folder/geodatabase -- which can be user-generated or hardcoded -- and then call it in the other tools by its name, surrounded by % symbols.


arcmap - Inserting newlines in Rectangle Text elements via ArcPy causes overlap?


I ran across an issue the other day when I tried to use ArcPy's mapping module to edit rectangle text elements with newlines (\n) in an ArcMap document. Here's what the output looked like:


enter image description here


Here's the code I used to generate that output. The first column are rectangle text elements Text1, Text2, Text3 going down; the second column are "plain" text elements Text4, Text5, and Text6 going down.


import os
import arcpy

HomeDir = r"C:\Desktop"
arcpy.env.workspace = HomeDir


CurrentMXD = arcpy.mapping.MapDocument(r"C:\Desktop\TextTest.mxd")
OutputFilename = r"C:\Desktop\TextTest.pdf"
if os.path.exists(OutputFilename):
os.remove(OutputFilename)

for TextElement in arcpy.mapping.ListLayoutElements(CurrentMXD, "TEXT_ELEMENT"):
TextElementName = TextElement.name

String1 = "The quick brown fox jumped over the lazy dog.\nShe sells sea shells by the sea shore."
String2 = "The quick brown fox \njumped over the lazy dog.\nShe sells sea shells by the sea shore."

String3 = "The quick brown fox jumped \nover the lazy dog.\nShe sells sea shells by the sea shore."

if TextElementName == "Text1":
TextElement.text = String1
if TextElementName == "Text2":
TextElement.text = String2
if TextElementName == "Text3":
TextElement.text = String3
if TextElementName == "Text4":
TextElement.text = String1

if TextElementName == "Text5":
TextElement.text = String2
if TextElementName == "Text6":
TextElement.text = String3

arcpy.mapping.ExportToPDF(CurrentMXD, OutputFilename)

So far, it looks like the presence of the messed up text depends on whether the line is longer enough to wrap, and whether the line before the newline is longer than the line after the newline.


Any ideas about what could be going wrong? Is there a workaround? I could use plain text elements and worry about wrapping lines using Python, but I'm hoping I can figure something out.



Answer




I ran into this as well. It's because ArcGIS requires Windows line endings which are both a carriage return and a line feed. Bit of a pain. Fortunately it's easy to get around. In Python by instead of just \n (which is linefeed - see the Python docs for more if you're keen), use \r\n.


Thursday 13 February 2020

dem - Calculating cells elevation (height) above nearest stream cell using ArcGIS Spatial Analyst


I have searched for some ideas on how to formulate this analysis, but have come up short. I imagine there is a pretty straight forward answer; any help would be greatly appreciates.


Data: DEM, 1/9th arc second (~10-meter), 239,891,470 cells, projected; Streams, National Hydrology Data set, high resolutions, 9,470 features vector; Streams, Same as above, but converted to a raster


Problem: Need to create a raster of the same resolution and extent as the DEM where each cells's value is its height (in meters) above the nearest stream (line or cell).


I am working with ArcGIS 10.1 with spatial analyst. I would prefer to run this analysis in Python, but can do it in R if it is easier.


Edit to clarify: each cell in my stream raster has the value of the elevation at that point. All cells that are not within a stream have a value of zero. For every cell in the extent of the DEM, I am looking to find the nearest stream cell and subtract the value of that cell (the elevation) from the elevation of the cell being analyzed. Calculate the difference in elevation between each cell and the nearest cell representing a stream.



Answer



@whuber has it right when they commented:



If "nearest" means in terms of distance on the map, then just subtract the Euclidean allocation of the stream elevations from the DEM. The procedure is illustrated in an answer to a related question.




First, get the euclidean allocation. Second, use the minus operation to get the difference between the two.


Can be easily done in python.


Topology in QGIS shapefile, GRASS and v.clean (nested polygons)


Following this: Nested polygons in QGIS overlap and don't show the real colour gradient


I applied the GRASS function v.clean -> break as suggested by @ahmadhanb and the shapefile now shows correctly all the nested basins without using the Feature blending mode = Multiply.


However I noticed that in the attribute table the v.clean function created ~2 times more variables than the original shapefile (e.g. original shapefile 250 polygons, after v.clean 500 polygons). There are then duplicated polygons that eventually create issues when doing a Join.


How can I remove the duplicated polygons created by v.clean and keep ONLY the 'good' ones that appear in the QGIS environment when selected?


I tried to apply rmdupl within v.clean but the shapefile did not change at all!




When to use ArcPy over ArcMap GUI?


I've recently began trying out basic geoprocessing tools by using arcpy (clipping layers, buffering, etc).


However, this seems slower than using the ArcMap GUI.



Are there any advantages to using ArcPy for simple tasks that can otherwise be done directly through the GUI?



Answer



As you have noted the geoprocessing tools are slower than their ArcMap GUI equivalents, if the latter exist.


My superficial understanding is that this is because ArcMap already has all the functions needed by the ArcMap tool loaded in memory, whereas geoprocessing tools load that functionality before starting to run.


The compelling use case for geoprocessing tools comes when you want to string automation operations together using ArcPy, ModelBuilder, etc.


If you are already in ArcMap, and if there is a tool available in its GUI that does what you want, then I would recommend using that.


PostGIS point distance from polygon not 0 when on edge



I have a pretty simple problem:


SELECT * FROM ST_Distance(ST_GeographyFromText('POLYGON((70 -40,70 -39,71 -39,71 -40,70 -40))'), ST_GeographyFromText('POINT(70.48 -39)'), False)

returns 118.45656881 Shouldn't it be 0? Isn't the point on the edge of the polygon? Why does


SELECT * FROM ST_Distance(ST_GeographyFromText('POLYGON((70 -40,70 -39,71 -39,71 -40,70 -40))'), ST_GeographyFromText('POINT(70.48 -40)'), False)

return 0?


SELECT * FROM ST_Distance(ST_GeographyFromText('POINT(70 -39.5)'), ST_GeographyFromText('POINT(70 -39)'));
SELECT * FROM ST_Distance(ST_GeographyFromText('POINT(70 -39.5)'), ST_GeographyFromText('POINT(70 -40)'));


Why does this return two different (although only slightly) numbers?



Answer



The video explains it very well. If you take the edge in question and segmentize it


 SELECT st_asewkt( st_segmentize(geog,10000)) FROM ST_GeographyFromText('LINESTRING (70 -39,71 -39)') geog

you going to have the following :


enter image description here


the point is the the point you are asking about. The line represents the edge. North is up.


(distance measurement done by hand,no snapping)


As for the 2nd query that returns 0? It does because hte point is actually inside your polygon.



To be more clear the intersection between the two geometries is empty because they simply don't overlap.


select st_asewkt(st_intersection(ST_GeographyFromText('POLYGON((70 -40,70 -39,71 -39,71 -40,70 -40))'), ST_GeographyFromText('POINT(70.48 -39)')))

st_asewkt
--------------
SRID=4326;GEOMETRYCOLLECTION EMPTY

arcgis desktop - ESRI Basemap No High Resolution



Did ESRI drop their high res. basemaps. I made a map in November and today I needed to make a small adjustment to it. The high resolution basemap is gone and my map looks bad. I have attached the 2 maps to show the vast difference. Has ESRI done something or is there a setting that I am missing? I am using ArcGIS 10.4.1 for desktop.


HighResMap


LowResMap




postgis - Finding all polygons inside of a polygon


How to find all the polygons inside of a polygon? I tried ST_Intersects and ST_Contains but they are booleans variables. I have a list of polygons in WKT format, so I am trying to pull all the polygons given a region in their WKT format.


Does anyone have an idea how to approach this?


select geom from polygon where 
ST_Intersects(ST_GeographyFromText

('POLYGON((210000 2400000, 300000 2300000, 330000 3708400, 210000 2400000))'),
polygon.geom);

Answer



You are missing ST_AsText, as @WKT mentions above. This should work for you:


select ST_AsText(geom) from polygon
where ST_Intersects(ST_GeographyFromText
('POLYGON((210000 2400000, 300000 2300000, 330000 3708400, 210000 2400000))'),
polygon.geom);

arcgis desktop - Setting in memory workspace in ArcPy?


Is it possible to set env.workspace = "in_memory" in ArcPy using both ArcGIS Pro and the ArcGIS 10.2.2 (or 10.3) architectures?


What I am trying to do is get the output of a snap pour point operation written to memory as opposed to disk.


I realize I could write to disk then bring it into memory but this would not help. I am trying to optimize a series of tasks minimizing the use of writing to physical media the output of a process that is merely required as an input to the next process.



Answer



I'm going to throw an answer on here because both answers thus far aren't 100% correct.



There are 2 items which can vary from tool to tool.



  1. if it honors the workspace environment (this item is always documented on the tool help page)

  2. if it can make use of the in_memory workspace (this item may not be explicitly documented. You're more likely to see a note if it DOES NOT support in_memory)


To simply answer the "can you set the environment workspace to in_memory". The answer is YES.


>>> import arcpy
>>> arcpy.env.workspace = r"in_memory"
>>> arcpy.CopyFeatures_management(r"c:\temp\foo.shp", "myinmemoutput")


>>> arcpy.Exists("myinmemoutput")
True

Snap Pour Point does honor the workspace environment per it's documentation and explained Python samples. And a test shows you can write output to in_memory and work with that variable reference...to put into another tool, or save the result


>>> import arcpy
>>> arcpy.env.workspace = r"in_memory"
>>> arcpy.CheckOutExtension("SPATIAL")
u'CheckedOut'
>>> snapOut = arcpy.sa.SnapPourPoint("e:/gpservices101/hydro/US30m/test.gdb/sourcepoint", "e:/gpservices101/hydro/US30m/Region08a/Input/elev_cm", 1,"PourPtID")
>>> snapOut

in_memory\SnapPou_sour1
>>> arcpy.Exists(snapOut)
True
>>> snapOut.save(r"c:\temp\todisk.tif")
>>> arcpy.Exists(r"c:\temp\todisk.tif")
True

python - Formula for coordinates (lat, lon), azimuth and distance.


I would like to know how to find a set of coordinates on a small circle some distance from another point. circles


For example, lets say A = (39.73, -104.98) #Denver



B = (39.83, -106.06) #point approximately 50 nautical miles away from A along the great circle track to C


C = (40.75, -111.88) #Salt Lake City, approximately 323 nautical miles away from A


How do I get the set of orange coordinates when I don't know the angle but I do know the cross track distance?


For concreteness, lets say the outer points near B have a cross track distance of +- 30 nautical miles and the inner points have a cross track of +- 15nm.


And the points near C have a cross track of +- 25 nm.


Also I have an initial bearing from A to C of -1.34 #radians



Answer



If the unknown points are on the specified circles, and if you assume a simple circular reference surface (a sphere, not an ellipsoid), then it is easy to get the angle at A from AC to the unknown points:


angle = ctd / scd, where


ctd is your so-called cross-track distance (what I'd call arc distance), and



scd is your small-circle distance away from A


Then you "just" calculate the new coordinates of unknown point, U, say:


coordsU = traverse (coordsA, azimAC + angle, scd), where


coordsA are coordinates at A,


azimAC is azimuth AC, (conveniently, you already have this in radians),


angle and scd are as above,


traverse is the appropriate function for the so-called "direct" problem of spherical trigonometry -- for now, left as a separate exercise for you to research, or already available from some library.


arcgis 10.1 - Zonal statistics for millions of overlapping polygons in ArcPy?


I am running ArcGIS 10.1 with 64-bit GP on Windows 8.1 x64.


I want to calculate the zonal statistics of buffer polygons based on parcels and a kernel density raster. The 4.7 million polygons overlap and are stored in a single feature class, making the problem particularly difficult. I've tried to use the supplemental Spatial Analyst Tools developed by ESRI in various ways, but the script isn't fast enough. I've tried to compute the zonal statistics piece by piece, using small census areas (with several hundred buffer polygons each) and large ones (thousands of buffer polygons). None of these implementations perform adequately. I've included the following features in a Python script:



  • loading buffer polygons into the in_memory workspace before loading into the ZS tool

  • creating feature layers where needed

  • writing the results to a single text file



The code for performing ZS one polygon at a time is posted below, however the tool only completes about one parcel every minute. I'm wondering if anyone has any insight on how to improve zonal statistics performance or the script in general.


import arcpy, os, string, csv
from arcpy import env

# arcpy.ImportToolbox("D:\ClusterAnalysis_Deployment\SpatialAnalystSupplementalTools\Spatial Analyst Supplemental Tools.pyt")
arcpy.CheckOutExtension("Spatial")

# LOCATIONS #
# Local machine paths

gdb = r"D:\path\to\geodatabase.gdb"
results = r"\path\to\results"

# Virtual machine paths

# INPUT DATA #
b = gdb + "\\" + "ParcelBuffers"
kd = gdb + "\\" + "kdensity750"

env.workspace = gdb

env.overwriteOutput = True

# TEMPORARY DATA #
buffers = arcpy.MakeFeatureLayer_management(b, "bdata")
kdensity = arcpy.Raster(kd)

parcelList = []

cursor = arcpy.da.SearchCursor(buffers, ("PIN"))
for row in cursor:

parcelList.append(row[0])
del cursor

countDict = {}
countDict["Count:"] = 0

print "Script setup... done"

# GEOPROCESSING #
for PIN in parcelList:

parcel_ram = "in_memory" + "\\" + "parcel"
zs_table = results + "\\" + "zs_table_" + str(PIN) + ".dbf"
solution = results + "\\" + "ZS_solutions.txt"

arcpy.SelectLayerByAttribute_management(buffers, "NEW_SELECTION", "\"PIN\" = \'" + str(PIN) + "\'")
count = int(arcpy.GetCount_management(buffers).getOutput(0))
arcpy.CopyFeatures_management(buffers, parcel_ram)
arcpy.SelectLayerByAttribute_management(buffers, "CLEAR_SELECTION")
arcpy.gp.ZonalStatisticsAsTable_sa(parcel_ram, "PIN", kdensity, zs_table, "DATA", "ALL")


table = arcpy.gp.MakeTableView_management(zs_table, "zs_table_lyr")

fields = arcpy.ListFields(table)
field_names = map(lambda n:n.name, fields)
header = string.join(field_names, "\t")

with open(solution, 'w') as x:
x.write(header + "\n")
cursor = arcpy.SearchCursor(table)
for row in cursor:

values = string.join(map(lambda n: str(row.getValue(n)), field_names), "\t")
x.write(values + "\n")
del row, cursor

countDict["Count:"] = countDict["Count:"] + count
print "Zonal statistics computed and saved for parcel: " + str(countDict["Count:"])

arcpy.Delete_management(parcel_ram)
arcpy.Delete_management(zs_table)
arcpy.Delete_management(table)

del count

After some further testing of Felix's script I think the way my data sets are set up is wrong, leading to the spatial join feature query to pull up polygons that do overlap. Here are the data sets and their relevant attributes:



  • fishnet feature (ID'd by long integer "FnID")

  • minBound feature (ID'd by long integer "FnID")

  • parentLR (Join FID is FnID and Target FID is parID, where parID is long integer field created for buffer polygons)


The script:


##  Zonal statistics on a set of overlapping

import arcpy, traceback, os, sys, time
from arcpy import env

arcpy.CheckOutExtension("Spatial")

env.workspace = "C:\Geoprocessing\ClusterAnalysis"
env.overwriteOutput = True

# Field name of fishnet ID
parID="FnID"

parID2="FnID_1"

# Location of input map document containing the minBound (minimum bounding geometry)
# layer and parentLR (spatial join of fishnet to polygons) layer

mapDocument = env.workspace + "\\" + "ClusterAnalysis.mxd"
raster = env.workspace + "\\" + "ClusterAnalysis.gdb\intKD750"
rast = env.workspace + "\\" + "ClusterAnalysis.gdb\parcelRaster"
kdensity = arcpy.Raster(raster)


dbf="stat.dbf" # output zonal statistics table
joinLR="SD.shp" # intermediate data set for determining non-overlapping polygons
subset="subset"

##try:
def showPyMessage():
arcpy.AddMessage(str(time.ctime()) + " - " + message)
def Get_V(aKey):
try:
return smallDict[aKey]

except:
return (-1)
def pgonsCount(aLayer):
result=arcpy.GetCount_management(aLayer)
return int(result.getOutput(0))

# Define variables based on layers in ArcMap document
mxd = arcpy.mapping.MapDocument(mapDocument)
layers = arcpy.mapping.ListLayers(mxd)
minBound,parentLR=layers[0],layers[1]


# Query parentLR polygons based on fishnet grid of minBound polygon
with arcpy.da.SearchCursor(minBound, ("SHAPE@","FnID")) as clipper:
for rcrd in clipper:
feat=rcrd[0]
env.extent=feat.extent
fp='"FnID"='+str(rcrd[1])
parentLR.definitionQuery=fp
nF=pgonsCount(parentLR)
arcpy.AddMessage("Processing subset %i containing %i polygons" %(rcrd[1],nF))

arcpy.AddMessage("Defining neighbours...")
arcpy.SpatialJoin_analysis(parentLR, parentLR, joinLR, "JOIN_ONE_TO_MANY")
arcpy.AddMessage("Creating empty dictionary")

dictFeatures = {}
with arcpy.da.SearchCursor(parentLR, parID) as cursor:
for row in cursor:
dictFeatures[row[0]]=()
del row, cursor


arcpy.AddMessage("Assigning neighbours...")
nF=pgonsCount(joinLR)
arcpy.SetProgressor("step", "", 0, nF,1)
with arcpy.da.SearchCursor(joinLR, (parID,parID2)) as cursor:
for row in cursor:
aKey=row[0]
aList=dictFeatures[aKey]
aList+=(row[1],)
dictFeatures[aKey]=aList
arcpy.SetProgressorPosition()

del row, cursor
arcpy.AddMessage("Defining non-overlapping subsets...")
runNo=0
while (True):
parentLR.definitionQuery=fp
toShow,toHide=(),()
nF=len(dictFeatures)
arcpy.SetProgressor("step", "", 0, nF,1)
for item in dictFeatures:
if item not in toShow and item not in toHide:

toShow+=(item,)
toHide+=(dictFeatures[item])
arcpy.SetProgressorPosition()
m=len(toShow)
quer='"parID" IN '+str(toShow)+ " AND "+fp
if m==1:
quer='"parID" = '+str(toShow[0])+ " AND "+fp
parentLR.definitionQuery=quer
runNo+=1
arcpy.AddMessage("Run %i, %i polygon(s) found" % (runNo,m))

arcpy.AddMessage("Running Statistics...")
arcpy.gp.ZonalStatisticsAsTable_sa(parentLR, "PIN", kdensity, dbf,"DATA", "MEAN")
arcpy.AddMessage("Data transfer...")
smallDict={}
with arcpy.da.SearchCursor(dbf, (parID,"MEAN")) as cursor:
for row in cursor:
smallDict[row[0]]=row[1]
del row, cursor
with arcpy.da.UpdateCursor(parentLR, (parID,"MEAN")) as cursor:
for row in cursor:

aKey=row[0]
row[1]=Get_V(aKey)
cursor.updateRow(row)
del row, cursor
for item in toShow:
del dictFeatures[item]
m=len(dictFeatures)
if m==0:
break
parentLR.definitionQuery=fp

del smallDict, dictFeatures
parentLR.definitionQuery=''

Produces the zonal statistics error: http://resources.arcgis.com/en/help/main/10.1/index.html#//00vq0000000s010423




Wednesday 12 February 2020

arcgis desktop - Symbolizing layer with graduated colors by date/time field?


I've got a layer of polygons where each has a date associated with it. I want to color the polygons with graduated colors, but I want the colors to start with the newest polygons and fade to the oldest polygons. So far I haven't been able to find a way to do this in ArcMap (<= 9.3.1) other than to create a new field and calculate some integer or float value based on the date. Is it possible to use a date or time field directly?



  • The data is natively a date-time value

  • There are far too many records to use any unique-value renderer

  • Converting the date to a numerical field is an answer for a different question.




qgis processing - Using pyqgis logic to stop process if output has no geometry


I have a pyqgis script which customises the functionality of the grass algorithm r.lake.coords:


##Seed Points Lake=name

##DEM_input=raster
##Seed_vector_points=vector point
##Overspill_vector_boundary=vector polygon
##output=folder


# Import classes
from qgis.core import QgsRasterLayer
from qgis.utils import iface
from PyQt4.QtCore import QFileInfo

import processing

# Define layers
point_layer = processing.getObject(Seed_vector_points)
layer = processing.getObject(DEM_input)
boundary_layer = processing.getObject(Overspill_vector_boundary)
fileName = layer.source()
fileInfo = QFileInfo(fileName)
baseName = fileInfo.baseName()
rlayer = QgsRasterLayer(fileName, baseName)


# Get extent of the region
extent = iface.mapCanvas().extent()
xmin = extent.xMinimum()
xmax = extent.xMaximum()
ymin = extent.yMinimum()
ymax = extent.yMaximum()

for feat in point_layer.getFeatures():
geom = feat.geometry()

attr = feat.attribute("Z_ADD_1")
i = 0
while i < 5:
output_1 = processing.runalg('grass:r.lake.coords', layer, "%f"%(attr), "%f,%f"% (geom.asPoint().x(), geom.asPoint().y()) , False, "%f,%f,%f,%f"% (xmin, xmax, ymin, ymax), 0.00, None)
output_2 = processing.runalg('saga:rastercalculator', output_1['lake'], '','a>0', False, 3, None)
output_3 = processing.runalg('gdalogr:polygonize', output_2['RESULT'], "DN", None)
output_4 = processing.runalg('qgis:dissolve',output_3['OUTPUT'], True, "DN", None)
output_5 = processing.runalg("qgis:extractbylocation", output_4['OUTPUT'], boundary_layer, u'within', output + '/' + baseName + "_%.6f_%.6f" %(geom.asPoint().x(), geom.asPoint().y()) + "_%.2f"%(attr) + 'm.shp')
i += 1
attr += 0.5


My problem is that the final algorithm which the script calls upon produces additional blank shapefiles where there is no geometry, this is technically correct if the outputs from output_4 are not within the overspill boundary.


However, I would ideally like to have the script set up using some sort of python logic where it avoids producing final shapefiles with no geometry.



Answer



Great answer by @LaughU! If you want to prevent processed layers with no geometry to have no output, you could use the None parameter to save the output as a temporary file initially, check if this has any features and if so, write to a shapefile otherwise skip it:


...
output_5 = processing.runalg("qgis:extractbylocation", output_4['OUTPUT'], boundary_layer, u'within', 0.00, None)
layer = processing.getObject(output_5['OUTPUT'])
if not layer.featureCount() == 0:
QgsVectorFileWriter.writeAsVectorFormat(layer, output + '/' + baseName + "_%.6f_%.6f" %(geom.asPoint().x(), geom.asPoint().y()) + "_%.2f"%(attr) + 'm.shp', "utf-8", None, "ESRI Shapefile")

i += 1
attr += 0.5

arcgis 10.2 - Downloading ArcObjects SDKs for .NET?


Can anyone can suggest where I can freely download ArcObjects SDKs for ArcGIS?


I am using version 10.2.



Answer



While Esri does have some freely available tools, the ArcObjects SDK (.NET or Java) is not one of them. ArcObjects is a licensed product, available from the same distribution mechanism as the rest of the ArcGIS platform.


In general, the procedure for gaining access to Esri software is:



  1. Contact an Esri representative (via web site, US regional office or international distributor, as appropriate)


  2. Register a login on the Esri global account system

  3. License the software (there are a multitude of options here, including non-commercial and educational licenses) -- the licensed software will be associated with a Global Account

  4. Connect to the My Esri Customer Service portal with your Esri login

  5. Navigate to "My Organizations" -> "Downloads" and identify the product and version you want to download, then press the "View Downloads" button (ArcObjects is actually accessible through both the "ArcGIS for Server" and "ArcGIS for Desktop" products)

  6. Review the list of files, and press "Download" for the software package(s) you need (don't forget any ancillary products, like 64-bit Background Geoprocessing, or Additional Products like database client drivers)

  7. Execute the Windows extraction utility, which will create a folder with a Windows "Setup.exe" executable (Linux and Unix media are often GNUzipped tarfiles and generic resources may be distributed as ISO files).

  8. Execute the setup utility (assuming you have local access rights to do so)


If you do not have Customer Service access rights, you'll need to find someone in your organization that does, and have them either download the media kit(s) for you (then start at step 7) or ask them to have Esri add your global login as an approved download user for the organization (and start at step 4).


Using Python script to control GRASS GIS from outside?


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


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




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


import os
import sys

gisbase = os.environ['GISBASE'] = 'C:\OSGeo4W64\apps\grass\grass-6.4.3' #GISBASE needs to point the root of the GRASS installation directory
gisrc = 'C:\Users\Heinz\Documents\grassdata'
gisdbase = 'C:\Users\Heinz\Documents\grassdata'
location = 'newLocation'
mapset = 'TC'
LD_LIBRARY_PATH = 'C:\OSGeo4W64\apps\grass\grass-6.4.3\lib'

PATH = 'C:\OSGeo4W64\apps\grass\grass-6.4.3\etc';'C:\OSGeo4W64\apps\grass\grass-6.4.3\etc\python';'C:\OSGeo4W64\apps\grass\grass-6.4.3\lib';'C:\OSGeo4W64\apps\grass\grass-6.4.3\bin';'C:\Python27';'C:\OSGeo4W64\apps\Python27';'C:\OSGeo4W64\apps\msys'
PYTHONLIB = 'C:\Python27'
PYTHONPATH = 'C:\OSGeo4W64\apps\grass\grass-6.4.3\etc\python'
GRASS_SH = 'C:\OSGeo4W64\apps\msys\bin\sh.exe'


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

import grass.script as grass


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


As I run module, I got the error:


enter image description here


enter image description here


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




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


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

But still got the same error.



Did I misunderstand the instruction in the post?




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


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


GISBASE = C:\OSGeo4W64\apps\grass\grass-6.4.3
GISRC = C:\Users\Heinz\Documents\grassdata
Path =
C:\Python27
C:\OSGeo4W64\apps\grass\grass-6.4.3
C:\Users\Heinz\Documents\grassdata

C:\OSGeo4W64\apps\grass\grass-6.4.3\lib
C:\OSGeo4W64\apps\grass\grass-6.4.3\etc\python
C:\OSGeo4W64\apps\grass\grass-6.4.3\etc
C:\OSGeo4W64\apps\grass\grass-6.4.3\bin
C:\OSGeo4W64\apps\Python27
C:\OSGeo4W64\apps\msys
C:\OSGeo4W64\apps\grass\grass-6.4.3\scripts
C:\OSGeo4W64\bin
C:\OSGeo4W64\apps\Python27\Lib\site-packages


then I ran my python script, and got the same error: enter image description here


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




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


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

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


I also found this python-related file:



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

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


I have tried this in cmd: enter image description here


I will keep trying to make my script work!



Answer



After some hard work, I got my GRASS working with Python using a .pth file in the sitepackages folder.


To try this, do to the following steps:



  1. Go to the folder C:\OSGeo4W\apps\python27\lib\site-packages


  2. Create a file called grass.pth and open it with an editor


  3. Enter the following two lines (assuming your GRASS was installed with OSGeo4W; check the paths to be sure):


    C:\OSGeo4W\apps\grass\grass-6.4.3\etc\python C:\OSGeo4W\bin




  4. Save and close the file




  5. Try the following statement in your Python Script



    import grass.script as grass




P.S.: To make that work, you have to set your enviroment variables as described in posts above.


Another option might be to add the C:\OSGeo4W\bin path to your script with sys.path.append as you did with other paths


Once it works, check out the following post for getting to run the grass tools: GRASS Geoprocessing in Python Script


Cannot Install QGIS Mapserver On Ubuntu 12.04



I'm trying to install qgis-mapserver on my ubuntu 12.04. I just installed qgis desktop 2.6.0-brighton successfully. Here is my steps when installed the desktop :


sudo apt-get install python-software-properties 
sudo add-apt-repository ppa:ubuntugis/ubuntugis-unstable
sudo apt-get update
sudo apt-get install qgis python-qgis qgis-plugin-grass

But when I tried to install the mapserver :


sudo apt-get install qgis-mapserver libapache2-mod-fcgid

I got this error :



The following packages have unmet dependencies:
qgis-mapserver : Depends: qgis-providers (= 2.4.0-0precise2) but 2.6.0-0precise1 is to be installed

Answer



The package is called "qgis-server" (no "map" part anymore) now.


Using GDAL tool rasterize with QGIS so that points represent center of grid


I want to convert point data into a grid and used the gdal tool rasterize. The tool works but unfortunately the result looks like this.



enter image description here


The converted point does not represent the center of the grid but the topleft corner. I need it to be in the center. I couldn't find any settings in the tool nor Explanation in the GDAL-Documentation to accomplish that.


Does anyone have suggestions?


I rasterrized the points into raster and NOT raster to point.


If I upload the data directly as raster I get them correct but that doesn't always work. In the picture the underlying raster image is the correct one.


enter image description here



Answer



Just modify the extent to so that xmin = xmin - 1/2 pixel, xmax = xmax + 1/2 pixel, ymin = min - 1/2 pixel, ymax = ymax + 1/2 pixel.


i.e assuming your desired pixel size is 1m and your point layer extent is xmin, xmax, ymin, ymax = 32346264.0,32346299.0,5631990.0,5631999.0


Then set your extent to 32346263.5,32346299.5,5631989.5,5631999.5



enter image description here


enter image description here


Tuesday 11 February 2020

arcobjects - Can you programmatically change the button image for a ESRI.ArcGIS.Desktop.AddIns.Button?


I'm new to the AddIn customization model in ArcGIS 10 and I'm trying to figure out how to change a button image for a button. I'm looking for the equivalent to the BaseCommand UpdateBitmap method and haven't had any luck. Is this a limitation of the AddIn model? Or did I miss something obvious?



Answer



Ok, I've worked something out using ICommandItem.FaceID. I wasn't initially thinking of using ICommandItem, but it gets the job done.


ICommandItem commandItem = ArcMap.Application.Document.CommandBars.Find(...;
commandItem.FaceID = ESRI.ArcGIS.ADF.COMSupport.OLE.GetIPictureDispFromBitmap(bitmap);

qgis - Delineate polygonal catchment areas from a raster using open-source GIS


I'm 90% sure this basic question has been asked before but looking through the other posts I simply do not see a useful answer.


I need to calculate polygonal(vector) catchment areas from a raster based on minimum driving distance/cost from a set of points. I've already got the raster, and for that matter, the nodes(intersections with cost attribute) used to interpolate it. What I need now are polygons for each point describing the area nearer it than any other point. This is basically like a voronoi diagram but based on network distance rather than euclidean space.


How exactly can I do this using only open-source tools? I'm comfortable with the command line if I just know what I'm looking for. I feel like GDAL is going to be the answer, but if so, what exactly is the command I'm looking for?


Thanks!




openstreetmap - Orphan road segments in QGIS?


I am working with OSM road data for a developing nation. Included in the dataset are various small networks of roads that are orphans or islands i.e. not connected to the main road network. For example, four islands / orphans are evident in the section of roads below.


enter image description here



My aim here is to have all the roads connected that I can run a Steiner analysis without getting errors from nodes that do not connect.


How can I get QGIS to select these so that I can either delete them or connect them to other sections of road network?




python 2.7 - OSgeo4W64 Windows 7 & ImportError: No Module Named Site


I realize this question has been asked a lot, but I was not able to resolve my issue through the answers posted elsewhere... I am new to Python, gdal, OSgeo4W, etc. but have experience with ArcGIS.


I'm receiving the dreaded "ImportError: No module named site" error from Windows Cmd after running the following batch file by clicking on it directly from windows explorer (it is save as a .bat file).



ECHO Off
C:\OSGeo4W64\bin\gdalbuildvrt.exe depth_max.vrt *.tif
python C:\OSGeo4W64\bin\gdal2tiles.py -s EPSG:2278 -z 12 -w none "depth_max.vrt" depth_max
pause




Couple observations/comments:



  • The file runs fine when running it from the OSgeo4W window

  • gdalbuildvrt.exe executes fine in either scenario, the error is associated with the python command

  • Python 2.7 is installed at C:\Python27\ArcGIS10.2\ and C:\OSGeo4W64\apps\Python27

  • Python 3.6 is installed at C:\OSGeo4W64\apps\Python36

  • Running Windows 7 (Note: This may be run on Windows 10 as well)

  • gdalbuildvrt does not have a .py file; only .exe, which I found odd but I guess it is not a python script but a standalone program?

  • I tried including the SET commands for PYTHONHOME and PYTHONPATH in the batch file and it either did not work or broke python (it would then say DLLs were not found)



My eventual goal is to run this from a SHELL command via another program.




Knowing acceleration, gyro and compass data how to calculate north direction



I have a small video camera with an attached sensor. I can get from the sensor his acceleration, gyro and compass data. Knowing this data I need to rotate the camera so that it will look strictly north.



I know how to rotate the camera, but how to calculate the right angle?


How can I do it? Or where can I read about it? Or, if this is not a right place to ask such a question, what is the right place?


EDIT:


All coordinates are in an Earth-fixed inertial frame.The compass and accelerometer readings are unit vectors, and the gyroscope readings are in radians per second




GeoServer CSS Styling variable


May I set variable in CSS, like as




I have this:


* {
fill: [recode(strTrim(name),
'Akmola Region', '#6495ED',
'Aktobe Region', '#B0C4DE')];
fill-opacity:0.4;

stroke: lightgrey;
}

but I want set variables for Fill and Stroke-width parameters. Like this:



Akmola Region


name
Akmola Region






#
Akmola_fill_color
6495ED



0.4


#605f5f


Akmola_stroke_width
0.1








Aktobe Region


name
Aktobe Region






#
Aktobe_fill_color
B0C4DE



0.4


#605f5f


Aktobe_stroke_width
0.1







I want emulate hover event when mouse hover on map object. Maybe it is possible with css capabilities?




qgis - How to add Direction and Distance to attribute table?


anyone who can help. i just want to add direction (Bearing: i.e N 25 35 E) and Distance (Length: 125 meters) as my new field in polyline/line data. is there a plug-in to generate these fields? i tried to used "export/add geometry columns" in my line data, but only "Length" value was added.




gdal - Creating compressed GeoTIFF with QGIS Clipper tool?


I have a geo-tiff created by Pix4D which is about 375 Mb in file size.


I have clipped it to a smaller area using Qgis 2.8(using Raster ‣ Extraction ‣ Clipper), and the output file is about 900 MB in size.


Why did the clipping process generate such a large file, when the source file was smaller?


How do I run the clipper tool so that it produces a smaller size?



Answer



The Clipper tool makes an uncompressed image by default. Read the GDAL manual of your format and add manually the compression options into the gdal_translate command that is shown in the lowest pane.


For example for GeoTIFF read http://gdal.org/frmt_gtiff.html and use for example-co COMPRESS=DEFLATE -co PREDICTOR=2 which gives a well compressed, lossless output for topographic maps.


enter image description here


If you think that it would be trivial to add compression options to the user interface of the Clipper tool it is not because more than 30 formats are supported for the output and all of them have different settings for compression if compression is supported at all.



Monday 10 February 2020

web mapping - Seeking tool(set) for creating community map?




When you move to a new place, it's very useful to have maps tailored to your interests available at the click of a mouse or tap of a smart phone. The following projects are great examples of more personal local maps:



These resources are fantastic, adding independence and providing a counter to the homogeneity of massive worldwide maps such as Google and Bing.



How can I best implement this on a small scale for a local organisation?


For example maps are needed to locate fruit trees for an urban agriculture project such as abundance.


The criteria are as follows:



  • Low price or free

  • User friendly, allowing more than 1 person to contribute

  • Scalable

  • (preferable but not essential) levels of security so sensitive information is only viewable by some people


The options I have considered for such a tasks, along with +s and -s are:




  • QGIS Cloud: + works on very accessible GIS, free. - seemed glitchy, not necessarily user friendly or easy to add new features for new users

  • Create geojson files and serve them from Github, using OpenLayers on the client side to process the files, allowing different layers to be added/removed (I have actually experimented with this on, adding a couple of shop polygons not yet on OSM or Google Maps) + simple, lightweight and free, - not user friendly - must edit geojson files externally

  • Geoserver: I've used this before, and enjoyed it until AWS, my remote server, started charging me + scalable, user interface, security levels - potential costs, - overkill for a small project.

  • MapQuest Open + uses pre-existing OSM data, - seemed limited in terms of functionality


Which of these would you recommend?


Have I missed any better options from the list?


And how should I organize different tools to work harmoniously together in the final solution?


Will try to keep people updated with progress on this as I think there is a real need for development of interactive local and community maps.





arcgis desktop - Exporting elevations from AutoCAD in format that ArcMap can read


I have elevations in Autocad, I can shift between a triangles, contours etc if it makes any difference.


How can I export these into something I can use in ArcMap?


I have tried the Make shape function with both contours and triangles, which work fine. I can display the contours or polygons. The only hickup is that the Z values are not transferred with it.


I do not know any CAD myself, and those in my organization that do haven't got any clue about ArcGIS. So I'm kind of stuck between not knowing what I can get myself and the ones who do doesn't know what I need.


Google found me some supposed workflows that I fail to complete.


How do I export from AutoCAD (2010, 2011, 2013) or Civil3D to shape/TIN/ESRI Grid?



Answer




If you have Civil 3D and your elevation data is in a TIN surface, you should be able to run SURFACEEXPORTTODEM, or, from the contextual menu that pops up when you have the surface selected:enter image description here


You can export to a DEM (geotiff), which is easily readable by other GIS software.


Releasing PyQGIS file locks?


I was wondering what triggers the release of file locks in pyQGIS?


I am trying to delete a few data sources (used temporarily) by calling QgsVectorFileWriter.deleteShapeFile, but I have to quit QGIS before I can do that. I have loaded the sources into QgsVectorLayer objects. Must all of these objects and references to them be garbage collected before I can delete the source? Is there a way to force this?




I've managed to create a minimal code sample that fails. Make sure temp dir is empty before running.


from qgis.core import *

import processing, os, gc

project_temp_dir = "C:/Path/To/My/Dir/"
layer1_path = project_temp_dir + "layer1.shp"
layer2_path = project_temp_dir + "layer2.shp"
input_layer = QgsMapLayerRegistry.instance().mapLayersByName('in_layer')[0]
if not input_layer.isValid(): raise Exception("Failed to grab input layer")

# Create layer 1
err = QgsVectorFileWriter.writeAsVectorFormat(input_layer, layer1_path, "utf-8", input_layer.crs())

if err != QgsVectorFileWriter.NoError: raise Exception("Failed to write layer 1")

# Load layer 1
layer1 = QgsVectorLayer(layer1_path, "lyr1", "ogr")
if not layer1.isValid(): raise Exception("Failed to load layer 1")

# Use layer 1 to create layer 2, read-only makes no difference
# if not layer1.setReadOnly(): raise Exception("Could not set layer 1 to read-only")
processing.runalg("qgis:reprojectlayer", layer1, "EPSG:54030", layer2_path)


# Load layer 2
layer2 = QgsVectorLayer(layer2_path, "lyr2", "ogr")
if not layer2.isValid(): raise Exception("Failed to load layer 2")

del layer1
del layer2
del input_layer
gc.collect()
print "Garbage: " + str(gc.garbage) # Empty


# Remove data sources for layers - FAILS!!
for f in os.listdir(project_temp_dir):
if f.endswith(".shp") and not os.path.isdir(project_temp_dir + f):
if not QgsVectorFileWriter.deleteShapeFile(project_temp_dir + f):
# F*%&ing locks.
print "Failed to clear project temp directory."



I found that it works if I use QgsVectorFileWriter to create layer2, instead of the processing algorithm. I get the same error if try the qgis:clip algorithm. So is this a bug in processing? Am I using it wrong?



Answer




Sorry to keep answering my own questions, but I think I found a solution.


As it turns out, it works well if you add the layer to the map registry, and then remove it again. The map registry takes ownership of the layer, so when it is deleted from the registry, the locks are freed. Note that you have to add the layer to the legend (.addMapLayer(layer, addToLegend = False) won't work).


Still not sure whether to call this a solution or a workaround, but it does the job.


# ...

# Replace the following code (note: should do error checking on map registry functions):

# Load layer 1
layer1 = QgsVectorLayer(layer1_path, "lyr1", "ogr")
if not layer1.isValid(): raise Exception("Failed to load layer 1")

QgsMapLayerRegistry.instance().addMapLayer(layer1) #!!!!

# Use layer 1 to create layer 2
processing.runalg("qgis:reprojectlayer", layer1, "EPSG:54030", layer2_path)

# Load layer 2
layer2 = QgsVectorLayer(layer2_path, "lyr2", "ogr")
if not layer2.isValid(): raise Exception("Failed to load layer 2")
QgsMapLayerRegistry.instance().addMapLayer(layer2) #!!!!


# Remove layer references
QgsMapLayerRegistry.instance().removeMapLayer(layer1.id()) #!!!!
QgsMapLayerRegistry.instance().removeMapLayer(layer2.id()) #!!!!

# Remove data sources for layers
for f in os.listdir(project_temp_dir):
if f.endswith(".shp") and not os.path.isdir(project_temp_dir + f):
# ...

If anyone has more info, I'd be happy to learn more about this.



Creating surface from points in ArcGIS for Desktop?


how to create surface from points in ArcGis? I have a points (centroids of the neighbourhood and number of people living in every neighbourhood) and i would like to create surface that would present density of population.



Answer



There are a wide variety of interpolation tools available for this type of analysis. Calculating Kernel Density will likely yield good results for you. You have the option to incorporate population values for each feature. If you specify population values each feature, these count or quantity values will be spread across the landscape to produce a continuous surface. Here are a few tips:


Run Kernel Density (Spatial Analyst) and choose a search radius and population field. By default, ArcGIS takes the shortest width or height of the output extent and divides by 30. This is supposedly a good rule rule-of-thumb that is robust to effects from outliers. However, you can change this value to make the results more meaningful for your analysis. In this example, I used a search radius of 50 which is greater than the default 34m search radius (Figure 1).


ArcGIS automatically categorizes the output into classes. You can display the default classes, change the classes to, for example, jenks, or use a continuous scheme. Finally, for cartographic purposes, I would recommend using bilinear resampling during display (Figure 2).


Figure 1


enter image description here


Figure 2



enter image description here


Lost access to Desktop Help of ArcMap 10.5?


Today I found that I can't aceess help on system tools through "Tool Help" button. It opens help window which says:



Can’t reach this page



  • Make sure the web address //ieframe.dll/dnserrordiagoff.htm# is correct

  • Search for this site on Bing

  • Refresh the page




enter image description here


I am afraid it is Windows October update in action. Does anybody know how to fix this?


I tried profile reset and repair installation, no luck.


It is Windows 10 I worry about. The October update was delayed, because it was deleting files on PCs.



Answer



Solution posted on GeoNet did the trick. Answer by cringo at the bottom of that post is the one.


Sunday 9 February 2020

javascript - Openlayers 3 equivalent of render intents of OpenLayers.Stylemap?


I want to convert my openlayers 2 features to openlayers 3. I have following part in openlayers 2. I need an equivalent of StyleMap in openlayers 3 so that I can distinguish the default action and select action.


var subMap = new OpenLayers.StyleMap({
"default": new OpenLayers.Style({

fillColor:"orange",
strokeColor:"orange"
}),
"select": new OpenLayers.Style({
fillColor:"#00FFFF",
strokeColor:"#00FFFF"
})
});

I can write the following for "default" case. But where to mention for the "select" case?



var subMap= new ol.style.Style({
fill: new ol.style.Fill({
color: "orange"
}),
stroke: new ol.style.Stroke({
color:"orange"
})
});

Can anyone please help in this regard?




Answer



In OpenLayers 3 (or 4), features selection is handled by an interaction. The relevant doc about selection is here: http://openlayers.org/en/latest/apidoc/ol.interaction.Select.html


Here's a piece of code that allows to select any features in your map, and to apply a special style to it (styleSelect):


  // render selected features with a red stroke
var styleSelect = new ol.style.Style({
stroke: new ol.style.Stroke({
color: 'red',
width: 2
})
});


// a normal select interaction to handle click
var select = new ol.interaction.Select({style: styleSelect});
map.addInteraction(select);

The default (not selected) style can be defined elsewhere, using another ol.style.Style object.


javascript - allow client to upload Layers to Openlayers 3


I want to allow the clients to add vector layers in form of a geojson to openlayers 3. It should look like in this fiddle . I know how to add multiple layers but I want the client to be able to add a Layer and see his added layer. The layer doesn't have to be saved for his next visit at the page.



Answer



there are two different way you can try,
1) you need to host your web server to receive the file that client uploaded. The file will be processed in your server and return the geojson file to the front-end.

here is an example (with php)



then you will get a geojson file to show in ol3.

2) the other one is using File API to process the file without write the server side code, and read the file content to show in ol3.
dev.w3.org/2006/webapi/FileAPI/


i also make a pure javascript lib to convert zip file with shp, dbf to geojson, just like user upload their file and show in ol3.




(this demo support .prj so no need to set EPSG and assign proper Encoding to show correct attribure in ol3)


arcpy - Python conditional on current and previous row?


I need to calculate a new field with a python conditional where some values only occur if a sequence in a separate attribute field (X) occurs.

For instance my data set looks like this. Y is the attribute I want to fill out:


ID   X  Y
1 1 -
2 1 -
3 0 -
4 0 -
5 0 -
6 1 -
7 1 -
8 1 -

9 0 -
10 0 -
11 1 -

Values for Y:
If X = 1, Y = 1
If X = 0 and previous X = 1, Y = 2
If X = 0 and previous x = 0, Y = 0


So my data would look like:


ID   X  Y

1 1 1
2 1 1
3 0 2
4 0 0
5 0 0
6 1 1
7 1 1
8 1 1
9 0 2
10 0 0

11 1 1

I have seen previous lines stored using:


prev_X = ""

But don't know how to incorporate it into my script which looks like:


fc = r'.....'

rows = arcpy.UpdateCursor(fc)


Y = 1
prev_X = "" #<=========== DOESN'T WORK
for row in rows:
if row.X == 1:
Y = 1
row.Y = Y
rows.updateRow(row)

while row.X == 0 and prev_X == 1:
Y = 2

row.Y = Y
rows.updateRow(row)

while row.X == 0 ande prev_X == 0:
Y = 0
row.Y = Y
rows.updateRow(row)
rows.next()

row.X = X

rows.updateRow(row)

del row, rows

Answer



Instead of using while statements, I'd go with your idea of using a prevX var. The key is that you declare the variable after you know the value of X is for each turn through the loop.


By declaring prevX in each of your control structures (if or elif statement), you know what the value of X was in that loop.


This is just a nested conditional (you could accomplish this by using and between conditionals too).


fc = r'.....'
rows = arcpy.UpdateCursor(fc)


for row in rows:

if row.X == 0:
row.Y = 1
rows.updateRow(row)
prevX = 0 #Now the script knows that the last value of X was 0

elif row.X == 1:
if prevX == 1:
row.Y = 2

rows.updateRow(row)
prevX = 1 #Now the script knows that the last value of X was 1

elif prevX == 0:
row.Y = 0
rows.updateRow(row)
prevX = 0 #Now the script knows that the last value of X was 0

In your question, you didn't state what the value of Y would be if X = 0 and there was no previous value of X (i.e. if this was the first row).


qgis - How to convert a shapefile into OSM .PBF format


I have a large .gdbtable file (1.3 GB) with a road network that I would like to convert to .PBF format keeping the topology information.


I can load the .gdbtable file in QGIS so I thought it should be straightforward to write/export the file in .PBF but I cannot find this option. I would rather use QGIS to do this but an ArcGIS solution would be welcomed as well.


Any ideas on how to do this?




Answer



QGIS uses the GDAL OSM driver, which is read-only. So no chance to write .pbf files directly.


You might look into http://wiki.openstreetmap.org/wiki/Ogr2osm to get .osm files out of OGR supported formats.


To get a .pbf file, you can use OSMconvert in a second step.


Saturday 8 February 2020

Using exitQgis() in PyQGIS?



After infinitely many attempts to fix this, I can't prevent python to crash (segmentation fault) when scripting in QGIS. I am using the qgis2 build from osgeo4mac, installed with homebrew on El Capitan. I am interested in running something like:


from qgis.core import *


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

# DO STUFF HERE

app.exitQgis()

The above script will work in the sense that all the STUFF I want it to do will be properly executed, and I believe that app.exitQgis() will also remove the data provider and layer registries from memory. However, python will crash unexpectedly on quit because of a SIGSEGV error.


Upon reading this answer from spacedman, I managed to avoid the crash with the following structure:



from qgis.core import *

def main(app):
# DO STUFF HERE
return

app = QgsApplication([], False)
app.initQgis()
main(app)
app.exitQgis()


Although I have no idea why, the above runs smoothly. I however encounter the SIGSEGV error again when I try to run something more involved, using the Processing toolbox:


from qgis.core import *
from processing.core.Processing import Processing

def main(app):
# DO STUFF HERE
Processing.runAlgorithm(...) # DO MORE INVOLVED STUFF
return


app = QgsApplication([], False)
app.initQgis()
Processing.initialize()
main(app)
app.exitQgis()

Again, the above script will properly execute everything, but python will crash on quit.


I am out of solutions and I was wondering what would happen if I simply omit the app.exitQgis() line from the above script. This would do all what needs to be done, avoid the crash, but would not "clean up".


Is this a terrible idea?


Will this use up all my RAM if I run it multiple times?



Update:


The segfault seems to be triggered by objects loaded in memory. For example this will crash:


from qgis.core import *
app = QgsApplication([], False)
app.initQgis()
layer = QgsVectorLayer('shapefile.shp', 'all', 'ogr')
app.exitQgis()

But this will not:


from qgis.core import *

app = QgsApplication([], False)
app.initQgis()
layer = QgsVectorLayer('shapefile.shp', 'all', 'ogr')
del layer # CLEANING UP HERE
app.exitQgis()

which is why a function structure also prevents the crash: loaded layers are local variables within the function environment. Regarding the use of the Processing toolbox, the number of objects in memory does indeed increase after the use of Processing.runAlgorithm(...). Using the garbage collector module, import gc, I find that gc.get_count() increases right after an algorithm from Processing is used.




some details on the crash. the python ouput produced by the app.exitQgis() line:


src/core/auth/qgsauthmethodregistry.cpp: 154: (~QgsAuthMethodRegistry)[936ms] cleanup: Basic

src/core/auth/qgsauthmethodregistry.cpp: 154: (~QgsAuthMethodRegistry) [0ms] cleanup: Identity-Cert
src/core/auth/qgsauthmethodregistry.cpp: 154: (~QgsAuthMethodRegistry) [0ms] cleanup: PKI-PKCS#12
src/core/auth/qgsauthmethodregistry.cpp: 154: (~QgsAuthMethodRegistry) [0ms] cleanup: PKI-Paths
src/core/qgsproviderregistry.cpp: 234: (clean) [3ms] cleanup:DB2
src/core/qgsproviderregistry.cpp: 234: (clean) [0ms] cleanup:WFS
src/core/qgsproviderregistry.cpp: 234: (clean) [1ms] cleanup:arcgisfeatureserver
src/core/qgsproviderregistry.cpp: 234: (clean) [0ms] cleanup:arcgismapserver
src/core/qgsproviderregistry.cpp: 234: (clean) [1ms] cleanup:delimitedtext
src/core/qgsproviderregistry.cpp: 234: (clean) [0ms] cleanup:gdal
src/core/qgsproviderregistry.cpp: 234: (clean) [0ms] cleanup:gpx

src/core/qgsproviderregistry.cpp: 234: (clean) [1ms] cleanup:memory
src/core/qgsproviderregistry.cpp: 234: (clean) [0ms] cleanup:mssql
src/core/qgsproviderregistry.cpp: 234: (clean) [1ms] cleanup:ogr
src/providers/ogr/qgsogrconnpool.cpp: 41: (~QgsOgrConnPool) [0ms] Entering.
src/providers/ogr/qgsogrconnpool.cpp: 41: (~QgsOgrConnPool) [0ms] Leaving.
src/core/qgsproviderregistry.cpp: 234: (clean) [0ms] cleanup:ows
src/core/qgsproviderregistry.cpp: 234: (clean) [0ms] cleanup:postgres
src/core/qgsproviderregistry.cpp: 234: (clean) [0ms] cleanup:spatialite
src/core/qgsproviderregistry.cpp: 234: (clean) [0ms] cleanup:virtual
src/core/qgsproviderregistry.cpp: 234: (clean) [0ms] cleanup:wcs

src/core/qgsproviderregistry.cpp: 234: (clean) [1ms] cleanup:wms
src/providers/ogr/qgsogrconnpool.cpp: 36: (QgsOgrConnPool) [78ms] Entering.
src/providers/ogr/qgsogrconnpool.cpp: 36: (QgsOgrConnPool) [0ms] Leaving.
[Finished in 7.8s with exit code -11]

the OSX crash report:


Date/Time:             2017-08-03 13:06:51.003 -0400
OS Version: Mac OS X 10.11.3 (15D21)
Report Version: 11
Anonymous UUID: 08F4C3E0-65A6-5136-C3BB-0D146F412EAC



Time Awake Since Boot: 47000 seconds

System Integrity Protection: enabled

Crashed Thread: 0 Dispatch queue: com.apple.main-thread

Exception Type: EXC_BAD_ACCESS (SIGSEGV)
Exception Codes: EXC_I386_GPFLT

Exception Note: EXC_CORPSE_NOTIFY

Thread 0 Crashed:: Dispatch queue: com.apple.main-thread
0 libgdal.20.dylib 0x00000001178771f4 OGR_DS_GetDriver + 14
1 libogrprovider.so 0x0000000126408000 QgsOgrProviderUtils::OGRDestroyWrapper(void*) + 40
2 libogrprovider.so 0x00000001263ee5e4 QgsOgrProvider::close() + 66
3 libogrprovider.so 0x00000001263ee348 QgsOgrProvider::~QgsOgrProvider() + 164
4 libogrprovider.so 0x00000001263ee694 QgsOgrProvider::~QgsOgrProvider() + 14
5 org.qgis.qgis2_core 0x0000000114ebdfc6 QgsVectorLayer::~QgsVectorLayer() + 58
6 _core.so 0x0000000114077a30 sipQgsVectorLayer::~sipQgsVectorLayer() + 14

7 _core.so 0x0000000114079c41 release_QgsVectorLayer(void*, int) + 66
8 sip.so 0x0000000108ac48c1 forgetObject + 114
9 sip.so 0x0000000108ac46b1 sipWrapper_dealloc + 14
10 org.python.python 0x00000001070dc257 subtype_dealloc + 667
11 org.python.python 0x00000001070bbee2 dict_dealloc + 129
12 org.python.python 0x00000001070bd13d insertdict_by_entry + 249
13 org.python.python 0x00000001070bba5a insertdict + 51
14 org.python.python 0x00000001070bb221 dict_set_item_by_hash_or_entry + 40
15 org.python.python 0x00000001070bfd78 _PyModule_Clear + 177
16 org.python.python 0x00000001071218ce PyImport_Cleanup + 958

17 org.python.python 0x000000010712c2ec Py_Finalize + 297
18 org.python.python 0x000000010713e153 Py_Main + 2455
19 libdyld.dylib 0x0000000107a595ad start + 1


qgis - Seeking detailed description of qgs file structure?


Is there a detailed description on the .qgs file structure somewhere?


I checked at http://mrcc.com/qgis.dtd and got




The requested URL /qgis.dtd was not found on this server.



and at https://svn.osgeo.org/qgis/trunk/qgis/qgis.dtd which wasn't very detailed.




How to install stable GDAL on Ubuntu 11.10?


I am trying to install the stable ubuntu package from here https://launchpad.net/~ubuntugis/+archive/ppa/+packages on my ubuntu 11.10 server



I have added the ubuntugis stable repository to server's list of software sources with these commands


sudo apt-get install python-software-properties  
sudo add-apt-repository ppa:ubuntugis/ppa

But, apt-get install libgdal1-1.8.0 command fails to find the package. Is the ubuntu 11.10 server not supported? I am newbie with this, so i am not sure if the missing oneric (ubuntu 11.10) entries is a sign that it is supported.



Answer



GDAL 1.9 is available for Oneiric 11.10 in ubuntugis unstable. https://launchpad.net/~ubuntugis/+archive/ubuntugis-unstable?field.series_filter=oneiric


According to this email exchange from 2 Feb 2012:



stable PPA is only upgraded twice a year, or for security releases. It does not contain any oneiric packages because oneiric is new. and we do not push the packages since we don't know if they are stable since no user tested them at the moment of the upgrade. In april, the unstable-oneiric packages will be moved to stable.




cartography - Good websites for map-lovers?





I often visit strange-maps, cool-maps, makingmaps and hipkiss.org.


Do you know of other good websites for map-lovers?



Answer






Websites:



David Rumsey Map Collection http://www.davidrumsey.com/



Map Gallery http://www.cartotalk.com/index.php?showforum=14


Ancient World Mapping Center http://www.unc.edu/awmc/mapsforstudents.html


Huge Resource - Berkeley Map Collection
http://www.lib.berkeley.edu/EART/MapCollections.html


These are good for OpenStreetMap-lovers:



The International Cartographic Association website has a section with some beautiful maps selected every month.


How to set the URL from an icon in OpenLayers 3?


I've tried to set the URL from an icon but without success. In fact, I've been searching in OpenLayers API but I haven't seen any method to set it. In OpenLayers 2 I've seen that It's possible.


I want to be able to change the URL at any time.



Do you have any idea?



Answer



To answer your question, now that I understand it, there is no way to set the src of an icon after instantiation. My guess would be because the creation of a style.Icon object is tied in to the image cache somehow. The developers obviously thought not to duplicate that logic: https://github.com/openlayers/ol3/blob/v3.7.0/src/ol/style/iconstyle.js#L589


This example here shows you how to create an icon by referencing an image file. Note the data/icon.png specified in the icon configuration is relative and resolves to http://openlayers.org/en/v3.7.0/examples/data/icon.png. You could either use an absolute or relative url depending on what you need.


var iconStyle = new ol.style.Style({
image: new ol.style.Icon(/** @type {olx.style.IconOptions} */ ({
anchor: [0.5, 46],
anchorXUnits: 'fraction',
anchorYUnits: 'pixels',
opacity: 0.75,

src: 'data/icon.png'
}))
});

As with all other style types in OpenLayers, it is better to create a cache of styles instead of one style per feature. An example on how to do this would be roughly as follows:


var styleCache = [];

for (var i = 0; i < dataSet.length; i++) {
var data = dataSet[i];
var style;

if (styleCache[data.styleAttribute]) {
styleCache[data.styleAttribute] = new ol.style.Style({
image: new ol.style.Icon({ etc.. })
});
}
style = styleCache[data.styleAttribute];
var feature = new ol.Feature({ etc.. });
feature.setStyle(feature);
}


If you need to change the style on the fly, you would just reference one of the styles in the style cache, or add a new style to the cache, then reference it. On an existing feature you can change the style using feature.setStyle(iconStyle).


I hope this helps.


QGIS opening file - very slow on PostGIS database with many layers


When opening a QGIS file with many layers (20 or so) it takes +5 minutes to load. All the layers come from a single PostGIS database.


When looking into the server logs, I see that there are many requests like:


LOG:  SELECT typname,typtype,typelem,typlen FROM pg_type WHERE oid=1700

LOG: SELECT pg_catalog.format_type(atttypid,atttypmod) FROM pg_attribute WHERE attrelid=34757 AND attnum=3
LOG: SELECT description FROM pg_description WHERE objoid=34757 AND objsubid=3

Since there are thousands (+3000) of these queries, I suppose the problem lies there...


Any idea?




color - QGIS coloring NULL values


I am using QGIS Version 2.14.3. I have a map with polygons and a lot of different attributes. Some polygons have NULL values in some columns. I've tried to color the polygons with properties->style->graduated. But in some columns half of the values are NULL and now they don't have any colour.


I know that i could color the NULL values when I chose "categorized" coloring, but I have way to many different entries to do that, because QGIS chooses a colour for every different entry.


I need to have intervals for coloring. Is there an easy way to do that?




Friday 7 February 2020

distance - Is there an R equivalent of NEAR in ESRI ArcGIS?


I have 54 points in one data set A I have 300+ data points in a dataset B.


I wish to find the distance for each point in dataset-set B, to points in dataset A.


Nearest Neighbor packages appear to only find neighbors within a given data-set, not a second data set.


rdist was close, but I didn't want a matrix, just nearest.--"Given two sets of locations computes the full Euclidean distance matrix among all pairings or a sparse version for points within a fixed threshhold distance."


This looked good, saying "Read two geocoded point sets from Comma Separated Value (CSV) files into R data objects. Assign to each member of first point set the geographically closest point from second set" but is so dated (2010), so something more recent must exist...


Same issue with this



Other packages find distances between points and lines, points and polygons, etc. viz: http://cran.r-project.org/web/packages/geosphere/geosphere.pdf


The omission of P2P suggests I am missing something simple...




Google Earth Engine - filter landsat by cloud cover over multiple polygons


I want to make a Landsat image collection that is totally cloud free over a SET of polygons. I have about 300 polygons (all within a single Landsat tile). I found this question, which is doing exactly what I want over a single polygon, but I can't figure out how to get it to work over a set of polygons (which are loaded as a feature class).


Here is what I tried to do to modify the solution from the question posted above, where imgcollection is my Landsat TOA Image Collection and fc is my feature class of polygons. I was trying to use map and clip together to select all the polygons in the feature class for each image (based on the answers in this related question).



  var combine = imgcollection.map(function(imgcollection) { return 
imgcollection.clip(fc); });

After this line, I tried to continue with the solution from the first question, but it doesn't work. The error I get is as follows, which confuses me. I'm trying to use the variable "combine" as the geometry - am I misunderstanding this and actually using something else?


  ImageCollection (Error)
Image.reduceRegion, argument 'geometry': Invalid type. Expected:
Geometry. Actual: ImageCollection.

EDIT: full code


  var fc = ee.FeatureCollection('ft:omit');


var imgcollection = ee.ImageCollection("LANDSAT/LC08/C01/T1_RT_TOA");

var combine = imgcollection.map(function(img) { return img.clip(fc); });

var withCloudiness = combine.map(function(image) {
var cloud = ee.Algorithms.Landsat.simpleCloudScore(image).select('cloud');
var cloudiness = cloud.reduceRegion({
reducer: 'mean',
geometry: combine,

scale: 30,
});
return image.set(cloudiness);
});

var filteredCollection = withCloudiness.filter(ee.Filter.lt('cloud', 10));
print(filteredCollection);

Answer



The error is quite clear. Parameter geometry of function reduceRegion must be a Geometry. In the documentation says:




geometry (Geometry, default: null): The region over which to reduce data. Defaults to the footprint of the image's first band.



So, if you want to reduce over fc:


var fc = ee.FeatureCollection('ft:omit');

var imgcollection = ee.ImageCollection("LANDSAT/LC08/C01/T1_RT_TOA");

var combine = imgcollection.map(function(img) { return img.clip(fc); });

var withCloudiness = combine.map(function(image) {

var cloud = ee.Algorithms.Landsat.simpleCloudScore(image).select('cloud');
var cloudiness = cloud.reduceRegion({
reducer: 'mean',
geometry: fc.geometry(),
scale: 30,
maxPixels: 1e13,
});
return image.set(cloudiness);
});


var filteredCollection = withCloudiness.filter(ee.Filter.lt('cloud', 10));
print(filteredCollection);

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...