# [R-sig-Geo] Principal Directions with surf.ls

Gabor Grothendieck ggrothendieck at gmail.com
Wed Mar 14 17:30:03 CET 2007

```Thanks.  In my case the region of interest is no larger
than a city block so I fall into your simpler case.

On 3/14/07, Roger Bivand <Roger.Bivand at nhh.no> wrote:
> 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
> > > 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