Python Swig Examples

From GRASS-Wiki
Jump to navigation Jump to search
  Important: SWIG has been removed from GRASS 6.4+ in favour of "ctypes" which 
             provide a more stable and portable environment. This page needs to be
             updated accordingly.

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

NumPtr in turn depends on the more common Numeric python (NumPy) module. See the example on the GRASS and Python wiki page.


Setting up the NumPtr module

(a copy now resides in swig/python/NumPtr/)

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/grass64/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.demo

Note: For the full version of m.distance have a look at swig/python/examples/m.distance in the GRASS source code. It allows multi-segment input from coordinate pairs piped from stdin and DDD:MM:SS.SSS formatted coordinates for Lat/Lon.

 #!/usr/bin/python
 ############################################################################
 #
 # MODULE:       m.distance.demo
 #
 # AUTHOR(S):    Hamish Bowman, Dunedin, New Zealand
 #
 # PURPOSE:      Find distance between two points
 #		 Uses geodetic distance calculations for Lat/Lon locations
 #		 Demonstrates GRASS SWIG-Python interface
 #
 # COPYRIGHT:    (c) 2008 The GRASS Development Team
 #
 #               This program is free software under the GNU General Public
 #               License (>=v2). Read the file COPYING that comes with GRASS
 #               for details.
 #
 ############################################################################
 
 #%Module
 #%  description: Find distance between two points.
 #%  keywords: miscellaneous, distance, measure
 #%End
 #%Option
 #%  key: coord
 #%  type: string
 #%  required: yes
 #%  multiple: no
 #%  key_desc: x1,y1,x2,y2
 #%  description: Starting and ending coordinates
 #%End
 
 import sys
 import os
 
 def main(): 
 
     #### add your code here ####
 
     # self calculated answers:
     #
     # 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
 
     #os.putenv("GIS_OPT_COORD", '609000,4913700,589980,4928010')
     coords = os.getenv("GIS_OPT_COORD").split(',')
 
     x1 = float(coords[0])
     y1 = float(coords[1])
     x2 = float(coords[2])
     y2 = float(coords[3])
 
     # run this before starting python to append module search path:
     #   export PYTHONPATH=/usr/src/grass63/swig/python
     #   check with "import sys; sys.path"
     # or:
     sys.path.append("/usr/src/grass63/swig/python")
     # put this inside main() so --help isn't slow.
     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
     
     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)
 
     #### end of your code ####
     return 
 
 if __name__ == "__main__":
     if ( len(sys.argv) <= 1 or sys.argv[1] != "@ARGS_PARSED@" ):
         os.execvp("g.parser", [sys.argv[0]] + sys.argv)
     else:
 	main();

See also