Sf package: Close a polygon fom complex shape

Question

Almost two weeks that my students and I are trying to find a solution on a trivial problem. We want to calculate the distance between two GPS points avoiding the coast according to this tutorial: https://www.r-bloggers.com/2020/02/three-ways-to-calculate-distances-in-r/

To this aim, we downloaded a shape file from the considered region, and as this shapefile only contained LINESTRING and that we have to turn our map into a grid, we decided to :

1/ Create a "corner point" to avoir any mistake

2/ Close the polygon But when the polygon is close, the shape is absolutely not what is expected. Lets see:

To download the same data we find shapefile here: https://www.marineregions.org/gazetteer.php?p=details&id=19888

aoi_boundary_HARV <- st_read("coasts_subnational/coasts_subnational.shp")

#close polygone
coord <- as.numeric(as.character(aoi_boundary_HARV$geometry[[1]]))
coord <- data.frame(matrix(coord, ncol = 2))
colnames(coord) <- c("x","y")
corner1 <- c(x = coord$x[nrow(coord)], y= coord$y[nrow(coord)])
corner2 <- c(x = coord$x[1], y= coord$y[1])
corner <- c(x = corner2[1], y = corner1[2])
names(corner) <- names(cord_Table)
cord_Table_plusCorner <- rbind(cord_Table, corner)
cord_Table_plusCorner <- as.data.frame(cord_Table_plusCorner)
colnames(cord_Table_plusCorner) <- c("x","y")
polygone_land <- data.frame(id = 1, 
                            x=cord_Table_plusCorner$x, 
                            y=cord_Table_plusCorner$y)

plot(polygone_land$x, polygone_land$y, type = "b")

So the shape seems good, we then needed to make it a polygon

CRS <- st_crs(aoi_boundary_HARV)
plot_aoi_boundary_HARV <- st_as_sf(polygone_land, 
                                   coords = c("x","y"), 
                                   crs = CRS)

polys = plot_aoi_boundary_HARV %>% 
  dplyr::group_by(id) %>% 
  dplyr::summarise() %>%
  st_cast("POLYGON")

plot(polys)

we get the following polygon:

We can clearly see that the external shape is more or less correct but the polygon is absolutely not well designed...

Thats where we would need a little help!!

We also used caveman, but did not work. We tried already those solutions: 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

But honestly it feels to bang our head against a wall... Would you have a tip ?

Thanks a lot in advance

Charlotte, Florian and Matthieu

Answer 1

I suggest avoiding splitting the coast linestring to individual coordinates; a possible alternative is intersecting it with its bounding box. A nice approach is using lwgeom::st_split() as a "blade" to cut the bounding box into pieces.

This will give you a number of polygons, and you need to choose one - which one will depend on what is your object of interest. Assuming you want the sea consider this piece of code:

library(sf)
library(dplyr)

coast <- st_read("coasts_subnational.shp") %>% 
  st_geometry() # o need for data now

bbox <- st_bbox(coast) %>% 
  st_as_sfc()

polygon <- bbox %>% 
  lwgeom::st_split(coast) %>% 
  st_collection_extract("POLYGON")

plot(polygon[3], col = 'steelblue')