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

[marcelino.delacruz at upm.es]

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ó:

>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:
>> 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.
>> 
>> _______________________________________________
>> R-sig-Geo mailing list
>> R-sig-Geo at r-project.org
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
>-- 
>Edzer Pebesma
>Institute for Geoinformatics (ifgi), University of Münster
>Weseler Straße 253, 48151 Münster, Germany. Phone: +49 251
>8333081, Fax: +49 251 8339763  http://ifgi.uni-muenster.de
>http://www.52north.org/geostatistics      e.pebesma at wwu.de
>
>_______________________________________________
>R-sig-Geo mailing list
>R-sig-Geo at r-project.org
>https://stat.ethz.ch/mailman/listinfo/r-sig-geo