-
-
Save BlaiseKelly/66db4312548f41f17366849a3229ce90 to your computer and use it in GitHub Desktop.
Run hysplit function
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
#' Function to command HYSPLIT (Hybrid Single Particle Lagrangian Integrated | |
#' Trajectory Model) from R. | |
#' | |
#' \code{run_hysplit} will return a data frame which is ready for use in | |
#' \strong{openair}'s \code{traj*} functions. | |
#' | |
#' @param latitude Latitude of a receptor location. | |
#' | |
#' @param longitude Longitude of a receptor location. | |
#' | |
#' @param start What is the start date of data to be returned? Ideally, the | |
#' date format should be \code{yyyy-mm-dd}, but the UK locale convention of | |
#' \code{dd/mm/yyyy} will also work. Years as strings or integers work too and | |
#' will floor-rounded. | |
#' | |
#' @param end What is the end date of data to be returned? Ideally, the | |
#' date format should be \code{yyyy-mm-dd}, but the UK locale convention of | |
#' \code{dd/mm/yyyy} will also work. Years as strings or integers work too and | |
#' will be ceiling-rounded. | |
#' | |
#' @param interval What interval should the trajectories be run at? Default is | |
#' \code{"3 hour"}, as in \strong{openair}. | |
#' | |
#' @param runtime How many hours from time zero should the trajectory be run? | |
#' Default is \code{-96}, i.e., a back trajectory of 96 hours, as in | |
#' \strong{openair}. | |
#' | |
#' @param start_height Height of receptor location. Default is \code{10} metres, | |
#' as in \strong{openair}. | |
#' | |
#' @param model_height Height of model domain. Default is \code{10000} metres, | |
#' as in \strong{openair}. | |
#' | |
#' @param hysplit_exec Location of HYSPLIT's executable files. Ensure executable | |
#' permisions have been set for the \code{hyts_std} application if on a Unix | |
#' system. HYSPLIT's \code{bdyfiles} directory with the \code{ASCDATA.CFG} will | |
#' also need to be in the correct location. | |
#' | |
#' @param hysplit_input Location of input meteorological files for HYSPLIT. | |
#' | |
#' @param hysplit_output Location of where HYSPLIT should write its trajectory | |
#' outputs while working. | |
#' | |
#' @param delete Should temporary files be silently deleted before and after the | |
#' model runs? Default is \code{FALSE}. | |
#' | |
#' @param drop Should the \code{"year"}, \code{"month"}, \code{"day"}, | |
#' \code{"hour"}, and \code{"receptor"} variables be dropped from the bound | |
#' object? Default is \code{TRUE}. | |
#' | |
#' @param site An optional site string to be added to the returned data frame. | |
#' | |
#' @param source An optional source string to be added to the returned data frame. | |
#' | |
#' @param verbose Should the function give messages on what trajectory is being | |
#' processed. Default is \code{FALSE} and will estimate the time until | |
#' completion. | |
#' | |
#' @seealso \code{\link{importTraj}}, \code{\link{trajPlot}}, | |
#' \url{http://ready.arl.noaa.gov/HYSPLIT.php} | |
#' | |
#' @author Stuart K. Grange and David Carslaw | |
#' | |
#' @import lubridate | |
#' @import stringr | |
#' | |
#' @examples | |
#' \dontrun{ | |
#' | |
#' | |
#' # Run back trajectories for first 10 days at 1500 metres and extended run time for | |
#' # analysis of air aloft or Gibraltar | |
#' data_gibraltar <- run_hysplit( | |
#' latitude = 36.134, | |
#' longitude = -5.347, | |
#' runtime = -96, | |
#' start_height = 1500, | |
#' model_height = 10000, | |
#' start = 2015, | |
#' end = "2015-01-10", | |
#' hysplit_exec = "~/Hysplit/exec", | |
#' hysplit_input = "~/TrajData", | |
#' hysplit_output = "~/temp", | |
#' site = "gibraltar") | |
#' | |
#' # On a Windows system with a default hysplit installation | |
#' data_gibraltar <- run_hysplit( | |
#' latitude = 36.134, | |
#' longitude = -5.347, | |
#' runtime = -96, | |
#' start_height = 10, | |
#' model_height = 10000, | |
#' start = 2015, | |
#' end = "2015-01-10", | |
#' hysplit_exec = "~/hysplit4/exec", | |
#' hysplit_input = "~/trajData", | |
#' hysplit_output = "~/temp", | |
#' site = "gibraltar") | |
#' | |
#' } | |
#' | |
#' @export | |
run_hysplit <- function (latitude = 51.5, | |
longitude = -0.1, | |
start = NA, | |
end = NA, | |
interval = "3 hour", | |
runtime = -96, | |
start_height = 10, | |
model_height = 10000, | |
hysplit_exec = "exec/", | |
hysplit_input, | |
hysplit_output = "hysplit_output/", | |
delete = TRUE, | |
drop = TRUE, | |
site = NA, | |
source = NA, | |
verbose = FALSE) { | |
require(stringr) | |
require(lubridate) | |
require(dplyr) | |
require(purrr) | |
require(readr) | |
# Parse arguments | |
# Expand paths | |
hysplit_exec <- path.expand(hysplit_exec) | |
hysplit_input <- path.expand(hysplit_input) | |
hysplit_output <- path.expand(hysplit_output) | |
# check these places exist | |
if (file.access(hysplit_exec, mode = 0) != 0) | |
stop(paste("File path", hysplit_exec, "does not exist")) | |
if (file.access(hysplit_input, mode = 0) != 0) | |
stop(paste("File path", hysplit_input, "does not exist")) | |
if (file.access(hysplit_output, mode = 0) != 0) | |
stop(paste("File path", hysplit_output, "does not exist")) | |
if (file.access(hysplit_exec, mode = 2) != 0) | |
stop(paste("File path", hysplit_exec, "is not writable")) | |
# write default config file that returns met data | |
set_config() %>% | |
write_config_list(dir = hysplit_exec) | |
# Add final path separator | |
hysplit_exec <- str_c(hysplit_exec, .Platform$file.sep) | |
hysplit_input <- str_c(hysplit_input, .Platform$file.sep) | |
hysplit_output <- str_c(hysplit_output, .Platform$file.sep) | |
# Selection after testing if directory is empty | |
files_old <- list.files(hysplit_output, full.names = TRUE, | |
pattern = "hysplit_output.txt") | |
# Delete old files | |
if (delete) { | |
# Delete | |
message("Deleting old files...") | |
file.remove(files_old) | |
} else { | |
# Selection control | |
while (length(files_old) > 1) { | |
# Get input | |
input <- readline("There are hysplit output files in the output directory. \nShould these be deleted (y/n)? \n") | |
# Parse | |
input <- str_to_upper(input) | |
input <- ifelse( | |
input %in% c("YES", "Y", "T", "TRUE", "TR", "TRU"), TRUE, FALSE) | |
if (input) { | |
message("Deleting old files...") | |
file.remove(files_old) | |
break | |
} else { | |
break | |
} | |
} | |
} | |
# Dates | |
# Start and end dates | |
start <- parse_date_arguments(start, "start") | |
end <- parse_date_arguments(end, "end") | |
# Receptor location and starting height | |
coordinates <- str_c(latitude, longitude, start_height, sep = " ") | |
# Set up where the control file is to be written | |
control_file <- file.path(hysplit_exec, "CONTROL") | |
# Store working directory because this will be changed | |
wd <- getwd() | |
# Create date sequence | |
date_sequence <- seq(start, end, interval) | |
# Apply function which runs the model multiple times | |
message(str_c("Running ", length(date_sequence), " HYSPLIT trajectories...")) | |
# For plyr's progress bar | |
#progress <- ifelse(verbose, "none", "time") | |
# Run model | |
purrr::walk(date_sequence, run_trajectory, | |
hysplit_exec = hysplit_exec, | |
hysplit_input = hysplit_input, | |
hysplit_output = hysplit_output, | |
control_file = control_file, | |
coordinates = coordinates, | |
runtime = runtime, | |
model_height = model_height) | |
# Change working directory back to original after system calls | |
setwd(wd) | |
# Bind output files | |
message("Binding HYSPLIT files...") | |
# Get file list | |
file_list <- list.files(hysplit_output, "hysplit_output.txt", full.name = TRUE) | |
# Load files as in openair, but drop some things usually | |
df <- purrr::map_df(file_list, read_hysplit_file, drop = drop) | |
# Add variables which are not in openair's files | |
df$start_height <- start_height | |
if (!is.na(site)) df$site <- site | |
if (!is.na(source)) df$source <- source | |
# Delete files | |
if (delete) file.remove(file_list) | |
# Return | |
df | |
} | |
# Define the function which creates a control file and calls the hy_std | |
# application. | |
# | |
# No export | |
run_trajectory <- function (date, | |
hysplit_exec, | |
hysplit_input, | |
hysplit_output, | |
control_file, | |
coordinates, | |
runtime, | |
model_height) { | |
# Get pieces of the date | |
date_year <- year(date) | |
date_month <- month(date) | |
date_day <- day(date) | |
date_hour <- hour(date) | |
# Pad zeros | |
date_month <- str_pad(date_month, width = 2, pad = "0") | |
date_day <- str_pad(date_day, width = 2, pad = "0") | |
date_hour <- str_pad(date_hour, width = 2, pad = "0") | |
# Format for control file | |
date_control <- str_c(date_year, date_month, date_day, date_hour, sep = " ") | |
# Use date to create a file name | |
file_name_export <- str_c( | |
str_replace_all(date_control, " ", ""), "_hysplit_output.txt") | |
# Get date string | |
current_month_pattern <- str_c(date_year, date_month) | |
past_month_pattern <- get_year_and_month(current_month_pattern, - 1) | |
future_month_pattern <- get_year_and_month(current_month_pattern, 1) | |
# Add other pieces of the file names | |
past_month_pattern <- str_c("RP", past_month_pattern , ".gbl") | |
current_month_pattern <- str_c("RP", current_month_pattern , ".gbl") | |
future_month_pattern <- str_c("RP", future_month_pattern , ".gbl") | |
# Create the file and directory list for the control file | |
file_list <- c(past_month_pattern, current_month_pattern, future_month_pattern) | |
# For back trajectories, do not use the future month | |
if (runtime < 0) { | |
# But only if not the final day of the month, otherwise a stop error occurs | |
if (!floor_date(date, "day") == ceiling_date(date, "month") - days(1)) | |
file_list <- file_list[1:2] | |
} | |
# Add directory to file names | |
file_dir_list <- str_c(hysplit_input, file_list, sep = "\n") | |
# Write control file | |
# This will replace the contents of the current file if it exists | |
write_to_control_file(date_control, control_file, append = FALSE) | |
# Starting locations | |
write_to_control_file("1", control_file) | |
# Write coordinates and starting height of model | |
write_to_control_file(coordinates, control_file) | |
# Write runtime of model, hours forward or backwards for trajectories | |
write_to_control_file(runtime, control_file) | |
# Vertical motion option, top of model, and input grids (number of files) | |
write_to_control_file(str_c("0\n", model_height, "\n", length(file_list)), | |
control_file) | |
# Write input directory and file names | |
write_to_control_file(file_dir_list, control_file) | |
# Output directory | |
write_to_control_file(hysplit_output, control_file) | |
# Output file | |
write_to_control_file(file_name_export, control_file) | |
# Change working directory to hysplit application | |
setwd(hysplit_exec) | |
# Message file name | |
# if (verbose) message(file_name_export) | |
# System call, quiet arguments are os specific | |
if (.Platform$OS.type == "windows") { | |
# Run | |
system("./hyts_std", show.output.on.console = FALSE) | |
} else { | |
# Run, ensure executable permissions are set | |
system("./hyts_std", ignore.stdout = TRUE) | |
} | |
} | |
# Three files are needed for the model to create a back trajectory, | |
# the current month, previous month, and future month | |
get_year_and_month <- function (pattern, difference = 1) { | |
# Do some date things | |
date <- ymd(str_c(pattern, "01")) + months(difference) | |
year <- year(date) | |
month <- month(date) | |
month <- str_pad(month, 2, pad = "0") | |
# Combine | |
pattern <- str_c(year, month) | |
# Return | |
pattern | |
} | |
# Write table function | |
write_to_control_file <- function(string, file, append = TRUE) { | |
# Write to file | |
write.table(string, file, col.names = FALSE, row.names = FALSE, quote = FALSE, | |
append = append) | |
} | |
# Function to read hysplit files | |
# | |
# Taken from openair manual to keep bound trajectory files in the same format | |
# as the openair package. | |
# | |
# No export | |
read_hysplit_file <- function (file, drop) { | |
# Load file, error catching is for when two or three input met files are used | |
# and results in a different length file header | |
tryCatch({ | |
if (isTRUE(grep("PRESSURE", read_lines(file)[6]) == 1)) { | |
df <- read.table(file, header = FALSE, skip = 6) | |
} | |
if (isTRUE(grep("PRESSURE", read_lines(file)[7]) == 1)) { | |
df <- read.table(file, header = FALSE, skip = 7) | |
} | |
}, error = function (e) { | |
return() | |
} | |
) | |
if (is.null(nrow(df))) return() | |
# Drop | |
# df <- subset(df, select = -c(V2, V7, V8)) | |
df <- df %>% | |
select(!any_of(c("V2", "V7", "V8"))) | |
# Rename | |
df <- dplyr::rename(df, c(receptor = V1, year = V3, month = V4, | |
day = V5, hour = V6, hour.inc = V9, | |
lat = V10, lon = V11, | |
height = V12, pressure = V13, | |
theta = V14, air_temp = V15, | |
rainfall = V16, mixdepth = V17, | |
rh = V18, sp_humidity = V19, | |
h2o_mixrate = V20, terr_msl = V21, | |
sun_flux = V22)) | |
# Clean two digit years | |
df$year <- ifelse(df$year < 50, df$year + 2000, df$year + 1900) | |
# Transform pieces of date to date | |
df$date2 <- with(df, ISOdatetime(year, month, day, | |
hour, min = 0, sec = 0, | |
tz = "GMT")) | |
# Drop variables | |
if (drop) df <- subset(df, select = -c(year, | |
month, | |
day, | |
hour, | |
receptor)) | |
# Transform arrival time, minus hours from hour.inc variable | |
df$date <- df$date2 - 3600 * df$hour.inc | |
# Return | |
df | |
} | |
# No export | |
parse_date_arguments <- function (date, what) { | |
if (what == "start") { | |
if (!is.na(date) & nchar(date) == 4) | |
date <- stringr::str_c(date, "-01-01") | |
date <- ifelse(is.na(date), as.character(lubridate::floor_date(Sys.Date(), | |
"year")), date) | |
} | |
if (what == "end") { | |
if (!is.na(date) & nchar(date) == 4) | |
date <- stringr::str_c(date, "-12-31 18:00") | |
date <- ifelse(is.na(date), as.character(lubridate::ceiling_date(Sys.Date(), | |
"year")), date) | |
} | |
date <- lubridate::parse_date_time(date, c("ymd", "dmy", "dmy_HM", "ymd_HM")) | |
date | |
} | |
# Hysplit config settings from https://github.com/rich-iannone/splitr, modified to return met variables by default | |
#' Create a configuration list for a trajectory or dispersion model | |
#' | |
#' @param tratio The advection stability ratio. By default, this is set to | |
#' `0.75`. | |
#' @param initd The initial distribution. Set to `0` by default. | |
#' @param kpuff The horizontal puff dispersion growth. Setting to `0` results in | |
#' linear growth whereas option `1` uses an empirical growth scheme. By | |
#' default, this is set to `0` (linear growth). | |
#' @param khmax The maximum duration (in hours) for a particle or trajectory. | |
#' This is set to the absolute maximum by default, which is `9999` hours. | |
#' @param kmixd The methodology for modeling the mixed layer depth. There are | |
#' three options: `0` for using the input (the default), `1` for using | |
#' temperature, and `2` for using the TKE. | |
#' @param kmix0 The minimum mixing depth in meters. This is `250` by default. | |
#' @param kzmix How to perform vertical mixing adjustments. There are three | |
#' options: `0` for not doing this at all (the default), `1` uses the PBL | |
#' average, and `3` scales with *TVMIX*. | |
#' @param kdef The modeling method for horizontal turbulence. We can use the | |
#' vertical with option `0` (the default), or, use the deformation method with | |
#' option `1`. | |
#' @param kbls How to derive the boundary layer stability. Two options are | |
#' available for this. We can use fluxes (option `1`, default) or use wind | |
#' temperature (option `2`). | |
#' @param kblt The boundary layer turbulence parameterization to use. This can | |
#' either be Beljaars (option `1`), Kanthar (option `2`, default), or TKE | |
#' (option `3`). | |
#' @param conage Particle to- or from-puff conversion at *CONAGE*. In units of | |
#' hours, with a default of `48`. | |
#' @param numpar The number of puffs or particles released per cycle. By default | |
#' this is set to `2500`. | |
#' @param qcycle An optional cycling of emissions, in units of hours. | |
#' @param efile An absolute path to an optional temporal emissions file. | |
#' @param tkerd The unstable turbulent kinetic energy ratio. This is set to | |
#' `0.18`. | |
#' @param tkern The stable turbulent kinetic energy ratio. This is set to | |
#' `0.18`. | |
#' @param ninit How to do particle initialization; with `0` no particle | |
#' initialization is done, with `1` this is done once, options `2` and `3` are | |
#' the *add* and *replace* methods. | |
#' @param ndump Should the particles be dumped to a file and, if so, how often? | |
#' Using `0` results in no writing particles to file (the default), and any | |
#' non-zero value performs this writing once per number of hours specified. | |
#' @param ncycl The *PARDUMP* output cycle time. | |
#' @param pinpf The particle input file name (default is `"PARINIT"`). This is | |
#' useful for initialization or boundary conditions. | |
#' @param poutf The particle output file name (default is `"PARDUMP"`). | |
#' @param mgmin The minimum meteorological subgrid size. The default value is | |
#' `10`. | |
#' @param kmsl The starting height reference. If it is to signify a distance | |
#' above ground level (AGL) then use option `0` (the default). If it is | |
#' instead relative to mean sea level, then option `1` should be used. | |
#' @param maxpar The maximum number of particles to be carried in simulation. By | |
#' default this is `10000`. | |
#' @param cpack The binary concentration packing. Here are the options: `0` for | |
#' none, `1` for nonzero, `2` for points, and `3` for polar. Option `1` is the | |
#' default option. | |
#' @param cmass Informs grid computation. The two options are to compute grid | |
#' concentrations (with `0`, the default) or to compute grid mass (option | |
#' `1`). | |
#' @param dxf,dyf The horizontal x- and y-grid adjustment factors for an | |
#' ensemble. By default, these are both set to `1.0`. | |
#' @param dzf The vertical factor for an ensemble. This is `0.01` by default. | |
#' @param ichem The chemistry conversion module to employ. Option `0` does no | |
#' chemistry (the default). With option `1` we use a matrix method, option `2` | |
#' does conversion, and option `3` works on dust. | |
#' @param maxdim The maximum number of pollutants to carry on one particle. By | |
#' default, this is `1`. | |
#' @param kspl The standard splitting interval in units of hours. By default, | |
#' this is `1` hour. | |
#' @param krnd The enhanced merge interval in hours. By default this is `6` | |
#' hours. | |
#' @param frhs The standard horizontal puff rounding fraction for the merge | |
#' process. This is `1.0` by default. | |
#' @param frvs The vertical puff rounding fraction, which is `0.01` by default. | |
#' @param frts The temporal puff rounding fraction, which is `0.10` by default. | |
#' @param frhmax The maximum value for the horizontal rounding parameter. This | |
#' is `3.0` by default. | |
#' @param splitf The automatic size adjustment factor for horizontal splitting. | |
#' By default, the value is `1.0`. | |
#' @param tm_pres,tm_tpot,tm_tamb,tm_rain,tm_mixd,tm_relh,tm_sphu,tm_mixr,tm_dswf,tm_terr | |
#' Options to include meteorology along trajectory points. These are the | |
#' pressure variable marker flag (`tm_pres`), the potential temperature | |
#' (`tm_tpot`), the ambient temperature (`tm_tamb`), the rainfall rate | |
#' (`tm_rain`), the mixed layer depth (`tm_mixd`), the relative humidity | |
#' (`tm_relh`), the specific humidity (`tm_sphu`), the mixing rate | |
#' (`tm_mixr`), the downward short-wave flux (`tm_dswf`), and the terrain | |
#' height (`tm_terr`). Setting any of these to `0` disables output, whereas | |
#' `1` will enable output of these data points. By default, all are set to | |
#' `0`. | |
#' | |
#' @export | |
set_config <- function(tratio = 0.75, | |
initd = 0, | |
kpuff = 0, | |
khmax = 9999, | |
kmixd = 0, | |
kmix0 = 25, | |
kzmix = 0, | |
kdef = 0, | |
kbls = 1, | |
kblt = 2, | |
conage = 48, | |
numpar = 2500, | |
qcycle = 0.0, | |
efile = NULL, | |
tkerd = 0.18, | |
tkern = 0.18, | |
ninit = 1, | |
ndump = 1, | |
ncycl = 1, | |
pinpf = "PARINIT", | |
poutf = "PARDUMP", | |
mgmin = 10, | |
kmsl = 0, | |
maxpar = 10000, | |
cpack = 1, | |
cmass = 0, | |
dxf = 1.0, | |
dyf = 1.0, | |
dzf = 0.01, | |
ichem = 0, | |
maxdim = 1, | |
kspl = 1, | |
krnd = 6, | |
frhs = 1.0, | |
frvs = 0.01, | |
frts = 0.10, | |
frhmax = 3.0, | |
splitf = 1.0, | |
tm_pres = 1, | |
tm_tpot = 1, | |
tm_tamb = 1, | |
tm_rain = 1, | |
tm_mixd = 1, | |
tm_relh = 1, | |
tm_sphu = 1, | |
tm_mixr = 1, | |
tm_dswf = 1, | |
tm_terr = 1) { | |
arg_names <- formals(set_config) %>% names() | |
arg_vals <- mget(arg_names) | |
if (is.null(arg_vals$efile)) { | |
arg_vals$efile <- "''" | |
} else if (!is.null(arg_vals$efile)) { | |
arg_vals$efile <- paste0("'", arg_vals$efile, "'") | |
} | |
if (is.null(arg_vals$pinpf)) { | |
arg_vals$pinpf <- "''" | |
} else if (!is.null(arg_vals$pinpf)) { | |
arg_vals$pinpf <- paste0("'", arg_vals$pinpf, "'") | |
} | |
if (is.null(arg_vals$poutf)) { | |
arg_vals$poutf <- "''" | |
} else if (!is.null(arg_vals$poutf)) { | |
arg_vals$poutf <- paste0("'", arg_vals$poutf, "'") | |
} | |
arg_vals[!vapply(arg_vals, FUN = is.null, FUN.VALUE = logical(1))] | |
} | |
write_config_list <- function(config_list, dir) { | |
paste0( | |
"&SETUP\n", | |
paste0(names(config_list), " = ", config_list, ",\n", collapse = ""), | |
"/\n" | |
) %>% | |
cat(file = file.path(dir, "SETUP.CFG")) | |
} | |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment