st_intersection only returning part of an intersecting geometry (line)

Potentially unnecessary background:

I am attempting to calculate the sinuosity of road within a 500m buffer of a point on that road.

Workflow is:

Buffer the point 500m
Use that buffer to clip the road layer to a 500m buffer around the point
Use st_intersection to pull out the specific road that the point intersects from the clipped road network from Step 2
Calculate the sinuosity of that specific intersecting road

Actual problem

I am having a problem where when I use st_intersection to pull the specific road segment that the point intersects, it is only returning PART of that road segment. I just can't figure out why it is not returning the full road segment geometry. Reproducible code highlighting the problem is below.

library(tidyverse)
library(sf)

#### create dummy point and line
point <- data.frame(
  utme = c(245862.2),
  utmn   = c(4891593),
  id = c (2)) %>%
  st_as_sf(coords = c("utme", "utmn"), na.fail = FALSE, crs = "epsg:26920") 
   
   
   line_df <- matrix(c( 
          245843.6, 4891449, 
          245844 ,4891505, 
          245845.7, 4891518, 
          245848 ,4891529, 
          246065.2, 4891972, 
          246066.5 ,4891989, 
          246067.8 ,4892006, 
          246073.1 ,4892075,
          246077 ,4892125, 
          246079.8, 4892169,
          246080.3, 4892174,
          246080.4 ,4892176, 
          246082.4 ,4892202, 
          246082.6 ,4892206, 
          245848 ,4891529, 
          245855.8, 4891568, 
          245862.2 ,4891593, 
          245871.7 ,4891618, 
          245882.1 ,4891642,
          245891.4 ,4891658,
          245908 ,4891688,
          245910.9, 4891691,
          245941.8 ,4891733,
          245967.4 ,4891766,
          246003.3 ,4891815,
          246023.9 ,4891847, 
          246025 ,4891849,
          246034.3, 4891867, 
          246042.6 ,4891886, 
          246049 ,4891906, 
          246050.1, 4891911,
          246057.6 ,4891939,
          246065.2 ,4891972),  byrow=TRUE, ncol = 2) %>%
     data.frame() %>%
     # rowwise() %>%
     mutate(id = 1 ) %>%
     rename(utme = X1, utmn = X2) %>%
     st_as_sf( coords = c("utme", "utmn"), crs = "epsg:26920") %>% 
     group_by(id) %>% 
     summarize() %>%
     st_cast("MULTILINESTRING") 
   
   plot(point)
   plot(line_df)
   
   sin_test <- st_intersection(line_df, st_buffer(point, .1)) %>%
     select(id)
   
   plot(sin_test)

#### As you can see from the plot(sin_test) and the plot(line_df), 
#### the sin_test object has had it's #### geometry clipped.....
#### Thanks in advance for the help!

Solution

The behaviour you observed is what st_intersection() is supposed to do. If you want to find roads that are intersecting with the buffer, you are after st_intersects() predicate function. Though youd probably want to use it through subsetting / filtering geometries by other geometries or through spatial joins.

library(tidyverse)
library(sf)
library(mapview)

# check input:
point
#> Simple feature collection with 1 feature and 1 field
#> Geometry type: POINT
#> Dimension:     XY
#> Bounding box:  xmin: 245862.2 ymin: 4891593 xmax: 245862.2 ymax: 4891593
#> Projected CRS: NAD83 / UTM zone 20N
#> # A tibble: 1 × 2
#>      id           geometry
#> * <dbl>        <POINT [m]>
#> 1     2 (245862.2 4891593)
line_df
#> Simple feature collection with 1 feature and 1 field
#> Geometry type: MULTILINESTRING
#> Dimension:     XY
#> Bounding box:  xmin: 245843.6 ymin: 4891449 xmax: 246082.6 ymax: 4892206
#> Projected CRS: NAD83 / UTM zone 20N
#> # A tibble: 1 × 2
#>      id                                                                 geometry
#>   <dbl>                                                    <MULTILINESTRING [m]>
#> 1     1 ((245843.6 4891449, 245844 4891505, 245845.7 4891518, 245848 4891529, 2…
mapview(line_df, label = "line_df", color = "darkgreen") +
mapview(point, label = "point", col.regions = "darkred")

# line 1 intersecting with point 1 buffer:
st_intersects(line_df, st_buffer(point, .1))
#> Sparse geometry binary predicate list of length 1, where the predicate
#> was `intersects'
#>  1: 1

st_intersects is a default predicate when subsetting geometries by other geometries:

# subset by intersecting geometries:
(line1 <- line_df[st_buffer(point, .1),])
#> Simple feature collection with 1 feature and 1 field
#> Geometry type: MULTILINESTRING
#> Dimension:     XY
#> Bounding box:  xmin: 245843.6 ymin: 4891449 xmax: 246082.6 ymax: 4892206
#> Projected CRS: NAD83 / UTM zone 20N
#> # A tibble: 1 × 2
#>      id                                                                 geometry
#>   <dbl>                                                    <MULTILINESTRING [m]>
#> 1     1 ((245843.6 4891449, 245844 4891505, 245845.7 4891518, 245848 4891529, 2…

Or when using a spatial join:

# join by buffer intersection, include all features from x (line_df) and 
# add fields from intersecting y features:
(line2 <- st_join(line_df, st_buffer(point, .1)))
#> Simple feature collection with 1 feature and 2 fields
#> Geometry type: MULTILINESTRING
#> Dimension:     XY
#> Bounding box:  xmin: 245843.6 ymin: 4891449 xmax: 246082.6 ymax: 4892206
#> Projected CRS: NAD83 / UTM zone 20N
#> # A tibble: 1 × 3
#>    id.x                                                           geometry  id.y
#> * <dbl>                                              <MULTILINESTRING [m]> <dbl>
#> 1     1 ((245843.6 4891449, 245844 4891505, 245845.7 4891518, 245848 4891…     2

For the task you describe, there's also st_is_within_distance (here the result is identical to the previous join):

# join using a different predicate, st_is_within_distance
line3 <- st_join(line_df, point, join = st_is_within_distance, dist = .1)

Resulting line geometries are identical to line_df geometry, i.e. not clipped by the buffer:

identical(st_geometry(line_df), st_geometry(line1))
#> [1] TRUE
identical(st_geometry(line_df), st_geometry(line2))
#> [1] TRUE
identical(st_geometry(line_df), st_geometry(line3))
#> [1] TRUE

^{Created on 2025-01-05 with reprex v2.1.1}