Skip to content

DEM interpolation with Delaunay triangulation

7 messages · gianni lavaredo, Rolf Turner, Edzer Pebesma +1 more

Original
Rolf Turner
#
On 21/11/11 04:30, gianni lavaredo wrote:
I have no idea what you're on about.  If you could make your question
clearer and more specific I might be able to make some suggestions
as to how you could use spatstat/deldir.

An explicit simple (toy) example might help.

     cheers,

         Rolf Turner
2 days later
Edzer Pebesma
#
On 11/20/2011 11:02 PM, Rolf Turner wrote:
I have wondered this myself too -- the issue is the so-called TIN,
triangular interpolation network. Say, we have (field) observations
(x,y,z), we form a delanay triangulation on all (x,y) pairs, and we want
to obtain the prediction of z for an arbitrary (x0,y0) pair, based on
the linear trend surface fit through the three (x1,y1), (x2,y2), (x3,y3)
points that form the triangle in which (x0,y0) falls.

You could do, of course, a linear trend surface using three nearest
neighbours, e.g. by

require(gstat)
loadMeuse()
x = krige(log(zinc)~x+y, meuse, meuse.grid, nmax=3)
spplot(x[1],col.regions=bpy.colors(),
	sp.layout=list("sp.points", meuse))


but that takes the three nearest observations, which don't need to be
the three points in the triangle (and often, they are not, because
meuse.grid contains points outside the convex hull of the meuse
observation locations).

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Marcelino de la Cruz
    

    

      
      #
    

  


  

      
Well, I think that this function (all with spatstat tools) could do the
job.

I included a weigthed linear interpolation but this could be easily
modified to fit other models.

It only takes as arguments a marked point pattern (xyz) and the number of
points in x and y dimensions to build the "arbitrary" x0 y0 points


tin <- function(X, nx, ny){
      require(spatstat)
      lltes<-delaunay(X)
      gri <-  gridcentres(X$window, nx=nx, ny=ny)
      gri.ppp <- ppp(gri$x,gri$y, window=X$window,
                           marks=rep(1,length(gri$x)))
        cat("\n","number of triangles =",
length(lltes[[3]]),"\n\n")
       for(i in 1:length(lltes[[3]])){
	     progressreport(i, length(lltes[[3]]))

	      #grid points within the triangulation
	     xoyo <- unmark(gri.ppp[lltes[[3]][[i]]])

	     # original points defining the triangle
	     xyz <- X[lltes[[3]][[i]]]
	       # z values of the three points
	       z<-xyz$marks

	     #distance from grid points to vertices
	      cdis <- crossdist(xoyo,xyz)

	      # weighted linear interpolation within triangle
	      grim <- apply(cdis, 1, function(f) sum(f*z)/sum(f))

	     #assign interpolated values
	       gri.ppp[lltes[[3]][[i]]]$marks <- grim
        }
   return(gri.ppp)
}


# See how does it work:

data(longleaf)
plot(longleaf)

cosa <- tin(longleaf, nx=100, ny=100)

plot(cosa)

image(x=1:100,y=1:100, z=t(matrix(cosa$marks, 100,100, byrow=T)))



Cheers, Marcelino



Con fecha 23/11/2011, "Edzer Pebesma" <edzer.pebesma at uni-muenster.de>
escribi?:

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

      
      
    

  
  


  


      

    

    
1 day later

    

      
      

        


  

        

      Marcelino de la Cruz
    

    

      
      #
    

  


  

      
Hi all,

Thanks to Rolf's suggestion I have incorporated the arguments x and y to
let the user select the precise set of points where to interpolate.

I have also corrected the formula (now is a  inverse distance weighted
interpolation) and have incorporated also the "power parameter" p as
an argument.

Cheers,

MArcelino






tin <- function(X, nx, ny,x=NULL, y=NULL,p=1){
     require(spatstat)
     lltes<-delaunay(X)

       if(is.null(x)){
            gri <-  gridcentres(X$window, nx=nx, ny=ny)
            gri.ppp <- ppp(gri$x,gri$y, window=X$window,
                                 marks=rep(0,length(gri$x)))
	}
	if(!is.null(x)){
	    gri.ppp<- ppp(x=x, y=y, window=X$window,
	                         marks=rep(0, length(x)))
        }

       cat("\n","number of triangles =",
length(lltes[[3]]),"\n\n")
      for(i in 1:length(lltes[[3]])){
            progressreport(i, length(lltes[[3]]))

             #grid points within the triangulation
            xoyo <- unmark(gri.ppp[lltes[[3]][[i]]])

            # original points defining the triangle
            xyz <- X[lltes[[3]][[i]]]
              # z values of the three points
              z<-xyz$marks

            #distance from grid points to vertices
             cdis <- crossdist(xoyo,xyz)

             # inverse distance weighted interpolation within triangle
             grim <- apply(cdis, 1, function(f) sum((1/f^p)*z)/sum(1/f^p))

            #assign interpolated values
              gri.ppp[lltes[[3]][[i]]]$marks <- grim
       }
  return(gri.ppp)
}

            

      
      
      
    

  
  


  


      

    

    

      
      

        


  

        

      Edzer Pebesma
    

    

      
      #
    

  


  

      
In that case, you might want to reconsider the name of this function, as
TIN is something different from inverse distance weighted; if you are on
an edge of the triangle, you need to get the straight line interpolation
from the two z values on that edge, the third needs to get zero weight...

        
On 11/24/2011 02:27 PM, marcelino.delacruz at upm.es wrote:

            

      
      
      
    

        

      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Marcelino de la Cruz
    

    

      
      #
    

  


  

      
Thank you for your suggestions Edzer.

I have rewritten the function (that now fits a linear trend on x and y),
and have reconsidered a more appropriate name.
I include a more ilustrative example with the Meuse data.

Cheers,

Marcelino


maybetin <- function(X, nx, ny, x=NULL, y=NULL, na.v=0){
     require(spatstat)
     lltes<-delaunay(X)

       if(is.null(x)){
            gri <-  gridcentres(X$window, nx=nx, ny=ny)
            gri.ppp <- ppp(gri$x,gri$y, window=X$window,
                                 marks=rep(na.v,length(gri$x)))
	}
	if(!is.null(x)){
	    gri.ppp<- ppp(x=x, y=y, window=X$window,
	                         marks=rep(na.v, length(x)))
        }

       cat("\n","number of triangles =",
length(lltes[[3]]),"\n\n")
      for(i in 1:length(lltes[[3]])){
            progressreport(i, length(lltes[[3]]))

             #grid points within the triangulation
            xoyo <- unmark(gri.ppp[lltes[[3]][[i]]])

            # original points defining the triangle
            xyz <- X[lltes[[3]][[i]]]
              # z values of the three points
              z<-xyz$marks
              mtrend <-with(xyz, lm(marks~x+y))

              grim <- predict(mtrend,
	             newdata=data.frame(x = xoyo$x, y=xoyo$y))

            #assign interpolated values
              gri.ppp[lltes[[3]][[i]]]$marks <- grim
       }
  return(gri.ppp)
}



# A very brute force TIN (maybe) of the log(zinc) data of meuse

require(gstat)
loadMeuse()
require(spatstat)
spatstat.options(gpclib=TRUE)

#data pre-processing
  library(maptools) # for the "as" method
  meuse.ppp <- as(meuse,"ppp")
  meusegrid.ppp <- as(meuse.grid,"ppp")

  #readjust the window for not missing points of meusegrid
   meuse.ppp$window <- meusegrid.ppp$window

  # use only log(zinc) as mark
      meuse.ppp$marks <-log( meuse at data$zinc)

#compute (maybe) TIN
  xgrid <-maybetin(meuse.ppp,x=meusegrid.ppp$x,y=meusegrid.ppp$y)

#A very brute force kind of representation
   meuse.grid at data$tin<-xgrid$marks

   spplot(meuse.grid["tin"],col.regions=bpy.colors(),
          sp.layout=list("sp.points", meuse))

   meuse.grid at data$tine<-exp(xgrid$marks)
   spplot(meuse.grid["tine"],col.regions=bpy.colors(),
      sp.layout=list("sp.points", meuse))









Con fecha 24/11/2011, "Edzer Pebesma" <edzer.pebesma at uni-muenster.de>
escribi?: