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>