如何按点分割“MULTILINESTRING”sf对象？

问题描述 投票：0回答：1

我已将连接线分组以形成一个连续的 

MULTILINESTRING
对象，现在想根据点将其分成更小的组。我找到了几个示例，说明如何使用 
LINESTRING
几何图形执行此操作，但不适用于 
MULTILINESTRING
几何图形。请参阅此处和此处。使用这两种方法中的任何一种，我都会得到一个包含 
LINESTRING
的 sf 对象。当然，我可以根据我的小组专栏将这些聚合回 
MULTILINESTRING
，但这只能让我回到开始的地方。我需要将相交点每一侧的所有线分组为具有唯一组号的
MULTILINESTRING
。我希望以最有效的方式进行这种分割，因为我的现实数据集非常大。请参阅下面我的示例数据集。

library(ggplot2)
library(sf)
        testArea<-structure(list(group = 3, geom = structure(list(structure(list(
    structure(c(6502476.72250405, 6502454.02996413, 1962622.20532215, 
    1962613.17515647), dim = c(2L, 2L)), structure(c(6502476.72250405, 
    6502453.90463631, 1962622.20532215, 1962632.57964523), dim = c(2L, 
    2L)), structure(c(6502366.41465381, 6502371.44253089, 6502473.30682847, 
    6502476.72250405, 1962490.68032272, 1962495.74363281, 1962498.12158081, 
    1962622.20532215), dim = c(4L, 2L)), structure(c(6502228.73645148, 
    6502211.24928172, 1962628.79028273, 1962643.49694623), dim = c(2L, 
    2L)), structure(c(6502213.34114106, 6502223.08324756, 1962491.00643755, 
    1962467.08522557), dim = c(2L, 2L)), structure(c(6502228.73645148, 
    6502250.32728755, 6502307.72317013, 6502325.46394831, 1962628.79028273, 
    1962625.21515864, 1962624.4385854, 1962633.67150655), dim = c(4L, 
    2L)), structure(c(6502671.5715238, 6502661.95999447, 1962487.84240189, 
    1962462.80997165), dim = c(2L, 2L)), structure(c(6502476.72250405, 
    6502494.79891147, 6502552.95102614, 6502581.00116688, 1962622.20532215, 
    1962617.1669464, 1962616.94680248, 1962626.03963207), dim = c(4L, 
    2L)), structure(c(6502228.73645148, 6502210.75486013, 1962628.79028273, 
    1962623.37789197), dim = c(2L, 2L)), structure(c(6502366.41465381, 
    6502376.15413564, 6502474.69757372, 6502493.52365156, 6502542.63214913, 
    6502548.52485389, 1962490.68032272, 1962480.89097223, 1962483.04615164, 
    1962483.45789623, 1962483.56846032, 1962488.53531389), dim = c(6L, 
    2L)), structure(c(6502366.41465381, 6502367.14956047, 6502405.43065189, 
    1962490.68032272, 1962469.13115323, 1962446.99602689), dim = 3:2), 
    structure(c(6502366.41465381, 6502306.02074572, 1962490.68032272, 
    1962490.85880005), dim = c(2L, 2L)), structure(c(6502306.02074572, 
    6502331.86321372, 1962490.85880005, 1962451.68335348), dim = c(2L, 
    2L)), structure(c(6502548.52485389, 6502671.5715238, 1962488.53531389, 
    1962487.84240189), dim = c(2L, 2L)), structure(c(6502476.72250405, 
    6502493.10337681, 6502550.52878688, 6502569.00578405, 1962622.20532215, 
    1962611.71616989, 1962610.93959664, 1962603.28344397), dim = c(4L, 
    2L)), structure(c(6502213.34114106, 6502180.82119296, 1962491.00643755, 
    1962465.7459894), dim = c(2L, 2L)), structure(c(6502228.73645148, 
    6502245.1176523, 6502303.41609214, 6502326.35075755, 1962628.79028273, 
    1962618.30080239, 1962617.95959572, 1962614.20074497), dim = c(4L, 
    2L)), structure(c(6502366.41465381, 6502361.39957197, 6502228.35587481, 
    6502228.73645148, 1962490.68032272, 1962495.69540456, 1962496.08812031, 
    1962628.79028273), dim = c(4L, 2L)), structure(c(6502548.52485389, 
    6502569.88406314, 1962488.53531389, 1962462.09999931), dim = c(2L, 
    2L)), structure(c(6502366.41465381, 6502376.00387347, 6502452.04538806, 
    1962490.68032272, 1962471.47498056, 1962457.93334098), dim = 3:2), 
    structure(c(6502306.02074572, 6502294.69202822, 1962490.85880005, 
    1962453.50651257), dim = c(2L, 2L)), structure(c(6502306.02074572, 
    6502213.34114106, 1962490.85880005, 1962491.00643755), dim = c(2L, 
    2L)), structure(c(6502671.5715238, 6502693.53867146, 1962487.84240189, 
    1962463.86377531), dim = c(2L, 2L)), structure(c(6502548.52485389, 
    6502538.80833788, 1962488.53531389, 1962464.75189689), dim = c(2L, 
    2L)), structure(c(6502693.53867146, 6502693.05278005, 1962463.86377531, 
    1962448.23880656), dim = c(2L, 2L)), structure(c(6502507.63714039, 
    6502497.52331547, 6502366.41465381, 1962494.67375307, 1962493.54777107, 
    1962490.68032272), dim = 3:2), structure(c(6502507.63714039, 
    6502517.90418021, 6502555.63835672, 6502736.0713948, 6502735.87224822, 
    6502735.84239264, 1962494.67375307, 1962493.70754765, 1962493.49527773, 
    1962492.47887556, 1962559.16378155, 1962569.20477198), dim = c(6L, 
    2L))), class = c("XY", "MULTILINESTRING", "sfg"))), class = c("sfc_MULTILINESTRING", 
"sfc"), precision = 0, bbox = structure(c(xmin = 6502180.82119296, 
ymin = 1962446.99602689, xmax = 6502736.0713948, ymax = 1962643.49694623
), class = "bbox"), crs = structure(list(input = "PROJCS[\"NAD_1983_StatePlane_California_III_FIPS_0403_Feet\",GEOGCS[\"GCS_North_American_1983\",DATUM[\"D_North_American_1983\",SPHEROID[\"GRS_1980\",6378137.0,298.257222101]],PRIMEM[\"Greenwich\",0.0],UNIT[\"Degree\",0.0174532925199433]],PROJECTION[\"Lambert_Conformal_Conic\"],PARAMETER[\"False_Easting\",6561666.666666666],PARAMETER[\"False_Northing\",1640416.666666667],PARAMETER[\"Central_Meridian\",-120.5],PARAMETER[\"Standard_Parallel_1\",37.06666666666667],PARAMETER[\"Standard_Parallel_2\",38.43333333333333],PARAMETER[\"Latitude_Of_Origin\",36.5],UNIT[\"Foot_US\",0.3048006096012192]]", 
    wkt = "PROJCRS[\"NAD83 / California zone 3 (ftUS)\",\n    BASEGEOGCRS[\"NAD83\",\n        DATUM[\"North American Datum 1983\",\n            ELLIPSOID[\"GRS 1980\",6378137,298.257222101,\n                LENGTHUNIT[\"metre\",1]],\n            ID[\"EPSG\",6269]],\n        PRIMEM[\"Greenwich\",0,\n            ANGLEUNIT[\"Degree\",0.0174532925199433]]],\n    CONVERSION[\"unnamed\",\n        METHOD[\"Lambert Conic Conformal (2SP)\",\n            ID[\"EPSG\",9802]],\n        PARAMETER[\"Latitude of false origin\",36.5,\n            ANGLEUNIT[\"Degree\",0.0174532925199433],\n            ID[\"EPSG\",8821]],\n        PARAMETER[\"Longitude of false origin\",-120.5,\n            ANGLEUNIT[\"Degree\",0.0174532925199433],\n            ID[\"EPSG\",8822]],\n        PARAMETER[\"Latitude of 1st standard parallel\",37.0666666666667,\n            ANGLEUNIT[\"Degree\",0.0174532925199433],\n            ID[\"EPSG\",8823]],\n        PARAMETER[\"Latitude of 2nd standard parallel\",38.4333333333333,\n            ANGLEUNIT[\"Degree\",0.0174532925199433],\n            ID[\"EPSG\",8824]],\n        PARAMETER[\"Easting at false origin\",6561666.66666667,\n            LENGTHUNIT[\"US survey foot\",0.304800609601219],\n            ID[\"EPSG\",8826]],\n        PARAMETER[\"Northing at false origin\",1640416.66666667,\n            LENGTHUNIT[\"US survey foot\",0.304800609601219],\n            ID[\"EPSG\",8827]]],\n    CS[Cartesian,2],\n        AXIS[\"(E)\",east,\n            ORDER[1],\n            LENGTHUNIT[\"US survey foot\",0.304800609601219,\n                ID[\"EPSG\",9003]]],\n        AXIS[\"(N)\",north,\n            ORDER[2],\n            LENGTHUNIT[\"US survey foot\",0.304800609601219,\n                ID[\"EPSG\",9003]]]]"), class = "crs"), n_empty = 0L)), row.names = c(NA, 
-1L), sf_column = "geom", agr = structure(c(group = NA_integer_), levels = c("constant", 
"aggregate", "identity"), class = "factor"), class = c("sf", 
"tbl_df", "tbl", "data.frame"))
    
