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

[marcelino.delacruz at upm.es]

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

>
>
>On 11/20/2011 11:02 PM, Rolf Turner wrote:
>> 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.
>> 
>> 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.
>
>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).
>
>> 
>>     cheers,
>> 
>>         Rolf Turner
>> 
>> _______________________________________________
>> 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