Last active
April 12, 2019 07:33
Star
You must be signed in to star a gist
Downloads and plots essential climate variable resembling the ECWMF plot
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(tidyverse) | |
library(raster) | |
library(rnaturalearth) | |
library(sf) | |
library(gridExtra) | |
library(grid) | |
# set coordinate systems | |
robinson <- CRS("+proj=robin +over") | |
EU <- CRS("+init=epsg:32661") | |
# set your user credientials | |
user <- wf_set_key(service = "cds") | |
# create a synoptic variable query | |
l <- list( | |
dataset = "ecv-for-climate-change", | |
format = "tgz", | |
variable = "surface_air_temperature", | |
year = "2019", | |
month = "03", | |
time_aggregation = "1_month", | |
product_type = "anomaly", | |
target = "download.tgz" | |
) | |
# download the data (file location is returned) | |
file <- wf_request(l, user = user) | |
# load the grid data using raster | |
g <- raster(file) | |
# rotate the data to -180 / 180 notation | |
g <- rotate(g) | |
# 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)) | |
coastline_eu <- st_as_sf(spTransform(coastline, EU)) | |
# convert gridded raster dato dataframe | |
g_df <- g %>% | |
projectRaster(., res=50000, crs = robinson) %>% | |
rasterToPoints %>% | |
as.data.frame() %>% | |
`colnames<-`(c("x", "y", "val")) | |
# crop and reproject to EU overview | |
g_c <- crop(g, extent(c(-180, 180, 0, 90))) | |
g_eu <- g_c %>% projectRaster(crs = EU) | |
g_eu <- crop(g_eu, extent(c(-1237268, 6357552, -5168481, 2360296))) | |
# create a bounding box for the robinson projection | |
bb_eu <- st_make_grid(st_bbox(g_eu), n = 1) | |
bb_eu <- st_transform(bb_eu, as.character(EU)) | |
# trim coastline | |
coastline_eu <- st_intersection(bb_eu, coastline_eu) | |
eu_df <- g_eu %>% | |
rasterToPoints %>% | |
as.data.frame() %>% | |
`colnames<-`(c("x", "y", "val")) | |
world_map <- ggplot()+ | |
geom_tile(data = g_df, aes(x=x,y=y,fill=val)) + | |
scale_fill_gradient2(low = "dodgerblue4", | |
high = "red", | |
limits = c(cellStats(g, min),cellStats(g, max))) + | |
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_void() + | |
theme(legend.position = 'none') | |
eu_map <- ggplot()+ | |
geom_tile(data = eu_df,aes(x=x,y=y,fill=val)) + | |
scale_fill_gradient2(low = "dodgerblue4", | |
high = "red", | |
limits = c(cellStats(g, min),cellStats(g, max)), | |
name = expression(degree*"C")) + | |
geom_sf(data=coastline_eu, | |
colour='grey25', | |
linetype='solid', | |
size=0.5) + | |
geom_sf(data=bb_eu, | |
colour='grey25', | |
linetype='solid', | |
fill = NA, | |
size=0.7) + | |
theme_void() + | |
coord_sf(datum=NA) | |
p <- arrangeGrob( | |
world_map, | |
eu_map, | |
layout_matrix = rbind(c(1,1,2)), | |
top = textGrob('Surface air temperature anomaly for March 2019 relative to 1981 - 2010', gp = gpar(fontsize = 18, font = 8)), | |
vp=viewport(width=0.98, height=1)) | |
ggsave(file="temp_anomalies.pdf", p, height = 5) |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment