Alternatives to using Arcpy

Question

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 1

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

Answer 2

You will find a number of other similar questions on this site that ask the same basic question and have very good references. The most similar (and detailed) is:

• What are the Python tools/modules/add-ins crucial in GIS?

Others include:

• Pure Python Library for Geometry Operations
• What tools in Python are available for doing great circle distance + line creation?
• Python module to delete SHP features (without Desktop GIS installed)

Answer 3

A good starting point would be the Geospatial Data Abstraction Library. It is actually made up oftwo libraries -- GDAL for manipulating geospatial raster data and OGR for manipulating geospatial vector data but people usually just call it GDAL.

There's a geoprocessing with Python using open source GIS course at the Utah State University. You might want to check it out, too.

Answer 4

In a lot of my academic research I work with LiDAR data doing surface analysis for geomorphology. I quickly found that performing a lot of operations using arcpy was very slow, especially on large datasets. As a result I began using:

• pyshp to manipulate shapefiles and update attribute tables
• numpy to manage ASCII rasters and perform kernel-based analysis such as curvature calculations
• scipy to perform statistical analysis on results and perform curve fitting for surfaces
• matplotlib to plot graphs and other graphical results, such as basic maps for quick visualizations

I would also recommend the book, Quantitative Modeling of Earth Surface Processes to anyone who wants to learn more about analyzing raster surfaces. The book comes with great code samples in C++, which are much more efficient than the ArcGIS tools. These algorithms can can also be ported to Python without needing anything more complex than numpy, although they run much faster in C++.

Answer 5

For people using ESRI I think GRASS would be a very similar environment with a GUI python environment and organized in separate 'toolkits' for different tasks (raster, vector, solar toolkits etc.). The scripting has other options besides Python but that is how I use it.

Definitely check out this great link which is up-to-date (I believe): http://grass.osgeo.org/wiki/GRASS_and_Python

EDIT: another link for those with background in ESRI: http://grass.osgeo.org/wiki/GRASS_migration_hints

I also second the motion of GDAL. It is invaluable and I would be lost without it.

Answer 6

I think the answers given so far cover basically all package out there worth mentioning (espically GDAL, OGR, pyshp, NumPy)

But there is also the GIS and Python Software Laboratory, that hosts a couple of interesting modules. They are:

• Fiona: OGR's neater API
• Rtree: spatial index for Python GIS
• Shapely: Python package for manipulation and analysis of features in the Cartesian plane

Personally I started to play around with GDAL/OGR lately and found them very impressive in respect to speed and coverage of analysis tools.

Here some examples of how to use the methods (taken from this excellent source which is a very good starting point):

# To select by attribute:
.SetAttributeFilter("soil = 'clay'")

# To select by location, either:
.SetSpatialFilter(<geom>)   

# or
.SetSpatialFilterRect(<minx>, <miny>, <maxx>, <maxy>)

# DataSource objects have a method `ExecuteSQL(<SQL>)`
.ExecuteSQL("SELECT* FROM sites WHERE soil = 'clay' ORDER BY id DESC")


# Plus all the well known tools, like:

# intersect
poly2.Intersect(<geom_1>)

# disjoint?
<geom>.Disjoint(geom_1)

# touches (on the edge?)
<geom>.Touches(geom_1)

# cross each other?
<geom>.Crosses(geom_1)

# within?
<geom>.Within(geom_1)

#contains?
<geom>.Contains(ptB)

# overlaps?
<geom>.Overlaps(geom_1)

## geoprecessing
<geom>.Union(<geom_1>)
<geom>.Intersection(<geom_1>)
<geom>.Difference(<geom_1>)
<geom>.SymmetricDifference(<geom_1>)

# Buffer (returns a new geometry)
<geom>.Buffer(<distance>)

# Are the geometries equal?
<geom1>.Equal(<geom2>)

# Returns the shortest distance between the two geometries
<geom1>.Distance(<geom2>)

# Returns the geometry's extent as a list (minx, maxx, miny, maxy)
<geom>.GetEnvelope()

The nice thing about these tools is that you are very flexible in how to implement them. I wrote for instance my own class CreateGeometry() to easily create vector files from the scratch. If your interested I can also post it here, even though I think it is beyond the scope of the question.

Answer 7

I know your question is Python-centric, but R has a wealth of value statistical analysis methods, some of which can be used for spatial analysis. @Whuber has a good answer here illustrating how to clip a raster to a box in two lines.

Answer 8

My solution, the quick solution, is to use GDAL with Python.

You need to

import subprocess

command = "gdalwarp -of GTiff -cutline clipArea.shp -cl area_of_interest -crop_to_cutline inData.asc outData.tiff"

subprocess.call(['C:\Temp\a b c\Notepad.exe'])

(From answer here: Clipping raster with vector layer using GDAL)

Of course, you should be able to achieve this using pure Python, but I have not needed to do it. And I almost always have GDAL around! The flexibility of GDAL is fantastic, especially in a linux environment. It handles huge rasters, it can be tied together with Python or Shell scripts and there are functions for many things. See also OGR for vector based tools.

Answer 9

If you don't mind running PostGIS it can do most spatial data processing for you.

PDF cheatsheet:

http://www.postgis.us/downloads/postgis20_cheatsheet.pdf

It integrates with python:

https://publicwiki.deltares.nl/display/OET/Accessing+PostgreSQL+PostGIS+with+Python

With supporting tools like SPIT within Quantum GIS or pgAdmin you are well equipped to set up PostGIS. You can then use python control the PostGIS operations on your spatial data.