Downloads and animates CDS essential climate variables

essential_climate_variables_animation.R

# Before you start install all the referenced packages below
# including the development release of ecmwfr you will also
# need the rnaturalearthdata in addition to the normal package
# To run the code copy and paste segments as this is still
# rough around the edges and requires user interaction to
# set your CDS API key.

# load libraries
if(!require(devtools)){install.packages("devtools")}
devtools::install_github("khufkens/ecmwfr")
library("ecmwfr")
library("raster")
library("tidyverse")
library("rnaturalearth")
library("sf")
library("gridExtra")
library("grid")
library("gganimate")

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

# download global coastline data from naturalearth
coastline <- ne_coastline(scale = 110, returnclass = c("sp"))

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

# transform the coastline to robinson
# roundabout through sp object, otherwise the date time line
# doesn't wrap nice (large landmass polygon)
coastline_robinson <- st_as_sf(spTransform(coastline,robinson))

# set your user credientials
user <- wf_set_key(service = "cds")

l <- list(
  format = "zip",
  year = c("1980", "1981", "1982", "1983", "1984", "1985", "1986", "1987",
           "1988", "1989", "1990", "1991", "1992", "1993", "1994", "1995",
           "1996", "1997", "1998", "1999", "2000", "2001", "2002", "2003",
           "2004", "2005", "2006", "2007", "2008", "2009", "2010", "2011",
           "2012", "2013", "2014", "2015", "2016", "2017", "2018", "2019"),
  variable = "surface_air_temperature",
  month = c("01", "02", "03", "04", "05", "06",
            "07", "08", "09", "10", "11", "12"),
  time_aggregation = "12_month",
  product_type = "anomaly",
  dataset = "ecv-for-climate-change",
  target = "download.zip"
)

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

# unzip zip file (when multiples are called this will be zipped)
unzip(file, exdir = tempdir())

files <- list.files(tempdir(), "*.grib", full.names = TRUE)

# load the grid data using raster
g <- stack(files)

# rotate the data to -180 / 180 notation
# this step take a real long time (hours) due to
# the many read and conversion operations
# be warned
g <- rotate(g)

# convert gridded raster data dataframe
g_df <- g %>%
  projectRaster(., res=50000, crs = robinson) %>%
  rasterToPoints %>%
  as.data.frame() %>%
  `colnames<-`(c("x", "y", names(g))) %>%
  gather(key = "layer", value = "val", -c("x","y")) %>%
  mutate(layer = paste0(substr(layer, 31, 36),"01")) %>%
  mutate(date = as.Date(layer, "%Y%m%d"))

# formulate graphing element
world_map <- ggplot() +
  theme_void() +
  geom_tile(data = g_df, aes(x = x, y = y, fill = val, group = date)) +
  scale_fill_gradient2(low = "dodgerblue4",
                       high = "red",
                       limits = c(min(g_df$val),
                                  max(g_df$val)),
                       name = expression(Delta*degree*"C")) +
  geom_sf(data=bb_robinson,
          colour='grey25',
          linetype='solid',
          fill = NA,
          size=0.7) +
  geom_sf(data=coastline_robinson,
            colour='grey25',
            linetype='solid',
            size=0.5) +
  theme(plot.title = element_text(hjust = 0.5),
        plot.margin = unit(c(2,2,2,2))) +
  labs(title = "12 month running mean temperature anomaly on {current_frame}") +
  transition_manual(date)

# animate the plot
gganimate::animate(world_map, width = 640, height = 400)

# save the animation
anim_save("~/ecv_anim.gif")