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

[Dylan Beaudette]

On Tuesday 10 July 2007 09:32, Helena Mitasova wrote:
> 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
>
>

While we are on this topic, it might be interesting to allow power users to 
specify the algorithm. I think that the paper that you are mentioning might 
be in my filing cabinet... just cannot quite recall the authors! I will check 
when back in town.

cheers,

dylan