Skip to content

extract point values from contour lines

2 messages · Barry Rowlingson

Original
Barry Rowlingson
# 
On Wed, Sep 4, 2013 at 10:09 AM, pierluigi de rosa
<pierluigi.derosa at gmail.com> wrote:
If you can lower the bar to helping people answer your question, you
will get more answers!

For example, giving some test data and code to generate it, and some
idea of what you expect to get out. We don't care too much about your
specifics (ie the 1:10000 and 1:25000 contours).

 What format is your contour data? Is it the output from the
'contourLines' function? Or is it a SpatialLinesDataFrame with a
height attribute? Or something else?

A good question would be something like:

"""
I have contour data in the same form as from contourLines, for example:

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

and I want to be able to get the height at any point in the area from
the contour line data (yes, I could just sample the underlying grid
but I don't have that data, I only have the contours).

"""

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

If all your contour lines are complete rings then you can do
point-in-polygon tests to see which rings the point is inside and
hence know its height is between that and the next ring, and then
maybe interpolate linearly by distance to ring, but there's some edge
cases when a point is inside a summit ring when you only know it is
bigger than X and less than X+10 (if 10 is your contour interval).

Searching for 'interpolation from contour lines' might be useful.

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

            

      
      
      
    

        
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