[GRASS-SVN] r45896 - grass/branches/develbranch_6/doc/python

svn_grass at osgeo.org svn_grass at osgeo.org
Sun Apr 10 17:55:22 EDT 2011


Author: martinl
Date: 2011-04-10 14:55:22 -0700 (Sun, 10 Apr 2011)
New Revision: 45896

Added:
   grass/branches/develbranch_6/doc/python/m.distance.py
   grass/branches/develbranch_6/doc/python/vector_example_ctypes.py
Modified:
   grass/branches/develbranch_6/doc/python/raster_example_ctypes.py
Log:
backport ctypes examples from trunk


Copied: grass/branches/develbranch_6/doc/python/m.distance.py (from rev 45895, grass/trunk/doc/python/m.distance.py)
===================================================================
--- grass/branches/develbranch_6/doc/python/m.distance.py	                        (rev 0)
+++ grass/branches/develbranch_6/doc/python/m.distance.py	2011-04-10 21:55:22 UTC (rev 45896)
@@ -0,0 +1,148 @@
+#!/usr/bin/env python
+############################################################################
+#
+# MODULE:       m.distance
+#
+# AUTHOR(S):    Hamish Bowman, Dunedin, New Zealand
+#               Updated for Ctypes by Martin Landa <landa.martin gmail.com>
+#
+# PURPOSE:      Find distance between two points
+#               If the projection is latitude-longitude, this distance
+#                 is measured along the geodesic.
+#               Demonstrates GRASS Python Ctypes interface
+#
+# COPYRIGHT:    (c) 2008, 2011 Hamish Bowman, and 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.
+#
+############################################################################
+#
+# Requires GRASS Python Ctypes interface
+# Requires Numeric module (NumPy) from  http://numpy.scipy.org/
+# Requires NumPrt module from  http://geosci.uchicago.edu/csc/numptr/
+#
+
+#%module
+#% label: Finds the distance between two or more points.
+#% description: If the projection is latitude-longitude, this distance is measured along the geodesic.
+#% keywords: miscellaneous, distance, measure
+#%end
+#%option
+#% key: coord
+#% type: string
+#% required: no
+#% multiple: yes
+#% key_desc: easting,northing
+#% description: Comma separated list of coordinate pairs
+#%end
+#%flag
+#% key: i
+#% description: Read coordinate pairs from stdin
+#%end
+
+import os, sys
+
+import grass.script as grass
+
+from grass.lib.gis import *
+
+def main(): 
+    G_gisinit('m.distance')
+
+    # calc distance
+    
+    proj_type = 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
+
+    # parser always creates at least an empty variable, and sys.argv is
+    # toast, so no way to check if option was given. So it hangs if
+    # --q was the only option given and there is no data from stdin.
+    coords = []
+    if flags['i']:
+        # read line by line from stdin
+        while True:
+            line = sys.stdin.readline().strip()
+            if not line: # EOF
+                break
+            else:
+                coords.append(line.split(','))
+    else:
+        # read from coord= command line option
+        p = None
+        for c in options['coord'].split(','):
+            if not p:
+                p = [c]
+            else:
+                p.append(c)
+                coords.append(p)
+                p = None
+    
+    if len(coords) < 2:
+       grass.fatal("A minimum of two input coordinate pairs are needed")
+    
+    # init variables
+    overall_distance = 0.0
+    coord_array = c_double * len(coords)
+    x = coord_array()
+    y = coord_array()
+    if proj_type == 2:
+        # lat/lon scan for DDD:MM:SS.SSSS
+        easting = c_double()
+        northing = c_double()
+        G_scan_easting(coords[0][0], byref(easting), PROJECTION_LL)
+        G_scan_northing(coords[0][1], byref(northing), PROJECTION_LL)
+        x[0] = float(easting.value)
+        y[0] = float(northing.value)
+    else:
+        # plain coordinates
+        x[0] = float(coords[0][0])
+        y[0] = float(coords[0][1])
+    
+    for i in range(1, len(coords)):
+        if proj_type == 2:
+            easting = c_double()
+            northing = c_double()
+            G_scan_easting(coords[i][0], byref(easting), PROJECTION_LL)
+            G_scan_northing(coords[i][1], byref(northing), PROJECTION_LL)
+            x[i] = float(easting.value)
+            y[i] = float(northing.value)
+        else:
+            x[i] = float(coords[i][0])
+            y[i] = float(coords[i][1])
+        
+        segment_distance = G_distance(x[i-1], y[i-1], x[i], y[i])
+        overall_distance += segment_distance
+        
+        print "segment %d distance is %.2f meters" % (i, segment_distance)
+        
+        # add to the area array
+    
+    print "\ntotal distance is %.2f meters\n" % overall_distance
+    
+    # calc area
+    if len(coords) < 3:
+       return 0
+    
+    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)
+
+    # do not need to close polygon (but it doesn't hurt if you do)
+    area = G_area_of_polygon(x, y, len(coords))
+    print "area is %.2f square meters\n" % area
+    
+    # we don't need this, but just to have a look
+    if proj_type == 1:
+        G_database_units_to_meters_factor()
+        grass.message("Location units are %s" % G_database_unit_name(True).lower())
+    
+    return 0
+
+if __name__ == "__main__":
+    options, flags = grass.parser()
+    sys.exit(main())

