rspatialterrapoint-in-polygon

Find which spatial polygon (object in SpatVector vect) covers a data point


I have created a SpatVector of several polygons using the terra package as

library(terra)
er <- rbind(c(6225, -6050), c(11250, -5400), c(10800, -3100), c(7560, -3545), c(5850, -3800))
re <- rbind(c(5850, -3800), c(7560, -3545), c(6900, -100), c(5250, -350))
wr <- rbind(c(10750, -3100), c(15600, -2450), c(14650, 2850), c(12650, 2575), c(11600, 2430), c(4950, 1500), c(5250, -350), c(5850, -3800), c(7560, -3545))
pa <- rbind(c(4950, 1500), c(5250, -350), c(-2500, -1500), c(-2100, -4100), c(-4700, -4400), c(-4930, -2800), c(-5300, 0))
tb <- rbind(c(-4700, -4400), c(-4425, -6350), c(-1775, -6100), c(-2100, -4100))
exl <- rbind(c(-4700, -4400), c(-4425, -6350), c(-6300, -6525), c(-6750, -3000), c(-4930, -2800))
exu <- rbind(c(11600, 2430), c(12650, 2575), c(12400, 4050), c(11350, 3925))
whole_space <- rbind(cbind(object = 1, part = 1, er, hole = 0),
            cbind(object = 2, part = 2, wr, hole = 0),
            cbind(object = 3, part = 3, pa, hole = 0),
            cbind(object = 4, part = 4, tb, hole = 0),
            cbind(object = 5, part = 5, re, hole = 0),
            cbind(object = 6, part = 6, exl, hole = 0),
            cbind(object = 7, part = 7, exu, hole = 0))
whole_space <- vect(whole_space, "polygons")

I now have a data point, say c(8506, -4010). How can I find the polygon/object which covers the data point (if any)?


Solution

  • If you have the point co-ordinates as a vector, you can turn it into a SpatVector first:

    point <- c(8506, -4010)
    
    vpoint <- vect(t(point), 'points')
    

    Then you can iterate through your polygons to see whether the point still exists after being cropped by each polygon. If the point still exists, it lies within said polygon:

    which(sapply(seq_along(whole_space), function(i) {
      length(crop(whole_space[i], vpoint)) > 0 }
    ))
    #> [1] 1
    

    We can confirm this is correct by drawing the point and the first member of whole_space:

    plot(whole_space[1])
    plot(vpoint, add = TRUE, col = 'red')
    

    enter image description here