Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

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?

share|improve this question
1  
Related: Clipping raster with vector layer using GDAL –  blah238 Sep 28 '12 at 16:52

10 Answers 10

up vote 22 down vote accepted
+100

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

share|improve this answer
1  
This is a great applied answer, thanks for putting this together. –  Aaron Oct 6 '12 at 15:39

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

share|improve this answer

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.

share|improve this answer

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.

share|improve this answer

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

share|improve this answer

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.

share|improve this answer

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.

share|improve this answer
    
New in GRASS GIS 7: pyGRASS, see ing.unitn.it/~zambelli/projects/pygrass –  markusN Apr 22 at 17:59
    

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.

share|improve this answer
    
Thanks for the USU link, it is extremely useful. –  Aaron Oct 7 '12 at 19:02

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.

share|improve this answer
6  
To bring it back to Python, you can use the RPy library. RPy is a very simple, yet robust, Python interface to the R Programming Language. It can manage all kinds of R objects and can execute arbitrary R functions (including the graphic functions). All errors from the R language are converted to Python exceptions. Any module installed for the R system can be used from within Python. –  RyanDalton Sep 28 '12 at 18:02
    
@Ryan Thanks for the excellent advice on the RPy library! R is another favorite programing language. –  Aaron Sep 28 '12 at 21:33

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.