Modified: grass/branches/develbranch_6/doc/python/raster_example_ctypes.py
===================================================================
--- grass/branches/develbranch_6/doc/python/raster_example_ctypes.py	2011-04-10 21:47:05 UTC (rev 45895)
+++ grass/branches/develbranch_6/doc/python/raster_example_ctypes.py	2011-04-10 21:55:22 UTC (rev 45896)
@@ -1,94 +1,80 @@
 #!/usr/bin/env python
 
 """
-This is an example of how to use "ctypes" to access GRASS's C API from
-within a Python script. See the GRASS Programmer's manual for details
-on the e.g. G_*() functions in the GRASS C libraries (libgis et al.).
+Sample Python script to access raster data using GRASS Ctypes
+interface
 
-USAGE:  example_ctypes.py [Raster map name]
-        If a raster map name is not given it will prompt you for one.
-	This raster map is then opened, read, and data values printed
-	to stdout. You may wish to use 'g.region' to set the rows x
-	columns to something small (say 10 x 5) to avoid overwhelming
-	yourself: `g.region rows=10 cols=5`
+You may wish to use 'g.region' to set the rows x columns to something
+small (say 10 x 5) to avoid overwhelming yourself: `g.region rows=10
+cols=5`
 """
 
-
 # FIXME: as an example it should make extensive use of code comments and document
 #  each and every step along the way.  (e.g. explain c_char_p().value memory pointer
 #  to string conversion for Python programmers not familar with C pointers)
 #
 #  FIXME: explain at a basic level what ctypes is & does.
 
-import os, sys, subprocess
-# actiavate ctypes
-from ctypes import *
+import os
+import sys
 
+from grass.lib.gis    import *
+from grass.lib.raster import *
+
 # check if GRASS is running or not
 if not os.environ.has_key("GISBASE"):
-    print "You must be in GRASS GIS to run this program."
-    sys.exit(1)
+    sys.exit("You must be in GRASS GIS to run this program")
 
-# load in the main GRASS C library
-grass = CDLL("libgrass_gis.so")
-
-# parse command line arguements, prompt user for a raster map name if one wasn't given
-if len(sys.argv)==2:
+# parse command line arguments, prompt user for a raster map name if one wasn't given
+if len(sys.argv) == 2:
   input = sys.argv[1]
 else:
-  input = raw_input("Raster Map Name? ")
+  input = raw_input("Name of raster map? ")
 
-# initialize the C library
-s = subprocess.Popen(['g.version','-r'], stdout=subprocess.PIPE).communicate()[0]
-for line in s.splitlines():
-    if line.startswith('Revision:'):
-        version = '$' + line + '$'
-grass.G__gisinit(version, '')
+# initialize GRASS library
+G_gisinit('')
 
-
 # find map in search path
-mapset = grass.G_find_cell2(input, '')
-mapset = c_char_p(mapset).value
-
+mapset = G_find_raster2(input, '')
 if not mapset:
-    print "Raster map <%s> not found." % input
-    sys.exit(1)
+    sys.exit("Raster map <%s> not found" % input)
 
-
 # determine the inputmap type (CELL/FCELL/DCELL)
-data_type = grass.G_raster_map_type(input, mapset)
+data_type = Rast_map_type(input, mapset)
 
