Add hexbin computation

NAMESPACE

@@ -2,6 +2,10 @@
 
 S3method("[",harp_fcst)
 S3method(arrange,harp_fcst)
+S3method(bin_fcst_obs,harp_det_point_df)
+S3method(bin_fcst_obs,harp_ens_point_df)
+S3method(bin_fcst_obs,harp_ens_point_df_long)
+S3method(bin_fcst_obs,harp_list)
 S3method(bind_point_verif,harp_verif)
 S3method(bind_point_verif,list)
 S3method(c,harp_fcst)
@@ -56,6 +60,7 @@ S3method(spread_members,default)
 S3method(spread_members,harp_fcst)
 S3method(transmute,harp_fcst)
 export("%>%")
+export(bin_fcst_obs)
 export(bind_bootstrap_score)
 export(bind_point_verif)
 export(bind_rows.harp_fcst)

R/det_verify.R

@@ -1,23 +1,6 @@
 #' Compute verification scores for deterministic forecasts.
 #'
-#' @param .fcst A \code{harp_fcst} object with tables that have a column for
-#'   observations, or a single forecast table.
-#' @param parameter The name of the column for the observed data.
-#' @param thresholds A numeric vector of thresholds for which to compute the
-#'   threshold based scores. Set to NULL (the default) to only compute summary
-#'   scores.
-#' @param groupings The groups for which to compute the scores. See
-#' \link[dplyr]{group_by} for more information of how grouping works.
-#' @param circle If set the parameter is assumed to be cyclic for bias
-#'   calculations. Should be this distance around a circle in the units of the
-#'   parameter, so would typically have a value of 360 for degrees or `2 * pi`
-#'   for radians.
-#' @param show_progress Logical - whether to show a progress bar. The default is
-#' `TRUE`
-#' @param ... Reserved for methods
-#'   `TRUE`.
-#' @return A list containting two data frames: \code{det_summary_scores} and
-#'   \code{det_threshold_scores}.
+#' @inheritParams ens_verify
 #' @export
 #'
 #' @examples
