Count number of times that sf linestrings intersect grid cell in r

Question

Consider a set of linestrings and a polygon grid (sf geometry):

library(sf)

#creating data example
id <- c("844", "844", "844", "844", "844","855", "855", "855", "855", "855")

lat <- c(-30.6456, -29.5648, -28.6667, -31.5587, -30.6934, -29.3147, -28.0538, 
         -26.5877, -26.6923, -27.40865)
long <- c(-50.4879, -49.8715, -51.8716, -50.4456, -50.9842, -51.9787, -47.2343, 
          -49.2859, -48.19599, -49.64302)

df <- data.frame(id = as.factor(id), lat, long)

#converting to sf
df.sf <- df %>% 
  sf::st_as_sf(coords = c("long", "lat"), crs = 4326)

#creating linestrings
df.line <- df.sf %>% 
  dplyr::group_by(id) %>%
  dplyr::summarize() %>%
  sf::st_cast("LINESTRING") 

#creating grid
xy <- sf::st_coordinates(df.sf)

grid <- sf::st_make_grid(sf::st_bbox(df.sf),
                        cellsize = 1, square = FALSE) %>%
  sf::st_as_sf() %>%
  dplyr::mutate(cell = 1:nrow(.))

intersection <- sf::st_intersection(grid, df.line)

If I use the st_intersection() function, it will count once per feature of linestrings that intersect each cell. But, the same feature may have crossed that cell twice or more, and I would like that number of times to be counted. It's possible?

I found this question with a answer that is very similar to what I need, but I would like to do the process with sf objects instead of sp.

Count lines crossing raster cells with R

Answer 1

You can use the rmapshaper package to 'erase' the part of the lines that are close to the grid. Then count the lines within each cell with st_intersection.

library(sf)
library(tidyverse)

########### from the question ############
#creating data example
id <- c("844", "844", "844", "844", "844","855", "855", "855", "855", "855")

lat <- c(-30.6456, -29.5648, -28.6667, -31.5587, -30.6934, -29.3147, -28.0538, 
         -26.5877, -26.6923, -27.40865)
long <- c(-50.4879, -49.8715, -51.8716, -50.4456, -50.9842, -51.9787, -47.2343, 
          -49.2859, -48.19599, -49.64302)

df <- data.frame(id = as.factor(id), lat, long)

#converting to sf
df.sf <- df %>% 
  sf::st_as_sf(coords = c("long", "lat"), crs = 4326)

#creating linestrings
df.line <- df.sf %>% 
  dplyr::group_by(id) %>%
  dplyr::summarize() %>%
  sf::st_cast("LINESTRING") 

#creating grid
xy <- sf::st_coordinates(df.sf)

grid <- sf::st_make_grid(sf::st_bbox(df.sf),
                         cellsize = 1, square = FALSE) %>%
  sf::st_as_sf() %>%
  dplyr::mutate(cell = 1:nrow(.))
###### End from question ##########


# Transform to a crs that uses meters for buffering
df.line <- st_transform(df.line, 3857)
grid <- st_transform(grid, 3857)

# Create a geometry to 'cut' the lines with by buffering the grid by 100m.
# You may want to change the buffer distance.
blade <- st_cast(grid, 'MULTILINESTRING') %>% 
  st_buffer(100) %>% 
  st_union() %>% 
  st_as_sf()

# erase the parts of the lines that cross the buffered grids
line_split <- rmapshaper::ms_erase(df.line, blade) %>%
  st_cast('LINESTRING')

split_count <- st_intersection(line_split, grid) %>%
                  group_by(id, cell) %>% 
                  count() %>%
                  arrange(desc(n))

head(split_count)
#> Simple feature collection with 6 features and 3 fields
#> Geometry type: GEOMETRY
#> Dimension:     XY
#> Bounding box:  xmin: -5780918 ymin: -3698418 xmax: -5619363 ymax: -3325137
#> Projected CRS: WGS 84 / Pseudo-Mercator
#> # A tibble: 6 × 4
#>   id     cell     n                                                     geometry
#>   <fct> <int> <int>                                               <GEOMETRY [m]>
#> 1 844      15     2 MULTILINESTRING ((-5644929 -3650436, -5674823 -3594332), (-…
#> 2 855       9     2 MULTILINESTRING ((-5688528 -3325137, -5724777 -3358756, -57…
#> 3 844       8     1 LINESTRING (-5675023 -3593956, -5675535 -3592995, -5675023 …
#> 4 844       9     1 LINESTRING (-5713732 -3433012, -5733856 -3399891, -5746569 …
#> 5 844      11     1            LINESTRING (-5619363 -3698418, -5644729 -3650812)
#> 6 844      12     1            LINESTRING (-5649610 -3538547, -5713628 -3433183)

In the example data, there are two grid cells that are crossed more than once by either of the lines, circled in red:

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