Zero-width buffer trick makes the Caspian Sea disappear (using R)

Question

I am new to plotting with coordinate systems. I am trying to generate a Pacific-centered map, and am using this wonderful code as a starting point. I have a bunch of data points from databases I'm generating and manipulating in R that I want to overlay on this map.

However, because I am using a different set of maps than the reference code (I do not want political boundaries), I am having to apply the zero-width gBuffer trick, which is causing the Caspian Sea to be filled in as land. Here is what I believe to be code for a minimal repro and the map:

library(data.table)
library(ggplot2)
library(rgdal)
library(rgeos)
library(maps)
library(maptools)
library(raster)
download.file(url = "https://www.naturalearthdata.com/http//www.naturalearthdata.com/download/110m/physical/ne_110m_land.zip", 
              destfile = "gallery/ne_110m_land.zip")
unzip(zipfile="gallery/ne_110m_land.zip", exdir = "gallery/ne_110m_land")
# delete the zip file
file.remove("gallery/ne_110m_land.zip")
# read the shapefile with readOGR from rgdal package
NE_land <- readOGR(dsn = "gallery/ne_110m_land", layer = "ne_110m_land")
class(NE_land) # is a SpatialPolygonsDataFrame object

shift <- 180+30

WGS84 <- CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")
split.line = SpatialLines(list(Lines(list(Line(cbind(180-shift,c(-90,90)))), ID="line")), proj4string=WGS84)

# apply the gBuffer trick
NE_land <- gBuffer(NE_land, byid=TRUE, width=0)

# intersecting line with country polygons
line.gInt <- gIntersection(split.line, NE_land)

# create a very thin polygon (buffer) out of the intersecting "split line"
bf <- gBuffer(line.gInt, byid=TRUE, width=0.000001)  

NE_land.split <- gDifference(NE_land, bf, byid=TRUE)
class(NE_land.split) # is a SpatialPolygons object

PROJ <- "+proj=eck4 +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs" 

Land.DT <- data.table(map_data(as(NE_land.split, "SpatialPolygonsDataFrame")))
# Shift coordinates
Land.DT[, long.new := long + shift]
Land.DT[, long.new := ifelse(long.new > 180, long.new-360, long.new)]
# project coordinates 
Land.DT[, c("X","Y") := data.table(project(cbind(long.new, lat), proj=PROJ))]

ggplot() + 
  geom_polygon(data = Land.DT, 
               aes(x = X, y = Y, group = group), 
               colour = "gray70", 
               fill = "gray90", 
               size = 0.25) +
  coord_equal() +
  theme_void()

Map:

Is there any way to avoid this issue? I'm guessing this kind of map technology does not play well with inland seas (I note that the Great Lakes are missing from North America, and Lake Victoria and friends are not present in Africa).

If the zero-buffer trick is removed, I get the following error:

NE_land is invalid Error in RGEOSBinTopoFunc(spgeom1, spgeom2, byid, id, drop_lower_td, unaryUnion_if_byid_false, : TopologyException: Input geom 0 is invalid: Ring Self-intersection at or near point -132.71000788000001 54.040009320000003 at -132.71000788000001 54.040009320000003

Answer 1

I think this is much easier with sf (see this question https://stackoverflow.com/questions/68278789/how-to-rotate-world-map-using-mollweide-projection-with-sf-rnaturalearth-ggplot/68313482#68313482).

Find here an example, if you need to ge back to sp use as(,"Spatial"):

library(sf)
f <- tempfile(fileext = ".zip")

download.file(
  url = "https://www.naturalearthdata.com/http//www.naturalearthdata.com/download/110m/physical/ne_110m_land.zip",
  destfile = f
)

unzip(zipfile = f, exdir = tempdir(), junkpaths = TRUE)

NE_land <- st_read(file.path(tempdir(), "ne_110m_land.shp")) %>%
  st_make_valid()


#> Simple feature collection with 127 features and 3 fields
#> Geometry type: POLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -180 ymin: -90 xmax: 180 ymax: 83.64513
#> Geodetic CRS:  WGS 84

# define a long & slim polygon that overlaps the meridian line & set its CRS to match
# that of world

# Centered in lon 150

lon_center <- 150

target_crs <- st_crs(paste0(
  "+proj=eck4 +x_0=0 +y_0=0 +lat_0=0 +lon_0=",
  lon_center
))
offset <- 180 - lon_center

polygon <- st_polygon(x = list(rbind(
  c(-0.0001 - offset, 90),
  c(0 - offset, 90),
  c(0 - offset, -90),
  c(-0.0001 - offset, -90),
  c(-0.0001 - offset, 90)
))) %>%
  st_sfc() %>%
  st_set_crs(4326)


# modify world dataset to remove overlapping portions with world's polygons
NE_land2 <- NE_land %>%
  st_difference(polygon) %>%
  st_transform(crs = target_crs)
#> Warning: attribute variables are assumed to be spatially constant throughout all
#> geometries

library(ggplot2)

ggplot(data = NE_land2) +
  geom_sf(
    colour = "gray70",
    fill = "gray90",
    size = 0.25
  ) +
  theme_void()

# And if you need to get back to Spatial*

tosp <- as(NE_land2, "Spatial")

class(tosp)
#> [1] "SpatialPolygonsDataFrame"
#> attr(,"package")
#> [1] "sp"

^{Created on 2022-06-17 by the reprex package (v2.0.1)}