0

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

2

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!

  • A quick note for posterity: In my own data, I first had to assign the right projection to the created polygon shapefile, before I could overlay w/ the point data. With that tweak, it works perfectly! – Leah Bevis Sep 5 '18 at 19:43

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.