问题描述
编辑:此问题现已在此处正确回答: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)
此代码更具逻辑性。投射线条以指向点,合并点,然后创建多边形,但是也不能产生正确的结果。
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)
解决方法
此问题现已在此处正确回答:https://gis.stackexchange.com/a/378085/171458,请参阅该帖子以获取出色的解决方案和解释(不由我提供)。