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

Question

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)?

score 2 · Accepted Answer · answered Feb 09 '23 at 12:55

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')

Wimpel · Answer 2 · 2023-02-09T13:08:05.063

I always return to the sf-package for this kind of operations. The sfheaders-package has convenience-functions for directly turning vector/matrices into simple feature objects.

library(sf)
library(sfheaders)
mypoint <- sfheaders::sf_point(c(8506, -4010))
sf::st_within(mypoint, st_as_sf(whole_space), sparse = FALSE)
#      [,1]  [,2]  [,3]  [,4]  [,5]  [,6]  [,7]
# [1,] TRUE FALSE FALSE FALSE FALSE FALSE FALSE

so, your (first) point is in polygon 1

If you do not need the complete matrix, drop the sparse-argument (default = TRUE)

sf::st_within(mypoint, st_as_sf(whole_space))
# Sparse geometry binary predicate list of length 1, where the predicate was `within'
#  1: 1

of you can subset your polygons which have points inside, using the logical output of the sparce matrix

st_as_sf(whole_space)[colSums(sf::st_within(mypoint, st_as_sf(whole_space), sparse = FALSE)) > 0, ]
# Simple feature collection with 1 feature and 0 fields
# Geometry type: POLYGON
# Dimension:     XY
# Bounding box:  xmin: 5850 ymin: -6050 xmax: 11250 ymax: -3100
# CRS:           NA
# geometry
# 1 POLYGON ((6225 -6050, 11250...

etc.. etc..

Nice, +1. I thought an sf approach might be a bit simpler, though stuck with terra as this seemed to be the framework within which the OP was working. — Allan Cameron, Feb 09 '23 at 12:58
I do not have much hands-on experience with terra... perhaps I should use it more :) — Wimpel, Feb 09 '23 at 13:00

Robert Hijmans · Answer 3 · 2023-02-09T22:48:53.180

You can use indexing like this

v <- vect(cbind(8506, -4010))
whole_space[v, ]
# class       : SpatVector 
# geometry    : polygons 
# dimensions  : 1, 0  (geometries, attributes)
# extent      : 5850, 11250, -6050, -3100  (xmin, xmax, ymin, ymax)
# coord. ref. :

You can also do

v <- vect(cbind(8506, -4010))
i <- is.related(whole_space, v, "intersects")
whole_space[i,]

Illustration

plot(whole_space)
text(whole_space, cex=.8)
points(v, col="red", cex=1.2, pch="x")
lines(whole_space[v], col="blue", lwd=2)

This also works for multiple points

v <- vect(cbind(c(8506,-3729.8,0), c(-4010, -1525.117, 2000)))
w <- whole_space[v]

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

3 Answers3