У меня есть простой набор данных с линиями, и я пытаюсь найти части линий, которые находятся на расстоянии более 30 м от соответствующих конечных точек с помощью пакета sf и R. Сначала я пробовал это:
library(tidyverse)
library(sf)
lines = lines |>
mutate(geom = st_cast(x = geometry, to = "LINESTRING"))|> #Make sure all geoms are Linestrings
mutate(edges = st_line_sample(x = geom, sample = c(0,1)))|> #Find Endpoints of every line
mutate(buff = st_buffer(edges, dist = 30)) |> #30m Buffer around Endpoints
mutate(geom = st_difference(x = geom, y = buff)) #Error! Creates more geoms that expected
Однако последняя строка выдает эту ошибку: geom
должен быть размером 897 или 1, а не 804206. Я ожидал, что он выдаст точно такое количество геометрий, так как в наборе данных есть строки. Однако, судя по сообщению об ошибке, st_difference() создает гораздо больше геометрий. Я думаю, что что-то неправильно понимаю в том, как mutate() и st_differnce() работают в комбинации.
Я использовал этот обходной путь:
lines = lines |>
mutate(geom = st_cast(x = geometry, to = "LINESTRING"))|> #Make sure all geoms are Linestrings
mutate(edges = st_line_sample(x = geom, sample = c(0,1)))|> #Find Endpoints of every line
mutate(buff = st_buffer(edges, dist = 30)) #30m Buffer around Endpoints
# mutate(geom = st_difference(x = geom, y = buff)) #Error! Creates more geoms that expected
#Workaround:
hed = st_difference(x = lines #Remove parts of hedges that are
y = st_union(lines$buff)) #within 30m of an Endpoint.
Но это не совсем то, что я ищу, так как он также удаляет части линий, которые находятся рядом с конечными точками других линий, но я хочу удалить только части линий, которые находятся в пределах 30 м от их собственных конечных точек.
Вот несколько примеров данных из моего набора данных для воспроизводимости (извините за беспорядок):
lines = structure(list(id = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1), geometry = structure(list(
structure(c(599534.62516477, 599554.318937559, 5764376.16380116,
5764315.9767656), dim = c(2L, 2L), class = c("XY", "LINESTRING",
"sfg")), structure(c(599590.720979866, 599682.072092792,
5764201.56638851, 5763913.81271676), dim = c(2L, 2L), class = c("XY",
"LINESTRING", "sfg")), structure(c(599685.38979868, 599731.651412595,
5763901.42515664, 5763751.07786744), dim = c(2L, 2L), class = c("XY",
"LINESTRING", "sfg")), structure(c(599464.01445783, 599384.800422916,
599356.793826265, 599342.811078276, 599325.378924622, 599310.023394395,
599290.690077839, 599271.914703978, 599252.532906822, 599229.797215977,
599202.097394155, 5763610.48439324, 5763616.7157656, 5763621.27416296,
5763624.43500734, 5763633.40050881, 5763639.38340642, 5763645.89049208,
5763653.4729713, 5763659.37588095, 5763666.60960097, 5763676.88897741
), dim = c(11L, 2L), class = c("XY", "LINESTRING", "sfg")),
structure(c(599153.617861475, 599146.140179938, 599135.468484187,
599123.38553011, 599112.149996558, 599090.341255877, 599071.816141346,
599049.86653477, 599041.201483651, 599022.964061505, 598980.441861548,
598941.588973157, 598919.370224986, 598897.926361804, 598885.228478353,
598871.067613897, 598859.525982862, 598849.894242691, 598839.611642872,
598822.882502728, 598804.930836576, 598784.897051154, 598768.484115503,
598766.328345305, 5763689.47173725, 5763690.27503453, 5763689.7135454,
5763688.17320693, 5763687.08021086, 5763686.89528389, 5763684.95513276,
5763683.21890848, 5763683.29902017, 5763684.50905542, 5763682.79160655,
5763681.49667385, 5763680.35723431, 5763684.29905714, 5763685.58823485,
5763686.09333662, 5763686.68206958, 5763686.66585263, 5763685.57204623,
5763678.77303556, 5763666.73451416, 5763652.53771155, 5763641.20759758,
5763639.99848025), dim = c(24L, 2L), class = c("XY", "LINESTRING",
"sfg")), structure(c(599227.531845866, 599202.859751073,
599179.521265781, 599155.284117019, 599129.715033211, 599107.414246027,
599092.184531882, 5764647.31843057, 5764734.50086842, 5764816.29444422,
5764899.29136926, 5764987.67428473, 5765065.81724858, 5765121.03303979
), dim = c(7L, 2L), class = c("XY", "LINESTRING", "sfg")),
structure(c(600594.087149005, 600588.916469261, 5764258.86810428,
5764248.29906133), dim = c(2L, 2L), class = c("XY", "LINESTRING",
"sfg")), structure(c(600579.203159313, 600595.259237174,
5764230.32853955, 5764201.37738517), dim = c(2L, 2L), class = c("XY",
"LINESTRING", "sfg")), structure(c(599238.932778222, 599290.882703812,
5764369.81007426, 5764379.22078406), dim = c(2L, 2L), class = c("XY",
"LINESTRING", "sfg")), structure(c(599687.588911078, 599711.323481742,
5763699.23995036, 5763709.54091522), dim = c(2L, 2L), class = c("XY",
"LINESTRING", "sfg"))), class = c("sfc_LINESTRING", "sfc"
), precision = 0, bbox = structure(c(xmin = 598766.328345305,
ymin = 5763610.48439324, xmax = 600595.259237174, ymax = 5765121.03303979
), class = "bbox"), crs = structure(list(input = "ETRS89 / UTM zone 32N",
wkt = "PROJCRS[\"ETRS89 / UTM zone 32N\",\n BASEGEOGCRS[\"ETRS89\",\n ENSEMBLE[\"European Terrestrial Reference System 1989 ensemble\",\n MEMBER[\"European Terrestrial Reference Frame 1989\"],\n MEMBER[\"European Terrestrial Reference Frame 1990\"],\n MEMBER[\"European Terrestrial Reference Frame 1991\"],\n MEMBER[\"European Terrestrial Reference Frame 1992\"],\n MEMBER[\"European Terrestrial Reference Frame 1993\"],\n MEMBER[\"European Terrestrial Reference Frame 1994\"],\n MEMBER[\"European Terrestrial Reference Frame 1996\"],\n MEMBER[\"European Terrestrial Reference Frame 1997\"],\n MEMBER[\"European Terrestrial Reference Frame 2000\"],\n MEMBER[\"European Terrestrial Reference Frame 2005\"],\n MEMBER[\"European Terrestrial Reference Frame 2014\"],\n ELLIPSOID[\"GRS 1980\",6378137,298.257222101,\n LENGTHUNIT[\"metre\",1]],\n ENSEMBLEACCURACY[0.1]],\n PRIMEM[\"Greenwich\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n ID[\"EPSG\",4258]],\n CONVERSION[\"UTM zone 32N\",\n METHOD[\"Transverse Mercator\",\n ID[\"EPSG\",9807]],\n PARAMETER[\"Latitude of natural origin\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433],\n ID[\"EPSG\",8801]],\n PARAMETER[\"Longitude of natural origin\",9,\n ANGLEUNIT[\"degree\",0.0174532925199433],\n ID[\"EPSG\",8802]],\n PARAMETER[\"Scale factor at natural origin\",0.9996,\n SCALEUNIT[\"unity\",1],\n ID[\"EPSG\",8805]],\n PARAMETER[\"False easting\",500000,\n LENGTHUNIT[\"metre\",1],\n ID[\"EPSG\",8806]],\n PARAMETER[\"False northing\",0,\n LENGTHUNIT[\"metre\",1],\n ID[\"EPSG\",8807]]],\n CS[Cartesian,2],\n AXIS[\"(E)\",east,\n ORDER[1],\n LENGTHUNIT[\"metre\",1]],\n AXIS[\"(N)\",north,\n ORDER[2],\n LENGTHUNIT[\"metre\",1]],\n USAGE[\n SCOPE[\"Engineering survey, topographic mapping.\"],\n AREA[\"Europe between 6°E and 12°E: Austria; Belgium; Denmark - onshore and offshore; Germany - onshore and offshore; Norway including - onshore and offshore; Spain - offshore.\"],\n BBOX[38.76,6,84.33,12]],\n ID[\"EPSG\",25832]]"), class = "crs"), n_empty = 0L)), row.names = c(NA,
-10L), sf_column = "geometry", agr = structure(c(id = NA_integer_), levels = c("constant",
"aggregate", "identity"), class = "factor"), class = c("sf",
"tbl_df", "tbl", "data.frame")) "tbl_df", "tbl", "data.frame"))
@AndrewChisholm Спасибо за совет. Я добавил несколько примеров данных, которые можно использовать для воспроизведения ошибки.
Не уверен в «правильности» или производительности, хотя обработка этого с помощью rowwise()
, похоже, работает. Разорвал трубу только для того, чтобы собрать набор данных для визуализации; последние строки на сером фоне.
library(sf)
#> Linking to GEOS 3.9.3, GDAL 3.5.2, PROJ 8.2.1; sf_use_s2() is TRUE
library(dplyr)
library(tidyr)
library(ggplot2)
lines2 <- lines |>
mutate(row_n = as.factor(row_number()), .before = 1) |> # add row number for visualization
mutate(geom = st_cast(x = geometry, to = "LINESTRING"))|> # Make sure all geoms are Linestrings
mutate(edges = st_line_sample(x = geom, sample = c(0,1)))|> # Find Endpoints of every line
mutate(buff = st_buffer(edges, dist = 30)) # 30m Buffer around Endpoints
lines3 <- lines2 |> rowwise() |>
# as some lines are completely covered by buffers there will be some 0-lenght results
# and mutate is not happy with those; creating a list column ...
mutate(geom = list(st_difference(x = geom, y = buff))) |>
ungroup() |>
# ... and unnesting it drops those records
unnest(geom) |>
mutate(length_1 = st_length(geometry),
length_2 = st_length(geom))
# leghts before and after st_difference()
lines3 |> select(row_n, geom, length_1, length_2)
#> Simple feature collection with 6 features and 3 fields
#> Active geometry column: geometry
#> Geometry type: LINESTRING
#> Dimension: XY
#> Bounding box: xmin: 598766.3 ymin: 5763610 xmax: 599731.7 ymax: 5765121
#> Projected CRS: ETRS89 / UTM zone 32N
#> # A tibble: 6 × 5
#> row_n geom lengt…¹ lengt…² geometry
#> <fct> <LINESTRING [m]> [m] [m] <LINESTRING [m]>
#> 1 1 (599544 5764348, 599545 57643… 63.3 3.33 (599534.6 5764376, 59955…
#> 2 2 (599599.8 5764173, 599673 576… 302. 242. (599590.7 5764202, 59968…
#> 3 3 (599694.2 5763873, 599722.8 5… 157. 97.3 (599685.4 5763901, 59973…
#> 4 4 (599434.1 5763613, 599384.8 5… 273. 213. (599464 5763610, 599384.…
#> 5 5 (599123.7 5763688, 599123.4 5… 402. 342. (599153.6 5763689, 59914…
#> 6 6 (599219.4 5764676, 599202.9 5… 493. 433. (599227.5 5764647, 59920…
#> # … with abbreviated variable names ¹length_1, ²length_2
ggplot(lines2, aes(color = row_n)) +
geom_sf(aes(geometry = buff), alpha = .2) +
geom_sf(aes(geometry = edges), alpha = .2) +
# resulting lines in grey
geom_sf(data = lines3, aes(geometry = geom), linewidth = 3, color = "grey", alpha = .7) +
geom_sf(aes(geometry = geom)) +
geom_sf_text(aes(geometry = geom, label = row_n), color = "black", nudge_x = -50, nudge_y = 40) +
theme_bw()
Входные данные:
lines <- structure(list(id = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1), geometry = structure(list(
structure(c(
599534.62516477, 599554.318937559, 5764376.16380116,
5764315.9767656
), dim = c(2L, 2L), class = c(
"XY", "LINESTRING",
"sfg"
)), structure(c(
599590.720979866, 599682.072092792,
5764201.56638851, 5763913.81271676
), dim = c(2L, 2L), class = c(
"XY",
"LINESTRING", "sfg"
)), structure(c(
599685.38979868, 599731.651412595,
5763901.42515664, 5763751.07786744
), dim = c(2L, 2L), class = c(
"XY",
"LINESTRING", "sfg"
)), structure(c(
599464.01445783, 599384.800422916,
599356.793826265, 599342.811078276, 599325.378924622, 599310.023394395,
599290.690077839, 599271.914703978, 599252.532906822, 599229.797215977,
599202.097394155, 5763610.48439324, 5763616.7157656, 5763621.27416296,
5763624.43500734, 5763633.40050881, 5763639.38340642, 5763645.89049208,
5763653.4729713, 5763659.37588095, 5763666.60960097, 5763676.88897741
), dim = c(11L, 2L), class = c("XY", "LINESTRING", "sfg")),
structure(c(
599153.617861475, 599146.140179938, 599135.468484187,
599123.38553011, 599112.149996558, 599090.341255877, 599071.816141346,
599049.86653477, 599041.201483651, 599022.964061505, 598980.441861548,
598941.588973157, 598919.370224986, 598897.926361804, 598885.228478353,
598871.067613897, 598859.525982862, 598849.894242691, 598839.611642872,
598822.882502728, 598804.930836576, 598784.897051154, 598768.484115503,
598766.328345305, 5763689.47173725, 5763690.27503453, 5763689.7135454,
5763688.17320693, 5763687.08021086, 5763686.89528389, 5763684.95513276,
5763683.21890848, 5763683.29902017, 5763684.50905542, 5763682.79160655,
5763681.49667385, 5763680.35723431, 5763684.29905714, 5763685.58823485,
5763686.09333662, 5763686.68206958, 5763686.66585263, 5763685.57204623,
5763678.77303556, 5763666.73451416, 5763652.53771155, 5763641.20759758,
5763639.99848025
), dim = c(24L, 2L), class = c(
"XY", "LINESTRING",
"sfg"
)), structure(c(
599227.531845866, 599202.859751073,
599179.521265781, 599155.284117019, 599129.715033211, 599107.414246027,
599092.184531882, 5764647.31843057, 5764734.50086842, 5764816.29444422,
5764899.29136926, 5764987.67428473, 5765065.81724858, 5765121.03303979
), dim = c(7L, 2L), class = c("XY", "LINESTRING", "sfg")),
structure(c(
600594.087149005, 600588.916469261, 5764258.86810428,
5764248.29906133
), dim = c(2L, 2L), class = c(
"XY", "LINESTRING",
"sfg"
)), structure(c(
600579.203159313, 600595.259237174,
5764230.32853955, 5764201.37738517
), dim = c(2L, 2L), class = c(
"XY",
"LINESTRING", "sfg"
)), structure(c(
599238.932778222, 599290.882703812,
5764369.81007426, 5764379.22078406
), dim = c(2L, 2L), class = c(
"XY",
"LINESTRING", "sfg"
)), structure(c(
599687.588911078, 599711.323481742,
5763699.23995036, 5763709.54091522
), dim = c(2L, 2L), class = c(
"XY",
"LINESTRING", "sfg"
))
), class = c("sfc_LINESTRING", "sfc"), precision = 0, bbox = structure(c(
xmin = 598766.328345305,
ymin = 5763610.48439324, xmax = 600595.259237174, ymax = 5765121.03303979
), class = "bbox"), crs = structure(list(
input = "ETRS89 / UTM zone 32N",
wkt = "PROJCRS[\"ETRS89 / UTM zone 32N\",\n BASEGEOGCRS[\"ETRS89\",\n ENSEMBLE[\"European Terrestrial Reference System 1989 ensemble\",\n MEMBER[\"European Terrestrial Reference Frame 1989\"],\n MEMBER[\"European Terrestrial Reference Frame 1990\"],\n MEMBER[\"European Terrestrial Reference Frame 1991\"],\n MEMBER[\"European Terrestrial Reference Frame 1992\"],\n MEMBER[\"European Terrestrial Reference Frame 1993\"],\n MEMBER[\"European Terrestrial Reference Frame 1994\"],\n MEMBER[\"European Terrestrial Reference Frame 1996\"],\n MEMBER[\"European Terrestrial Reference Frame 1997\"],\n MEMBER[\"European Terrestrial Reference Frame 2000\"],\n MEMBER[\"European Terrestrial Reference Frame 2005\"],\n MEMBER[\"European Terrestrial Reference Frame 2014\"],\n ELLIPSOID[\"GRS 1980\",6378137,298.257222101,\n LENGTHUNIT[\"metre\",1]],\n ENSEMBLEACCURACY[0.1]],\n PRIMEM[\"Greenwich\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n ID[\"EPSG\",4258]],\n CONVERSION[\"UTM zone 32N\",\n METHOD[\"Transverse Mercator\",\n ID[\"EPSG\",9807]],\n PARAMETER[\"Latitude of natural origin\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433],\n ID[\"EPSG\",8801]],\n PARAMETER[\"Longitude of natural origin\",9,\n ANGLEUNIT[\"degree\",0.0174532925199433],\n ID[\"EPSG\",8802]],\n PARAMETER[\"Scale factor at natural origin\",0.9996,\n SCALEUNIT[\"unity\",1],\n ID[\"EPSG\",8805]],\n PARAMETER[\"False easting\",500000,\n LENGTHUNIT[\"metre\",1],\n ID[\"EPSG\",8806]],\n PARAMETER[\"False northing\",0,\n LENGTHUNIT[\"metre\",1],\n ID[\"EPSG\",8807]]],\n CS[Cartesian,2],\n AXIS[\"(E)\",east,\n ORDER[1],\n LENGTHUNIT[\"metre\",1]],\n AXIS[\"(N)\",north,\n ORDER[2],\n LENGTHUNIT[\"metre\",1]],\n USAGE[\n SCOPE[\"Engineering survey, topographic mapping.\"],\n AREA[\"Europe between 6°E and 12°E: Austria; Belgium; Denmark - onshore and offshore; Germany - onshore and offshore; Norway including - onshore and offshore; Spain - offshore.\"],\n BBOX[38.76,6,84.33,12]],\n ID[\"EPSG\",25832]]"
), class = "crs"), n_empty = 0L)), row.names = c(
NA,
-10L
), sf_column = "geometry", agr = structure(c(id = NA_integer_), levels = c(
"constant",
"aggregate", "identity"
), class = "factor"), class = c(
"sf",
"tbl_df", "tbl", "data.frame"
))
Created on 2023-02-21 with reprex v2.0.2
Найдите способ поделиться некоторыми данными, которые воспроизводят проблему. Используйте команду
R
dput
, чтобы создать код, который можно вставить в исходные вопросы.