如何正确使用st_difference()和mutate()?

1mrurvl1  于 2023-02-27  发布在  其他
关注(0)|答案(1)|浏览(118)

我有一个包含线的简单要素数据集,我尝试使用sf包和R查找距离各自端点超过30m的线部分。这是我最初尝试的:

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"))
kgsdhlau

kgsdhlau1#

虽然通过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"
))

创建于2023年2月21日,使用reprex v2.0.2

相关问题