[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