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

[Dylan Beaudette]

On Tuesday 10 July 2007 11:11, Brad Douglas wrote:
> 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).

I would be interested in these, and have a rather high-density slope and 
aspect field-data set which could be used for some comparisons.

Dylan