# Combining a feature collection of MULTILINESTRINGs into one closed polygon in R

I have a shapefile consisting of 262 rows of MULTILINESTRINGs. Combined this shapefile (the political boundary of the arctic region) wraps the world and I would like to check for a large set of coordinates if they fall within this shapefile. I know how to do that, but that does require having a closed polygon, instead of a combination of lines. Some of my solutions have come pretty far, but none actually work.

I have looked at this: Convert a linestring into a closed polygon when the points are not in order and this page: Converting points to polygons by group. I only have access to R to do this task, but none of these solutions work.

``````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
``````

This almost works, but is just an approximation and especially at the bottom of the plot a big chunk is incorrect.

``````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)
``````

This code feels more logical. Cast lines to point, combine points and then create a polygon, but does not yield the correct result either.

``````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)

``````

• in the zip file there's "amaplim_lam_poly.shp" which is a polygon like the one you are expecting, the trouble here is that your lines are not ordered "by the geometry", and, as far as I know, `sf` has not implemented and order by geometry as PostGIS Commented Oct 30, 2020 at 17:50
• You are right. I was too eager in loading in this data to have a look at other available shapefiles. I am going to leave this question open though, as there might be a solution possible that is viable for other situations like this one. Commented Oct 31, 2020 at 7:38

There was a lot of trouble trying to cast the file to polygon using the order of the lines or vertices; first I gave it a try with PostGIS ordering by the geometry, after converting to `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 = "")
``````

with the above code you get:

I also tried with `order by geom` in the sql query, with a worse result

Then I decided to inspect the lines, as to find a way to order them; `st_touches(amap, amap)` returned this:

``````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
``````

As you may see there was a linestring touching three other linestrings; in fact, there were up to four linestrings touching a single one, which makes it really complicated to make a line out of them; further inspecting, there were several linestrings which were closed, such as:

``````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)
``````

In this case, line 79, the blue one, is touching 3 other lines.

WKT shows that the first and last points of the purple linestring are the same, this would yield an invalid polygon geometry:

"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))"

BUT

There's this function: `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")
``````

• Stellar answer in both explanation and result. Commented Nov 1, 2020 at 19:37