m.distance: major update for ctypes

git-svn-id: https://svn.osgeo.org/grass/grass/trunk@45895 15284696-431f-4ddb-bdfa-cd5b030d7da7

doc/python/m.distance.py

@@ -1,206 +1,149 @@
-#!/usr/bin/python
+#!/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 SWIG-Python interface
+#               Demonstrates GRASS Python Ctypes interface
 #
-# COPYRIGHT:    (c) 2008 Hamish Bowman, and The GRASS Development Team
+# 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.
+#               This program is free software under the GNU General
+#               Public License (>=v2). Read the file COPYING that
+#               comes with GRASS for details.
 #
 ############################################################################
 #
-# Requires GRASS SWIG-Python interface
+# Requires GRASS Python Ctypes interface
 # Requires Numeric module (NumPy) from  http://numpy.scipy.org/
 # Requires NumPrt module from  http://geosci.uchicago.edu/csc/numptr/
 #
 
-"""
-     ****  FIXME: needs to be ported from SWIG to Ctypes!  ****
-"""
-
-
-#%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: x,y
-#%  description: Comma separated list of coordinate pairs
-#%End
-#%Flag
-#%  key: i
-#%  description: Read coordinate pairs from stdin
-#%End
-
+#%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: yes
+#% multiple: yes
+#% key_desc: easting,northing
+#% description: Comma separated list of coordinate pairs
+#%end
+#%flag
+#% key: i
+#% description: Read coordinate pairs from stdin
+#% suppress_required: yes
+#%end
 
 import os, sys
 
-if not os.environ.has_key("GISBASE"):
-    print "You must be in GRASS GIS to run this program."
-    sys.exit(1)
+import grass.script as grass
 
+from grass.lib.gis import *
 
 def main(): 
+    G_gisinit('m.distance')
 
-    #### add your code here ####
-
-    # run this before starting python to append module search path:
-    #   export PYTHONPATH=$PYTHONPATH:/usr/local/grass-7.0.svn/etc/python
-    #   check with "import sys; sys.path"
-    # or:
-    sys.path.append("/usr/local/grass-7.0.svn/etc/python")
-    from grass.lib import gis as g7lib
- 
-    # for passing pointers
-    import Numeric
-    import NumPtr
-
-    g7lib.G_gisinit('m.distance')
-    # returns 0 on success
-
-
-    ### calc distance ###
+    # calc distance
     
-    proj_type = g7lib.G_begin_distance_calculations()
+    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.
-    coord_ans  = os.getenv("GIS_OPT_COORD")
-    stdin_flag = bool(int(os.getenv("GIS_FLAG_I")))
-
-    if stdin_flag is True:
-        coords = []
+    # 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 1:
+        while True:
             line = sys.stdin.readline().strip()
-            if not line:   # EOF
+            if not line: # EOF
                 break
             else:
-                coords += line.split(',')
+                coords.append(line.split(','))
     else:
         # read from coord= command line option
-        coords = coord_ans.split(',')
-
-
-    if len(coords) < 4:
-       print "A minimum of two input coordinate pairs are needed"
-       return
-
-
+        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 = Numeric.array(0., Numeric.Float64)
-        eastPtr = NumPtr.getpointer(easting)
-        northing = Numeric.array(0., Numeric.Float64)
-        northPtr = NumPtr.getpointer(northing)
-
-        # TODO: for clarity, figure out how to replace "3" with
-        #       the defined LOCATION_LL constant from gis.i
-        g7lib.G_scan_easting(coords[0], eastPtr, 3)
-        g7lib.G_scan_northing(coords[1], northPtr, 3)
-        x1 = float(easting)
-        y1 = float(northing)
+        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 old coordinates
-        x1 = float(coords[0])
-        y1 = float(coords[1])
-
-    x = [x1]
-    y = [y1]
-
-    for i in range(1, (len(coords) / 2)):
-
+        # 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:
-            g7lib.G_scan_easting (coords[ i*2 + 0 ], eastPtr, 3)
-            g7lib.G_scan_northing(coords[ i*2 + 1 ], northPtr, 3)
-            x2 = float(easting)
-            y2 = float(northing)
+            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:
-            x2 = float(coords[ i*2 + 0 ])
-            y2 = float(coords[ i*2 + 1 ])
-
-        segment_distance = g7lib.G_distance(x1, y1, x2, y2)
+            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
-
-        # setup for the next loop
-        x1 = x2
-        y1 = y2
-
-        x += [x2]
-        y += [y2]
-  
-    print
-    print " total distance is %.2f meters" % overall_distance
-    print
-
-
-    ### calc area ###
-    if len(coords) < 6:
-        return
- 
-    g7lib.G_begin_polygon_area_calculations()
+    
+    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)
-    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 = G_area_of_polygon(x, y, len(coords))
+    print "area is %.2f square meters\n" % area
     
-    area = g7lib.G_area_of_polygon(Xptr, Yptr, npoints)
-    print "AREA:  %10.2f square meters" % area
-    print
-
-
     # we don't need this, but just to have a look
-    if False:
-        if proj_type == 1:
-            g7lib.G_database_units_to_meters_factor()
-            # 1.0
-        print "Location units are", g7lib.G_database_unit_name(True)
-
-
-    #### end of your code ####
-    return 
+    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__":
-    if ( len(sys.argv) <= 1 or sys.argv[1] != "@ARGS_PARSED@" ):
-        os.execvp("g.parser", [sys.argv[0]] + sys.argv)
-    else:
-        main();
-
+    options, flags = grass.parser()
+    sys.exit(main())