Forked from khufkens/essential_climate_variables_animation.R
Created
April 3, 2020 12:43
-
-
Save matteodefelice/50018448c1dab3fddbfc895ed2780d85 to your computer and use it in GitHub Desktop.
Downloads and animates CDS essential climate variables
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
# 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