Geografic-Information-System/scripts/r.plane/r.plane.py

#!/usr/bin/env python
# -*- coding: utf-8 -*-
############################################################################
#
# MODULE:	r.plane for GRASS 5.7; based on r.plane for GRASS 5
# AUTHOR(S):	1994, Stefan Jäger, University of Heidelberg/USGS
#               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-2012 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.
#% keyword: raster
#% keyword: elevation
#%end
#%option G_OPT_R_OUTPUT
#%end
#%option
#% key: dip
#% type: double
#% gisprompt: -90-90
#% answer: 0.0
#% description: Dip of plane in degrees
#% required : yes
#%end
#%option
#% key: azimuth
#% type: double
#% gisprompt: 0-360
#% answer: 0.0
#% description: Azimuth of the plane in 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: CELL,FCELL,DCELL
#% description: Type of raster map to be created
#% answer: FCELL
#%end

import sys
import os
import math
import string
import grass.script as grass

def main():
    name = options['output']
    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 == "CELL":
	round = "round"
	dtype = "int"
    elif type == "FCELL":
	round = ""
	dtype = "float"
    else:
	round = ""
	dtype = "double"

    grass.mapcalc("$name = $type($round(x() * $kx + y() * $ky + $kz))",
		  name = name, type = dtype, round = round, kx = kx, ky = ky, kz = kz)

    grass.run_command('r.support', map = name, history = '')
    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()