Python Swig Examples: Difference between revisions

From GRASS-Wiki
Jump to navigation Jump to search
m (→‎Overview: change tone)
Line 7: Line 7:


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)
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 to GRASS's C libraries is to ''"Let the experts solve the problem for you."''
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 ==
== Passing arrays and specific data types ==

Revision as of 09:37, 3 March 2008

Advanced SWIG Python Interface Examples


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? "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/grass63/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:

#!/usr/bin/python
# m.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/grass63/swig/python
#   check with "import sys; sys.path"
# or:
import sys
sys.path.append("/usr/src/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 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)