point in polygon using terra package in R

Question

I want to do a point in polygon analysis using the terra package. I have set of points and I want to extract which polygons do they fall in. The below sample data shows it for a single polygon

library(terra)
crdref <- "+proj=longlat +datum=WGS84"
  
longitude <- c(-116.7, -120.4, -116.7, -113.5, -115.5, -120.8, -119.5, -113.7, -113.7, -110.7)
latitude <- c(45.3, 42.6, 38.9, 42.1, 35.7, 38.9, 36.2, 39, 41.6, 36.9)
lonlat <- cbind(longitude, latitude)
pts <- vect(lonlat, crs=crdref)
  
lon <- c(-116.8, -114.2, -112.9, -111.9, -114.2, -115.4, -117.7)
lat <- c(41.3, 42.9, 42.4, 39.8, 37.6, 38.3, 37.6)
lonlat <- cbind(id=1, part=1, lon, lat)
pols <- vect(lonlat, type="polygons", crs=crdref)
  
plot(pols, border='blue', col='yellow', lwd=3)
points(pts, col='red', pch=20, cex=3)

terra::extract(pts, pols)
       id.y id.x
[1,]    1   NA

But the output is not clear to me. I simply need which polygon does each point falls into.

Answer 1

Example data

library(terra)
longitude <- c(-116.7, -120.4, -116.7, -113.5, -115.5, -120.8, -119.5, -113.7, -113.7, -110.7)
latitude <- c(45.3, 42.6, 38.9, 42.1, 35.7, 38.9, 36.2, 39, 41.6, 36.9)
pts <- vect(cbind(longitude, latitude), crs="+proj=longlat")
  
lon <- c(-116.8, -114.2, -112.9, -111.9, -114.2, -115.4, -117.7)
lat <- c(41.3, 42.9, 42.4, 39.8, 37.6, 38.3, 37.6)
lonlat <- cbind(id=1, part=1, lon, lat)
pols <- vect(lonlat, type="polygons", crs="+proj=longlat")

Here are four approaches

## 1
e <- extract(pols, pts) 
e[!is.na(e[,2]), 1]
#[1] 3 4 8 9
 
## 2
relate(pols, pts, "contains") |> which()
#[1] 3 4 8 9

## 3 
is.related(pts, pols, "intersects") |> which()
#[1] 3 4 8 9

## 4
pts$id <- 1:nrow(pts)
intersect(pts, pols)$id 
#[1] 3 4 8 9

So you were on the right track with extract but you had the arguments in the wrong order. extract may be the easiest if you have multiple polygons.

More examples here

Answer 2

I think you are looking for terra::intersect.

points(intersect(pols, pts), col = 'blue')

Answer 3

There's a relate function that allows specification of "contains" as the desired relationship:

relate(pols, pts, "contains")
      [,1]  [,2] [,3] [,4]  [,5]  [,6]  [,7] [,8] [,9] [,10]
[1,] FALSE FALSE TRUE TRUE FALSE FALSE FALSE TRUE TRUE FALSE