Mapping in ggplot2 and R - animated raster maps

animate raster maps

Although maps are generally meant to summarize multi-dimensional data in a simple a format as possible animation can provide additional insight into the data at hand. Here, I’m building upon a previous tutorial to illustrate how to animate climatological raster data.

In this example I’m going to use data from the ECMWF Atmosphere Data Store (ADS) service and their forecasts of black carbon transport throughout the atmosphere of the northern hemisphere in response to heavy forest fires in the west of the US and over eastern Siberia.

So, let’s get started!

First we setup our working environment ! You will need our ecmwfr package for downloading climate data, the tidyverse for processing and plotting, raster to deal with gridded climate data and sf + rnaturalearth for vector features, showtext for fancy fonts, finally both gganimate and gifski to render the final animated gif file.

# load libraries
library("ecmwfr")
library("raster")
library("tidyverse")
library("rnaturalearth")
library("sf")
library("gganimate")
library("gifski")

# custom fonts
library(showtext)
font_add_google("Prata", regular.wt = 400)
showtext_auto()

Downloading the data using our ecmwfr package is easy (see links below). For any given date you can download forecasts up to 120 hours into the future. I’ll use data downloaded earlier this season, but you can always use more current data in your own analysis / map.

# formulate a ECMWF API request
request <- list(
  date = "2021-07-21/2021-07-21",
  type = "forecast",
  format = "netcdf_zip",
  variable = "black_carbon_aerosol_optical_depth_550nm",
  time = "00:00",
  leadtime_hour = as.character(1:120),
  area = c(90, -180, 0, 180),
  dataset_short_name = "cams-global-atmospheric-composition-forecasts",
  target = "download.zip"
)

# download the data (file location is returned)
file <- wf_request(
  request, 
  user = "xyz"
  )

# unzip zip file (when multiples are called this will be zipped)
unzip(file, exdir = tempdir())
files <- list.files(tempdir(), "*.nc", full.names = TRUE)

# copy files to the data directory
file.copy(files, "data/carbon.nc")

With all data and libraries loaded we can start composing our figure. As with the previous global map I’ll use the Robinson projection. I set this as a variable and reproject all required layers.

To provide context to the map we download a simple features (sf) vector of country outlines using the rnaturalearth package. In addition, we define a bounding box outlining the total extent of the map. Both are reprojected to the Robinson projection.

# set coordinate systems
robinson <- CRS("+proj=robin +over")

# create a bounding box for the robinson projection
bb <- sf::st_union(sf::st_make_grid(
  st_bbox(c(xmin = -179.999,
            xmax = 179.999,
            ymax = 90,
            ymin = -1), crs = st_crs(4326)),
  n = 100))
bb_robinson <- st_transform(bb, as.character(robinson))

# download global coastline data from naturalearth
countries <- ne_countries(scale = 110, returnclass = c("sf"))

# clip countries to bounding box
# and transform
countries_robinson <- countries %>%
  st_buffer(0) %>%
  st_intersection(st_union(bb)) %>%
  st_transform(robinson)

For rapid rendering we reproject and resample the original map. Using pipes, and in one movement, we convert this reprojected raster file to a dataframe. Note that you can change the resolution of the output by adjusting the resampling value. For a full resolution map you can disable resampling.

# load the grid data using raster
g <- stack("data/carbon.nc")

# convert gridded raster data dataframe
g_df <- g %>%
  projectRaster(., res=50000, crs = robinson) %>%
  rasterToPoints %>%
  as.data.frame() %>%
  `colnames<-`(c("x", "y", names(g))) %>%
  pivot_longer(cols = starts_with("X20"),
               names_to = "layer",
               values_to = "val") %>%
  mutate(layer = substr(layer, 2, 14)) %>%
  mutate(date = as.POSIXct(layer, "%Y.%m.%d.%H", tz = "UTC")
         )

With that all the required elements are preformated and can be used within a ggplot routine. We can now build the map, only taking care of small changes particular to the gganimate package in order to animate the map over time (or any variable for that matter).

# formulate graphing element
world_map <- ggplot() +
  theme_void() +
  geom_tile(
    data = g_df,
    aes(
      x = x,
      y = y,
      fill = log(val),
      group = date)) +
  scale_fill_viridis_c(
    option = "B"
  ) +
  geom_sf(data=countries_robinson,
          colour='white',
          linetype='solid',
          fill = NA,
          size=0.2) +
  geom_sf(data=bb_robinson,
          colour='white',
          linetype='solid',
          fill = NA,
          size= 1) +
  coord_sf(ylim = c(-17309.98, 8582690)) +
  theme(
    plot.title = element_text(
      family = "Prata",
      size = 50,
      hjust = 0.5),
    plot.subtitle = element_text(
      hjust = 0.5,
      size = 25),
    plot.caption = element_text(
      color = "grey50",
      size = 12,
      hjust = 0.9),
    legend.position = "none"
    ) +
  labs(
    title = "Fire season",
    subtitle = "",
    caption = "Data source: Copernicus ADS at {current_frame}") +
  transition_manual(date)

Note that in the above map we define an graphing element called transition_manual(date) which defines that if rendered the routine should step through all dates of the input data frame. It is therefore important to format your input data well and provide it as a consistent data frame (tibble).

With all this specified you can then ask gganimate to render the final plot.

#animate the plot
gganimate::animate(
  world_map,
  width = 1000,
  height = 400,
  renderer=gifski_renderer("map.gif")
  )

And that’s it! This will give you the animated map as shown below (to save space the animated gif was converted to a mp4 movie). You see dominant wind and weather systems carry the smoke from the western US to the east coast and beyond. Smoke from Siberia almost makes it to the west coast of the US.

As always, keep an eye on our blog for new entries in this series on map building, and tips and tricks to improve your R mapping skills. The code for this project can be found on our github R tutorials page. If unfamiliar with the package, the ecmwfr package can be found here.

Avatar
Koen Hufkens, PhD
Partner, Researcher

As an earth system scientist and ecologist I model ecosystem processes.

Related

Previous