[GRASS-dev] using slope/aspect functions inside our module

[Helena Mitasova]

On Tue, 2007-07-10 at 17:03 +0100, Glynn Clements wrote:
> Dylan Beaudette wrote:
> 
> > > > Interesting.
> > > > Why r.slope.aspect uses
> > > >
> > > > dx = ... / (4 * ewres);
> > > >
> > > > while r.shaded.relief uses
> > > >
> > > > dx = ... / (8 * ewres);
> > > >
> > > > ??
> > >
> > > Oops. It should be 8; the ewres/nsres values are actually the
> > > distances between the centres of the top/bottom rows and left/right
> > > columns in the 3x3 window, which are two rows/columns apart. IOW, the
> > > 4*ewres should have been 4*(2*ewres) = 8*ewres, ditto for nsres.
> > 
> > Yikes, does this mean that previous calculations are now wrong?
> 
> No. I just made a mistake in "paraphrasing" the code.
> 
> The actual code for r.slope.aspect has:
> 
> 	V = G_distance(east, north, east, south) * 4 / zfactor;
> 	H = G_distance(east, ns_med, west, ns_med) * 4 / zfactor;
> 	...
> 	dx = ((*c1 + *c4 + *c4 + *c7) - (*c3 + *c6 + *c6 + *c9)) / H;
> 	dy = ((*c7 + *c8 + *c8 + *c9) - (*c1 + *c2 + *c2 + *c3)) / V;
> 
> I paraphrased the result of G_distance nsres/ewres, but it's actually
> twice that.
> 
> The factor of 4 corresponds to the weights of 1+2+1 = 4, while the
> factor of 2 corresponds to the distance in cells.
> 
> It might be clearer to write the calculation as:
> 
> 	dx = ((c1/4 + c4/2 + c7/4) - (c3/4 + c6/2 + c9/4)) / (2 * ewres);
> 	dy = ((c7/4 + c8/2 + c9/4) - (c1/4 + c2/2 + c3/4)) / (2 * nsres);
> 
> r.slope.aspect premultiplies the 4 into the divisors to eliminate the
> divisions from the weighting.
> 
> I'm not sure of the theoretical basis behind the weights used. A
> simpler approach would be to ignore the corners and use:
> 
> 	dx = (c4 - c6) / (2 * ewres);
> 	dy = (c8 - c2) / (2 * nsres);
> 
> This still produces the correct result for a sampled planar surface
> (it's also the approach used by r.bilinear and r.resamp.interp).
> 
> Ultimately, computing the slope of a DEM requires some form of
> interpolation, and the choice of algorithm is largely arbitrary (if
> you don't interpolate, the surface is a grid of horizontal rectangles,
> so the slope is always either zero or infinite).

just a note for record - the choice of algorithm does make a difference
(I think it was Brad who sent me an excellent paper on it some time
ago),the equation for the polynomial that is used for approximation here
is in the GRAS book appendix and also in the upcoming Geomorphometry
book (it has a chapter about Geomorphometry in GRASS). As long as the
3x3 neighborhood and a bivariate 2nd order polynomial is enough, this is
a very simple and accurate algorithm (introduced by Horn in 1981 for a
pioneering work for computing illuminated terrain maps). Some simpler
algorithms can produce pretty bad results.

But the way how it is written in r.slope.aspect is quite confusing - I
think it is mostly due to the fact that it needs to call G_distance
rather than using plain ewres to get the latlong converted to actual
distance. We already had a discussion on that topic too,

Helena



 
>