rgeospatialspatialrgeo

How to transform a linestring geometry type to a SpatialPolygon


I am having troubles to transform a coastline of the east Atlantic area (west coast of Africa) into a SpatialPolygon to worh with it later. The problem is that when I transform it into a SpatialPolygon a line that merges the north of africa with the south is inclined and crosses the map as you can see in the picture I attached. One of the issues I think that can be causing this problem is the UTM tranformation maybe. I am working with several areas including the north and the south hemisphere, so may be this can be influencing the final product?The idea behing is to used for the creation of a mesh with R-INLA in order to project some sampling points in the ocean.

African west coast

library(sp)
library(sf)
library(rnaturalearth)
library(rnaturalearthdata)
library(FRK)

coastline <- ne_download(scale = "medium", type = "coastline", category = "physical", returnclass = "sf") #Download the data
bbox_west_africa <- st_bbox(c(xmin = -43.18, xmax = 19.34, ymin = -26.57, ymax = 29.49), crs = st_crs(coastline)) #to crop the data later

coastline_west_africa <- st_crop(coastline, bbox_west_africa) #The study area I need

coordinates <- st_coordinates(coastline) #Extract the coordinates to transform it

coordinates_df <- as.data.frame(coordinates) #Transform it to a Data.Frame


AfricanPolygon <- df_to_SpatialPolygons(coordinates_df,"L1",c("X","Y"), proj = CRS("+proj=longlat +ellps=sphere")) #tranform a data.frame to a SpatialPolygon
AfricanPolygon
plot(AfricanPolygon)

AfricanPolygon class : SpatialPolygons features : 34 extent : -43.18, 15.39318, -27.46624, 33.43104 (xmin, xmax, ymin, ymax) crs : +proj=longlat +ellps=sphere +no_defs


Solution

  • You are cutting a line from a set of lines and then trying to make a polygon - so sf joins the ends of the lines.

    Get a polygon of Africa and crop that instead. Either find one online or build it by dissolving African countries from Natural Earth.

    mapunits = ne_download(scale = "medium", type = "map_units", returnclass = "sf") 
    africa = mapunits[mapunits$CONTINENT=="Africa",]
    sf_use_s2(FALSE) # pretend lat-long is cartesian
    coast = st_crop(st_union(africa), bbox_west_africa)
    plot(coast)
    

    enter image description here

    You've got a bit of South America in your image, which I assume you don't want.

    Note the use of sf_use_s2(FALSE) so sf doesn't use spherical geometry, in which bounding boxes aren't boxes. If you want to transform everything to UTM and work in that, do that.