sf MULTIPOLYGONs not mergeable across countries

[courtiol]

Dear maintainer, I am not quite sure whether this is a sf or rnaturalearth issue (or just me), but when trying to combine polygons of all countries using sf output from rnaturalearth and sf function st_union, the results I obtain are imperfect irrespective of the scale:

library(sf)
#> Linking to GEOS 3.10.2, GDAL 3.4.1, PROJ 8.2.1; sf_use_s2() is TRUE
packageVersion("sf")
#> [1] '1.0.14'
packageVersion("rnaturalearth")
#> [1] '0.3.3.9000'
borders <- rnaturalearth::ne_countries(scale = 'large', returnclass = 'sf')
#> The legacy packages maptools, rgdal, and rgeos, underpinning the sp package,
#> which was just loaded, were retired in October 2023.
#> Please refer to R-spatial evolution reports for details, especially
#> https://r-spatial.org/r/2023/05/15/evolution4.html.
#> It may be desirable to make the sf package available;
#> package maintainers should consider adding sf to Suggests:.
borders <- borders[, 0]
plot(sf::st_union(borders), bg = 2, col = 3)

borders <- rnaturalearth::ne_countries(scale = 'medium', returnclass = 'sf')
borders <- borders[, 0]
plot(sf::st_union(borders), bg = 2, col = 3)

borders <- rnaturalearth::ne_countries(scale = 'small', returnclass = 'sf')
borders <- borders[, 0]
plot(sf::st_union(borders), bg = 2, col = 3)
#> Error in wk_handle.wk_wkb(wkb, s2_geography_writer(oriented = oriented, : Loop 0 is not valid: Edge 0 crosses edge 78

^{Created on 2023-10-05 with reprex v2.0.2}

Note the with the scale set to small, the plotting does not work at all due to a polygon error.

[courtiol]

For info: calling sf::sf_use_s2(FALSE) circumvents the issue, but my understanding is that it should not be needed. (Am I wrong @edzer?)

[edzer]

We try to explain here why many geometries cannot be valid on the sphere (S2) and on the flat surface (R2, what you are looking at) at the same time. After manipulations assuming validity on S2, you'd have to un-S2 them to get back things that "plot nicely" on R2.

[PMassicotte]

Thank you for the comment @edzer. Bought the book recently and halfway through. Nice reading.

[edzer]

Thanks, @PMassicotte !

[courtiol]

Thanks to you both and to @liamdbailey who suggested a projection. Here is a solution:

library(sf)
#> Linking to GEOS 3.10.2, GDAL 3.4.1, PROJ 8.2.1; sf_use_s2() is TRUE
packageVersion("sf")
#> [1] '1.0.14'
packageVersion("rnaturalearth")
#> [1] '0.3.3.9000'
sf_use_s2(TRUE)
borders <- rnaturalearth::ne_countries(scale = 'large', returnclass = 'sf')
#> The legacy packages maptools, rgdal, and rgeos, underpinning the sp package,
#> which was just loaded, were retired in October 2023.
#> Please refer to R-spatial evolution reports for details, especially
#> https://r-spatial.org/r/2023/05/15/evolution4.html.
#> It may be desirable to make the sf package available;
#> package maintainers should consider adding sf to Suggests:.
borders <- borders[, 0]
borders <- st_transform(borders, "ESRI:54030")
borders <- st_union(borders)
borders <- st_transform(borders, 4326)
plot(borders, bg = 2, col = 3)

borders <- rnaturalearth::ne_countries(scale = 'medium', returnclass = 'sf')
borders <- borders[, 0]
borders <- st_transform(borders, "ESRI:54030")
borders <- st_union(borders)
borders <- st_transform(borders, 4326)
plot(borders, bg = 2, col = 3)

borders <- rnaturalearth::ne_countries(scale = 'small', returnclass = 'sf')
borders <- borders[, 0]
borders <- st_transform(borders, "ESRI:54030")
borders <- st_make_valid(borders) ## fix validity issue
borders <- st_union(borders)
borders <- st_transform(borders, 4326)
plot(borders, bg = 2, col = 3)

^{Created on 2023-10-05 with reprex v2.0.2}

[edzer]

Yes, that is effectively the same as using sf_use_s2(FALSE): you instruct the geom engine to treat the Earth as a geometry on a flat space. Without reprojecting, doing this on R2 is analogue to assuming (implicitly) the data are projected on plate carree.