Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

point in polygon using terra package in R

Tags:

r

terra

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.

like image 889
89_Simple Avatar asked Nov 01 '25 07:11

89_Simple


2 Answers

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

like image 114
Robert Hijmans Avatar answered Nov 02 '25 20:11

Robert Hijmans


I think you are looking for terra::intersect.

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

enter image description here

like image 39
hrvg Avatar answered Nov 02 '25 22:11

hrvg