[R-sig-Geo] DEM interpolation with Delaunay triangulation

[marcelino.delacruz at upm.es]

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)
}


>
># See how does it work:
>
>data(longleaf)
>plot(longleaf)
>
>cosa <- tin(longleaf, nx=100, ny=100)
>
>cosa2 <- tin(longleaf, x=1:100, y=1:100)


>plot(cosa)

>plot(cosa2)


>
>image(x=1:100,y=1:100, z=t(matrix(cosa$marks, 100,100, byrow=T)))
>
>
>
>Cheers, Marcelino
>
>
>>> On 21/11/11 04:30, gianni lavaredo wrote:
>>>> Dear Researchers,
>>>>
>>>> I am looking a package o method to interpolate a point data.frame (x,y,z)
>>>> with Delaunay triangulation  (Trianguled interpoaltion network) and
>>>> create
>>>> a polygons data,frame o a raster. I know there are several fuction in
>>>> library RSAGA, spatstat (with "deldir"), library geometry.