Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Downloads and animates CDS essential climate variables
# 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")
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment