I have a 1km buffer shapefile and a roads shapefile and I'm trying to figure out the total length of roads that intersect with each buffer so I can calculate road density within each buffer. There are a few buffers that do not intersect with any roads so they are cut from the result. I still need the values for those buffers even if they are equal to zero. Any insight would be appreciated thank you. My code is below:
### Loop to calculate each buffer and erase the original point, one point/buffer at a time
#START OF LOOP
for (i in 1:nrow(p_rep)){
buff1km <- outerBuffer(p_rep[i,], dist)
st_write(buff1km, "1kmbuff.shp", driver = "ESRI Shapefile", sep="", append = TRUE)
} #END OF LOOP
### Read in the buffer shapefile
buffer_only <- st_read("1kmbuff.shp")
### read in roads shapefile:
## Road shapefile accessed from state of michigan GIS opendata
## https://gis-michigan.opendata.arcgis.com/datasets/Michigan::all-roads-v17a/about
roads <- read_sf(dsn = "All_Roads_(v17a).shp")
# Distance (length) unit of measure is meters.
roads
roads.prep <- st_transform(roads, crs(buffer_only))
roads.buff <- sf::st_intersection(roads.prep, buffer_only) %>%
mutate(length_m = st_length(geometry))
I tried using st_intersection but it cut out the rows that didn't intersect.
Example data for roads and center points for buffers can be found at: https://drive.google.com/drive/folders/1ll-VnNsGv29PHgAnZShNHU9QcvgX0v7F?usp=sharing
You can join buffer_only
to the resulting frame. If spatial st_join()
is too expensive and you can identify your buffers by some unique id, you can use regular dplyr::right_join()
or left_join
, just make sure not to pass 2 sf objects:
library(sf)
#> Linking to GEOS 3.12.1, GDAL 3.8.4, PROJ 9.3.1; sf_use_s2() is TRUE
library(dplyr)
st_intersection(roads, buffer_only) |>
mutate(length = st_length(geometry)) |>
right_join(st_drop_geometry(buffer_only)) |>
st_sf() |>
tidyr::replace_na(list(length = 0))
#> Warning: attribute variables are assumed to be spatially constant throughout
#> all geometries
#> Joining with `by = join_by(id_buff)`
#> Simple feature collection with 3 features and 3 fields (with 1 geometry empty)
#> Geometry type: LINESTRING
#> Dimension: XY
#> Bounding box: xmin: 1 ymin: 1 xmax: 3 ymax: 2
#> CRS: NA
#> id_road id_buff length geometry
#> 1 1 2 2 LINESTRING (1 1, 3 1)
#> 2 2 2 1 LINESTRING (2 2, 2 1)
#> 3 NA 1 0 LINESTRING EMPTY
Example data:
buffer_only <-
st_as_sfc(c("POINT (1 4)", "POINT (2 1)")) |>
st_buffer(1) |>
st_sf(geometry = _, id_buff = 1:2)
buffer_only
#> Simple feature collection with 2 features and 1 field
#> Geometry type: POLYGON
#> Dimension: XY
#> Bounding box: xmin: 0 ymin: 0 xmax: 3 ymax: 5
#> CRS: NA
#> id_buff geometry
#> 1 1 POLYGON ((2 4, 1.99863 3.94...
#> 2 2 POLYGON ((3 1, 2.99863 0.94...
roads <-
st_as_sfc(c("LINESTRING (0 1, 4 1)", "LINESTRING (2 3, 2 1)")) |>
st_sf(geometry = _, id_road = 1:2)
roads
#> Simple feature collection with 2 features and 1 field
#> Geometry type: LINESTRING
#> Dimension: XY
#> Bounding box: xmin: 0 ymin: 1 xmax: 4 ymax: 3
#> CRS: NA
#> id_road geometry
#> 1 1 LINESTRING (0 1, 4 1)
#> 2 2 LINESTRING (2 3, 2 1)
plot(buffer_only$geometry, col = "gold", axes = T)
plot(roads$geometry, lwd = 3, add = T)
Created on 2024-10-31 with reprex v2.1.1