-if data_type == 0:
+if data_type == CELL_TYPE:
     ptype = POINTER(c_int)
     type_name = 'CELL'
-elif data_type == 1:
+elif data_type == FCELL_TYPE:
     ptype = POINTER(c_float)
     type_name = 'FCELL'
-elif data_type == 2:
+elif data_type == DCELL_TYPR:
     ptype = POINTER(c_double)
     type_name = 'DCELL'
 
 print "Raster map <%s> contains data type %s." % (input, type_name)
 
-in_fd = grass.G_open_cell_old(input, mapset)
-in_rast = grass.G_allocate_raster_buf(data_type)
+in_fd   = Rast_open_old(input, mapset)
+in_rast = Rast_allocate_buf(data_type)
 in_rast = cast(c_void_p(in_rast), ptype)
 
-rows = grass.G_window_rows()
-cols = grass.G_window_cols()
-print "Current region is %d rows x %d columns.\n" % (rows, cols)
+rows = Rast_window_rows()
+cols = Rast_window_cols()
+print "Current region is %d rows x %d columns" % (rows, cols)
 
 # iterate through map rows
 print "Map data:"
 for row_n in xrange(rows):
     # read a row of raster data into memory, then print it
-    grass.G_get_raster_row(in_fd, in_rast, row_n, data_type)
+    Rast_get_row(in_fd, in_rast, row_n, data_type)
     print row_n, in_rast[0:cols]
+    # TODO check for NULL
 
-#TODO: NULL -> NaN
-
-
 # closed map and cleanup memory allocation
-grass.G_close_cell(in_fd)
-grass.G_free(in_rast)
+Rast_close(in_fd)
+G_free(in_rast)
 
+def check_null(value):
+    if math.isnan(value):
+        return 'null'
+    return value

Copied: grass/branches/develbranch_6/doc/python/vector_example_ctypes.py (from rev 45892, grass/trunk/doc/python/vector_example_ctypes.py)
===================================================================
--- grass/branches/develbranch_6/doc/python/vector_example_ctypes.py	                        (rev 0)
+++ grass/branches/develbranch_6/doc/python/vector_example_ctypes.py	2011-04-10 21:55:22 UTC (rev 45896)
@@ -0,0 +1,60 @@
+#!/usr/bin/python
+
+"""
+Sample Python script to access vector data using GRASS Ctypes
+interface
+"""
+
+import os, sys
+
+from grass.lib.gis    import *
+from grass.lib.vector import *
+
+if not os.environ.has_key("GISBASE"):
+    sys.exit("You must be in GRASS GIS to run this program.")
+
+if len(sys.argv) == 2:
+  input = sys.argv[1]
+else:
+  input = raw_input("Name of vector map? ")
+
+# initialize GRASS library
+G_gisinit('')
+
+# find vector map in the search path
+mapset = G_find_vector2(input, "")
+if not mapset:
+    sys.exit("Vector map <%s> not found" % input)
+
+# define map structure
+map_info = pointer(Map_info())
+
+# define open level (level 2: topology)
+Vect_set_open_level(2)
+
+# open existing vector map
+Vect_open_old(map_info, input, mapset)
+
+# query
+print 'Vector map        :', Vect_get_full_name(map_info)
+print 'Vector is 3D      :', Vect_is_3d(map_info)
+print 'Vector DB links   :', Vect_get_num_dblinks(map_info)
+print 'Map Scale         : 1:%d' % Vect_get_scale(map_info)
+
+# vector box tests
+box = bound_box()
+c_easting1  = 599505.0
+c_northing  = 4921010.0
+c_easting2  = 4599505.0
+
+Vect_get_map_box(map_info, byref(box))
+print 'Position 1 in box ?', Vect_point_in_box(c_easting1, c_northing, 0, byref(box))
+print 'Position 2 in box ?', Vect_point_in_box(c_easting2, c_northing, 0, byref(box))
+
+print 'Number of features:', Vect_get_num_lines(map_info)
+print 'Number of points  :', Vect_get_num_primitives(map_info, GV_POINT)
+print 'Number of lines   :', Vect_get_num_primitives(map_info, GV_LINE)
+print 'Number of areas   :', Vect_get_num_areas(map_info)
+
+# close map
+Vect_close(map_info)



More information about the grass-commit mailing list