Mapping in ggplot2 and R - visualizing dynamic features

plotting wind direction and intensity

Plotting dynamic weather data such as wind, which have both a magnitude and direction, can be a challenge. How do you visually communicate direction and magnitude? In these cases streamlines might help. Here the length of the streamline scales proportionally to the magnitude of the wind force, while the direction of the wind is represented by the shape of the streamline (both in width and curvature).

In this example I’m going to use data from the ECMWF Climate Data Store (CDS) service, in particular 10m u and v components of wind, and mean sea level pressure to visualize an intense wind field combined with the pressure field. In this example I use the wind and pressure data at noon around the maximum intensity of hurricane Katrina on August the 28th, 2005.

Let’s get started!

First of all a set of libraries are required.

# read in libraries
library(tidyverse)
library(raster)
library(metR)
library(rnaturalearth)
library(ecmwfr)

Downloading the data using the BlueGreen Labs ecmwfr package is easy (see links below).

# formulate a ecmwfr query
request <- list(
  product_type = "reanalysis",
  format = "netcdf",
  variable = c(
    "10m_u_component_of_wind",
    "10m_v_component_of_wind",
    "mean_sea_level_pressure"),
  year = "2005",
  day = "28",
  month = "08",
  time = "12:00",
  area = c(38, -104, 20, -73),
  dataset_short_name = "reanalysis-era5-single-levels",
  target = "wind.nc"
)

# requires preset credentials using
# wf_set_key(), see ecmwfr docs
wf_request(
  request = request,
  user = "2088",
  path = "data/"
  )

With all data and libraries loaded we can start composing our figure. To provide context to the map we download a simple features (sf) vector of country outlines using the rnaturalearth package.

# coastline sf object
land <- ne_coastline(
  scale = 50,
  returnclass = "sf")

In addition you have to wrangle the netcdf raster data into a long format compatible with ggplot and the metR package I use to plot the wind streamlines.

# convert pressure
pressure <- raster::raster(
  "data/wind.nc",
  varname = "msl"
  ) %>%
  rasterToPoints() %>%
  as.data.frame() %>%
  rename(
    'lat' = 'y',
    'lon' = 'x',
    'pressure' = 'Mean.sea.level.pressure'
  ) %>%
  pivot_longer(
    cols = "pressure",
    names_to = "layer",
    values_to = "pressure"
  ) %>%
  mutate(
    pressure = pressure / 100
  )

# convert wind
u <- raster::raster("data/wind.nc",
                   varname = "u10") %>%
  rasterToPoints() %>%
  as.data.frame() %>%
  rename(
    'lat' = 'y',
    'lon' = 'x',
    'u' = starts_with("X10")
  ) %>%
  pivot_longer(
    cols = "u",
    names_to = "layer",
    values_to = "u"
  ) %>%
  dplyr::select(-layer)

v <- raster::raster(
  "data/wind.nc",
  varname = "v10") %>%
  rasterToPoints() %>%
  as.data.frame() %>%
  rename(
    'lat' = 'y',
    'lon' = 'x',
    'v' = starts_with("X10")
  ) %>%
  pivot_longer(
    cols = "v",
    names_to = "layer",
    values_to = "v"
  ) %>%
  dplyr::select(-layer)

wind <- left_join(u, v)

With that all the required elements formatted we can now build the map!

# NOTE: make sure to load the theme settings
# below to beautify the plot !!
p <- ggplot() +
  geom_contour_fill(
    data = pressure,
    aes(
      lon,
      lat,
      z = pressure,
      fill = stat(level)
    ),
    breaks = MakeBreaks(3)
  ) +
  geom_sf(data = land,
          color = "white",
          fill = NA,
          size = 1.2
  ) +
  geom_sf(data = land,
          color = "grey50",
          fill = NA,
          size = 0.4
  ) +
  # using the streamline geom from the awesome
  # metR package
  metR::geom_streamline(
    data = wind,
    aes(
      x = lon,
      y = lat,
      dx = u, # define east west displacement
      dy = v, # define north south displacement
      alpha = ..step..
    ),
    color = "grey",
    size = 0.4,
    L = 2,
    res = 2,
    n = 30,
    arrow = NULL,
    lineend = "round",
    inherit.aes = FALSE) +
  scale_alpha(guide = "none") +
  scale_fill_divergent_discretised(
    high = "#7f3b08",
    mid = "#f7f7f7",
    low = "#2d004b",
    name = "mean sea leavel pressure (hPa) \n \n ",
    midpoint = 984
  ) +
  labs(
    title = "Hurricane Katrina",
    subtitle = "peak intensity at August 28, 2005",
    x = "",
    y = ""
  ) +
  coord_sf(xlim = c(-95, -75),
           ylim = c(36, 23)) +
  theme_map()

# print the plot
print(p)

That’s it! Much of the credit goes to Elio Campitelli’s awesome metR package for providing streamline and other meteorological visualization aids within the ggplot plotting environment. As always, please provide credit where credit is due if your work relies heavily on this package.

Keep an eye on our blog for new entries in this series on map building, and tips and tricks to improve your R mapping skills. The code for this project can be found on our github R tutorials page. If unfamiliar with the package, the ecmwfr package can be found here.

# cleaned up theme settings
theme_map <- function(...) {
  theme_minimal() +
    theme(
      legend.key = element_rect(fill = "white", size = 1),
      legend.key.size = unit(0.5, "cm"),
      legend.key.width = unit(2.5, "cm"),
      legend.margin = margin(10, 10, 10, 10),
      legend.text = element_text(
        size = 10,
        angle = 45,
        vjust = 1,
        hjust = 1),
      legend.position = "bottom",
      axis.line = element_blank(),
      panel.grid.major = element_blank(),
      panel.grid.minor = element_blank(),
      plot.background = element_rect(fill = "#ffffff", color = NA),
      panel.background = element_rect(fill = "#ffffff", color = NA),
      legend.background = element_rect(fill = "#ffffff", color = NA),
      text = element_text(family = "Montserrat", color = "grey30", size=12),
      panel.border = element_blank(),
      plot.margin = margin(5, 0, 5, 0),
      axis.title.y = element_text(face = "bold", vjust = 10),
      axis.title.x = element_text(face = "bold", vjust = -2),
      panel.spacing = unit(1, "lines"),
      ...
    )
}
Avatar
Koen Hufkens, PhD
Partner, Researcher

As an earth system scientist and ecologist I model ecosystem processes.

Related

Next
Previous