Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Downloads and plots essential climate variable resembling the ECWMF plot
# 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
You can’t perform that action at this time.