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 _______________________________________________ grass-dev mailing list [email protected] http://grass.itc.it/mailman/listinfo/grass-dev
