Re-implement r.plane in Python and r.mapcalc

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

scripts/r.plane/r.plane.py

@@ -0,0 +1,127 @@
+#!/usr/bin/env python
+#
+############################################################################
+#
+# MODULE:	r.plane for GRASS 5.7; based on r.plane for GRASS 5
+# AUTHOR(S):	CERL?; updated to GRASS 5.7 by Michael Barton
+#               Dec 2004: Alessandro Frigeri & Ivan Marchesini
+#               Modified to produce floating and double values maps
+#               Converted to Python by Glynn Clements
+# PURPOSE:	Creates a raster plane map from user specified inclination and azimuth
+# COPYRIGHT:	(C) 2004,2008 by 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.
+#
+#############################################################################
+
+#%Module
+#%  description: Creates raster plane map given dip (inclination), aspect (azimuth) and one point.
+#%  keywords: raster, elevation
+#%End
+#%option
+#% key: name
+#% type: string
+#% gisprompt: new,cell,raster
+#% description: Name of raster plane to be created
+#% required : yes
+#%end
+#%option
+#% key: dip
+#% type: double
+#% gisprompt: -90-90
+#% answer: 0.0
+#% description: Dip of plane. Value must be between -90 and 90 degrees
+#% required : yes
+#%end
+#%option
+#% key: azimuth
+#% type: double
+#% gisprompt: 0-360
+#% answer: 0.0
+#% description: Azimuth of the plane. Value must be between 0 and 360 degrees
+#% required : yes
+#%end
+#%option
+#% key: easting
+#% type: double
+#% description: Easting coordinate of a point on the plane
+#% required : yes
+#%end
+#%option
+#% key: northing
+#% type: double
+#% description: Northing coordinate of a point on the plane
+#% required : yes
+#%end
+#%option
+#% key: elevation
+#% type: double
+#% description: Elevation coordinate of a point on the plane
+#% required : yes
+#%end
+#%option
+#% key: type
+#% type: string 
+#% options: int,float,double
+#% description: Type of the raster map to be created
+#% required : yes
+#%end
+
+import sys
+import os
+import math
+import string
+import grass
+
+def main():
+    name = options['name']
+    type = options['type']
+    dip = float(options['dip'])
+    az  = float(options['azimuth'])
+    ea  = float(options['easting'])
+    no  = float(options['northing'])
+    el  = float(options['elevation'])
+
+    reg = grass.region()
+
+    ### test input values ###
+    if abs(dip) >= 90:
+	grass.fatal("dip must be between -90 and 90.")
+
+    if az < 0 or az >= 360:
+	grass.fatal("azimuth must be between 0 and 360")
+
+    ### now the actual algorithm
+    az_r  = math.radians(-az)
+    sinaz = math.sin(az_r)
+    cosaz = math.cos(az_r)
+
+    dip_r = math.radians(-dip)
+    tandip = math.tan(dip_r)
+
+    kx = sinaz * tandip
+    ky = cosaz * tandip
+    kz = el - ea * sinaz * tandip - no * cosaz * tandip
+
+    if type == "int":
+	round = "round"
+    else:
+	round = ""
+
+    t = string.Template('$name = $type($round(x() * $kx + y() * $ky + $kz))')
+    e = t.substitute(name = name, type = type, round = round, kx = kx, ky = ky, kz = kz)
+    grass.run_command('r.mapcalc', expression = e)
+
+    grass.raster_history(name)
+
+    grass.message("Done.")
+    t = string.Template("Raster map <$name> generated by r.plane " +
+			"at point $ea E, $no N, elevation $el with dip = $dip degrees and " +
+			"aspect = $az degrees ccw from north.")
+    grass.message(t.substitute(name = name, ea = ea, no = no, el = el, dip = dip, az = az))
+
+if __name__ == "__main__":
+    options, flags = grass.parser()
+    main()