Python Swig Examples
More Complicated SWIG Python Interface Examples
Distance and area calculations
GRASS's libgis (C API) distance and area functions automatically switch to using geodetic calculations when in a Lat/Lon location. As you might expect, they preform many common GIS functions in very well tested and mature code. Reimplementing them in python is a waste of time and potentially buggy. The idea of using the SWIG interface to GRASS's C libraries is to "Let the experts solve the problem for you."
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? "del SomePtr"?)
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/grass/grass63/swig/python/
Once NumPtr is installed with can run our script:
#!/usr/bin/python
# g.distance -- demo SWIG interface
# (c) 2008 Hamish Bowman, and the GRASS Development Team
# Licensed as GPL >=2
# run this before starting python to append module search path:
# export PYTHONPATH=/usr/src/grass/svn/grass63/swig/python
# check with "import sys; sys.path"
# or:
import sys
sys.path.append("/usr/src/grass/svn/grass63/swig/python")
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
# 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
x1 = 609000
x2 = 589980
y1 = 4913700
y2 = 4928010
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 dont' 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
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)