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

[Brad Douglas]

On Tue, 2007-07-10 at 12:32 -0400, 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.

Yeah, that was a paper from PE&RS.  Recently, there was a new paper with
a "new" algorithm, but it wasn't compared to many others.

Paper I sent Helena:
Q. Zhou, et. al., 2004, Error Analysis on Grid-Based Slope and Aspect
Algorithms, PE&RS (Photogrammetric Engineering & Remote Sensing),
70(8):957-962

> 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,

I implemented several of the algorithms in r.slope.aspect, but I never
committed them to CVS.  I didn't think others would be terribly
interested because the default one used is acceptable for most everyone
and it might confuse users with so many (3-4) options.

If there is interest, I'll find some time to commit it.

Interpolation algorithms included:
3rd Order Finite Distace
3rd Order Finite Distance Weighted by Reciprocal of Distance
3rd Order Finite Distance Weighted by Reciprocal of Squared Distance

The algorithms are similar, but can improve results depending on the
situation.  The last one is the overall best performer (and I believe
the algorithm already used in r.slope.aspect).


-- 
73, de Brad KB8UYR/6 <rez touchofmadness com>