[R-sig-Geo] Principal Directions with surf.ls
Roger Bivand
Roger.Bivand at nhh.no
Wed Mar 14 15:47:38 CET 2007
On Wed, 14 Mar 2007, Gabor Grothendieck wrote:
> Thanks. I was hoping that either there were an existing
> function or the components of the output of surf.ls were
> sufficiently documented to allow one to do it oneself.
>
> surf.ls seems to do all the real work in a C function
> and I don't relish having to go through that but
> since its just a regression I guess I can ditch surf.ls
> and use lm in which case I can get the coefficients
> explicitly and figure it out from there.
There is one extra step that is necessary. It has been known since the
1970's (D. Unwin, also Brian Ripley's 1981 book) that if the northings are
in metres north of the equator, numerical problems cut in very fast when
using trend surface, because the powers become very large numbers.
Typically, for a small site, there will be no difference between y^2 and
y^3, say, from the point of view of the computer.
So the coordinates need to be scaled so that they do not explode, and both
the functions in spatial and in gstat for trend surface do this, in
surf.ls() by transforming both x and y coordinates to lie between +/-1.
Other than that, lm() ought to be fine. For blocks, it may be sufficient
just to go from m to km as the unit.
Roger
>
> On 3/14/07, Roger Bivand <Roger.Bivand at nhh.no> wrote:
> > On Wed, 14 Mar 2007, Gabor Grothendieck wrote:
> >
> > > Suppose we plot a quadratic surface as in:
> > >
> > > library(spatial)
> > > example(surf.ls)
> > >
> > > How can we plot lines on the graph to show the
> > > principal directions/axes of the elliptical contours
> > > and how do we get them?
> >
> > In the example there are no closed ellipses, with an open, north-facing
> > depression with flow from the west, south-west, south, south-east and
> > east towards the north at the centre of the northern edge.
> >
> > The contour lines can be retrieved with contourLines() on the same object
> > as contour() is used on. However, they are just a generalised
> > visualisation of the surface represented by the matrix. The underlying
> > question might then be: what is the aspect and slope of the surface, in
> > this case a smoothed surface, what are the plan and especially profile
> > curvatures, and what are their length attributes (ie. patches of similar
> > aspect along a profile curvature.
> >
> > In GIS terms this would involve geomorphometrics, often computed using a 3
> > x 3 filter, and then flow analysis to find the directions (on a digital
> > elevation model, these would be similar to predicting drainage channel
> > networks. Landslip and avalanche studies often need slope length too,
> > which involves measuring patch characteristics in aspect and profile
> > curvature - that is extracting patches with interesting pattern
> > signatures, and checking whether the patch is "straight" in some sense,
> > and longer than a minimum accumulation threshold to feed mass movement.
> >
> > These are general observations in the case that we don't have a smoothed
> > surface. But here we have a model, and I'd guess that the underlying model
> > could be used, rather than its predictions.
> >
> > I couldn't see a package with a vector field plot, there was a suggestion
> > in:
> >
> > http://finzi.psych.upenn.edu/R/Rhelp02a/archive/35933.html
> >
> > but you'd still need the slope and aspect to get the length and direction
> > of the arrows (for display only). A classic reference is:
> >
> > Horn, B. K. P. (1981). Hill Shading and the Reflectance Map, Proceedings
> > of the IEEE, 69(1): 14-47
> >
> > and Horn's algorithm has stood the test of time, as has Zevenbergen and
> > Thorne:
> >
> > Zevenbergen, L. W., Thorne, C. R. 1987. Quantitative analysis of land
> > surface topography. Earth Surface Processes and Landforms, 12 (1), 4756.
> >
> > It could well be that the surface you are considering is something other
> > than elevation (temperature, etc.), but getting further would depend on
> > the specific characteristics of the study, and whether the surface is a
> > fitted trend surface or not. With a fitted trend surface model, it ought
> > to be possible to use the model. If thinking in terms of geomorphometrics
> > is unhelpful, please accept that geographers "see" things in those kinds
> > of ways!
> >
> > Hope this helps,
> >
> > Roger
> >
> > >
> > > Thanks.
> > >
> > > _______________________________________________
> > > R-sig-Geo mailing list
> > > R-sig-Geo at stat.math.ethz.ch
> > > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> > >
> >
> > --
> > Roger Bivand
> > Economic Geography Section, Department of Economics, Norwegian School of
> > Economics and Business Administration, Helleveien 30, N-5045 Bergen,
> > Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43
> > e-mail: Roger.Bivand at nhh.no
> >
> >
>
--
Roger Bivand
Economic Geography Section, Department of Economics, Norwegian School of
Economics and Business Administration, Helleveien 30, N-5045 Bergen,
Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43
e-mail: Roger.Bivand at nhh.no
More information about the R-sig-Geo
mailing list