|
@@ -120,7 +120,7 @@ def main():
|
|
|
units = options['units']
|
|
|
|
|
|
if not grass.find_file(input)['file']:
|
|
|
- grass.fatal("Map <%s> not found." % elev)
|
|
|
+ grass.fatal("Map <%s> not found." % input)
|
|
|
|
|
|
if input == output:
|
|
|
grass.fatal("Input elevation map and output relief map must have different names")
|
|
@@ -152,18 +152,19 @@ def main():
|
|
|
|
|
|
grass.message("Calculating shading, please stand by.")
|
|
|
|
|
|
- t = string.Template(r'''$elev_out = eval( \
|
|
|
-x=($zmult*$elev[-1,-1] + 2*$zmult*$elev[0,-1] + $zmult*$elev[1,-1] - \
|
|
|
- $zmult*$elev[-1,1] - 2*$zmult*$elev[0,1] - $zmult*$elev[1,1])/(8.*ewres()*$scale), \
|
|
|
-y=($zmult*$elev[-1,-1] + 2*$zmult*$elev[-1,0] + $zmult*$elev[-1,1] - \
|
|
|
- $zmult*$elev[1,-1] - 2*$zmult*$elev[1,0] - $zmult*$elev[1,1])/(8.*nsres()*$scale), \
|
|
|
-slope=90.-atan(sqrt(x*x + y*y)), \
|
|
|
-a=round(atan(x,y)), \
|
|
|
-a=if(isnull(a),1,a), \
|
|
|
-aspect=if(x!=0||y!=0,if(a,a,360.)), \
|
|
|
-cang = sin($alt)*sin(slope) + cos($alt)*cos(slope) * cos($az-aspect), \
|
|
|
-if(cang < 0.,0.,100.*cang), \
|
|
|
-if(isnull(cang), null(), 100.*cang))''')
|
|
|
+ t = string.Template(
|
|
|
+ r'''eval( \
|
|
|
+ x=($zmult*$elev[-1,-1] + 2*$zmult*$elev[0,-1] + $zmult*$elev[1,-1] - \
|
|
|
+ $zmult*$elev[-1,1] - 2*$zmult*$elev[0,1] - $zmult*$elev[1,1])/(8.*ewres()*$scale), \
|
|
|
+ y=($zmult*$elev[-1,-1] + 2*$zmult*$elev[-1,0] + $zmult*$elev[-1,1] - \
|
|
|
+ $zmult*$elev[1,-1] - 2*$zmult*$elev[1,0] - $zmult*$elev[1,1])/(8.*nsres()*$scale), \
|
|
|
+ slope=90.-atan(sqrt(x*x + y*y)), \
|
|
|
+ a=round(atan(x,y)), \
|
|
|
+ a=if(isnull(a),1,a), \
|
|
|
+ aspect=if(x!=0||y!=0,if(a,a,360.)), \
|
|
|
+ cang = sin($alt)*sin(slope) + cos($alt)*cos(slope) * cos($az-aspect), \
|
|
|
+ )
|
|
|
+ $elev_out = if(isnull(cang), null(), 100.*cang)''')
|
|
|
expr = t.substitute(alt = alt, az = az, elev = elev, elev_out = elev_out, scale = scale, zmult = zmult)
|
|
|
p = grass.feed_command('r.mapcalc')
|
|
|
p.stdin.write(expr)
|