将MULTILINESTRING的要素集合组合到R

问题描述

编辑:此问题现已在此处正确回答:https://gis.stackexchange.com/a/378085/171458

我有一个shapefile,由262行MULTILInesTRING组成。结合使用这个shapefile(北极地区的政治边界)包裹了整个世界,我想检查是否有大量坐标位于该shapefile中。我知道该怎么做,但这确实需要有一个封闭的多边形,而不是线条的组合。我的一些解决方案已经走得很远,但实际上没有任何作用。

我看过以下内容https://gis.stackexchange.com/questions/290170/convert-a-linestring-into-a-closed-polygon-when-the-points-are-not-in-order和本页:https://gis.stackexchange.com/questions/332427/converting-points-to-polygons-by-group

可以在此处下载shapefile:https://www.amap.no/work-area/document/868

我也在这里发布了它:https://gis.stackexchange.com/questions/378010/combining-a-feature-collection-of-multilinestrings-into-one-closed-polygon-in-r,但是我也在这里问它,因为我只有R可用于此任务。

library(sf)
library(tidyverse)
library(concaveman)

polar_projection <- "+proj=stere +lat_0=90 +lat_ts=71 +lon_0=0 +k=1 +x_0=0 +y_0=0 +datum=wgs84 +units=m +no_defs +ellps=wgs84 +towgs84=0,0"
amap <- read_sf("amaplim_geo_nw.shp",layer = "amaplim_geo_nw")
amap_transformed <- st_transform(amap,crs = polar_projection)

amap_transformed
Simple feature collection with 262 features and 8 fields
geometry type:  MULTILInesTRING
dimension:      XY
bBox:           xmin: -4358644 ymin: -3324819 xmax: 2893297 ymax: 4379941
CRS:            +proj=stere +lat_0=90 +lat_ts=71 +lon_0=0 +k=1 +x_0=0 +y_0=0 +datum=wgs84 +units=m +no_defs +ellps=wgs84 +towgs84=0,0
# A tibble: 262 x 9
   FNODE_ TNODE_ Lpoly_ Rpoly_  LENGTH AMAPLIM3G_ AMAPLIM3G1 STATUS                                                                         geometry
 *  <dbl>  <dbl>  <dbl>  <dbl>   <dbl>      <dbl>      <dbl>  <int>                                                            <MULTILInesTRING [m]>
 1     34     37      0      0 242366.          1      16682      1                                       ((-538630.6 -3054726,-463882.6 -3300696))
 2     38     36      0      0 215303.          2      16683      1                                       ((-235293.1 -3324819,-218966.3 -3094112))
 3     37     38      0      0 227471.          3          5      1 ((-463882.6 -3300696,-442961.7 -3303569,-417752.4 -3306851,-391767.9 -331003~
 4    386    386      0      0   4822.          4       4382      1 ((-3397144 278632.8,-3398147 277487.9,-3396967 276827.8,-3396465 277740.4,-~
 5    387    388      0      0  13976.          5       4499      1 ((-3526696 268595.1,-3523913 268160.4,-3522760 267883.9,-3520248 266594.5,-~
 6    388    389      0      0  20817.          6       4499      1 ((-3511875 266079.8,-3509774 265407.4,-3507925 263556.7,-3506566 262993.4,-~
 7    389    386      0      0 168362.          7       4382      1 ((-3490612 262851.4,-3488809 262562.4,-3485696 261529.9,-3482610 260126.9,-~
 8    391    387      0      0 568338.          8       4499      1 ((-3665413 25038.86,-3663223 25042.43,-3662560 26120.65,-3662054 27060.08,-~
 9    392    391      0      0  29253.          9       4499      1 ((-3682065 -2.25454e-10,-3679368 2426.171,-3677501 4582.511,-3676795 5716.63~
10    393    392      0      0 259118.         10       4499      1 ((-3805480 -196918.1,-3805796 -195642.3,-3805842 -193614.8,-3805680 -192121.~
# ... with 252 more rows

这几乎可以用,但是只是一个近似值,尤其是在图的底部,一大块是不正确的。

amap_transformed_poly <- concaveman(amap_transformed,concavity = 1)

ggplot() +
  geom_sf(data = amap_transformed,fill = NA,col = "red",size = 1) +
  geom_sf(data = amap_transformed_poly,fill = "blue",alpha = .3,col = NA,size = 1)

enter image description here

代码更具逻辑性。投射线条以指向点,合并点,然后创建多边形,但是也不能产生正确的结果。

amap_transformed_poly_2 <- st_cast(st_combine(st_cast(amap_transformed,to = "POINT")),to = "polyGON")

ggplot() +
  geom_sf(data = amap_transformed,size = 1) +
  geom_sf(data = amap_transformed_poly_2,fill = "green",size = 10)

enter image description here

解决方法

此问题现已在此处正确回答:https://gis.stackexchange.com/a/378085/171458,请参阅该帖子以获取出色的解决方案和解释(不由我提供)。