Issue with graticule across 180° for several country/territory EEZs

Question

I am helping a set of data managers across the Pacific plot fisheries data within their EEZ. In trying to make a standard example that works for everyone, I have an issue with the latitude/longitude graticule for several countries that either cross or are near the dateline (180°). It's difficult to provide a MWE given the size of map shapefiles (EEZ, 12nm, 24nm shapefiles, etc.). So I will see if by showing my code and an example of what happens, someone can suggest a general fix. As others in the past have sought help with Fiji, I'll use that as an example.

The map shapefiles are all downloaded from marineregions.org, using the 0-360 versions, so that the Pacific is centered. The world's EEZs are read in, and Fiji extracted to its own object. It's then plotted using reasonable lat/lon limits

library(ggplot2)
library(sf)
library(maps)

ALL_EEZ <- sf::st_read(paste0(map.dir,"EEZ/eez_v12_0_360.shp")) #this is the World's EEZ shapefile to download from marineregions.org
FJ_EEZ<- ALL_EEZ[ALL_EEZ$TERRITORY1=="Fiji",]
geo_bounds<- c(170, 185,-25,-10) #west long, east lon, south lat, north lat for Fiji
ww=abs(geo_bounds[2]-geo_bounds[1])*100 # set window width
wh=abs(geo_bounds[4]-geo_bounds[3])*100 # set window height

if (dev.cur()!=1) dev.off() #close other windows
windows(ww,wh) # open window of correct dimensions (for mercator)

gg1<-ggplot() + 
  geom_sf(data = FJ_EEZ, colour = "black", fill=alpha("steelblue1",0.6), linewidth = 0.75) +
  coord_sf(xlim = c(geo_bounds[1],geo_bounds[2]), ylim = c(geo_bounds[3],geo_bounds[4])) +
  theme(panel.grid.major = element_line(color = gray(.5), linetype = "dotted", linewidth = 0.3), panel.background = element_rect(fill = NA, colour="black"))

print(gg1)

This produces this figure:

Adding scale_y_continuous and scale_x_continuous such as

scale_x_continuous(breaks = seq(geo_bounds[1], geo_bounds[2], 5)) +
scale_y_continuous(breaks = seq(geo_bounds[3], geo_bounds[4], 5)) +

doesn't result in a graticule extending past 180°.

however if I change the geobounds to

geo_bounds<- c(170, 195,-25,-10) #west long, east lon, south lat, north lat

The graticule draws correctly but I have a lot of area plotted east of Fiji that I don't want in the plot.

I have the same plotting issue with the EEZs of Gilbert Islands, New Zealand, Tonga, Tuvalu, and Wallis and Futuna.

I have searched the web and stackoverflow extensively for a simple solution to this but have not found one. I note again I'm trying to make a standard example for all the countries and territories that works the same and gives a centered map of the EEZ, on which they later plot data. The geobounds for each is read on from a .csv file I prepare for all 37 EEZs we are working with.

Without resorting to new map projections, is there a way to keep sensible geobounds and have the graticule drawn for each map?

Answer 1

As noted in the comments, this is a known issue and the solutions offered here involve compromise e.g. they do not not avoid resorting to new projections. However, by creating helper functions, it will help streamline the process somewhat. They may look daunting, particularly if you have less experienced R users in your group, but distributing the functions as-is will allow plotting to be standardised.

UPDATED: the way the xlim/ylim values were calculated in the original functions were not nuanced enough and could produce 'unbalanced' plots, see Bahrain as an example. The updated functions ensure the sf is centred in the plot. Custom breaks and labels for the y-axis are also now included. Original answer can be viewed in edit history.

The following demonstrates two options:

• plotting just the EEZ
• plotting the EEZ and adding subsequent data

The helper functions:

1. take a valid ALL_EEZ$TERRITORY1 value and desired number of x-axis breaks (4 seems to be appropriate for all the countries you listed)
2. derive a longitude centroid xcent value for centering map and creating custom PROJ string prj
3. to ensure axis labels and breaks can round to 2dp and are symmetrical relative to the sf centroid, derive axis range xrange and adjust xlim values xlims accordingly
4. derive x-axis limits xlims using xrange, calculate 2dp steps xsteps and adjust xlims
5. use xlims and xsteps to calculate x-axis breaks xbreaks
6. create x-axis labels xlabels based on WGS84/EPSG:4326 where range is -180 to 180
7. recalculate xcent based on xbreaks and use to create custom PROJ4 string prj, project sf using prj, and update xbreaks to new CRS
8. repeat steps 3 - 6 for y-axis
9. plot result

Option 1: Plot EEZ only

library(sf)
library(ggplot2)

# Set data path
map.dir <- "C:/path/to/data/"

# Load full EEZ data previously unzipped, downloaded from
# https://marineregions.org/downloads.php
ALL_EEZ <- st_read(paste0(map.dir,"EEZ/eez_v12_0_360.shp")) 

# Helper function to manage antimeridian/>180 longitude coordinates
# EEZ only
am_plot <- \(country, breaks) {
  
  if (!(country %in% ALL_EEZ$TERRITORY1)) {
    stop(paste("Error: Country", country, "not found in ALL_EEZ$TERRITORY1.",
               "Please check spelling and try again"))
  }
  
  sf <- ALL_EEZ[ALL_EEZ$TERRITORY1 == country,]
  lims <- unname(st_bbox(sf))
  xcent <- mean(lims[c(1, 3)])
  xrange <- ceiling((diff(lims[c(1, 3)]) / 2) * 100) / 100
  xlims <- round(c(xcent - xrange, xcent + xrange), 2)
  xsteps <- round((xlims[2] - xlims[1]) / (breaks - 1), 2)
  xlims[2] <- xlims[1] + (breaks - 1) * xsteps
  xbreaks <- seq(xlims[1], xlims[2], by = xsteps)
  xlabels <- ifelse(xbreaks > 180, 
                    paste0(sprintf("%.2f", abs(xbreaks - 360)), "°E"), 
                    paste0(sprintf("%.2f", xbreaks), "°W"))
  xcent <- mean(xbreaks)
  
  prj <- paste0("+proj=longlat +datum=WGS84 +lon_0=", xcent, " +no_defs")
  sf <- st_transform(sf, prj)
  xlims <- unname(st_bbox(sf)[c(1,3)])
  xbreaks <- seq(xlims[1], xlims[2], length.out = breaks)
  
  ycent <- mean(lims[c(2, 4)])
  yrange <- ceiling((diff(lims[c(2, 4)]) / 2) * 100) / 100
  ylims <- round(c(ycent - yrange, ycent + yrange), 2)
  ysteps <- round((ylims[2] - ylims[1]) / (breaks - 1), 2)
  ylims[2] <- ylims[1] + (breaks - 1) * ysteps
  ybreaks <- seq(ylims[1], ylims[2], by = ysteps)
  ylabels <- ifelse(ybreaks > 0, 
                    paste0(sprintf("%.2f", ybreaks), "°N"), 
                    paste0(sprintf("%.2f", abs(ybreaks)), "°S"))
  
  
  ggplot() +
    geom_sf(data = sf,
            colour = "black",
            fill = alpha("steelblue1", 0.6),
            linewidth = 0.75) +
    scale_x_continuous(breaks = xbreaks,
                       labels = xlabels) +
    scale_y_continuous(breaks = ybreaks,
                       labels = ylabels) +
    coord_sf(datum = prj) +
    theme(panel.grid.major = element_line(color = grey(0.5),
                                          linetype = "dotted",
                                          linewidth = 0.3),
          panel.background = element_rect(fill = NA, colour = "black"))
  
}


if (dev.cur()!=1) dev.off()

windows()

