Python Swig Examples

More Complicated SWIG Python Interface Examples

Overview

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

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? "del SomePtr"?)

Setting up the NumPtr module

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/

Examples

Distance and area calculations

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)