Plotting with rnaturalearth, terra, and tidyterra

library(terra) # rasts and vects
library(tidyterra) # terra + ggplot
library(sf) # shapefiles
library(rnaturalearth) # boundary data
library(tidyverse) # everything else

Data

Here I’m just using an example raster of vertebrate richness in Brazil

rast <- rast("verts_richness_brazil.tif")
plot(rast)

Get boundary data w/ rnaturalearth

We can use this as a proxy for a basemap. Here we’ll use ne_countries and specify the continent. Note within that function you can alternatively specify individual countries. Another option to explore is ne_coastline, which gives the full coastline, but can be difficult to work with because it’s not enclosed polygons.

basemap_sf <- ne_countries(
  continent = "South America",
  scale = "medium", 
  returnclass = "sf") %>% 
  dplyr::select(1) %>% 
  st_transform(crs = crs(rast))  # make sure crs is same as raster
  
plot(basemap_sf)

# make vect to play more eaily with terra
# but we'll keep the sf for ggplot
basemap_vect <- basemap_sf %>% 
  vect()

plot(basemap_vect)

Simple plot with terra

plot(basemap_vect, col = "grey90")
plot(rast, add = TRUE)

Or can customize

Here I’ll use a color ramp from https://developers.arcgis.com/javascript/latest/visualization/symbols-color-ramps/esri-color-ramps/

colors <- colorRampPalette(
  colors = c("#85c1c8", "#af4980", "#c0182a", "#d33300", "#e99900", 
             "#ffff00"),
  bias = 0.25)

plot(basemap_vect, # grey fill behind raster
     col = "grey90", 
     border = "grey90", 
     axes = FALSE,
     mar = c(3.1, 3.1, 2.1, 7.1)) # to make room for legend
plot(rast, 
     add = TRUE, 
     col = colors(100), 
     maxcell = 5e+06, 
     axes = FALSE)
polys(basemap_vect, # lines on top of raster
      border = "grey20", 
      lwd = 0.25, 
      axes = FALSE)

Or can use tidyterra

ggplot() + 
  geom_sf(data = basemap_sf, # grey fill behind raster
          fill = "grey90",
          col = NA) + 
  geom_spatraster(
    data = rast,
    aes(fill = template_raster_10km_eckIV)) + # just the name of the layer
  geom_sf(data = basemap_sf, # lines on top of raster
          fill = NA,
          col = "grey20",
          lwd = 0.25) +
  scale_fill_gradientn(
    colors = colors(100),
    na.value = "transparent") +
  theme_void() + 
  labs(fill = NULL)