[GRASS-user] MDOW relief shading in GRASS

[Micha Silver]

On 09/30/2010 05:55 AM, maning sambale wrote:
> Hi,
>
> I am revisiting a previous script I made in producing
> "Multidirectional, oblique-weighted, shaded-relief".  ArcGis has this
> script:
> http://blogs.esri.com/Support/blogs/mappingcenter/archive/2008/10/07/updated-hillshade-toolbox.aspx
>
> The commads in ARC 6.0.1 GRID by Rober Mark (USGS) are as follows:
>
> shade225 = hillshade(hawaii, 225, 30, shade, 5)
> shade270 = hillshade(hawaii, 270, 30, shade, 5)
> shade315 = hillshade(hawaii, 315, 30, shade, 5)
> shade360 = hillshade(hawaii, 360, 30, shade, 5)
> h00 = resample(hawaii, 1000)
> h01 = focalmean(h00)
> h02 = focalmean(h01)
> h03 = focalmean(h02)
> asp = aspect(h03)
> asp1 = con(isnull(asp), 293, asp)
> w225 = sqr(sin((asp1 - 225) div deg ))
> w270 = sqr(sin((asp1 - 270) div deg))
> w315 = sqr(sin((asp1 - 315) div deg ))
> w360 = sqr(sin(asp1 div deg ))
> setcell minof
> temp = w225 * shade225 + w270 * shade270 + w315 *
> shade315 + w360 * shade360
> shade4 = int(temp div 2)
>
> I tried porting in GRASS shell:
>
> g.region rast=$GIS_OPT_input -p
>
> # Compute hillshade at azimuth 225, 270, 315 and 360 at 30 degrees sun
> illumination angle
> r.shaded.relief map=$GIS_OPT_input shadedmap=shade225 altitude=30
> azimuth=225 zmult=1 scale=1
> r.shaded.relief map=$GIS_OPT_input shadedmap=shade270 altitude=30
> azimuth=270 zmult=1 scale=1
> r.shaded.relief map=$GIS_OPT_input shadedmap=shade315 altitude=30
> azimuth=315 zmult=1 scale=1
> r.shaded.relief map=$GIS_OPT_input shadedmap=shade360 altitude=30
> azimuth=360 zmult=1 scale=1
>
> # Resample DEM map
> r.neighbors in=$GIS_OPT_input out=res1 method=average size="$GIS_OPT_window"
>
> # compute aspect map
> r.slope.aspect elevation=res1 format=degrees prec=float aspect=aspect
> zfactor=1.0 min_slp_allowed=0.0
>
>    

Hi Maning:
Markus M reminded me recently that the GRASS convention for aspect is 
cartesian direction. i.e 90° is north, 180° is west. How is this handled 
in Arc? Maybe that's the difference?


> # compute weights 225, 270, 315 and 360
> r.mapcalc "w225 = sin((aspect - 225)*sin(aspect - 225))"
> r.mapcalc "w270 = sin((aspect - 270)*sin(aspect - 270))"
> r.mapcalc "w315 = sin((aspect - 315)*sin(aspect - 315))"
> r.mapcalc "w360 = sin((aspect)*sin(aspect))"
>
>
> #compute weighted
> r.mapcalc "weightedshade = (((w225 * shade225) + (w270 * shade270) +
> (w315 * shade315) + (w360 * shade360))/2)"
> r.colors map=weightedshade color=grey
>
> But I get not so good results, here's the output using nc_sp_08's
> elev_ned_30m (X0=MDOW, X1= default r.shaded.relief module)
> http://farm5.static.flickr.com/4148/5037622535_113924312b_b.jpg
>
> Anything wrong with my commands?
>
>    


-- 
Micha Silver
Arava Development Co. +972-52-3665918
http://www.surfaces.co.il