[GRASS-SVN] r55932 - grass/branches/releasebranch_6_4/doc/python
svn_grass at osgeo.org
svn_grass at osgeo.org
Sat Apr 20 07:56:11 PDT 2013
Author: neteler
Date: 2013-04-20 07:56:11 -0700 (Sat, 20 Apr 2013)
New Revision: 55932
Added:
grass/branches/releasebranch_6_4/doc/python/m.distance.py
grass/branches/releasebranch_6_4/doc/python/raster_example_ctypes.py
grass/branches/releasebranch_6_4/doc/python/vector_example_ctypes.py
Removed:
grass/branches/releasebranch_6_4/doc/python/example_ctypes.py
Log:
backport of Python ctypes usage examples
Deleted: grass/branches/releasebranch_6_4/doc/python/example_ctypes.py
===================================================================
--- grass/branches/releasebranch_6_4/doc/python/example_ctypes.py 2013-04-20 14:28:03 UTC (rev 55931)
+++ grass/branches/releasebranch_6_4/doc/python/example_ctypes.py 2013-04-20 14:56:11 UTC (rev 55932)
@@ -1,49 +0,0 @@
-#!/usr/bin/env python
-import os, sys, subprocess
-from ctypes import *
-grass = CDLL("libgrass_gis.so")
-
-if not os.environ.has_key("GISBASE"):
- print "You must be in GRASS GIS to run this program."
- sys.exit(1)
-
-if len(sys.argv)==2:
- input = sys.argv[1]
-else:
- input = raw_input("Raster Map Name? ")
-
-# initialize
-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, '')
-
-# find map in search path
-mapset = grass.G_find_cell2(input, '')
-mapset = c_char_p(mapset).value
-
-# determine the inputmap type (CELL/FCELL/DCELL) */
-data_type = grass.G_raster_map_type(input, mapset)
-
-if data_type == 0:
- ptype = POINTER(c_int)
-elif data_type == 1:
- ptype = POINTER(c_float)
-elif data_type == 2:
- ptype = POINTER(c_double)
-
-infd = grass.G_open_cell_old(input, mapset)
-inrast = grass.G_allocate_raster_buf(data_type)
-inrast = cast(c_void_p(inrast), ptype)
-
-rows = grass.G_window_rows()
-cols = grass.G_window_cols()
-
-for rown in xrange(rows):
- grass.G_get_raster_row(infd, inrast, rown, data_type)
- print rown, inrast[0:cols]
-
-grass.G_close_cell(infd)
-grass.G_free(inrast)
-
Copied: grass/branches/releasebranch_6_4/doc/python/m.distance.py (from rev 55931, grass/branches/develbranch_6/doc/python/m.distance.py)
===================================================================
--- grass/branches/releasebranch_6_4/doc/python/m.distance.py (rev 0)
+++ grass/branches/releasebranch_6_4/doc/python/m.distance.py 2013-04-20 14:56:11 UTC (rev 55932)
@@ -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())
Copied: grass/branches/releasebranch_6_4/doc/python/raster_example_ctypes.py (from rev 55931, grass/branches/develbranch_6/doc/python/raster_example_ctypes.py)
===================================================================
--- grass/branches/releasebranch_6_4/doc/python/raster_example_ctypes.py (rev 0)
+++ grass/branches/releasebranch_6_4/doc/python/raster_example_ctypes.py 2013-04-20 14:56:11 UTC (rev 55932)
@@ -0,0 +1,79 @@
+#!/usr/bin/env python
+
+"""
+Sample Python script to access raster data using GRASS Ctypes
+interface
+
+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
+import sys
+
+from grass.lib.gis import *
+
+# check if GRASS is running or not
+if not os.environ.has_key("GISBASE"):
+ sys.exit("You must be in GRASS GIS to run this program")
+
+# 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("Name of raster map? ")
+
+# initialize GRASS library
+G_gisinit('')
+
+# find map in search path
+mapset = G_find_cell2(input, '')
+if not mapset:
+ sys.exit("Raster map <%s> not found" % input)
+
+# determine the inputmap type (CELL/FCELL/DCELL)
+data_type = G_raster_map_type(input, mapset)
+
+if data_type == CELL_TYPE:
+ ptype = POINTER(c_int)
+ type_name = 'CELL'
+elif data_type == FCELL_TYPE:
+ ptype = POINTER(c_float)
+ type_name = 'FCELL'
+elif data_type == DCELL_TYPE:
+ ptype = POINTER(c_double)
+ type_name = 'DCELL'
+
+print "Raster map <%s> contains data type %s." % (input, type_name)
+
+in_fd = G_open_cell_old(input, mapset)
+in_rast = G_allocate_raster_buf(data_type)
+in_rast = cast(c_void_p(in_rast), ptype)
+
+rows = G_window_rows()
+cols = G_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
+ G_get_raster_row(in_fd, in_rast, row_n, data_type)
+ print row_n, in_rast[0:cols]
+ # TODO check for NULL
+
+# closed map and cleanup memory allocation
+G_close_cell(in_fd)
+G_free(in_rast)
+
+def check_null(value):
+ if math.isnan(value):
+ return 'null'
+ return value
Copied: grass/branches/releasebranch_6_4/doc/python/vector_example_ctypes.py (from rev 55931, grass/branches/develbranch_6/doc/python/vector_example_ctypes.py)
===================================================================
--- grass/branches/releasebranch_6_4/doc/python/vector_example_ctypes.py (rev 0)
+++ grass/branches/releasebranch_6_4/doc/python/vector_example_ctypes.py 2013-04-20 14:56:11 UTC (rev 55932)
@@ -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