0

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:

enter image description here

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

5
  • Why do you need to do a zero buffer?
    – Ian Turton
    Jun 17, 2022 at 11:59
  • If that line 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 Apparently a known issue. Jun 17, 2022 at 12:28
  • then open the file in QGIS and fix it properly, see gis.stackexchange.com/questions/432792/…
    – Ian Turton
    Jun 17, 2022 at 12:44
  • "(I note that the Great Lakes are missing from North America, and Lake Victoria and friends are not present in Africa" this depends of the source shp that you use. In this case naturalearthdata with resolution 110m is a low-res shapefile, so some sort of simplification has been done.
    – dieghernan
    Jun 17, 2022 at 13:16
  • "Fix it properly" is not very helpful here. I am not all that familiar with QGIS and I do not even know what the problem is that is leading to this. Jun 17, 2022 at 13:20

1 Answer 1

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()

enter image description here


# 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)

1
  • This works wonderfully and I am still able to add all my points. Thank you so much. Jun 17, 2022 at 13:46

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.