计算 sf 线串与 r 中的网格单元相交的次数
Count number of times that sf linestrings intersect grid cell in r
考虑一组线串和一个多边形网格（sf 几何）：
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)
如果我使用 st_intersection() 函数，它将对与每个单元格相交的线串的每个特征计数一次。但是，相同的特征可能已经穿过该单元格两次或更多次，我希望计算该次数。可能吗？
我发现这个问题的答案与我需要的非常相似，但我想用 sf 对象而不是 sp 来完成这个过程。
Count lines crossing raster cells with R
我建议先将交点转换为 MULTILINESTRING，然后再转换为 LINESTRING（因为您从 GEOMETRY 开始，无法轻松转换为 LINESTRING，因此需要中间步骤）。
这样，您将从 intersection 数据集中的 22 行变为我建议的 grid_lines 数据集中的 24 行，并且两个单元格的 ID 和同一行交叉两次重复。
那么这是一个简单的 {dplyr} 聚合和排序案例/我建议在此阶段删除几何图形，因为它不再增加任何值，并在单元 ID 和行 ID 级别工作。
# 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
<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
您可以使用rmapshaper包来'erase'靠近网格的部分线。然后使用 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)
在示例数据中，有两个网格单元被任何一条线交叉多次，用红色圈出：
由 reprex package (v2.0.1)
于 2022-04-22 创建
考虑一组线串和一个多边形网格（sf 几何）：
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)
如果我使用 st_intersection() 函数，它将对与每个单元格相交的线串的每个特征计数一次。但是，相同的特征可能已经穿过该单元格两次或更多次，我希望计算该次数。可能吗？
我发现这个问题的答案与我需要的非常相似，但我想用 sf 对象而不是 sp 来完成这个过程。
Count lines crossing raster cells with R
我建议先将交点转换为 MULTILINESTRING，然后再转换为 LINESTRING（因为您从 GEOMETRY 开始，无法轻松转换为 LINESTRING，因此需要中间步骤）。
这样，您将从 intersection 数据集中的 22 行变为我建议的 grid_lines 数据集中的 24 行，并且两个单元格的 ID 和同一行交叉两次重复。
那么这是一个简单的 {dplyr} 聚合和排序案例/我建议在此阶段删除几何图形，因为它不再增加任何值，并在单元 ID 和行 ID 级别工作。
# 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
<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
您可以使用rmapshaper包来'erase'靠近网格的部分线。然后使用 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)
在示例数据中，有两个网格单元被任何一条线交叉多次，用红色圈出：
由 reprex package (v2.0.1)
于 2022-04-22 创建