split_point<-structure(list(PhasingCode = 4L, geom = structure(list(structure(c(6502366.41465381, 
1962490.68032272), class = c("XY", "POINT", "sfg"))), class = c("sfc_POINT", 
"sfc"), precision = 0, bbox = structure(c(xmin = 6502366.41465381, 
ymin = 1962490.68032272, xmax = 6502366.41465381, ymax = 1962490.68032272
), class = "bbox"), crs = structure(list(input = "PROJCS[\"NAD_1983_StatePlane_California_III_FIPS_0403_Feet\",GEOGCS[\"GCS_North_American_1983\",DATUM[\"D_North_American_1983\",SPHEROID[\"GRS_1980\",6378137.0,298.257222101]],PRIMEM[\"Greenwich\",0.0],UNIT[\"Degree\",0.0174532925199433]],PROJECTION[\"Lambert_Conformal_Conic\"],PARAMETER[\"False_Easting\",6561666.666666666],PARAMETER[\"False_Northing\",1640416.666666667],PARAMETER[\"Central_Meridian\",-120.5],PARAMETER[\"Standard_Parallel_1\",37.06666666666667],PARAMETER[\"Standard_Parallel_2\",38.43333333333333],PARAMETER[\"Latitude_Of_Origin\",36.5],UNIT[\"Foot_US\",0.3048006096012192]]", 
    wkt = "PROJCRS[\"NAD83 / California zone 3 (ftUS)\",\n    BASEGEOGCRS[\"NAD83\",\n        DATUM[\"North American Datum 1983\",\n            ELLIPSOID[\"GRS 1980\",6378137,298.257222101,\n                LENGTHUNIT[\"metre\",1]],\n            ID[\"EPSG\",6269]],\n        PRIMEM[\"Greenwich\",0,\n            ANGLEUNIT[\"Degree\",0.0174532925199433]]],\n    CONVERSION[\"unnamed\",\n        METHOD[\"Lambert Conic Conformal (2SP)\",\n            ID[\"EPSG\",9802]],\n        PARAMETER[\"Latitude of false origin\",36.5,\n            ANGLEUNIT[\"Degree\",0.0174532925199433],\n            ID[\"EPSG\",8821]],\n        PARAMETER[\"Longitude of false origin\",-120.5,\n            ANGLEUNIT[\"Degree\",0.0174532925199433],\n            ID[\"EPSG\",8822]],\n        PARAMETER[\"Latitude of 1st standard parallel\",37.0666666666667,\n            ANGLEUNIT[\"Degree\",0.0174532925199433],\n            ID[\"EPSG\",8823]],\n        PARAMETER[\"Latitude of 2nd standard parallel\",38.4333333333333,\n            ANGLEUNIT[\"Degree\",0.0174532925199433],\n            ID[\"EPSG\",8824]],\n        PARAMETER[\"Easting at false origin\",6561666.66666667,\n            LENGTHUNIT[\"US survey foot\",0.304800609601219],\n            ID[\"EPSG\",8826]],\n        PARAMETER[\"Northing at false origin\",1640416.66666667,\n            LENGTHUNIT[\"US survey foot\",0.304800609601219],\n            ID[\"EPSG\",8827]]],\n    CS[Cartesian,2],\n        AXIS[\"(E)\",east,\n            ORDER[1],\n            LENGTHUNIT[\"US survey foot\",0.304800609601219,\n                ID[\"EPSG\",9003]]],\n        AXIS[\"(N)\",north,\n            ORDER[2],\n            LENGTHUNIT[\"US survey foot\",0.304800609601219,\n                ID[\"EPSG\",9003]]]]"), class = "crs"), n_empty = 0L)), row.names = 23L, class = c("sf", 
"data.frame"), sf_column = "geom", agr = structure(c(PhasingCode = NA_integer_), levels = c("constant", 
"aggregate", "identity"), class = "factor"))

