Point and Polygon Intersect

Source: R/point.in.poly.R

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, ...)

Arguments

x

sp SpatialPointsDataFrame or SpatialPoints or sf point object
y

sp SpatialPolygonsDataFrame or sf polygon object
sp

(TRUE/FALSE) Return an sp class object, else returns sf class object
duplicate

(TRUE/FALSE) Return duplicated features with more than one polygon intersection
...

Additional arguments passed to sf::st_join

Value

A SpatialPointsDataFrame or sf

Author

Jeffrey S. Evans <jeffrey_evans@tnc.org>

Examples

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

#>   cadmium copper lead zinc  elev       dist   om ffreq soil lime landuse dist.m
#> 1    11.7     85  299 1022 7.909 0.00135803 13.6     1    1    1      Ah     50
#> 2     8.6     81  277 1141 6.983 0.01222430 14.0     1    1    1      Ah     30
#> 3     6.5     68  199  640 7.800 0.10302900 13.0     1    1    1      Ah    150
#> 4     2.6     81  116  257 7.655 0.19009400  8.0     1    2    0      Ga    270
#> 5     2.8     48  117  269 7.480 0.27709000  8.7     1    2    0      Ah    380
#> 6     3.0     61  137  281 7.791 0.36406700  7.8     1    2    0      Ga    470
#>   test.na PIDS         y pid
#> 1      NA    1 0.6715948 101
#> 2      NA    1 0.6715948 101
#> 3      NA    1 0.6715948 101
#> 4      NA    1 0.6715948 101
#> 5      NA    1 0.6715948 101
#> 6      NA    1 0.6715948 101
 
 # Count points in each polygon
 tapply(pts.poly$cadmium, pts.poly$pid, FUN=length)    

#> 101 102 103 
#>  50  57  48 
   
 #### 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)

#>     pt.ids poly.ids
#> 1        1        1
#> 1.1      1        3
#> 2        2        1
#> 3        3        1
#> 3.1      3        2
#> 4        4        1
   
if (FALSE) {
 # **** Should throw error due to lack of attributes ****  
   pts.poly <- point.in.poly(pts, polys, duplicate = FALSE) 
}
 
 # 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)  

#>   IDS        pty    p pid1 pid2 pid3
#> 1   1 0.99692867    1    1    3 <NA>
#> 2   2 0.02802265    2    1 <NA> <NA>
#> 3   3 0.02596293    3    1    2 <NA>
#> 4   4 0.02129025 <NA>    1 <NA> <NA>
#> 5   5 0.83191578 <NA>    2 <NA> <NA>
#> 6   6 0.95174343 <NA>    1    3 <NA>
   
 # Duplicates points for each new polygon, joins all attributes  
 pts.poly.dup <-  point.in.poly(x, y)
   head(pts.poly.dup@data)  

#>     IDS.x        pty IDS.y        py
#> 1       1 0.99692867     1 0.8987087
#> 1.1     1 0.99692867     3 0.9064305
#> 2       2 0.02802265     1 0.8987087
#> 3       3 0.02596293     1 0.8987087
#> 3.1     3 0.02596293     2 0.9470662
#> 4       4 0.02129025     1 0.8987087
 
 # Count points in each polygon
 tapply(pts.poly.dup$IDS.x, pts.poly.dup$IDS.y, FUN=length)      

#>  1  2  3 
#> 46 43 34