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?

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:

• Python Script examples for geoprocessing shapefiles without using arcpy
• 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

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 6

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

I have been working on an open-source geoprocessing library called WhiteboxTools that can be used in place of ArcPy in many applications. Currently there are nearly 300 tools available for processing raster, vector, and LiDAR (LAS) data, although the plan is to eventually port over all of the 400+ tools available in Whitebox GAT. Although the tools are developed using the Rust programming language (for efficiency), each tool is callable from Python, as in the following example:

from whitebox_tools import WhiteboxTools

wbt = WhiteboxTools()

# Set the working directory. This is the path to the folder containing the data,
# i.e. files sent to tools as input/output parameters. You don't need to set
# the working directory if you specify full path names as tool parameters.
wbt.work_dir = "/path/to/data/"

# The most convenient way to run a tool is to use its associated method, e.g.:
wbt.elev_percentile("DEM.tif", "output.tif", 15, 15)

# You may also provide an optional custom callback for processing output from the
# tool. If you don't provide a callback, and verbose is set to True, tool output
# will simply be printed to the standard output.
def my_callback(value):
    if user_selected_cancel_btn: # Assumes a 'Cancel' button on a GUI
        print('Cancelling operation...')
        wbt.cancel_op = True
    else:
        print(value)

wbt.breach_depressions('DEM.flt', 'DEM_breached.flt', callback=my_callback)

# List all available tools in WhiteboxTools
print(wbt.list_tools())

# Lists tools with 'lidar' or 'LAS' in tool name or description.
print(wbt.list_tools(['lidar', 'LAS']))

# Print the help for a specific tool.
print(wbt.tool_help("ElevPercentile"))

# Want to read the source code for a tool?
# 'view_code' opens a browser and navigates to a tool's  
# source code in the WhiteboxTools GitHub repository
wbt.view_code('watershed')

More detailed information can be found provided in the WhiteboxTools user manual. The library is stand-alone and does not have any other dependencies. You simply need to download the small (< 5Mb) file located here. The download file contains the WhiteboxTools exe, the whitebox_tools.py script, which provides the Python API for the library (imported on the top line of the above script), and the user manual. There is also a very basic tkinter GUI (wb_runner.py) for interfacing with the library.

The permissive MIT licence is intended to allow WhiteboxTools to be integrated as a back-end with other open-source GIS; Alexander Bruy has developed a QGIS plugin for the WhiteboxTools back-end. You may also mix and match tools from WhiteboxTools and ArcPy in a single script as needed. The library is still somewhat experimental, developed out of the University of Guelph's Geomorphometry and Hydrogeomatics Research Group, and is currently pre-1.0 release, which should be taken into account in usage.

Answer 10

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.

Answer 11

Using Python to clip a raster to a shapefile without ArcPy: http://geospatialpython.com/2011/02/clip-raster-using-shapefile.html

Answer 12

Two amazing tools that aren't mentioned here are:

GeoPandas. This makes so much GIS processing in Python more efficient, especially if you're interested in Pandas-like processing and visualization.

PYSAL. There are a lot of great geostatistical/geospatial data science tools in this package, and I find that it accomplishes a lot of what ArcPy does, but in an open-source framework. Well worth checking out.