The HYSPLIT model is quite useful if you want to understand where air has been (backtrajectory) and where it is going (forward trajectory). Backtrajectory analyses, for instance, use interpolated meteorological data to estimate the most likely central path over geographical areas that provided air to a receptor at a given time. Essentially, you’re following a parcel of air backward in hourly steps for a specified length of time.
The R package SplitR makes trajectory modelling with HYSPLIT pretty easy: use one function to run the model (hysplit_trajectory), and use another function to visualize the output (trajectory_plot). To get started, get the most recent version of SplitR from GitHub using devtools::install_github("rich-iannone/SplitR") in R.
The model requires input meteorological data files to determine wind trajectories during a given time and within a given domain. SplitR will download the necessary files (i.e., you won’t have to manually download those files from the NOAA FTP). Furthermore, while performing the HYSPLIT model runs, output files will be generated. Because of all this file management, it’s useful and important to set a working directory to hold input and output files.
# Load in the `SplitR` package
library(SplitR)
# Set the working directory
setwd("~/Documents/wd_splitR")
# Run the Hysplit model, determining the locations
# of air masses going back in time 14 days (336 h) from a
# site located in Vancouver, BC
wind_trajectory_back_14d <-
hysplit_trajectory(
start_lat_deg = 49.263, # the latitude in decimal degrees
start_long_deg = -123.250, # the longitude in decimal degrees
start_height_m_AGL = 50, # site is 50 m above ground level
simulation_duration_h = 336, # model run is 336 h
backtrajectory = TRUE, # the trajectory goes back in time
met_type = "reanalysis", # the meteorology data to use
run_type = "day", # the run is on a single day...
run_day = "2015-07-01", # ...this is the date
daily_hours_to_start = 0, # the run will start at 00:00 UTC
return_met_along_traj = TRUE) # provides extra met data to df
The working directory will now be full of files but, importantly, there will be a data frame called wind_trajectory_back_14d that contains a lot of useful of data. Here is a summary of all the columns and what they signify:
receptor: a numeric label for the receptoryear, month, day, hour: integer values for date/time componentshour.inc: the integer hour difference compared to the run starting timelat, lon, height: the latitude, longitude, and height (meters above ground level) of the air mass along the trajectorypressure: the air pressure along the trajectory (in hPa)date2: a POSIXct date-time value (in UTC) for the air mass along the trajectorydate: a POSIXct date-time value (in UTC) for the time of release or time of incidence at the receptor sitetheta: the potential temperature (in K) along the trajectoryair_temp: the ambient air temperature (in K) along the trajectoryrainfall: the rate of rainfall (in mm/h) along the trajectorymixdepth: the mixing depth (or mixing height, in meters) along the trajectoryrh: the relative humidity (%) along the trajectorysp_humidity: the specific humidity (in g/kg) along the trajectoryh2o_mixrate: the mixed layer depth (in meters) along the trajectoryterr_msl: the terrain height (meters above mean sea level) at the location defined by lat and longsun_flux: the downward solar radiation flux (in watts) along the trajectoryMany of these variables are useful for ad-hoc analyses of the data. The data frame object can also be used directly with the trajectory_plot function for plotting directly on an interactive map:
trajectory_plot(
traj_df = wind_trajectory_back_14d)
Notice that each of the hourly timesteps is rendered as a circle that presents an informative popup when clicked (try zooming in a bit for more accuracy in clicking a specific hourly marker). The path for each trajectory also has its own simplified popup. Both the points and the paths can have their visibility toggled using the controls Layers at the top-right of the display.
For the next example we will run the HYSPLIT model to determine the locations of several air masses going back in time (this time 2 days, or, 48 h) from a site located in Wichita, KN. Each trajectory is initiated every 2 h for 2 days starting at 00:00 UTC on July 1, 2015.
# Run the Hysplit model, determining the locations
# of air masses going back in time (48 h) from a
# site in Wichita
wind_trajectories_back <-
hysplit_trajectory(
start_lat_deg = 37.687959,
start_long_deg = -97.343538,
start_height_m_AGL = 50,
simulation_duration_h = 48,
backtrajectory = TRUE,
met_type = "reanalysis",
run_type = "range",
run_range = c("2015-07-15", "2015-07-16"),
daily_hours_to_start = seq(0, 22, 2),
return_met_along_traj = TRUE)
Plotting the trajectories onto a map:
trajectory_plot(
traj_df = wind_trajectories_back)
With many trajectories, it can be difficult to ascertain their temporal ordering. You can plot with a different color scheme (increasingly_gray), where lighter grays are for trajectories alighting to the orgin further back in time.
trajectory_plot(
traj_df = wind_trajectories_back,
color_scheme = "increasingly_gray")
The HYSPLIT model can be run in the forward direction to determine the locations of air masses leaving the same site located in Wichita. Temporal options are also the same as in the backward-running model. The only change to the function is the use of backtrajectory = FALSE.
# Run the Hysplit model, determining the locations
# of air masses going forward in time (48 h) from a
# site in Wichita
wind_trajectories_forward <-
hysplit_trajectory(
start_lat_deg = 37.687959,
start_long_deg = -97.343538,
start_height_m_AGL = 50,
simulation_duration_h = 48,
backtrajectory = FALSE, # the trajectory now goes forward in time
met_type = "reanalysis",
run_type = "range",
run_range = c("2015-07-15", "2015-07-16"),
daily_hours_to_start = seq(0, 22, 2),
return_met_along_traj = TRUE)
Plotting this with gray colors as before:
trajectory_plot(
traj_df = wind_trajectories_forward,
color_scheme = "increasingly_gray")
Finally, you can combine the forward and backward trajectories and get a map of winds passing through the receptor site at 37.687959ºN, -97.343538ºW, 50 m AGL:
trajectory_plot(
traj_df = rbind(wind_trajectories_back,
wind_trajectories_forward),
color_scheme = "increasingly_gray")
There can certainly be many improvements but the existing functions work pretty well! If you have any comments or suggestions, please do let me know!