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

#!/usr/bin/env python3

############################################################################
#
# 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 G_OPT_R_TYPE
# % answer: FCELL
# % required: no
# %end

import math
import string
import grass.script as gscript


def main():
    name = options["output"]
    type = options["type"]
    dip = float(options["dip"])
    az = float(options["azimuth"])
    try:
        ea = float(options["easting"])
        no = float(options["northing"])
    except ValueError:
        try:
            ea = float(gscript.utils.float_or_dms(options["easting"]))
            no = float(gscript.utils.float_or_dms(options["northing"]))
        except:
            gscript.fatal(_("Input coordinates seems to be invalid"))
    el = float(options["elevation"])

    # reg = gscript.region()

    # Test input values
    if abs(dip) >= 90:
        gscript.fatal(_("dip must be between -90 and 90."))

    if az < 0 or az >= 360:
        gscript.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"

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

    gscript.run_command("r.support", map=name, history="")
    gscript.raster_history(name)

    gscript.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."
    )
    gscript.message(t.substitute(name=name, ea=ea, no=no, el=el, dip=dip, az=az))


if __name__ == "__main__":
    options, flags = gscript.parser()
    main()