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