rgispolygonintersectionr-sf

Is there a way for st_intersection to return a value of 0 instead of cutting non intersecting features?

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

Solution

  • 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