[R-sig-Geo] extract point values from contour lines

Barry Rowlingson b.rowlingson at lancaster.ac.uk
Wed Sep 4 13:06:39 CEST 2013


On Wed, Sep 4, 2013 at 11:12 AM, Barry Rowlingson
<b.rowlingson at lancaster.ac.uk> wrote:

> One approach might be to consider the points of the contour line as
> samples from the underlying data and do a smoothing such as kriging or
> inverse distance weighting, but I suspect that having a lot of samples
> with the same height at a range of distances is going to mess up the
> variogram somewhat...

    x <- 10*1:nrow(volcano)
     y <- 10*1:ncol(volcano)
     image(x,y,volcano)
     lines = contourLines(x, y, volcano)

Following on from that, the Kriging approach doesn't look so bad:

require(plyr)
# convert output of contourLines to a spatial points data frame with
height in 'h':
xyhd=ldply(l,function(e){data.frame(x=e$x,y=e$y,h=e$level)})
coordinates(xyhd)=~x+y
require(automap)
a = autoKrige(h~1,xyhd)
plot(a)

 that shows you how the points of the contours have been interpolated
over the space. If we want to just get predictions at locations we
construct a new points data set and feed that in:

 pnew = data.frame(x=runif(100,100,800),y=runif(100,50,500))
 coordinates(pnew)=~x+y
 anew = autoKrige(h~1,xyhd,pnew)
 plot(anew)

and you get point estimates and kriging standard errors - these are in
anew$krige_output as a spatial points data frame.

Now I suppose I should compare those estimates with direct samples
from the original volcano DEM to see how well we've done.... Left as
an exercise...

There seems to be a bit of literature on interpolating from contour
lines via various other methods too.

Barry



More information about the R-sig-Geo mailing list