m.distance.py 5.6 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203
  1. #!/usr/bin/python
  2. ############################################################################
  3. #
  4. # MODULE: m.distance
  5. #
  6. # AUTHOR(S): Hamish Bowman, Dunedin, New Zealand
  7. #
  8. # PURPOSE: Find distance between two points
  9. # If the projection is latitude-longitude, this distance
  10. # is measured along the geodesic.
  11. # Demonstrates GRASS SWIG-Python interface
  12. #
  13. # COPYRIGHT: (c) 2008 Hamish Bowman, and The GRASS Development Team
  14. #
  15. # This program is free software under the GNU General Public
  16. # License (>=v2). Read the file COPYING that comes with GRASS
  17. # for details.
  18. #
  19. ############################################################################
  20. #
  21. # Requires GRASS SWIG-Python interface
  22. # Requires Numeric module (NumPy) from http://numpy.scipy.org/
  23. # Requires NumPrt module from http://geosci.uchicago.edu/csc/numptr/
  24. #
  25. #%Module
  26. #% label: Finds the distance between two or more points.
  27. #% description: If the projection is latitude-longitude, this distance is measured along the geodesic.
  28. #% keywords: miscellaneous, distance, measure
  29. #%End
  30. #%Option
  31. #% key: coord
  32. #% type: string
  33. #% required: no
  34. #% multiple: yes
  35. #% key_desc: x,y
  36. #% description: Comma separated list of coordinate pairs
  37. #%End
  38. #%Flag
  39. #% key: i
  40. #% description: Read coordinate pairs from stdin
  41. #%End
  42. import os, sys
  43. if not os.environ.has_key("GISBASE"):
  44. print "You must be in GRASS GIS to run this program."
  45. sys.exit(1)
  46. def main():
  47. #### add your code here ####
  48. # run this before starting python to append module search path:
  49. # export PYTHONPATH=/usr/src/grass70/swig/python
  50. # check with "import sys; sys.path"
  51. # or:
  52. sys.path.append("/usr/src/grass70/swig/python")
  53. # FIXME: install the g7lib.py bindings in $GISBASE/lib/ ?
  54. import swig.grass as g7lib
  55. # for passing pointers
  56. import Numeric
  57. import NumPtr
  58. g7lib.G_gisinit('m.distance')
  59. # returns 0 on success
  60. ### calc distance ###
  61. proj_type = g7lib.G_begin_distance_calculations()
  62. # returns 0 if projection has no metrix (ie. imagery)
  63. # returns 1 if projection is planimetric
  64. # returns 2 if projection is latitude-longitude
  65. # parser always creates at least an empty variable, and sys.argv is
  66. # toast, so no way to check if option was given. So it hangs if
  67. # --q was the only option given and there is no data from stdin.
  68. coord_ans = os.getenv("GIS_OPT_COORD")
  69. stdin_flag = bool(int(os.getenv("GIS_FLAG_I")))
  70. if stdin_flag is True:
  71. coords = []
  72. # read line by line from stdin
  73. while 1:
  74. line = sys.stdin.readline().strip()
  75. if not line: # EOF
  76. break
  77. else:
  78. coords += line.split(',')
  79. else:
  80. # read from coord= command line option
  81. coords = coord_ans.split(',')
  82. if len(coords) < 4:
  83. print "A minimum of two input coordinate pairs are needed"
  84. return
  85. # init variables
  86. overall_distance = 0.0
  87. if proj_type == 2:
  88. # lat/lon scan for DDD:MM:SS.SSSS
  89. easting = Numeric.array(0., Numeric.Float64)
  90. eastPtr = NumPtr.getpointer(easting)
  91. northing = Numeric.array(0., Numeric.Float64)
  92. northPtr = NumPtr.getpointer(northing)
  93. # TODO: for clarity, figure out how to replace "3" with
  94. # the defined LOCATION_LL constant from gis.i
  95. g7lib.G_scan_easting(coords[0], eastPtr, 3)
  96. g7lib.G_scan_northing(coords[1], northPtr, 3)
  97. x1 = float(easting)
  98. y1 = float(northing)
  99. else:
  100. # plain old coordinates
  101. x1 = float(coords[0])
  102. y1 = float(coords[1])
  103. x = [x1]
  104. y = [y1]
  105. for i in range(1, (len(coords) / 2)):
  106. if proj_type == 2:
  107. g7lib.G_scan_easting (coords[ i*2 + 0 ], eastPtr, 3)
  108. g7lib.G_scan_northing(coords[ i*2 + 1 ], northPtr, 3)
  109. x2 = float(easting)
  110. y2 = float(northing)
  111. else:
  112. x2 = float(coords[ i*2 + 0 ])
  113. y2 = float(coords[ i*2 + 1 ])
  114. segment_distance = g7lib.G_distance(x1, y1, x2, y2)
  115. overall_distance += segment_distance
  116. print "segment %d distance is %.2f meters" % (i, segment_distance)
  117. # add to the area array
  118. # setup for the next loop
  119. x1 = x2
  120. y1 = y2
  121. x += [x2]
  122. y += [y2]
  123. print
  124. print " total distance is %.2f meters" % overall_distance
  125. print
  126. ### calc area ###
  127. if len(coords) < 6:
  128. return
  129. g7lib.G_begin_polygon_area_calculations()
  130. # returns 0 if the projection is not measurable (ie. imagery or xy)
  131. # returns 1 if the projection is planimetric (ie. UTM or SP)
  132. # returns 2 if the projection is non-planimetric (ie. latitude-longitude)
  133. # do not need to close polygon (but it doesn't hurt if you do)
  134. npoints = len(x)
  135. # unset variables:
  136. #del [Xs, Xptr, Ys, Yptr]
  137. # or
  138. #Xs = Xptr = Ys = Yptr = None
  139. Xs = Numeric.array(x, Numeric.Float64)
  140. Xptr = NumPtr.getpointer(Xs)
  141. Ys = Numeric.array(y, Numeric.Float64)
  142. Yptr = NumPtr.getpointer(Ys)
  143. area = g7lib.G_area_of_polygon(Xptr, Yptr, npoints)
  144. print "AREA: %10.2f square meters" % area
  145. print
  146. # we don't need this, but just to have a look
  147. if False:
  148. if proj_type == 1:
  149. g7lib.G_database_units_to_meters_factor()
  150. # 1.0
  151. print "Location units are", g7lib.G_database_unit_name(True)
  152. #### end of your code ####
  153. return
  154. if __name__ == "__main__":
  155. if ( len(sys.argv) <= 1 or sys.argv[1] != "@ARGS_PARSED@" ):
  156. os.execvp("g.parser", [sys.argv[0]] + sys.argv)
  157. else:
  158. main();