style the package

R/cal_heat_stress.R

@@ -0,0 +1,107 @@
+#' cal_heat_stress
+#'
+#' Calculate Heat Stress Levels Based on Climate Projections
+#'
+#' The `cal_heat_stress` function computes daily or monthly grid-level heat stress levels for specified
+#' sectors, using a user-defined heat stress function and climate data. The function
+#' supports both daily and monthly climate inputs and evaluates heat stress over a specific year.
+#'
+#' @param TempRes temporal resolution of input, "day" or "month"
+#' @param SECTOR choice of sector
+#' @param HS choice of heat stress function
+#' @param YEAR_INPUT a vector of years of interest
+#' @param ... individual .nc file name of climate variables in order
+#'
+#' @importFrom raster stack overlay calc stackApply as.matrix ncell nlayers
+#' @importFrom terra rast
+#' @importFrom dplyr %>%
+#' @import assertthat
+#'
+#' @return A raster stack of daily grid-level heat stress level
+#' @export
+#'
+#' @examples
+#' \dontrun{
+#' GD_HS <- cal_heat_stress(
+#'   TempRes = "day", SECTOR = "MAIZ_I", HS = WBGT1, YEAR_INPUT = 2027,
+#'   "hurs_day_GFDL-ESM2M_rcp60_r1i1p1_EWEMBI_20210101-20301231.nc4",
+#'   "tas_day_GFDL-ESM2M_rcp60_r1i1p1_EWEMBI_20210101-20301231.nc4",
+#'   "ps_day_GFDL-ESM2M_rcp60_r1i1p1_EWEMBI_20210101-20301231.nc4"
+#' )
+#' }
+cal_heat_stress <- function(TempRes, SECTOR, HS, YEAR_INPUT, ...) {
+  # check the input temporal resolution
+  if (missing(TempRes) || TempRes == "" || !TempRes %in% c("day", "month")) {
+    stop("Error: Please specify TempRes: 'day' or 'month'.")
+  }
+
+  climate_vars <- list(...)
+  HS_arg_count <- length(formals(HS))
+  if (length(climate_vars) != HS_arg_count) {
+    stop("The number of climate variables does not match the HS function's requirements.")
+  }
+
+  # read in individual files of climate variables
+  for (i in seq_along(climate_vars)) {
+    assign(paste0("var", i, ".stack"), terra::rast(climate_vars[[i]]))
+  }
+
+  SECTOR_INDEX <- which(SECTOR_ALL == SECTOR)
+  sector_filter <- SECTOR_FLAG[[SECTOR_INDEX]]
+
+  # -----
+  # subset layers for year of interest
+  layer_dates <- terra::time(var1.stack)
+  YEAR <- format(layer_dates, "%Y")
+
+
+  assertthat::assert_that(
+    YEAR_INPUT %in% unique(YEAR),
+    msg = paste0("invalid YEAR_INPUT, make sure input files include data for ", YEAR_INPUT)
+  )
+
+  # subset layers falls into the year of interest
+  year_layer_index <- which(YEAR == YEAR_INPUT)
+
+  # check number of layers
+  if (TempRes == "day") {
+    assertthat::assert_that(
+      length(year_layer_index) %in% c(365, 366),
+      msg = paste0("For daily input, expect 365 or 366 layers, but got ", length(year_layer_index), " layers.")
+    )
+  } else if (TempRes == "month") {
+    assertthat::assert_that(
+      length(year_layer_index) == 12,
+      msg = paste0("For monthly input, expect 12 layers, but got ", length(year_layer_index), " layers.")
+    )
+  }
+
+  variable_list <- list()
+  for (i in 1:length(formals(HS))) {
+    # assign SECTOR grid filter raster
+    assign(paste0("var", i, ".SECTOR.y"), var.i.SECTOR.y <- get(paste0("var", i, ".stack"))[[year_layer_index]] * terra::rast(sector_filter))
+    rm(list = paste0("var", i, ".stack"))
+    layer_dates.y <- terra::time(get(paste0("var", i, ".SECTOR.y")))
+    assign(paste0("stack", i, ".SECTOR.y"), raster::brick(get(paste0("var", i, ".SECTOR.y"))))
+    rack <- get(paste0("stack", i, ".SECTOR.y"))
+    rm(list = paste0("stack", i, ".SECTOR.y"))
+    # names(rack) <- layer_dates.y
+    variable_list[[i]] <- rack
+  }
+  HS.stack.y <- do.call(overlay, c(variable_list, fun = HS))
+  HS.stack.y <- stack(HS.stack.y)
+  names(HS.stack.y) <- layer_dates.y # name it so that DAY2MON function can work as expected
+  # check number of layers
+  if (TempRes == "day") {
+    assert_that(
+      nlayers(HS.stack.y) %in% c(365, 366),
+      msg = paste0("Error: For daily input, expect 365 or 366 layers, but got ", nlayers(HS.stack.y), " layers.")
+    )
+  } else if (TempRes == "month") {
+    assert_that(
+      nlayers(HS.stack.y) == 12,
+      msg = paste0("Error: For monthly input, expect 12 layers, but got ", nlayers(HS.stack.y), " layers.")
+    )
+  }
+  return(HS.stack.y)
+}

R/PWC.R → R/cal_pwc.R

@@ -1,4 +1,10 @@
-#' PWC
+#' cal_pwc
+#'
+#' Calculate Physical Work Capacity (PWC)
+#'
+#' The `cal_pwc` function computes daily grid-level physical work capacity (PWC),
+#' which range from 0 to 1, based on heat stress levels, a chosen labor heat response function,
+#' and workload intensity.
 #'
 #' @param WBGT output from HeatStress function
 #' @param LHR choice of labor heat response function
@@ -12,10 +18,9 @@
 #'
 #' @examples
 #' \dontrun{
-#' GD_PWC <- PWC(HeatStress = HS.stack.y,  LHR = LHR_Foster, workload = "high")
+#' GD_PWC <- cal_pwc(HeatStress = HS.stack.y, LHR = LHR_Foster, workload = "high")
 #' }
-
-PWC <- function(WBGT, LHR, workload){
+cal_pwc <- function(WBGT, LHR, workload) {
   PWC.stack.y <- stack(lapply(1:nlayers(WBGT), function(nlay) {
     calc(WBGT[[nlay]], function(x) LHR(x, workload))
   }))

R/daily_to_monthly.R

@@ -0,0 +1,29 @@
+#' daily_to_monthly
+#'
+#' Convert Daily Raster Data to Monthly Averages
+#'
+#' This function takes a raster stack of daily values and computes the monthly averages for each pixel.
+#'
+#' @param input_rack Raster rack of daily value
+#'
+#' @importFrom raster stack overlay calc stackApply as.matrix ncell nlayers
+#' @importFrom dplyr %>%
+#'
+#' @return A raster rack of monthly value
+#' @export
+#'
+#' @examples
+#' \dontrun{
+#' month_PWC <- daily_to_monthly(input_rack = PWC.stack.y)
+#' month_HS <- daily_to_monthly(input_rack = HS.stack.y)
+#' }
+daily_to_monthly <- function(input_rack) {
+  # Create a vector of months corresponding to each layer
+  layer_names <- names(input_rack)
+  layer_dates <- as.Date(gsub("X", "", layer_names), format = "%Y.%m.%d")
+  MONTHS <- format(layer_dates, "%m")
+  # Calculate monthly mean using stackApply
+  monthly_mean <- stackApply(input_rack, indices = MONTHS, fun = mean)
+  monthly_mean <- stack(monthly_mean)
+  return(monthly_mean)
+}