Skip to content

Instantly share code, notes, and snippets.

@khufkens
Last active April 19, 2019 14:27
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
Star You must be signed in to star a gist
Save khufkens/8c7b82fffdf4c69ad068d39587f69943 to your computer and use it in GitHub Desktop.
Cold front transect along geopotential heights
# 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