您可以在下面看到线条图。中心点应将 

MULTILINESTRING
分成几个较小的 
MULTILINSTRING
。最后，我应该注意我的真实世界数据集由一个 sf 对象中的多个 
MULTILINESTRING
组成。理想情况下，我想将每个 
MULTILINESTRING
分成更小的组。每个较小的组都应该有一个彼此唯一的 ID。这可能可以在最后一步实现，并根据几何列的行号进行分配。

ggplot() + geom_sf(data = testArea, color = "red") + geom_sf(data = split_point, color = "blue")
r split multilinestring sf-package
1个回答
0
投票

嗯..恐怕这两种数据、两种方法都有一些奇怪的地方：

让我们看看绘制的数据：

library(ggplot2)
ggplot() + 
  geom_sf(data = testArea, color = "red") + 
  geom_sf(data = split_point, color = "blue")

我们只考虑那些与 

split_point

相交/接触的部分 
testArea |>
  sf::st_cast(to = "LINESTRING") |>
  sf::st_intersects(split_point, sparse = FALSE)
#> Warning in st_cast.sf(testArea, to = "LINESTRING"): repeating attributes for
#> all sub-geometries for which they may not be constant
#>        [,1]
#>  [1,] FALSE
#>  [2,] FALSE
#>  [3,]  TRUE
#>  [4,] FALSE
#>  [5,] FALSE
#>  [6,] FALSE
#>  [7,] FALSE
#>  [8,] FALSE
#>  [9,] FALSE
#> [10,]  TRUE
#> [11,]  TRUE
#> [12,]  TRUE
#> [13,] FALSE
#> [14,] FALSE
#> [15,] FALSE
#> [16,] FALSE
#> [17,] FALSE
#> [18,]  TRUE
#> [19,] FALSE
#> [20,]  TRUE
#> [21,] FALSE
#> [22,] FALSE
#> [23,] FALSE
#> [24,] FALSE
#> [25,] FALSE
#> [26,]  TRUE
#> [27,] FALSE

