Intersects point and polygon feature classes and adds polygon attributes to points
If duplicate argument is TRUE and more than one polygon intersection occurs, points will be duplicated (new row added) and all attributes joined. However, if duplicate is FALSE, with duplicate intersections, a new column for each unique intersecting polygon will be returned and the points will not be duplicated. For example, if a point intersect three polygons, three new columns will be added representing the polygons ID.
point.in.poly(x, y, sp = TRUE, duplicate = TRUE, ...)
sp SpatialPointsDataFrame or SpatialPoints or sf point object
sp SpatialPolygonsDataFrame or sf polygon object
(TRUE/FALSE) Return an sp class object, else returns sf class object
(TRUE/FALSE) Return duplicated features with more than one polygon intersection
Additional arguments passed to sf::st_join
A SpatialPointsDataFrame or sf
# NOT RUN {
#### Simple one-to-one feature overlay.
require(sp)
data(meuse)
coordinates(meuse) = ~x+y
meuse@data$test.na <- NA
sr1=Polygons(list(Polygon(cbind(c(180114, 180553, 181127, 181477, 181294,
181007, 180409, 180162, 180114), c(332349, 332057, 332342, 333250, 333558,
333676, 332618, 332413, 332349)))),'10')
sr2=Polygons(list(Polygon(cbind(c(180042, 180545, 180553, 180314, 179955, 179142,
179437, 179524, 179979, 180042), c(332373, 332026, 331426, 330889, 330683,
331133, 331623, 332152, 332357, 332373)))),'20')
sr3=Polygons(list(Polygon(cbind(c(179110, 179907, 180433, 180712, 180752, 180329,
179875, 179668, 179572, 179269, 178879, 178600, 178544, 179046, 179110),
c(331086, 330620, 330494, 330265, 330075, 330233, 330336, 330004,
329783, 329665, 329720, 329933, 330478, 331062, 331086)))),'30')
sr4=Polygons(list(Polygon(cbind(c(180304, 180403,179632,179420,180304),
c(332791, 333204, 333635, 333058, 332791)))),'40')
sr=SpatialPolygons(list(sr1,sr2,sr3,sr4))
polys=SpatialPolygonsDataFrame(sr, data.frame(row.names=c('10','20','30','40'),
PIDS=1:4, y=runif(4)))
polys@data$pid <- polys@data$PIDS + 100
plot(polys)
plot(meuse, pch=19, add=TRUE)
# Point in polygon overlay
pts.poly <- point.in.poly(meuse, polys)
head(pts.poly@data)
# Count points in each polygon
tapply(pts.poly$cadmium, pts.poly$pid, FUN=length)
#### Complex many-to-one feature overlay.
require(sf)
p <- sf::st_polygon(list(rbind(c(0,0), c(1,0), c(1,1), c(0,1), c(0,0))))
polys <- sf::st_sf(sf::st_sfc(p, p + c(.8, .2), p + c(.2, .8)))
pts <- sf::st_sf(sf::st_sample(polys, size=100))
# Duplicates points for each new polygon, no attributes so returns IDs for features
pts.poly.dup <- point.in.poly(pts, polys)
head(pts.poly.dup@data)
# }
# NOT RUN {
# **** Should throw error due to lack of attributes ****
pts.poly <- point.in.poly(pts, polys, duplicate = FALSE)
# }
# NOT RUN {
# Coerce to sp class objects
x <- as(pts, "Spatial")
x <- SpatialPointsDataFrame(x, data.frame(IDS=1:nrow(x), pty=runif(nrow(x))))
y <- as(polys, "Spatial")
y <- SpatialPolygonsDataFrame(y, data.frame(IDS=1:nrow(y), py=runif(nrow(y))))
# Returns point attributes with column for each unique polygon
pts.poly <- point.in.poly(x, y, duplicate = FALSE)
head(pts.poly@data)
# Duplicates points for each new polygon, joins all attributes
pts.poly.dup <- point.in.poly(x, y)
head(pts.poly.dup@data)
# Count points in each polygon
tapply(pts.poly.dup$IDS.x, pts.poly.dup$IDS.y, FUN=length)
# }
Run the code above in your browser using DataLab