sf: How to get back to MULTIPOLYGON from GEOMETRYCOLLECTION?

Question

I have a world country dataset, and would like to split it on the prime meridian, and re-center the data to focus on the Pacific.

I am trying to do this using Simple Features (sf), but am coming across an object-type issue I can't solve.

In order to split the data I tried the following:

   st_wg84 <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"

   # world country layer
   sfpolys <- rnaturalearth::ne_countries(scale = "medium", returnclass = "sf") 
   %>% st_sfc(crs = st_wg84 )

   # shift central/prime meridian towards west 
   shift <- 152 

   # create "split line" to split worldmap (split at Prime Meridian)
   split.line <- st_linestring(
     x = cbind(matrix(shift-180, 181, 1), matrix(-90:90,181,1))
    ) %>% 
     st_sfc(crs=st_wg84)

   # split country polygons along prime meridian
   sfpolys.split <- lwgeom::st_split(sfpolys, split.line)

Which works, resulting in a GEOMETRYCOLLECTION object, split along the desired line, containing the same number of features as the ingoing MULTIPOLYGON.

Next, I need to shift the coordinates in order to re-center the map, and to do this I must convert the polygon coordinates into a data frame.

    countries <- data.table(map_data(as(sfpolys.split, "Spatial")))

    # Shift coordinates to fall correctly on shifted map
    countries$long.shift <- countries$long + shift
    countries$long.shift <- ifelse(countries$long.shift > 180, 
    countries$long.shift - 360, countries$long.shift)

    # plot shifted map
    ggplot() + 
      geom_polygon(data=countries, 
        aes(x=long.shift, y=lat, group=group), 
        colour="black", fill="gray80", size = 0.25) +
      coord_equal()
  

However this conversion does not work with a GEOMETRYCOLLECTION, but it does with a MULTIPOLYGON.

So in order to get back to a MULTIPOLYGON I tried the following first:

sfpolys.split <- sfpolys.split %>% st_cast("MULTIPOLYGON")

But this results in the following error: "Error in m[1, ] : incorrect number of dimensions"

then I tried:

sfpolys.split <- sfpolys.split %>% st_collection_extract(type="POLYGON")

But this gives a POLYGON object, which I can't figure out how to group correctly into a MULTIPOLYGON.

Does anyone know either a better way of conducting this split and shift, or a simple way to get from a GEOMETRYCOLLECTION to a MULTIPOLYGON?

This is my desired result:

Dodgem

sf also has a function st_collection_extract() which might help. – Lithotrity

@Lithotrity st_collection_extract will only extract one of "POLYGON", "POINT", "LINESTRING". Infuriatingly, it doesn't extract any other type of geometry. – Rainproof

Answer 1

The GEOMETRYCOLLECTION is a list of geometriers, so we can extract the individual geometries.

Fortunately each of your GEOMETRYCOLLECTION geometries are POLYGONS, so we can wrap these up into MULTIPOLYGONS nicely

geoms <- lapply( sfpolys.split$geometry, `[` )
mp <- lapply( geoms, function(x) sf::st_multipolygon( x = x ) )

Then create an sfc

sfc_mp <- sf::st_sfc( mp )

and attach it to your object

sfpolys.split$mp <- sfc_mp
sfpolys.split <- sf::st_set_geometry( sfpolys.split, sfc_mp )

Here's a plot to check Greenland has been split. I've added a white border around each separate polygon

library(mapdeck)

sf_line <- sf::st_sf( geometry = split.line )

mapdeck() %>%
    add_path(
        data = sf_line
    ) %>%
    add_polygon(
        data = sfpolys.split
        , fill_colour = "name_pl"
        , stroke_colour = "#FFFFFF"
        , stroke_width = 50000
    )

The rest of your plotting code isn't reproducible, so I'm leaving that for you to sort out.

Aubreir