所以，只有七个与它相交。

t <- testArea |>
  sf::st_cast(to = "LINESTRING") |>
  sf::st_filter(split_point)

t
#> Simple feature collection with 7 features and 1 field
#> Geometry type: LINESTRING
#> Dimension:     XY
#> Bounding box:  xmin: 6502228 ymin: 1962447 xmax: 6502549 ymax: 1962629
#> Projected CRS: NAD83 / California zone 3 (ftUS)
#> # A tibble: 7 × 2
#>   group                                                                     geom
#> * <dbl>                                            <LINESTRING [US_survey_foot]>
#> 1     3     (6502366 1962491, 6502371 1962496, 6502473 1962498, 6502477 1962622)
#> 2     3 (6502366 1962491, 6502376 1962481, 6502475 1962483, 6502494 1962483, 65…
#> 3     3                      (6502366 1962491, 6502367 1962469, 6502405 1962447)
#> 4     3                                       (6502366 1962491, 6502306 1962491)
#> 5     3     (6502366 1962491, 6502361 1962496, 6502228 1962496, 6502229 1962629)
#> 6     3                      (6502366 1962491, 6502376 1962471, 6502452 1962458)
#> 7     3                      (6502508 1962495, 6502498 1962494, 6502366 1962491)

但是，要按点分割线（字符串），我们必须确保它既不是线（字符串）的起点，也不是终点。我们来看看：

t <- t  |>
  subset(lwgeom::st_startpoint(t) != sf::st_as_sfc(split_point))
#> Linking to GEOS 3.11.1, GDAL 3.6.2, PROJ 9.1.1; sf_use_s2() is TRUE

t
#> Simple feature collection with 1 feature and 1 field
#> Geometry type: LINESTRING
#> Dimension:     XY
#> Bounding box:  xmin: 6502366 ymin: 1962491 xmax: 6502508 ymax: 1962495
#> Projected CRS: NAD83 / California zone 3 (ftUS)
#> # A tibble: 1 × 2
#>   group                                                geom
#>   <dbl>                       <LINESTRING [US_survey_foot]>
#> 1     3 (6502508 1962495, 6502498 1962494, 6502366 1962491)

只剩下一行了。让我们检查一下端点：

t |>
  subset(lwgeom::st_endpoint(t) != sf::st_as_sfc(split_point))
#> Simple feature collection with 0 features and 1 field
#> Bounding box:  xmin: NA ymin: NA xmax: NA ymax: NA
#> Projected CRS: NAD83 / California zone 3 (ftUS)
#> # A tibble: 0 × 2
#> # ℹ 2 variables: group <dbl>, geom <GEOMETRY [US_survey_foot]>

Sooo..没有一根线可以拼接。

也许您应该看看在您的点上开始/结束的线？并让它们的组（多线串）来检查哪些与这些字符串相交但不与 split_point 相交？这就是我的路线。

创建于 2024-04-12，使用 reprex v2.1.0

最新问题
© www.soinside.com 2019 - 2024. All rights reserved.