-
-
Notifications
You must be signed in to change notification settings - Fork 291
/
r.plane.py
executable file
·128 lines (113 loc) · 3.38 KB
/
r.plane.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
#!/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-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 math
import string
import grass.script as gscript
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 = 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()