rggplot2gisr-sf

whole earth polygon for world map in ggplot2 (and sf)


When plotting world maps there is a problem with ggplot2: it colours the whole background with the same colour, including the corners of the plot which aren't actually part of the globe, see the snapshot below produced by the following code (it uses bleading edge sf abd ggplot2 versions but the issue is generic, see the blog post mentioned below):

    #install.packages("devtools")
    #devtools::install_github("tidyverse/ggplot2")
    #devtools::install_github("edzer/sfr")

    library(ggplot2)
    library(sf)
    library(rnaturalearth)
    library(dplyr)

    theme_map <- function(...) {
      theme_minimal() +
      theme(
        text = element_text(family = "Ubuntu Regular", color = "#22211d"),
        axis.line = element_blank(),
        axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks = element_blank(),
        axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        panel.grid.minor = element_line(color = "#ebebe5", size = 0.2),
        panel.grid.major = element_line(color = "#ebebe5", size = 0.2),
        plot.background = element_rect(fill = "#f5f5f2", color = NA),
        panel.background = element_rect(fill = "#f5f5f2", color = NA),
        legend.background = element_rect(fill = "#f5f5f2", color = NA),
        panel.border = element_blank(),
        ...
      )
    }

    crs <- "+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +datum=WGS84 +units=m +no_defs"
    ctrys50m <- ne_countries(scale = 50, type = "countries", returnclass = "sf") %>%
       select(iso_a3, iso_n3, admin)

    ggplot() + 
      geom_sf(data = ctrys50m, alpha = 0.15, fill="grey") +
      coord_map() +
      coord_sf(crs = crs) +
      theme_map()

enter image description here

In order to be able to nicely plot Earth contour, in D3.js a special GeoJSON type, {type: "Sphere"} has been added, see this thread which can be seen in action here: it is the exterior whole Earth black border in the following snapshot:

enter image description here

The only trick I found in R/ggplot2 is the one published by Matt Strimas-Mackey, in his blog entry Mapping the Longest Commericial Flights in R, see the Bounding box and graticules section and the make_bbox and project_recenter functions.

These are quite a lot of code and I was wondering whether some sf or geom_sf code would make to a cleaner/simpler code, so I tried:

    # whole world WSG84 bounding box
    sphere <- ne_download(category = "physical", type = "wgs84_bounding_box", returnclass = "sf")
    sphere_laea <- st_transform(sphere, crs)
    ggplot() + 
      geom_sf(data = sphere, fill = "#D8F4FF") +
      coord_sf(crs = crs) +
      geom_sf(data = ctrys50m, alpha = 0.15, fill="grey") +
      coord_map() +
      coord_sf(crs = crs) +
      theme_map()

What I get is just an extra "anti-meridian" (note the line from North pole...) and no oceans filled with #D8F4FF... And the polygon is quite irregular at the bottom (the D3.js gurus did some smart adaptive resampling to increase the accuracy of projected lines...)

enter image description here

Any ideas about what's wrong in my attempt to get a whole world polygon for ggplot2 world maps? (Thanks for reading this far!)


Solution

  • New solution

    I recently found a piece of code from Dominic Royé which I think solves the problem remarkably well, here for orthographic projection:

    library(ggplot2)
    library(sf)
    library(rnaturalearth)
    library(dplyr)
    
    
    ortho_crs <-'+proj=ortho +lat_0=30 +lon_0=0.5 +x_0=0 +y_0=0 +R=6371000 +units=m +no_defs +type=crs'
    ocean <- st_point(x = c(0,0)) |>
      st_buffer(dist = 6371000) |> # Earth's radius
      st_sfc(crs = ortho_crs)
    
    ctrys50m <- ne_countries(scale = 50, type = "countries", returnclass = "sf") %>%
      select(iso_a3, iso_n3, admin) |> 
      st_intersection(st_transform(ocean, 4326)) |>
      st_transform(crs = ortho_crs) 
    
    ggplot()  +
      geom_sf(data = ocean, fill = "#D8F4FF", alpha = 0.7) +
      geom_sf(data = ctrys50m, fill="grey") +
      theme_bw()
    

    enter image description here

    Old solution

    Following Dave's suggestion, I created a new graticule small enough so that the convex hull operation will result in a good enough approximation of the border of the globe. The following code gives very good results (see image after):

    library(ggplot2)
    library(sf)
    library(rnaturalearth)
    library(dplyr)
    
    crs <- "+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +datum=WGS84 +units=m +no_defs"
    
    ctrys50m <- ne_countries(scale = 50, type = "countries", returnclass = "sf") %>%
      select(iso_a3, iso_n3, admin)
    
    sphere <- st_graticule(ndiscr = 10000, margin = 10e-6) %>%
      st_transform(crs = crs) %>%
      st_convex_hull() %>%
      summarise(geometry = st_union(geometry))
    
    ggplot()  +
      geom_sf(data = sphere, fill = "#D8F4FF", alpha = 0.7) +
      geom_sf(data = ctrys50m, fill="grey") +
      theme_bw()
    

    enter image description here