# Plot
fiji_plot <- am_plot("Fiji", 4)

print(fiji_plot)

Option 2: Plot EEZ with other data

This is essentially the same as above with two main modifications:

1. the prj, xbreaks/ybreaks, and xlabels/ylabels variables are written to the global environment using <<- so they can be called outside the function
2. scale_x_continuous()/scale_y_continuous() and coord_sf() are called outside the function to sidestep a message warning that a coordinate system and scale is already present in the plot (in order to avoid potentially confusing your stakeholders)

# Helper function to manage antimeridian/>180 longitude coordinates
# Add subsequent data
am_plot <- \(country, breaks) {
  
  if (!(country %in% ALL_EEZ$TERRITORY1)) {
    stop(paste("Error: Country", country, "not found in ALL_EEZ$TERRITORY1.",
               "Please check spelling and try again"))
  }
  
  sf <- ALL_EEZ[ALL_EEZ$TERRITORY1 == country,]
  lims <- unname(st_bbox(sf))
  xcent <- mean(lims[c(1, 3)])
  xrange <- ceiling((diff(lims[c(1, 3)]) / 2) * 100) / 100
  xlims <- round(c(xcent - xrange, xcent + xrange), 2)
  xsteps <- round((xlims[2] - xlims[1]) / (breaks - 1), 2)
  xlims[2] <- xlims[1] + (breaks - 1) * xsteps
  xbreaks <- seq(xlims[1], xlims[2], by = xsteps)
  xlabels <<- ifelse(xbreaks > 180, 
                    paste0(sprintf("%.2f", abs(xbreaks - 360)), "°E"), 
                    paste0(sprintf("%.2f", xbreaks), "°W"))
  xcent <- mean(xbreaks)
  
  prj <<- paste0("+proj=longlat +datum=WGS84 +lon_0=", xcent, " +no_defs")
  sf <- st_transform(sf, prj)
  xlims <- unname(st_bbox(sf)[c(1,3)])
  xbreaks <<- seq(xlims[1], xlims[2], length.out = breaks)
  
  ycent <- mean(lims[c(2, 4)])
  yrange <- ceiling((diff(lims[c(2, 4)]) / 2) * 100) / 100
  ylims <- round(c(ycent - yrange, ycent + yrange), 2)
  ysteps <- round((ylims[2] - ylims[1]) / (breaks - 1), 2)
  ylims[2] <- ylims[1] + (breaks - 1) * ysteps
  ybreaks <<- seq(ylims[1], ylims[2], by = ysteps)
  ylabels <<- ifelse(ybreaks > 0, 
                     paste0(sprintf("%.2f", ybreaks), "°N"), 
                     paste0(sprintf("%.2f", abs(ybreaks)), "°S"))
  
  ggplot() +
    geom_sf(data = sf,
            colour = "black",
            fill = alpha("steelblue1", 0.6),
            linewidth = 0.75) +
    theme(panel.grid.major = element_line(color = grey(0.5),
                                          linetype = "dotted",
                                          linewidth = 0.3),
          panel.background = element_rect(fill = NA, colour = "black"))
  
}

if(dev.cur()!=1) dev.off()

windows()

# Plot
tuvalu_plot <- am_plot("Tuvalu", 4)

print(tuvalu_plot)

# Create sample data
set.seed(42)

sf_points <- ALL_EEZ[ALL_EEZ$TERRITORY1 == "Tuvalu",] %>%
  st_sample(50)

# OPTIONAL: Project to current plot data CRS
# sf_points <- st_transform(sf_points, prj)

# Add data and plot parameters to current plot
tuvalu_plot +
  geom_sf(data = sf_points, colour = "#F8766D", size = 4) +
  scale_x_continuous(breaks = xbreaks,
                     labels = xlabels) +
  scale_y_continuous(breaks = ybreaks,
                     labels = ylabels) +
  coord_sf(datum = prj)