rggplot2plotmapping

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

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 Fiji map with incomplete graticule

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. Fiji map with correct (but extended) graticule

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?

Solution

  • 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.

    The following demonstrates two options:

    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 cent value for centering map and creating custom PROJ string prj
    3. derive x-axis limits xlims and breaks xbreaks from st_bbox(sf) based on original CRS, and creates equivalent WGS84/EPSG:4326 x-axis labels xlabels. This replaces the need for your geo_bounds values
    4. project sf to prj CRS and recreate xlims and xbreaks
    5. 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,]
  cent <- suppressWarnings(round(st_coordinates(st_centroid(sf))[1]))
  xlims <- unname(st_bbox(sf)[c(1,3)])
  xbreaks <- round(seq(floor(xlims[1]), ceiling(xlims[2]), length.out = breaks), 2)
  xlabels <- ifelse(xbreaks > 180, 
                    paste0(sprintf("%.2f", xbreaks - 360), "°E"), 
                    paste0(sprintf("%.2f", xbreaks), "°W"))
  
  prj <- paste0("+proj=longlat +datum=WGS84 +lon_0=", cent, " +no_defs")
  sf <- st_transform(sf, prj)
  xlims <- unname(st_bbox(sf)[c(1,3)])
  xbreaks <- round(seq(floor(xlims[1]), ceiling(xlims[2]), length.out = breaks), 2)
  
  ggplot() +
    geom_sf(data = sf,
            colour = "black",
            fill = alpha("steelblue1", 0.6),
            linewidth = 0.75) +
    scale_x_continuous(breaks = xbreaks,
                       labels = xlabels) +
    coord_sf(datum = prj,
             xlim = range(xbreaks)) +
    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)

    1

    Option 2: Plot EEZ with other data

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

    1. the xlims, xbreaks, and xlabels variables are written to the global environment using <<- so they can be called outside the function
    2. scale_x_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
# and 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,]
  cent <- suppressWarnings(round(st_coordinates(st_centroid(sf))[1]))
  xlims <- unname(st_bbox(sf)[c(1,3)])
  xbreaks <- round(seq(floor(xlims[1]), ceiling(xlims[2]), length.out = breaks), 2)
  xlabels <<- ifelse(xbreaks > 180, 
                     paste0(sprintf("%.2f", xbreaks - 360), "°E"), 
                     paste0(sprintf("%.2f", xbreaks), "°W"))
  
  prj <<- paste0("+proj=longlat +datum=WGS84 +lon_0=", cent, " +no_defs")
  sf <- st_transform(sf, prj)
  xlims <- unname(st_bbox(sf)[c(1,3)])
  xbreaks <<- round(seq(floor(xlims[1]), ceiling(xlims[2]), length.out = breaks), 2)
  
  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) +
  coord_sf(datum = prj,
           xlim = range(xbreaks))

    2