Match SpatialPoints to Delaunay Triangles in R

Question

I have 2 SpatialPointsDataFrame files. I compute the Delaunay Triangles from the first shapefile (shp1) as in the code below. Then I wish to assign each point from my second shapefile (shp2) to one of the Delaunay Triangles held within the object vtess.

vtess <- deldir(shp1$coords.x1,shp1$coords.x2) 

I can assign each shp2 point to a tessellation tile in the following way:

tl <- tile.list(vtess)
shp2$tile <- 99
for (i in 1:length(shp2$tile)){
   shp2$tile[i] <- which.tile(shp2$coords.x1[i], shp2$coords.x2[i], tl)
}

However, while I can find the list of triangles using triang.list(vtess), there is no which.triang command, so it is not possible to similarly assign shp2 points to triangles.

I'm guessing I therefore need to (1) create a SpatialPolygon of the triangles, and then (2) overlap the shp2 points to the new polygon file. Is this the most efficient way? And if so, how do I create a polygon shapefile from the triangles held within vtess?

Cross-posted, with no reply, here: https://stackoverflow.com/questions/51700872/match-spatialpoints-to-delaunay-triangles

Answer 1

Here's how - use triang.list to get the points that make the triangles, then wrap it all up in some code that constructs SpatialPolygons from the sp package. Then do an overlay test.

Sample data:

library(sp)
library(deldir)
pts1 = data.frame(ID=LETTERS, x=runif(26),y=runif(26))
pts2 = data.frame(ID=LETTERS, x=runif(26),y=runif(26))
coordinates(pts1)=~x+y
coordinates(pts2)=~x+y

make the tesselation:

vtess=deldir(pts1@coords[,1],pts1@coords[,2])

build triangles:

tl = triang.list(vtess)

then this monster builds the triangles into spatial polygons:

polys = SpatialPolygons(
  lapply(1:length(tl),
         function(i){
           Polygons(
             list(
               Polygon(tl[[i]][c(1:3,1),c("x","y")])
                 ),ID=i)
           }
          )
        )

Then you can overlay your test points on those polygons:

> over(pts2,polys)
 1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 
 6 34 NA 32  8 NA 26 NA NA  9 30 NA 36 NA NA 10 NA 34 39 22 34 21 34 10 NA NA 

which gives the index number of the triangle. So the first point is over triangle 6, which is:

> tl[[6]]
  ptNum        x        y
1     2 0.404418 0.625027
2    25 0.678718 0.612262
3     9 0.487887 0.742302

made up of points 2, 25, and 9. You can verify this by plotting those three points, then the rest of them, then the first test point:

> plot(pts1[c(2,25,9),],pch=19)
> plot(pts1,add=TRUE,pch=3)
> plot(pts2[1,],add=TRUE,pch=19,cex=2,col="blue")

and you should see a blue dot in a triangle of black dots - with no plus signs in the triangle!