@@ -27,6 +10,7 @@ det_verify <- function(
   thresholds    = NULL,
   groupings     = "lead_time",
   circle        = NULL,
+  num_bins      = 30,
   show_progress = TRUE,
   ...
 ) {
@@ -48,6 +32,7 @@ det_verify.harp_ens_point_df <- function(
   thresholds    = NULL,
   groupings     = "lead_time",
   circle        = NULL,
+  num_bins      = 30,
   show_progress = TRUE,
   fcst_model    = NULL,
   ...
@@ -62,22 +47,22 @@ det_verify.harp_ens_point_df <- function(
   groupings <- purrr::map(groupings, ~union(c("sub_model", "member"), .x))
 
   det_verify(
-    .fcst, {{parameter}}, thresholds, groupings, circle,
+    .fcst, {{parameter}}, thresholds, groupings, circle, num_bins,
     show_progress, fcst_model
   )
 }
 
 
-#' @param fcst_model The name of the forecast model to use in the `fcst_model`
-#'  column of the output. If the function is dispatched on a `harp_list`
-#'  object, the names of the `harp_list` are automatically used.
+#' @rdname det_verify
+#' @inheritParams ens_verify.harp_ens_point_df
 #' @export
 det_verify.harp_det_point_df <- function(
   .fcst,
   parameter,
   thresholds    = NULL,
   groupings     = "lead_time",
   circle        = NULL,
+  num_bins      = 30,
   show_progress = TRUE,
   fcst_model    = NULL,
   ...
@@ -184,15 +169,25 @@ det_verify.harp_det_point_df <- function(
     fcst_df
   }
 
-  res <- list()
-
-  res[["det_summary_scores"]] <- purrr::map(
+  det_summary_scores <- list()
+  det_summary_scores[["basic"]] <- purrr::map(
     groupings, compute_summary_scores, .fcst
   ) %>%
     purrr::list_rbind() %>%
     fill_group_na(groupings) %>%
     dplyr::mutate(fcst_model = fcst_model, .before = dplyr::everything())
 
+  det_summary_scores[["hexbin"]] <- bin_fcst_obs(
+    .fcst, !!parameter, groupings, num_bins, show_progress
+  )[["det_summary_scores"]]
+
+  res <- list()
+
+  res[["det_summary_scores"]] <- Reduce(
+    function(x, y) suppressMessages(dplyr::inner_join(x, y)),
+    det_summary_scores
+  )
+
   if (is.numeric(thresholds)) {
 
     meta_cols     <- grep(
@@ -307,6 +302,7 @@ det_verify.harp_list <- function(
   thresholds    = NULL,
   groupings     = "lead_time",
   circle        = NULL,
+  num_bins      = 30,
   show_progress = TRUE
 ) {
 
@@ -322,7 +318,7 @@ det_verify.harp_list <- function(
     purrr::imap(
       .fcst,
       ~det_verify(
-        .x, {{parameter}}, thresholds, groupings, circle,
+        .x, {{parameter}}, thresholds, groupings, circle, num_bins,
         show_progress, fcst_model = .y
       )
     )

R/ens_verify.R

@@ -1,7 +1,7 @@
 #' Compute all verification scores for an ensemble.
 #'
-#' @param .fcst A \code{harp_fcst} object with tables that have a column for
-#'   observations, or a single forecast table.
+#' @param .fcst A `harp_df` or `harp_list` object with tables that have a column
+#'   for observations, or a single forecast table.
 #' @param parameter The name of the column for the observed data.
 #' @param verify_members Whether to verify the individual members of the
 #'   ensemble. Even if thresholds are supplied, only summary scores are
@@ -35,12 +35,16 @@
 #'   elements eps_model and member to use a member of an eps model in the
 #'   harp_fcst object for the climatology, or a data frame with columns for
 #'   threshold and climatology and also optionally lead_time.
+#' @param hexbin Logical. Whether to compute hexbins for forecast, observation
+#'   pairs. Defaults to `TRUE`. See \link{bin_fcst_obs} for more details.
+#' @param num_bins The number of bins into which to partition observations for
+#'   the hexbin computation.
 #' @param rank_hist Logical. Whether to compute the rank histogram. Defaults to
 #'   `TRUE`. Note that the computation of the rank histogram can be slow if
 #'   there is a large number (> 1000) of groups.
 #' @param crps Logical. Whether to compute the CRPS. Defaults to `TRUE`.
-#' @param brier Logical. Whether to compute the Brier score. Defaults to
-#'   `TRUE`. Will be ignored if no thresholds are set.
+#' @param brier Logical. Whether to compute the Brier score. Defaults to `TRUE`.
+#'   Will be ignored if no thresholds are set.
 #' @param roc Logical. Whether to compute the Relative Operating Characteristic
 #'   (ROC). Defaults to `TRUE`. Will be ignored if no thresholds are set.
 #' @param econ_val Logical. Whether to compute the economic value. Defaults to
@@ -66,6 +70,8 @@ ens_verify <- function(
   spread_drop_member = NULL,
   jitter_fcst        = NULL,
   climatology        = "sample",
+  hexbin             = TRUE,
+  num_bins           = 30,
   rank_hist          = TRUE,
   crps               = TRUE,
   brier              = TRUE,
@@ -101,6 +107,8 @@ ens_verify.harp_ens_point_df <- function(
   spread_drop_member = NULL,
   jitter_fcst        = NULL,
   climatology        = "sample",
+  hexbin             = TRUE,
+  num_bins           = 30,
   rank_hist          = TRUE,
   crps               = TRUE,
   show_progress      = TRUE,
@@ -170,9 +178,20 @@ ens_verify.harp_ens_point_df <- function(
       spread_drop_member = spread_drop_member
     )[["ens_summary_scores"]]
 
+    if (hexbin) {
+      ens_summary_scores[["hexbin"]] <- dplyr::select(
+        bin_fcst_obs(
+          .fcst, !!parameter, groupings = groupings, num_bins = num_bins,
+          show_progress = show_progress
+        )[["ens_summary_scores"]],
+        -dplyr::all_of("num_cases")
+      )
+    }
+
     if (rank_hist) {
       ens_summary_scores[["rh"]] <- ens_rank_histogram(
-        .fcst, !! parameter, groupings = groupings
+        .fcst, !! parameter, groupings = groupings,
+        show_progress = show_progress
       )[["ens_summary_scores"]]
     }
 
@@ -284,6 +303,8 @@ ens_verify.harp_list <- function(
   spread_drop_member = NULL,
   jitter_fcst        = NULL,
   climatology        = "sample",
+  hexbin             = TRUE,
+  num_bins           = 30,
   rank_hist          = TRUE,
   crps               = TRUE,
   brier              = TRUE,
@@ -312,7 +333,7 @@ ens_verify.harp_list <- function(
       function(x, y, z) ens_verify(
         x, !!parameter, verify_members, thresholds, groupings, circle,
         rel_probs, num_ref_members, z, jitter_fcst, climatology,
-        rank_hist, crps, brier, roc, econ_val,
+        hexbin, num_bins, rank_hist, crps, brier, roc, econ_val,
         show_progress, fcst_model = y
       )
     )