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

[Glynn Clements]

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).

-- 
Glynn Clements <glynn at gclements.plus.com>