Как правильно использовать st_difference() в сочетании с mutate()?

У меня есть простой набор данных с линиями, и я пытаюсь найти части линий, которые находятся на расстоянии более 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"))

Найдите способ поделиться некоторыми данными, которые воспроизводят проблему. Используйте команду Rdput, чтобы создать код, который можно вставить в исходные вопросы.

Andrew Chisholm 21.02.2023 09:34

@AndrewChisholm Спасибо за совет. Я добавил несколько примеров данных, которые можно использовать для воспроизведения ошибки.

Jelto 21.02.2023 10:22
Стоит ли изучать PHP в 2023-2024 годах?
Стоит ли изучать PHP в 2023-2024 годах?
Привет всем, сегодня я хочу высказать свои соображения по поводу вопроса, который я уже много раз получал в своем сообществе: "Стоит ли изучать PHP в...
Поведение ключевого слова "this" в стрелочной функции в сравнении с нормальной функцией
Поведение ключевого слова "this" в стрелочной функции в сравнении с нормальной функцией
В JavaScript одним из самых запутанных понятий является поведение ключевого слова "this" в стрелочной и обычной функциях.
Приемы CSS-макетирования - floats и Flexbox
Приемы CSS-макетирования - floats и Flexbox
Здравствуйте, друзья-студенты! Готовы совершенствовать свои навыки веб-дизайна? Сегодня в нашем путешествии мы рассмотрим приемы CSS-верстки - в...
Тестирование функциональных ngrx-эффектов в Angular 16 с помощью Jest
В системе управления состояниями ngrx, совместимой с Angular 16, появились функциональные эффекты. Это здорово и делает код определенно легче для...
Концепция локализации и ее применение в приложениях React ⚡️
Концепция локализации и ее применение в приложениях React ⚡️
Локализация - это процесс адаптации приложения к различным языкам и культурным требованиям. Это позволяет пользователям получить опыт, соответствующий...
Пользовательский скаляр GraphQL
Пользовательский скаляр GraphQL
Листовые узлы системы типов GraphQL называются скалярами. Достигнув скалярного типа, невозможно спуститься дальше по иерархии типов. Скалярный тип...
0
2
50
1
Перейти к ответу Данный вопрос помечен как решенный

Ответы 1

Ответ принят как подходящий

Не уверен в «правильности» или производительности, хотя обработка этого с помощью 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

Другие вопросы по теме