Python Swig Examples

From GRASS-Wiki
Revision as of 06:25, 3 March 2008 by ⚠️HamishBowman (talk | contribs) (init)
(diff) ← Older revision | Latest revision (diff) | Newer revision → (diff)
Jump to navigation Jump to search

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)