Hello,
In a given chain processing module (DN->high-level product) in development,
there is a need of calculating slope and aspect, i would like to use the
existing code as a library or something to the similar effect.
Anyone could explain if it is possible and how it could be possible?
thank you
Yann
On Jul 3, 2007, at 2:29 AM, Yann wrote:
Hello,
In a given chain processing module (DN->high-level product) in development,
there is a need of calculating slope and aspect, i would like to use the
existing code as a library or something to the similar effect.Anyone could explain if it is possible and how it could be possible?
it would be useful to have such library functions but you would have write them
based on r.slope.aspect. You can also look at
aspect.c in raster/r.flow
it should be easy to add slope to it or modify it as a separate function for slope.
It should use the same algorithm as r.slope.aspect but double check that.
If you need to handle lat/long, look at r.slope.aspect how to do it.
Helena
thank you
Yann_______________________________________________
grass-dev mailing list
grass-dev@grass.itc.it
http://grass.itc.it/mailman/listinfo/grass-dev
Yann wrote on 07/03/2007 08:29 AM:
Hello,
In a given chain processing module (DN->high-level product) in development,
there is a need of calculating slope and aspect, i would like to use the
existing code as a library or something to the similar effect.Anyone could explain if it is possible and how it could be possible?
A fancy solution would be to wrap r.slope.aspect in SWIG. See
swig/perl/R_slope_aspect/
But a good solution would not make a local copy of the code.
FYI: I have just updated the PERL part of the SWIG interface to
the current 6.3 state.
Markus
------------------
ITC -> dall'1 marzo 2007 Fondazione Bruno Kessler
ITC -> since 1 March 2007 Fondazione Bruno Kessler
------------------
Yann wrote:
In a given chain processing module (DN->high-level product) in development,
there is a need of calculating slope and aspect, i would like to use the
existing code as a library or something to the similar effect.Anyone could explain if it is possible and how it could be possible?
It isn't practical to use code from a module in a library.
Even if it was, calculating slope/aspect is so trivial that it's
easier to just write it from scratch than to extract the code from
r.slope.aspect.
FWIW, the calculation which r.slope.aspect uses is essentially:
dx = ((c1 + 2*c4 + c7) - (c3 + 2*c6 + c9)) / (4 * ewres);
dy = ((c7 + 2*c8 + c9) - (c1 + 2*c2 + c3)) / (4 * nsres);
key = dx*dx + dy*dy;
slp_in_perc = 100*sqrt(key);
slp_in_deg = atan(sqrt(key)) * radians_to_degrees;
aspect = (atan2(dy,dx)/degrees_to_radians);
where the cells are labelled:
c1 c2 c3
c4 c5 c6
c7 c8 c9
--
Glynn Clements <glynn@gclements.plus.com>
Interesting.
Why r.slope.aspect uses
dx = ... / (4 * ewres);
while r.shaded.relief uses
dx = ... / (8 * ewres);
??
BTW there is a typo in r.shaded.relief in line 137 ("Nap"):
g.message -e "Nap <$elev> not found."
Carlos
On 7/4/07, Glynn Clements <glynn@gclements.plus.com> wrote:
Yann wrote:
> In a given chain processing module (DN->high-level product) in development,
> there is a need of calculating slope and aspect, i would like to use the
> existing code as a library or something to the similar effect.
>
> Anyone could explain if it is possible and how it could be possible?It isn't practical to use code from a module in a library.
Even if it was, calculating slope/aspect is so trivial that it's
easier to just write it from scratch than to extract the code from
r.slope.aspect.FWIW, the calculation which r.slope.aspect uses is essentially:
dx = ((c1 + 2*c4 + c7) - (c3 + 2*c6 + c9)) / (4 * ewres);
dy = ((c7 + 2*c8 + c9) - (c1 + 2*c2 + c3)) / (4 * nsres);
key = dx*dx + dy*dy;
slp_in_perc = 100*sqrt(key);
slp_in_deg = atan(sqrt(key)) * radians_to_degrees;
aspect = (atan2(dy,dx)/degrees_to_radians);where the cells are labelled:
c1 c2 c3
c4 c5 c6
c7 c8 c9--
Glynn Clements <glynn@gclements.plus.com>_______________________________________________
grass-dev mailing list
grass-dev@grass.itc.it
http://grass.itc.it/mailman/listinfo/grass-dev
--
+-----------------------------------------------------------+
Carlos Henrique Grohmann - Guano
Visiting Researcher at Kingston University London - UK
Geologist M.Sc - Doctorate Student at IGc-USP - Brazil
Linux User #89721 - carlos dot grohmann at gmail dot com
+-----------------------------------------------------------+
_________________
"Good morning, doctors. I have taken the liberty of removing Windows
95 from my hard drive."
--The winning entry in a "What were HAL's first words" contest judged
by 2001: A SPACE ODYSSEY creator Arthur C. Clarke
"Carlos \"Guâno\" Grohmann" wrote:
Interesting.
Why r.slope.aspect usesdx = ... / (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.
--
Glynn Clements <glynn@gclements.plus.com>
Carlos wrote:
BTW there is a typo in r.shaded.relief in line 137 ("Nap"):
g.message -e "Nap <$elev> not found."
thanks, fixed.
Hamish
On Wednesday 04 July 2007 08:41, Glynn Clements wrote:
"Carlos \"Guâno\" Grohmann" 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?
dylan
On Jul 9, 2007, at 3:47 PM, Dylan Beaudette wrote:
On Wednesday 04 July 2007 08:41, Glynn Clements wrote:
"Carlos \"Guâno\" Grohmann" wrote:
Interesting.
Why r.slope.aspect usesdx = ... / (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 - it is correctly implemented - you would have found that it is wrong by now.
It is just that the variable that you divide with in the r.slope.aspect
uses 2*ewres (twice the resolution, not just one as is usual in
the published equations) - I got caught on that many times when checking
the r.slope.aspect code.
Helena
dylan
_______________________________________________
grass-dev mailing list
grass-dev@grass.itc.it
http://grass.itc.it/mailman/listinfo/grass-dev
On Monday 09 July 2007 13:04, Helena Mitasova wrote:
On Jul 9, 2007, at 3:47 PM, Dylan Beaudette wrote:
> On Wednesday 04 July 2007 08:41, Glynn Clements wrote:
>> "Carlos \"Guâno\" Grohmann" 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 - it is correctly implemented - you would have found that it is
wrong by now.
It is just that the variable that you divide with in the r.slope.aspect
uses 2*ewres (twice the resolution, not just one as is usual in
the published equations) - I got caught on that many times when checking
the r.slope.aspect code.Helena
Thanks for the quick reply! That is a relief (! no pun intended there).
Cheers,
Dylan
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@gclements.plus.com>
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
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
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>
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 DistanceThe 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