Python Swig Examples

From GRASS-Wiki
Revision as of 23:00, 3 March 2008 by ⚠️HamishBowman (talk | contribs) (→‎m.distance: link to full m.distance module)
Jump to navigation Jump to search

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?)

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

  • 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();