harp Point verification script

Copy - Paste the script below. You will need to set the dates, fcst_model, paths etc. Save it and run with source("<script_file_name.R>") or click on Source in RStudio.

library(harp)
library(purrr)

# The following variables are expected to change for each user / use case

dttm       <- seq_dttm(2018030800, 2018033100, "1d")
fcst_model <- c("REF", "SPP")
lead_time  <- seq(0, 48, 3)
fcst_type  <- "eps"
lags       <- "0s"
fcst_path  <- "/path/to/FCTABLE"
obs_path   <- "/path/to/OBSTABLE"
verif_path <- "/path/to/verif"
grps       <- "lead_time"

# Add more parameters below using the same format. For upper air parameters,
# don't forget the vc = "pressure"
params <- list(
  
  # Parameter that won't work!
  apple = list(
    thresholds = seq(5, 50, 5)
  ),
  
  T2m = list(
    thresholds = c(-20, -10, seq(-5, 25, 5)),
    scale_fcst = list(scaling = -273.15, new_units = "degC"),
    scale_obs  = list(scaling = -273.15, new_units = "degC")
  ), 
  
  S10m = list(
    thresholds = seq(0, 25, 5)
  ), 
  
  T = list(
    scale_fcst = list(scaling = -273.15, new_units = "degC"),
    scale_obs  = list(scaling = -273.15, new_units = "degC"),
    vc         = "pressure"
  )
  
)

# Function that runs the verification
run_verif <- function(prm_info, prm_name) {
  
  if (!is.null(prm_info$vc)) {
    vertical_coordinate <- prm_info$vc
  } else {
    vertical_coordinate <- NA_character_
  }
  
  # Read the forecast
  fcst <- read_point_forecast(
    dttm          = dttm,
    fcst_model    = fcst_model,
    fcst_type     = fcst_type,
    parameter     = prm_name,
    lead_time     = lead_time,
    lags          = lags,
    file_path     = fcst_path,
    vertical_coordinate = vertical_coordinate
  )
  
  # Find the common cases - for upper air parmeters we need to ensure 
  # that the level column  is included in the check
  fcst <- switch(
    vertical_coordinate,
    "pressure" = common_cases(fcst, p),
    "height"   = common_cases(fcst, z),
    common_cases(fcst)
  )
  
  # optional rescaling of forecasts using the scale_fcst part of the
  # params list. We use do.call to call the scale_point_forecast 
  # function with a named list containing the arguments. 
  if (!is.null(prm_info$scale_fcst)) {
    fcst <- do.call(
      scale_param, 
      c(list(x = fcst), prm_info$scale_fcst))
  }
  
  # Read the observations getting the dates and stations from 
  # the forecast
  obs <- read_point_obs(
    dttm       = unique_valid_dttm(fcst),
    parameter  = prm_name,
    obs_path   = obs_path,
    stations   = unique_stations(fcst),
    vertical_coordinate = vertical_coordinate
  )
  
  # optional rescaling of observations using the scale_obs part of the
  # params list. We use do.call to call the scale_point_forecast 
  # function with a named list containing the arguments.
  if (!is.null(prm_info$scale_obs)) {
    obs <- do.call(
      scale_point_obs, 
      c(list(x = obs, col = {{prm_name}}), prm_info$scale_obs)
    )
  }
  
  # Join observations to the forecast
  fcst <- join_to_fcst(fcst, obs)
  
  # Check for errors removing obs that are more than a certain number 
  # of standard deviations from the forecast. You could add a number 
  # of standard deviations to use in the params list 
  fcst <- check_obs_against_fcst(fcst, {{prm_name}})
  
  # Make sure that grps is a list so that it adds on the vertical 
  # coordinate group correctly
  if (!is.list(grps)) {
    grps <- list(grps)
  }
  
  grps <- switch(
    vertical_coordinate,
    "pressure" = map(grps, ~c(.x, "p")),
    "height"   = map(grps, ~c(.x, "z")),
    grps
  )
  
  # Do the verification
  if (fcst_type == "eps") {
    verif <- ens_verify(
      fcst, prm_name, thresholds = prm_info$thresholds, groupings = grps
    )
  } else {
    verif <- det_verify(
      fcst, prm_name, thresholds = prm_info$thresholds, groupings = grps
    )
  }
  
  # Save the scores
  save_point_verif(verif, verif_path = verif_path)
  
  # Return the data to the calling environment (normally global)
  verif
  
}

# Use possibly from the purrr package to allow the script to continue
# if it fails for a parameter - it returns NULL if it fails. See
# ?safely and ?quietly if you want to retain the errors.
possible_run_verif <- possibly(run_verif, otherwise = NULL)

# Use imap from the purrr package to map each element of the params list
# to the possible_run_verif function. imap passes the element of the list
# as the first argument and the name of the element as the second.
verif <- imap(params, possible_run_verif)

# You can open the results in a shiny app using 
# shiny_plot_point_verif(verif_path)