将 MULTILINESTRINGs 的特征集合组合成 R 中的一个封闭多边形
Combining a feature collection of MULTILINESTRINGs into one closed polygon in R
编辑： 这个问题现在在这里得到了正确的回答：https://gis.stackexchange.com/a/378085/171458
我有一个包含 262 行 MULTILINESTRING 的 shapefile。结合这个 shapefile（北极地区的政治边界）包裹了世界，我想检查一大组坐标是否落在这个 shapefile 中。我知道该怎么做，但这确实需要一个封闭的多边形，而不是线条的组合。我的一些解决方案已经取得了很大进展，但 none 确实有效。
我看过这个：https://gis.stackexchange.com/questions/290170/convert-a-linestring-into-a-closed-polygon-when-the-points-are-not-in-order and this page: 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,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,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)
这段代码感觉更符合逻辑。投线点，合并点，然后创建多边形，但也没有产生正确的结果。
amap_transformed_poly_2 <- st_cast(st_combine(st_cast(amap_transformed, to = "POINT")), to = "POLYGON")
ggplot() +
geom_sf(data = amap_transformed,
fill = NA,
col = "red",
size = 1) +
geom_sf(data = amap_transformed_poly_2,
fill = "green",
alpha = .3,
col = NA,
size = 10)
此问题现已在此处得到正确回答：https://gis.stackexchange.com/a/378085/171458，请参阅 post 以获得出色的解决方案和解释（不是我做的）。
编辑： 这个问题现在在这里得到了正确的回答：https://gis.stackexchange.com/a/378085/171458
我有一个包含 262 行 MULTILINESTRING 的 shapefile。结合这个 shapefile（北极地区的政治边界）包裹了世界，我想检查一大组坐标是否落在这个 shapefile 中。我知道该怎么做，但这确实需要一个封闭的多边形，而不是线条的组合。我的一些解决方案已经取得了很大进展，但 none 确实有效。
我看过这个：https://gis.stackexchange.com/questions/290170/convert-a-linestring-into-a-closed-polygon-when-the-points-are-not-in-order and this page: 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,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,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)
这段代码感觉更符合逻辑。投线点，合并点，然后创建多边形，但也没有产生正确的结果。
amap_transformed_poly_2 <- st_cast(st_combine(st_cast(amap_transformed, to = "POINT")), to = "POLYGON")
ggplot() +
geom_sf(data = amap_transformed,
fill = NA,
col = "red",
size = 1) +
geom_sf(data = amap_transformed_poly_2,
fill = "green",
alpha = .3,
col = NA,
size = 10)
此问题现已在此处得到正确回答：https://gis.stackexchange.com/a/378085/171458，请参阅 post 以获得出色的解决方案和解释（不是我做的）。