Intersecting Points and Polygons in R

Intersecting Points and Polygons in R

If you do overlay(pts, polys) where pts is a SpatialPointsDataFrame object and polys is a SpatialPolygonsDataFrame object then you get back a vector the same length as the points giving the row of the polygons data frame. So all you then need to do to combine the polygon data onto the points data frame is:

 o = overlay(pts, polys)
 pts@data = cbind(pts@data, polys[o,])

HOWEVER! If any of your points fall outside all your polygons, then overlay returns an NA, which will cause polys[o,] to fail, so either make sure all your points are inside polygons or you'll have to think of another way to assign values for points outside the polygon...

Determine lines that intersect a polygon in R

If you have your polygons and lines in a standard format you can use st_intersects
from the sf package. Below is an example using three counties of North Carolina as
polygons.

Load sf:

library(sf)
#> Linking to GEOS 3.8.0, GDAL 3.0.2, PROJ 6.2.1

Make polygonal data (disclaimer: for quick illustration geographical longitude,latitude
coordinates are used here, but really planar projected coordinates should be used):

nc <- st_read(system.file("shape/nc.shp", package="sf"), quiet = TRUE)
n <- 3
poly_dat <- st_sf(poly_id = LETTERS[1:n], nc$geometry[1:n])

Make line data:

line1 <- st_linestring(rbind(c(-81.8,36.5),c(-80.5,36.5)))
line2 <- st_linestring(rbind(c(-81,36.3),c(-80.5,36.3)))
line_dat <- st_sf(linename = c("1", "2"), geometry = st_sfc(line1, line2, crs = st_crs(poly_dat)))

Plot:

plot(poly_dat, reset = FALSE)
plot(line_dat, add = TRUE, col = 1:2, lty = 2)
legend("topright", legend = paste("Line", 1:2), col = 1:2, lty = 2)

Find intersections:

st_intersects(line_dat, poly_dat)
#> although coordinates are longitude/latitude, st_intersects assumes that they are planar
#> Sparse geometry binary predicate list of length 2, where the predicate was `intersects'
#>  1: 1, 2, 3
#>  2: 3

Intersection of polygons in R using sf

How can I : asess the degree of spatial proximity of each point to other equivalent points by looking at the number of others within 400m (5 minute walk).

Here is how to add a column to your initial sfc of pollings stations that tells you how many polling stations are within 400m of each feature in that sfc.

Note that the minimum value is 1 because a polling station is always within 400m of itself.

# n_neighbors shows how many polling stations are within 400m
polls %>% 
  mutate(n_neighbors = lengths(st_is_within_distance(polls, dist = 400)))

Similarly, for your sfc collection of intersecting polygons, you could add a column that counts the number of buffer polygons that contain each intersection polygon:

polls_intersection %>% 
  mutate(n_overlaps = lengths(st_within(geometry, polls_buffer_400)))

And this is the bit I'm not sure about, to get to the output I want (which will show "Hotspots" of polling stations in this case) how do I colour things?

If you want to plot these things I highly recommend using ggplot2. It makes it very clear how you associate an attribute like colour with a specific variable.

For example, here is an example mapping the alpha (transparency) of each polygon to a scaled version of the n_overlaps column:

library(ggplot2)
polls_intersection %>% 
  mutate(n_overlaps = lengths(st_covered_by(geometry, polls_buffer_400))) %>% 
  ggplot() + 
  geom_sf(aes(alpha = 0.2*n_overlaps), fill = "red") 

Lastly, there should be a better way to generate your intersecting polygons that already counts overlaps. This is built in to the st_intersection function for finding intersections of sfc objects with themselves.

However, your data in particular generates an error when you try to do this:

st_intersection(polls_buffer_400)

# > Error in CPL_nary_intersection(x) : 
#>  Evaluation error: TopologyException: side location conflict at 315321.69159061194 199694.6971799387.

I don't know what a "side location conflict" is. Maybe @edzer could help with that. However, most subsets of your data do not contain that conflict. For example:

# this version adds an n.overlaps column automatically:
st_intersection(polls_buffer_400[1:10,]) %>% 
  ggplot() + geom_sf(aes(alpha = 0.2*n.overlaps), fill = "red") 

