Python Swig Examples
Dilemma Content to be removed or conserved as a historical reference
Important: SWIG has been removed from GRASS 6.4+ in favour of "ctypes" which
provide a more stable and portable environment. This page needs to be
updated accordingly.
Advanced SWIG Python Interface Examples
- Simpler examples are given in the SWIG-Python README file and the GRASS and Python wiki page.
Overview
As you might expect, GRASS's C library functions preform many common GIS tasks using very well tested and mature code. Reimplementing them in Python would thus be a waste of programming effort and would introduce potentially buggy and slower code. (Well written C is hard to beat for speed) A major advantage of using the SWIG interface can be summarized as "the experts have already solved this problem for you, take advantage of that."
Passing arrays and specific data types
Some of GRASS's functions want to be passed arrays of numbers or of a specific data type (e.g. integer, double). To do this from Python the 3rd party NumPtr module is used to access a memory pointer object which can be passed to the function.
- (any way to then free() the memory? Does "del SomePtr" do it?)
NumPtr in turn depends on the more common Numeric python (NumPy) module. See the example on the GRASS and Python wiki page.
Setting up the NumPtr module
(a copy now resides in swig/python/NumPtr/)
NumPtr setup:
# NumPtr - Numeric Pointer Module for Python (GPL2) # http://geosci.uchicago.edu/csc/numptr/ # 23k .tgz ; 100k installed wget http://geosci.uchicago.edu/csc/numptr/NumPtr-1.1.tar.gz tar xzf NumPtr-1.1.tar.gz cd NumPtr-1.1 python setup.py build #python setup.py install --prefix=/home/user cp build/lib.linux-i686-2.4/*NumPtr.* /usr/src/grass64/swig/python/
Examples
Distance and area calculations
Once NumPtr is installed we can run our script.
GRASS's libgis (C API) distance and area functions automatically switch to using geodetic calculations when in a Lat/Lon location.
The following calculates the area of the default Spearfish region bounds and the distance between the region's far corners.
Spearfish uses a UTM (planimetric) projection.
m.distance.demo
Note: For the full version of m.distance have a look at swig/python/examples/m.distance in the GRASS source code. It allows multi-segment input from coordinate pairs piped from stdin and DDD:MM:SS.SSS formatted coordinates for Lat/Lon.
#!/usr/bin/python
############################################################################
#
# MODULE: m.distance.demo
#
# AUTHOR(S): Hamish Bowman, Dunedin, New Zealand
#
# PURPOSE: Find distance between two points
# Uses geodetic distance calculations for Lat/Lon locations
# Demonstrates GRASS SWIG-Python interface
#
# COPYRIGHT: (c) 2008 The GRASS Development Team
#
# This program is free software under the GNU General Public
# License (>=v2). Read the file COPYING that comes with GRASS
# for details.
#
############################################################################
#%Module
#% description: Find distance between two points.
#% keywords: miscellaneous, distance, measure
#%End
#%Option
#% key: coord
#% type: string
#% required: yes
#% multiple: no
#% key_desc: x1,y1,x2,y2
#% description: Starting and ending coordinates
#%End
import sys
import os
def main():
#### add your code here ####
# self calculated answers:
#
# G63> g.region -d && g.region -g
# G63> g.region -e
# north-south extent: 14310.000000
# east-west extent: 19020.000000
# calc length of hypotenuse:
# >>> pow( pow(14310,2) + pow(19020,2), 0.5 )
# 23802.027224587404
# calc area of current region
# >>> 14310 * 19020
# 272176200
#os.putenv("GIS_OPT_COORD", '609000,4913700,589980,4928010')
coords = os.getenv("GIS_OPT_COORD").split(',')
x1 = float(coords[0])
y1 = float(coords[1])
x2 = float(coords[2])
y2 = float(coords[3])
# run this before starting python to append module search path:
# export PYTHONPATH=/usr/src/grass63/swig/python
# check with "import sys; sys.path"
# or:
sys.path.append("/usr/src/grass63/swig/python")
# put this inside main() so --help isn't slow.
import python_grass6 as g6lib
g6lib.G_gisinit('m.distance')
# returns 0 on success
### calc distance ###
g6lib.G_begin_distance_calculations()
# returns 0 if projection has no metrix (ie. imagery)
# returns 1 if projection is planimetric
# returns 2 if projection is latitude-longitude
distance = g6lib.G_distance(x1, y1, x2, y2)
print "distance is", distance
# 23802.0272246 (ok, matches above calc.)
### calc area ###
g6lib.G_begin_polygon_area_calculations()
# returns 0 if the projection is not measurable (ie. imagery or xy)
# returns 1 if the projection is planimetric (ie. UTM or SP)
# returns 2 if the projection is non-planimetric (ie. latitude-longitude)
# we don't need this, but just to have a look
g6lib.G_database_units_to_meters_factor()
# 1.0
# passing an array of values
import Numeric
import NumPtr
# do not need to close polygon (but it doesn't hurt if you do)
x = [ x1, x2, x2, x1 ]
y = [ y1, y1, y2, y2 ]
npoints = len(x)
# unset variables:
#del [Xs, Xptr, Ys, Yptr]
# or
#Xs = Xptr = Ys = Yptr = None
Xs = Numeric.array(x, Numeric.Float64)
Xptr = NumPtr.getpointer(Xs)
Ys = Numeric.array(y, Numeric.Float64)
Yptr = NumPtr.getpointer(Ys)
area = g6lib.G_area_of_polygon(Xptr, Yptr, npoints)
print "area is", area
# 272176200.0 (ok, matches above calc)
#### end of your code ####
return
if __name__ == "__main__":
if ( len(sys.argv) <= 1 or sys.argv[1] != "@ARGS_PARSED@" ):
os.execvp("g.parser", [sys.argv[0]] + sys.argv)
else:
main();