Есть ли способ полностью привязать строку к многоугольнику в R?
У меня есть два файла SHP: один со множеством линий (железные дороги) и один со множеством полигонов (округа в США). Между этими двумя файлами вообще нет пересечений. Моя цель — привязать линии к ближайшей границе многоугольника, чтобы они пересекались и вычисляли длину линий в каждом многоугольнике.
Я пытался сделать это, используя пакет st_snap из sf, но не получил ожидаемых результатов.
Это побудило меня изучить эту функцию на некоторых простых данных и получить следующие результаты:
library(dplyr)
library(tidyverse)
library(sf)
library(units)
library(ggplot2)
poly = st_polygon(list(cbind(c(0, 0, 1, 1, 0), c(0, 1, 1, 0, 0))))
lines = st_multilinestring(list(
cbind(c(0, 1), c(1.01, 1.05)),
cbind(c(0, 1), c(-0.01, -.05)),
cbind(c(1.025, 1.05, 1.025), c(1.05, .5, -.05))
))
t <- st_intersection(poly, lines)
st_is_empty(t)
# [1] TRUE
ggplot()+
geom_sf(data = lines, col = "green")+
geom_sf(data = poly, col = "red", fill = NA)
Мой ожидаемый результат после привязки будет:
Но что я получаю:
для tolerance=0.5
snapped1 = st_snap(lines, poly, tolerance=0.5)
print(snapped1)
# MULTILINESTRING ((0 1, 1 1), (0 0, 1 0), (1 1, 1.05 0.5, 1 0))
ggplot()+
geom_sf(data = lines, col = "green")+
geom_sf(data = poly, col = "red", fill = NA)+
geom_sf(data = snapped1, col = "blue", alpha = 0.5)
для tolerance=1.005
snapped2 = st_snap(lines, poly, tolerance=1.005)
print(snapped2)
# MULTILINESTRING ((1 1, 0 1, 0 0, 1 0), (0 1, 1 1, 1 0, 0 0), (0 1, 1 1, 1.05 0.5, 1 0, 0 0))
Который содержит ожидаемый результат: (0 1, 1 1, 1 0, 0 0), но и нежелательные линии.
ggplot()+
geom_sf(data = lines, col = "green")+
geom_sf(data = poly, col = "red", fill = NA)+
geom_sf(data = snapped2, col = "blue", alpha = 0.5)
для tolerance=1.1
snapped3 = st_snap(lines, poly, tolerance=1.1)
print(snapped3)
# MULTILINESTRING ((1 1, 0 1, 0 0, 1 0), (0 1, 1 1, 1 0, 0 0), (1 1, 1 0, 0 0, 0 1))
Который также содержит ожидаемый результат: (0 1, 1 1, 1 0, 0 0), но также и нежелательные линии.
ggplot()+
geom_sf(data = lines, col = "green")+
geom_sf(data = poly, col = "red", fill = NA)+
geom_sf(data = snapped3, col = "blue", alpha = 0.5)
При любом большем допуске результат тот же, что и в третьем. Очевидно, я чего-то не понимаю в том, как работает эта функция. Я попытался просмотреть документацию, но не смог найти ничего полезного. Я был бы признателен за любую информацию об этом и о том, как получить ожидаемые результаты с помощью этой функции или какой-либо другой функции.
Спасибо.
@margusl Спасибо за помощь! Извините за отсутствие ясности. Мой первоначальный набор данных включал все железные дороги США и всех округов, в результате чего они пересекались. Однако из-за различных проблем, особенно связанных с прибрежными границами, некоторые железные дороги остались не связанными ни с одним округом. Следовательно, я создал шейп-файл, содержащий эти изолированные железные дороги. Мой подход к решению этой проблемы заключался в том, чтобы соединить эти железные дороги с ближайшим округом. Есть ли способ добиться этого без изменения длины?
Спасибо, получил его! Я предполагаю, что вы хотите, чтобы линии пересекались с многоугольниками, чтобы вы могли использовать пространственное соединение для присвоения названий округов линиям/сегментам, хотя st_join() также работает с довольно большим количеством других предикатов, помимо st_intersects по умолчанию, я бы, по крайней мере, проверил st_nearest_feature. Добавил еще один пример в свой ответ.
Ответы 1
st_snapпривязывает вершины и сегменты геометрии к вершинам другой геометрии.
[?sf::st_snap]
Если ваши полигоны имеют длинные (почти) прямые участки, вам, скорее всего, потребуется добавить больше точек (вершин), чтобы вершинам и сегментам линейных объектов было к чему привязываться. И может быть довольно сложно найти оптимальный баланс между допуском привязки и требуемой плотностью точек.
С помощью этого репрекса st_segmentize(..., dfMaxLength = .2) и st_snap(..., tolerance = .18) (т.е. чуть меньше dfMaxLength, чтобы избежать неправильных снимков) должны делать:
library(sf)
#> Linking to GEOS 3.11.2, GDAL 3.6.2, PROJ 9.2.0; sf_use_s2() is TRUE
library(ggplot2)
library(dplyr, warn.conflicts = FALSE)
poly <- st_polygon(list(cbind(c(0, 0, 1, 1, 0), c(0, 1, 1, 0, 0))))
st_as_text(poly)
#> [1] "POLYGON ((0 0, 0 1, 1 1, 1 0, 0 0))"
# add points to straight lines
poly_2 <- st_segmentize(poly, dfMaxLength = .2)
st_as_text(poly_2)
#> [1] "POLYGON ((0 0, 0 0.2, 0 0.4, 0 0.6, 0 0.8, 0 1, 0.2 1, 0.4 1, 0.6 1, 0.8 1, 1 1, 1 0.8, 1 0.6, 1 0.4, 1 0.2, 1 0, 0.8 0, 0.6 0, 0.4 0, 0.2 0, 0 0))"
lines <- st_multilinestring(list(
cbind(c(0, 1), c(1.01, 1.05)),
cbind(c(0, 1), c(-0.01, -.05)),
cbind(c(1.025, 1.05, 1.025), c(1.05, .5, -.05))
))
snapped <- st_snap(lines, poly_2, tolerance = .18)
st_as_text(snapped)
#> [1] "MULTILINESTRING ((0 1, 0.2 1, 0.4 1, 0.6 1, 0.8 1, 1 1), (0 0, 0.2 0, 0.4 0, 0.6 0, 0.8 0, 1 0), (1 1, 1 0.8, 1 0.6, 1 0.4, 1 0.2, 1 0))"
ggplot() +
geom_sf(data = st_cast(poly_2,"LINESTRING"), aes(color = "poly"), fill = NA, linewidth = 3) +
geom_sf(data = st_cast(poly_2, "MULTIPOINT"), aes(color = "poly"), size = 6) +
geom_sf(data = lines, aes(color = "lines"), linewidth = 2, alpha = .8) +
geom_sf(data = st_cast(lines, "MULTIPOINT"), aes(color = "lines"), size = 4, alpha = .8) +
geom_sf(data = snapped, aes(color = "snapped"), linewidth = 1, alpha = .8) +
geom_sf(data = st_cast(snapped, "MULTIPOINT"), aes(color = "snapped"), size = 2, alpha = .8) +
scale_color_viridis_d() +
theme_minimal()
-e- Пространственное соединение с предикатом ближайшего объекта
В зависимости от того, как организованы объекты в наборе данных железнодорожной сети, пространственное соединение с предикатом ближайшего объекта ( st_join(... , join = st_nearest_feature) ) может быть лучшим подходом к решению этой исходной проблемы. Однако обратите внимание на несколько неоднозначные случаи, такие как B&D в следующем примере.
poly_sf <-
st_make_grid(poly, n = 2) |>
st_sf(geometry = _) |>
mutate(name = LETTERS[row_number()] |> as.factor())
poly_sf
#> Simple feature collection with 4 features and 1 field
#> Geometry type: POLYGON
#> Dimension: XY
#> Bounding box: xmin: 0 ymin: 0 xmax: 1 ymax: 1
#> CRS: NA
#> geometry name
#> 1 POLYGON ((0 0, 0.5 0, 0.5 0... A
#> 2 POLYGON ((0.5 0, 1 0, 1 0.5... B
#> 3 POLYGON ((0 0.5, 0.5 0.5, 0... C
#> 4 POLYGON ((0.5 0.5, 1 0.5, 1... D
lines_sf <-
st_sfc(lines) |>
st_sf(geometry = _) |>
st_cast("LINESTRING")
lines_sf
#> Simple feature collection with 3 features and 0 fields
#> Geometry type: LINESTRING
#> Dimension: XY
#> Bounding box: xmin: 0 ymin: -0.05 xmax: 1.05 ymax: 1.05
#> CRS: NA
#> geometry
#> 1 LINESTRING (0 1.01, 1 1.05)
#> 2 LINESTRING (0 -0.01, 1 -0.05)
#> 3 LINESTRING (1.025 1.05, 1.0...
lines_w_names <-
st_join(lines_sf, poly_sf, join = st_nearest_feature)
lines_w_names
#> Simple feature collection with 3 features and 1 field
#> Geometry type: LINESTRING
#> Dimension: XY
#> Bounding box: xmin: 0 ymin: -0.05 xmax: 1.05 ymax: 1.05
#> CRS: NA
#> name geometry
#> 1 C LINESTRING (0 1.01, 1 1.05)
#> 2 A LINESTRING (0 -0.01, 1 -0.05)
#> 3 D LINESTRING (1.025 1.05, 1.0...
ggplot() +
geom_sf(data = poly_sf, aes(fill = name), show.legend = FALSE) +
geom_sf_label(data = poly_sf, aes(label = name), alpha = .8) +
geom_sf(data = lines_w_names, aes(color = name), show.legend = FALSE, linewidth = 2) +
geom_sf_label(data = lines_w_names, aes(label = name), alpha = .8) +
scale_fill_viridis_d() +
scale_color_viridis_d() +
theme_minimal()
Created on 2024-04-21 with reprex v2.1.0
Кстати, я не совсем уверен, поможет ли ответ, который вы ищете, в решении вашей первоначальной проблемы (длина железных дорог в каждом округе). Я также не до конца это понимаю - я бы предположил, что топологически действительные полигоны округов не изолированы, не имеют перекрытий и полос, поэтому часть «без пересечений» (между фигурами округов и железнодорожной сети?) немного сбивает с толку; при привязке железной дороги к границе между округами A и B она будет пересекаться как с A, так и с B; Если длина сегмента рельса не взята из таблицы атрибутов или не измерена перед изменением, изменение формы линии также приводит к изменению длины.