Intersection keeping non-intersecting polygons in sf - R

Please consider this approach: Once that you have your intersection, you can remove the intersecting parts with st_difference. That would effectively split the intersecting SA2 in zones based on ARKS, and leave the rest as they are originally. After that, you can rejoin the dataset with dplyr::bind_rows, having the ARKS layer, the SA2 intersected split and the SA2 non-intersected as they are originally:

library(sf)
library(ggplot2)
library(dplyr)

SA2 <- st_read("./SA2_2021_AUST_SHP_GDA94/SA2_2021_AUST_GDA94.shp")
ARKS <-
  st_read(
    "./Biodiversity_KoalaPrioritisationProjectNSW_ARKS/KoalaPrioritisationProjectNSW_ARKS.shp"
  )
ARKS <- st_transform(ARKS, 3577)
SA2 <- st_transform(SA2, 3577)
SA2 <- SA2[SA2$STE_CODE21 == "1", ]

intersection <- st_intersection(SA2, ARKS)

# Control var, you can remove this
intersection$koala <- "Yes"

# Difference
SA2_diff <-
  st_difference(SA2, st_union(st_geometry(intersection)))

# Check visually, we should see the gaps
ggplot(SA2_diff) +
  geom_sf(fill = "green")

# Join all back
SA2_end <- dplyr::bind_rows(
  SA2_diff,
  intersection
)

ggplot(SA2_end) +
  geom_sf(aes(fill = koala))

R intersecting multiple (2+) overlapping polygons

It sounds like you just want the buffer(s), but without having to deal with all of the overlapping sections. It doesn't matter if a person is within 400ft of one bus-stop or three, right?

If so, you can use the st_union function to "blend" the buffers together.

library(tidytransit)
library(sf)
library(mapview)
library(ggplot2)

# s2 true allows buffering in meters, s2 off later speeds things up
sf::sf_use_s2(TRUE)
ART2019Path <- file.path("/your/file/path/")
ART2019GTFS <- read_gtfs(ART2019Path)
ART2019StopLoc <- stops_as_sf(ART2019GTFS$stops) ### Make a spatial file for stops
ART2019Buffer <- st_buffer(ART2019StopLoc, dist = 121.92) ### Make buffer with 400ft (121.92m) radius

# might be needed due to some strange geometries in buffer, and increase speed
sf::sf_use_s2(FALSE)
#> Spherical geometry (s2) switched off

# MULTIPOLYGON  sfc object covering only the buffered areas,
#  there are no 'overlaps'.
buff_union <- st_union(st_geometry(ART2019Buffer))
#> although coordinates are longitude/latitude, st_union assumes that they are planar

buff_union
#> Geometry set for 1 feature 
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -77.16368 ymin: 38.83828 xmax: -77.04768 ymax: 38.9263
#> Geodetic CRS:  WGS 84
#> MULTIPOLYGON (((-77.08604 38.83897, -77.08604 3...

# Non-overlapping buffer & stops
ggplot() + 
  geom_sf(data = buff_union, fill = 'blue', alpha = .4) +
  geom_sf(data = ART2019StopLoc, color = 'black') + 
  coord_sf(xlim = c(-77.09, -77.07),
           ylim = c(38.885, 38.9)) 

# Overlapping buffer & stops
ggplot() +
  geom_sf(data = ART2019Buffer, fill = 'blue', alpha = .4) +
  geom_sf(data = ART2019StopLoc, color = 'black') +
  coord_sf(xlim = c(-77.09, -77.07),
           ylim = c(38.885, 38.9)) 

# Back to original settings
sf::sf_use_s2(TRUE)
#> Spherical geometry (s2) switched on

^{Created on 2022-04-18 by the reprex package (v2.0.1)}

polygon intersection in R

You may convert your data.frame to a SpatialPolygons.

library(sp)
polys <- lapply(unique(coords$ID), function(i) {
    Polygons(list(Polygon(coords[coords$ID==i, 1:2])), ID=i)
})
spa_polys <- SpatialPolygons(polys)

Also, I have found that gIntersects doesn't work for a point (one-point spatial polygon) or a line (two-point spatial polygon).

