R GIS相交线和多边形（大数据）答案

【问题标题】：R GIS intersecting lines and polygons (large data)R GIS相交线和多边形（大数据）
【发布时间】：2021-11-17 20:04:26
【问题描述】：

GIS 和 R 问题：

我需要将许多线与许多多边形相交。

理想情况下，我想在 R 中执行此操作。

当一条线穿过多边形边界时，我希望分割那条线（制作两条特征/线）。

然后，我想将新的（通常更小的）线与它们现在所在的多边形的属性相结合。

如果线条位于多边形之外，我希望从多边形附加到线条的新属性将为 NULL 或 NA 或空。

我很高兴将数据加载为 shapefiles 或 geopackages。

在 R 中执行此操作的最佳工具是什么？

（我很乐意使用 Raster、rgeos 和 SF 等包。但是，理想情况下，如果我们可以避免调用QGIS 或其他软件，除非它非常稳定，否则会很棒。我没有 Arc 的许可证）

非常感谢。

【问题讨论】：

您已经有了答案 - 使用 sf 包及其各种几何函数。制作一个可能有十个不同特征的小数据集，然后弄清楚如何去做。如果您仍然不能，请让我们访问该数据集并向我们展示您迄今为止的工作。
好的，谢谢，我试试看。我希望它比需要数小时（甚至数天）处理的 QGIS 更快。

标签： r vector gis polygon polyline

【解决方案1】：

让我们收集一些在线数据，一个多边形和几条线。我们将使用osmdata 为我的村庄获取OSM 数据。

library(osmdata)
#> Data (c) OpenStreetMap contributors, ODbL 1.0. https://www.openstreetmap.org/copyright
library(sf)
#> Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1; sf_use_s2() is TRUE
library(tidyverse)

让我们找到一个边界框并为其添加一点空间

Lbb <- getbb("Lubnów trzebnicki Poland")
addM <- matrix(data = c(-0.01, -0.01, 0.01, 0.01), nrow = 2, ncol = 2)
newBB <- Lbb + addM

让我们获取我们将在计算中使用的高速公路（线、多线）和边界（多边形、多多边形）

highways <- opq (newBB) |> 
  add_osm_feature (key = "highway") |>
  osmdata_sf()

highways <- highways$osm_lines |>
  select(osm_id, name, geometry) |>
  st_as_sf()

boundaries <- opq(newBB) |>
  add_osm_feature (key = "boundary", value = "administrative") |>
  osmdata_sf()

boundariesWsi <- boundaries$osm_multipolygons |>
  filter(boundary == "administrative" & admin_level == 8 & name == "Lubnów")

让我们为可视化绘制它

plot(boundariesWsi$geometry, col = "gray90")
plot(highways$geometry, lwd = 0.6, lty = 4, add = TRUE)

我们将创建一个对应于村庄边界的支持MULTILINESTRING，我们将使用它来查找那些跨越边界的路径

boundariesWsiAsLines <- st_cast(boundariesWsi, "MULTILINESTRING")

现在，让我们找到那些跨越村庄边界的interestingHighways：

interestingHighways <- highways |> 
  filter(st_intersects(boundariesWsiAsLines$geometry, highways$geometry, sparse = FALSE))

最后，找到这些高速公路的 IN 和 OUT 部分：

a<- st_as_sf(st_intersection(interestingHighways$geometry, boundariesWsi$geometry))
b<- st_as_sf(st_difference(interestingHighways$geometry, boundariesWsi$geometry))

plot(a, lwd = 4, col = "darkblue", add = TRUE)
plot(b, lwd = 4, col = "darkgreen", add = TRUE)

您可以拥有一个多面体列表，而不是一个多面体。要将属性分配给方式，您可以使用dplyr::mutate()。

^{由reprex package 创建于 2022-01-31 (v2.0.1)}

【讨论】：