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

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.

# 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
   <fct> <int> <int>
 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

CodePudding user response：

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