library(rgeos)
> gIntersects(spa_polys[13], spa_polys[1])
Error in RGEOSBinPredFunc(spgeom1, spgeom2, byid, func) : 
  IllegalArgumentException: Invalid number of points in LinearRing found 3 - must be 0 or >= 4
> gIntersects(spa_polys[13], spa_polys[2])
[1] FALSE
> gIntersects(spa_polys[13], spa_polys[3])
Error in RGEOSBinPredFunc(spgeom1, spgeom2, byid, func) : 
  IllegalArgumentException: point array must contain 0 or >1 elements
> gIntersects(spa_polys[13], spa_polys[13])
[1] TRUE

R Overlay points and polygons with a certain degree of tolerance

The example data -

set.seed(1)
library(raster)
library(rgdal)
library(sp)
p <- shapefile(system.file("external/lux.shp", package="raster"))
p2 <- as(0.30*extent(p), "SpatialPolygons")
proj4string(p2) <- proj4string(p)
pts1 <- spsample(p2-p, n=3, type="random")
pts2<- spsample(p, n=10, type="random")
pts<-rbind(pts1, pts2)

## Plot to visualize
plot(p, col=colorRampPalette(blues9)(12))
plot(pts, pch=16, cex=.5,col="red", add=TRUE)

Solution using sf and nngeo packages -

library(nngeo)

# Convert to 'sf'
pts = st_as_sf(pts)
p = st_as_sf(p)

# Spatial join
p1 = st_join(pts, p, join = st_nn)
p1

## Simple feature collection with 13 features and 5 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: 5.795068 ymin: 49.54622 xmax: 6.518138 ymax: 50.1426
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs
## First 10 features:
##   ID_1       NAME_1 ID_2           NAME_2 AREA                  geometry
## 1     1     Diekirch    4          Vianden   76 POINT (6.235953 49.91801)
## 2     1     Diekirch    4          Vianden   76 POINT (6.251893 49.92177)
## 3     1     Diekirch    4          Vianden   76  POINT (6.236712 49.9023)
## 4     1     Diekirch    1         Clervaux  312  POINT (6.090294 50.1426)
## 5     1     Diekirch    5            Wiltz  263  POINT (5.948738 49.8796)
## 6     2 Grevenmacher   12     Grevenmacher  210 POINT (6.302851 49.66278)
## 7     2 Grevenmacher    6       Echternach  188 POINT (6.518138 49.76773)
## 8     3   Luxembourg    9 Esch-sur-Alzette  251 POINT (6.116905 49.56184)
## 9     1     Diekirch    3          Redange  259 POINT (5.932418 49.78505)
## 10    2 Grevenmacher    7           Remich  129 POINT (6.285379 49.54622)

Plot showing which polygons and points are joined -

# Visuzlize join
l = st_connect(pts, p, dist = 1)
plot(st_geometry(p))
plot(st_geometry(pts), add = TRUE)
plot(st_geometry(l), col = "red", lwd = 2, add = TRUE)

EDIT:

# Spatial join with 100 meters threshold
p2 = st_join(pts, p, join = st_nn, maxdist = 100)
p2
## Simple feature collection with 13 features and 5 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: 5.795068 ymin: 49.54622 xmax: 6.518138 ymax: 50.1426
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs
## First 10 features:
##   ID_1       NAME_1 ID_2           NAME_2 AREA                  geometry
## 1    NA         <NA> <NA>             <NA>   NA POINT (6.235953 49.91801)
## 2    NA         <NA> <NA>             <NA>   NA POINT (6.251893 49.92177)
## 3     1     Diekirch    4          Vianden   76  POINT (6.236712 49.9023)
## 4     1     Diekirch    1         Clervaux  312  POINT (6.090294 50.1426)
## 5     1     Diekirch    5            Wiltz  263  POINT (5.948738 49.8796)
## 6     2 Grevenmacher   12     Grevenmacher  210 POINT (6.302851 49.66278)
## 7     2 Grevenmacher    6       Echternach  188 POINT (6.518138 49.76773)
## 8     3   Luxembourg    9 Esch-sur-Alzette  251 POINT (6.116905 49.56184)
## 9     1     Diekirch    3          Redange  259 POINT (5.932418 49.78505)
## 10    2 Grevenmacher    7           Remich  129 POINT (6.285379 49.54622)

---

Related Topics