triangular_interpolation: Function for various types of triangular ['linear']...

In GeoscienceAustralia/unstructInterp: Routines for interpolation from unstructured point data

Description

Delanay triangulation is supported, as well as a method based on linear interpolation from the 3 nearest neighbour of the interpolation point, with limiting.

Usage

triangular_interpolation(xy, vals, newPts, useNearestNeighbour = TRUE)

Arguments

xy	= point locations associated with the values to interpolate from (typically nx2 matrix)
vals	= values at xy (can be a matix with 1 or more colums and the same number of rows as xy)
newPts	= points where we want interpolated values (typically mx2 matrix)
useNearestNeighbour	= TRUE/FALSE (effect described above)

Details

If useNearestNeighbour = FALSE then it provides a wrapper around the delanay triangulation used in the 'geometry' package. Unfortunately the look-up can be slow with this method for large point clouds.
If useNearestNeighbour=TRUE, we find the 3 nearest xy neighbours of each point to interpolate to, and interpolate using the plane defined by those 3 neighbours. Limiting is used to ensure the interpolated value does not exceed the range of the xy neighbours. This method is fast since it relies only an a fast nearest neighbours implementation (via FNN)

Value

matrix/vector with as many columns as 'vals' and as many rows as 'newPts', containing the 'vals' interpolated to 'newPts'

Examples

# Make a single triangle in the plane z=x+y, and interpolate from it
   xy=matrix(c(0,0,0,1,1,1),ncol=2,byrow=TRUE)
   vals=c(0, 1, 2) # z=x+y
   newPts=matrix(c(0.5, 0.5, 0.3, 0.3), ncol=2, byrow=TRUE)

   out=triangular_interpolation(xy, vals, newPts)
   stopifnot(all.equal(out, c(1.0,0.6)))

   # Re-order triangle
   xy=xy[3:1,]
   vals=vals[3:1]
   out=triangular_interpolation(xy, vals, newPts)
   stopifnot(all.equal(out,c(1.0,0.6)))

   #another one, with formula z=0.5*x+0.2*y+7
   xy=matrix(c(-1, -1, 1, -0.5, 0.5, 1), ncol=2,byrow=2)
   vals=0.5*xy[,1]+0.2*xy[,2]+7
   newPts=matrix(c(0,0, 0.5, 0.3),ncol=2,byrow=TRUE)
   expectedVals=0.5*newPts[,1]+0.2*newPts[,2]+7
   out=triangular_interpolation(xy,vals,newPts)
   stopifnot(all.equal(out,expectedVals))

   # A point outside the triangle
   newPts=matrix(c(-1,0, -1, 1),ncol=2,byrow=TRUE)
   out=triangular_interpolation(xy,vals,newPts,useNearestNeighbour=FALSE)
   stopifnot(all(is.na(out)))
   # Failure is expected here if using approximate triangulation based on nearest neighbour methods

   # A single point
   newPts=matrix(c(0,0),ncol=2)
   out=triangular_interpolation(xy,vals,newPts)
   stopifnot(out==7)

   # Points on the triangle
   newPts=xy
   out=triangular_interpolation(xy,vals,newPts)
   stopifnot(all(out==vals))

   # Point on an edge
   newPts=matrix(0.5*(xy[1,]+xy[2,]),ncol=2)
   out=triangular_interpolation(xy,vals,newPts)
   stopifnot(all(out==0.5*(vals[1]+vals[2])))

GeoscienceAustralia/unstructInterp documentation built on May 28, 2019, 12:38 p.m.