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

#!/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)