Thursday, 27 July 2017

python - Alternatives to using Arcpy


I seem to use ESRI's Arcpy site package for virtually all of my python geoprocessing. To ESRI's credit, these are an incredible suite of tools that can help accomplish a great deal. However, I would also like to create geoprocessing scripts outside of the ESRI Arcpy domain. For example, if I want to clip a raster to a polygon, I would start with the following script from ESRI:


# Import system modules
import arcpy
from arcpy import env

from arcpy.sa import *

# Set environment settings
env.workspace = "C:/sapyexamples/data"

# Set local variables
inRaster = "elevation"
inMaskData = "mask.shp"

# Check out the ArcGIS Spatial Analyst extension license

arcpy.CheckOutExtension("Spatial")

# Execute ExtractByMask
outExtractByMask = ExtractByMask(inRaster, inMaskData)

# Save the output
outExtractByMask.save("C:/sapyexamples/output/extractmask")

I'm not sure how I would accomplish the same task programmatically without Arcpy. My questions for the serious programmers out there: What collection of python tools do you use to accomplish tasks that ESRI users would accomplish with the Arcpy site package? Where do I begin?



Answer




GDAL is the tool to use. In fact that entire call is one line for gdal_rasterize:


gdal_rasterize -l mask -i -burn -9999 mask.shp elevation.tif

if you knew the no data value of the dem


For some python control:


lyr = 'mask'
shp = 'mask.shp'
dem = 'elevation.tif'
ndv = -9999
p = os.Popen('gdal_rasterize -l %s -i -burn %d %s %s' % (lyr,ndv,shp,dem)


where your variables could be set in python


For full python:


from osgeo import gdal, ogr
from osgeo.gdalconst import *
shp = ogr.Open('mask.shp')
lyr = shp.GetLayer('mask')
dem = gdal.Open('elevation.tif', GA_Update)
ndv = dem.GetRasterBand(1).GetNoDataValue()
gdal.RasterizeLayer(dem, 1, lyr, None, ndv) # other options, such as transformer func, creation options...

dem = None

I just took a quick peek at the syntax for the C API, so my syntax for python is probably off a little. See gdal_alg.h: http://gdal.org/gdal__alg_8h.html


No comments:

Post a Comment

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...