Combinación de una colección de características de MULTILINESTRING en un polígono cerrado en R

Question

Tengo un shapefile que consta de 262 filas de MULTILINESTRINGs. Combinado este shapefile (el límite político de la región ártica) envuelve el mundo y me gustaría comprobar para un gran conjunto de coordenadas si caen dentro de este shapefile. Sé cómo hacerlo, pero eso requiere tener un polígono cerrado, en lugar de una combinación de líneas. Algunas de mis soluciones han llegado bastante lejos, pero ninguna funciona realmente.

He mirado esto: Convertir una cadena de líneas en un polígono cerrado cuando los puntos no están en orden y esta página: Convertir puntos en polígonos por grupo . Sólo tengo acceso a R para realizar esta tarea, pero ninguna de estas soluciones funciona.

El shapefile puede descargarse aquí: https://www.amap.no/work-area/document/868

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

Esto casi funciona, pero es sólo una aproximación y, especialmente en la parte inferior del gráfico, una gran parte es incorrecta.

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)

Este código parece más lógico. Convierte las líneas en puntos, combina los puntos y luego crea un polígono, pero tampoco da el resultado correcto.

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)

Answer 1

Hubo muchos problemas al intentar convertir el archivo en polígono utilizando el orden de las líneas o los vértices; primero lo intenté con PostGIS ordenando por la geometría, después de convertir a LINESTRING

pgsql2shp -f amap_order -u db -P *** db "select 'a',ST_MakeLine(geom order by st_x(st_startpoint(geom))) as geom from amap_geo

read_sf("amap_order.shp") %>% st_cast("POLYGON") %>% 
    st_transform( crs = polar_projection) %>%
  plot(graticule = T, axes = T, main = "")

con el código anterior se obtiene:

También probé con order by geom en la consulta sql, con un resultado peor

Entonces decidí inspeccionar las líneas, como para encontrar una manera de ordenarlas; st_touches(amap, amap) devolvió esto:

Sparse geometry binary predicate list of length 262, where the predicate was `touches'
first 10 elements:
 1: 3, 210
 2: 3, 208
 3: 1, 2
 4: 7, 133
 5: 6, 8
 6: 5, 7
 7: 4, 6, 133
 8: 5, 9
 9: 8, 10
 10: 9, 11

Como puedes ver, había una línea que tocaba a otras tres; de hecho, había hasta cuatro líneas que tocaban a una sola, lo que hace realmente complicado hacer una línea con ellas; si sigues inspeccionando, había varias líneas que estaban cerradas, como por ejemplo:

bbx = st_bbox(amap[c(79,82),])
plot(amap[c(79,80,81,82),5], xlim = c(bbx[1],bbx[3]), ylim = c(bbx[2], bbx[4]),
     main = "Rows 79,80,81,82", graticule = T)

En este caso, la línea 79, la azul, está tocando otras 3 líneas.

WKT muestra que el primer y el último punto de la línea púrpura son iguales, esto produciría una geometría de polígono inválida:

"MULTILINESTRING ((-60.41898 55.21146, -60.41361 55.24055, -60.39611 55.25361, -60.37389 55.26055, -60.35528 55.25944, -60.34805 55.25055, -60.39084 55.22249, -60.41898 55.21146))"

PERO

Hay una función: st_split :

library(lwgeom)
ch = amap_transformed %>% st_combine() %>% st_convex_hull()
ch_split = st_split(ch, st_union(amap_transformed))
ch_coll = st_collection_extract(ch_split) 
plot(ch_coll[1], col = "blue")