Search code examples
rterra

point in polygon using terra package in R

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)

enter image description here

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.

Solution

  • 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