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

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