Last active
April 19, 2019 14:27
-
-
Save khufkens/8c7b82fffdf4c69ad068d39587f69943 to your computer and use it in GitHub Desktop.
Cold front transect along geopotential heights
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) | |
request <- list( | |
product_type = "reanalysis", | |
format = "netcdf", | |
day = "13", | |
month = "11", | |
year = "2013", | |
pressure_level = c( | |
"300","350","400","450","500","550","600","650","700", | |
"750","775","800","825","850","875","900","925", | |
"950","975","1000"), | |
variable = "temperature", | |
time = "00:00", | |
area = "40/-87/24/-80", # top left (lat/lon) / bottom right (lat/lon) | |
dataset = "reanalysis-era5-pressure-levels", | |
target = "download.nc" | |
) | |
user <- wf_set_key(service = "cds") | |
# download the data from CDS | |
file <- wf_request(request, user = user) | |
# read in data into raster stack | |
s <- stack(file) | |
# convert gridded raster data dataframe | |
s_df <- s[[37]] %>% | |
rasterToPoints %>% | |
as.data.frame() %>% | |
`colnames<-`(c("x", "y", "val")) | |
# create a bounding box from the stack extent | |
bb <- sf::st_union(sf::st_make_grid( | |
st_bbox(extent(s), | |
crs = st_crs(4326)), | |
n = 100)) | |
# read in coastline and trim it to size using the bounding box | |
coastline <- ne_coastline(scale = 110, returnclass = c("sf")) | |
coastline <- sf::st_intersection(bb, coastline) | |
# create a SpatialLine along which to sample | |
x <- c(-86,-81) | |
y <- c(39,25) | |
l <- SpatialLines(list(Lines(Line(cbind(x,y)), ID="transect"))) | |
projection(l) <- projection(s) | |
# extract values for all layers along transect | |
v <- as.data.frame(extract(s, l)) | |
# convert these to a dataframe | |
colnames(v) <- request$pressure_level | |
v$row <- as.numeric(1:nrow(v)) | |
v_df <- v %>% | |
gather(key = "y", value = "val", -row) %>% | |
mutate(y = as.numeric(y)) %>% | |
`colnames<-`(c("x", "y", "val")) | |
# add height parameter | |
# to fill all tiles | |
v_df$h <- 50 | |
p1 <- ggplot()+ | |
geom_tile(data = s_df, aes(x=x,y=y,fill=val)) + | |
scale_fill_viridis_c(name = expression(degree*"K")) + | |
geom_sf(data=bb, | |
colour='grey25', | |
linetype='solid', | |
fill = NA, | |
size=0.7) + | |
geom_sf(data=coastline, | |
colour='grey25', | |
linetype='solid', | |
size=0.5) + | |
geom_line(data = data.frame(x,y), aes(x, y), | |
arrow = arrow(ends = "last", | |
type = "closed"), | |
col = "white") + | |
theme_minimal() + | |
theme(axis.text.x = element_text(angle = 45, hjust = 1)) + | |
labs(x = "", | |
y = "", | |
title = "Florida cold front", | |
subtitle = "13 Nov. 2013") | |
p2 <- ggplot(data = v_df)+ | |
geom_tile(aes(x = x,y = y, height = h, fill = val)) + | |
scale_fill_viridis_c(name = expression(degree*"K")) + | |
labs(y = "geopotential height (hPa)", | |
x = "distance (pixels)", | |
title = "Geopotential height section", | |
subtitle = "section along arrow") + | |
theme_minimal() | |
gridExtra::grid.arrange(p1, p2, ncol = 2) |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment