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

user_7547327 · Answer

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
#>                                                    
#> 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)}

user_7756889 · Answer

I suggest casting your intersection first to MULTILINESTRING and then to LINESTRING (since you start with GEOMETRY which can't be easily converted to LINESTRING, you need the intermediate step).

This way you will go from 22 rows in your intersection dataset to 24 in my proposed grid_lines dataset, with the id's of cells & lines of the two cells crossed twice by the same line duplicated.

Then it is a simple case of {dplyr} aggreggation & sorting / I suggest dropping geometry at this stage, as it no longer adds any value, and work at the level of cell id's and line id's.

# a visual overview
plot(st_geometry(grid))
plot(intersection["id"], add = T)

# split lines to two /or more, if applicable
grid_lines <- st_cast(intersection, "MULTILINESTRING") %>% 
  st_cast("LINESTRING")

library(dplyr)

grid_lines %>% 
  st_drop_geometry() %>% # no longer relevant...
  group_by(id, cell) %>% 
  tally() %>% 
  arrange(desc(n))

# A tibble: 22 × 3
# Groups:   id [2]
   id     cell     n
     
 1 844      15     2
 2 855       9     2
 3 844       8     1
 4 844       9     1
 5 844      11     1
 6 844      12     1
 7 844      16     1
 8 844      18     1
 9 844      19     1
10 855       5     1
# … with 12 more rows

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

2 Answers