tests/testthat/test-temperature.R

In climatekit: Unified Climate Indices for Temperature, Precipitation, and Drought

test_that("ck_frost_days counts days below 0", {
  dates <- as.Date("2024-01-01") + 0:9
  tmin <- c(-2, 3, -1, 5, -3, 0, 2, -4, 1, -1)
  result <- ck_frost_days(tmin, dates)
  expect_s3_class(result, "data.frame")
  expect_equal(result$value, 5)
  expect_equal(result$index, "frost_days")
  expect_equal(result$unit, "days")
})

test_that("ck_frost_days monthly aggregation works", {
  dates <- c(as.Date("2024-01-01") + 0:4, as.Date("2024-02-01") + 0:4)
  tmin <- c(-2, 3, -1, 5, -3, 0, 2, -4, 1, -1)
  result <- ck_frost_days(tmin, dates, period = "monthly")
  expect_equal(nrow(result), 2)
  expect_equal(result$value, c(3, 2))
})

test_that("ck_frost_days rejects non-numeric", {
  dates <- as.Date("2024-01-01") + 0:2
  expect_error(ck_frost_days("a", dates), "numeric")
})

test_that("ck_frost_days rejects mismatched lengths", {
  expect_error(ck_frost_days(c(1, 2), as.Date("2024-01-01")), "same length")
})

test_that("ck_ice_days counts days where tmax < 0", {
  dates <- as.Date("2024-01-01") + 0:4
  tmax <- c(-2, 3, -1, 5, -3)
  result <- ck_ice_days(tmax, dates)
  expect_equal(result$value, 3)
})

test_that("ck_summer_days counts days where tmax > 25", {
  dates <- as.Date("2024-07-01") + 0:4
  tmax <- c(22, 26, 28, 24, 30)
  result <- ck_summer_days(tmax, dates)
  expect_equal(result$value, 3)
})

test_that("ck_tropical_nights counts days where tmin > 20", {
  dates <- as.Date("2024-07-01") + 0:4
  tmin <- c(18, 21, 22, 19, 25)
  result <- ck_tropical_nights(tmin, dates)
  expect_equal(result$value, 3)
})

test_that("ck_growing_season returns days using 6-day spells and tavg", {
  # Create a year of data with clear warm season
  dates <- as.Date("2024-01-01") + 0:364
  tavg <- sin(seq(0, 2 * pi, length.out = 365)) * 15 + 5
  result <- ck_growing_season(tavg, dates)
  expect_s3_class(result, "data.frame")
  expect_true(result$value > 0)
  expect_equal(result$index, "growing_season")
})

test_that("ck_growing_season returns 0 when no 6-day spell", {
  dates <- as.Date("2024-01-01") + 0:29
  tavg <- rep(-10, 30)
  result <- ck_growing_season(tavg, dates)
  expect_equal(result$value, 0)
})

test_that("ck_growing_season uses 6-day spell not 5-day", {
  # Exactly 5 warm days should NOT trigger GSL
  dates <- as.Date("2024-01-01") + 0:364
  tavg <- rep(0, 365)
  tavg[100:104] <- 10  # only 5 warm days
  result <- ck_growing_season(tavg, dates)
  expect_equal(result$value, 0)

  # 6 warm days should trigger
  tavg[100:105] <- 10
  result <- ck_growing_season(tavg, dates)
  expect_true(result$value > 0)
})

test_that("ck_growing_season ends at first cold spell after July 1", {
  dates <- as.Date("2024-01-01") + 0:364
  # Warm from day 60 to day 250, then cold
  tavg <- rep(0, 365)
  tavg[60:250] <- 10
  result <- ck_growing_season(tavg, dates, lat = 50)
  # Season starts at day 60, ends when cold spell starts after July 1
  # Cold spell of 6+ days starts at day 251
  # End of season = day 250
  expect_equal(result$value, 250 - 60 + 1)
})

test_that("ck_heating_degree_days sums correctly", {
  dates <- as.Date("2024-01-01") + 0:4
  tavg <- c(10, 15, 20, 5, 18)
  result <- ck_heating_degree_days(tavg, dates, base = 18)
  # (18-10) + (18-15) + 0 + (18-5) + 0 = 8 + 3 + 13 = 24
  expect_equal(result$value, 24)
})

test_that("ck_cooling_degree_days sums correctly", {
  dates <- as.Date("2024-07-01") + 0:4
  tavg <- c(20, 25, 15, 30, 18)
  result <- ck_cooling_degree_days(tavg, dates, base = 18)
  # (20-18) + (25-18) + 0 + (30-18) + 0 = 2 + 7 + 12 = 21
  expect_equal(result$value, 21)
})

test_that("ck_growing_degree_days sums correctly", {
  dates <- as.Date("2024-07-01") + 0:4
  tavg <- c(15, 20, 8, 12, 25)
  result <- ck_growing_degree_days(tavg, dates, base = 10)
  # (15-10) + (20-10) + 0 + (12-10) + (25-10) = 5 + 10 + 2 + 15 = 32
  expect_equal(result$value, 32)
})

test_that("ck_diurnal_range computes mean range", {
  dates <- as.Date("2024-01-01") + 0:4
  tmin <- c(0, 5, 2, 3, 1)
  tmax <- c(10, 15, 12, 13, 11)
  result <- ck_diurnal_range(tmin, tmax, dates)
  expect_equal(result$value, 10)
  expect_equal(result$unit, "\u00b0C")
})

test_that("ck_diurnal_range rejects mismatched lengths", {
  dates <- as.Date("2024-01-01") + 0:4
  expect_error(ck_diurnal_range(1:5, 1:3, dates), "same length")
})

test_that("ck_warm_spell detects spells >= 6 days", {
  dates <- as.Date("2024-01-01") + 0:19
  # 90th percentile will be high; create a clear 7-day spell above it
  tmax <- c(rep(10, 13), rep(50, 7))
  result <- ck_warm_spell(tmax, dates, threshold = 0.5)
  expect_true(result$value >= 7)
})

test_that("ck_warm_spell returns 0 with no long spells", {
  dates <- as.Date("2024-01-01") + 0:9
  tmax <- c(30, 10, 30, 10, 30, 10, 30, 10, 30, 10)
  result <- ck_warm_spell(tmax, dates)
  expect_equal(result$value, 0)
})

test_that("all temperature functions handle NA values", {
  dates <- as.Date("2024-01-01") + 0:4
  tmin <- c(-2, NA, -1, 5, -3)
  result <- ck_frost_days(tmin, dates)
  expect_s3_class(result, "data.frame")
})

test_that("temperature functions reject non-Date dates", {
  expect_error(ck_frost_days(c(1, 2), c("a", "b")), "Date")
})

# --- Reference value tests ---

test_that("frost days at exactly 0C returns 0 (strict < 0)", {
  dates <- as.Date("2024-01-01")
  result <- ck_frost_days(c(0), dates)
  expect_equal(result$value, 0)
})

test_that("HDD with base=0 returns 0 when all temps >= 0", {
  dates <- as.Date("2024-01-01") + 0:4
  tavg <- c(0, 5, 10, 15, 20)
  result <- ck_heating_degree_days(tavg, dates, base = 0)
  expect_equal(result$value, 0)
})

test_that("CDD with base=0 returns sum of all temps", {
  dates <- as.Date("2024-01-01") + 0:4
  tavg <- c(5, 10, 15, 20, 25)
  result <- ck_cooling_degree_days(tavg, dates, base = 0)
  expect_equal(result$value, 75)
})

# ETCCDI extreme-value indices (TXx / TNx / TXn / TNn) ------------------------

test_that("ck_txx returns max of daily Tmax", {
  dates <- as.Date("2024-01-01") + 0:9
  tmax <- c(5, 10, 18, 12, 4, 8, 22, 3, 7, 6)
  result <- ck_txx(tmax, dates)
  expect_s3_class(result, "data.frame")
  expect_equal(result$value, 22)
  expect_equal(result$index, "txx")
  expect_equal(result$unit, "\u00b0C")
})

test_that("ck_txx monthly aggregation works", {
  dates <- c(as.Date("2024-01-01") + 0:4, as.Date("2024-02-01") + 0:4)
  tmax <- c(5, 10, 18, 12, 4, 8, 22, 3, 7, 6)
  result <- ck_txx(tmax, dates, period = "monthly")
  expect_equal(nrow(result), 2)
  expect_equal(result$value, c(18, 22))
})

test_that("ck_tnx returns max of daily Tmin (warmest night)", {
  dates <- as.Date("2024-07-01") + 0:9
  tmin <- c(15, 18, 22, 19, 14, 21, 23, 17, 20, 19)
  result <- ck_tnx(tmin, dates)
  expect_equal(result$value, 23)
  expect_equal(result$index, "tnx")
})

test_that("ck_txn returns min of daily Tmax (coldest day)", {
  dates <- as.Date("2024-01-01") + 0:9
  tmax <- c(5, 10, -3, 12, 4, 8, 22, -8, 7, 6)
  result <- ck_txn(tmax, dates)
  expect_equal(result$value, -8)
  expect_equal(result$index, "txn")
})

test_that("ck_tnn returns min of daily Tmin (coldest night)", {
  dates <- as.Date("2024-01-01") + 0:9
  tmin <- c(-2, 3, -1, 5, -8, 0, 2, -12, 1, -1)
  result <- ck_tnn(tmin, dates)
  expect_equal(result$value, -12)
  expect_equal(result$index, "tnn")
})

test_that("extreme-value functions reject non-numeric", {
  dates <- as.Date("2024-01-01") + 0:2
  expect_error(ck_txx("a", dates), "numeric")
  expect_error(ck_tnx("a", dates), "numeric")
  expect_error(ck_txn("a", dates), "numeric")
  expect_error(ck_tnn("a", dates), "numeric")
})

test_that("extreme-value functions reject mismatched lengths", {
  expect_error(ck_txx(c(1, 2), as.Date("2024-01-01")), "same length")
  expect_error(ck_tnn(c(1, 2), as.Date("2024-01-01")), "same length")
})

test_that("ETCCDI: txx >= tnx >= tnn and txx >= txn", {
  set.seed(42)
  dates <- as.Date("2024-01-01") + 0:364
  tmin <- rnorm(365, mean = 5, sd = 8)
  tmax <- tmin + abs(rnorm(365, mean = 8, sd = 3))
  expect_gte(ck_txx(tmax, dates)$value, ck_tnx(tmin, dates)$value)
  expect_gte(ck_tnx(tmin, dates)$value, ck_tnn(tmin, dates)$value)
  expect_gte(ck_txx(tmax, dates)$value, ck_txn(tmax, dates)$value)
})

# ETCCDI percentile-based indices (TX10p / TN10p / TX90p / TN90p) ------------

test_that("ck_tx10p returns approx 10% on uniform synthetic ref data", {
  set.seed(1)
  dates <- seq(as.Date("1961-01-01"), as.Date("1962-12-31"), by = "day")
  tmax <- 15 + 10 * sin(2 * pi * as.integer(format(dates, "%j")) / 365) +
          rnorm(length(dates))
  result <- ck_tx10p(tmax, dates, ref_start = 1961L, ref_end = 1962L)
  expect_s3_class(result, "data.frame")
  expect_equal(result$index[1], "tx10p")
  expect_equal(result$unit[1], "%")
  # Within the base period itself the long-run average should be near 10%.
  expect_true(all(result$value >= 0 & result$value <= 100))
  expect_true(abs(mean(result$value) - 10) < 5)
})

test_that("ck_tx90p returns approx 10% on uniform synthetic ref data", {
  set.seed(2)
  dates <- seq(as.Date("1961-01-01"), as.Date("1962-12-31"), by = "day")
  tmax <- 15 + 10 * sin(2 * pi * as.integer(format(dates, "%j")) / 365) +
          rnorm(length(dates))
  result <- ck_tx90p(tmax, dates, ref_start = 1961L, ref_end = 1962L)
  expect_true(abs(mean(result$value) - 10) < 5)
  expect_equal(result$index[1], "tx90p")
})

test_that("ck_tn10p and ck_tn90p run on Tmin and dispatch on percentile", {
  set.seed(3)
  dates <- seq(as.Date("1961-01-01"), as.Date("1962-12-31"), by = "day")
  tmin <- 5 + 8 * sin(2 * pi * as.integer(format(dates, "%j")) / 365) +
          rnorm(length(dates))
  r10 <- ck_tn10p(tmin, dates, ref_start = 1961L, ref_end = 1962L)
  r90 <- ck_tn90p(tmin, dates, ref_start = 1961L, ref_end = 1962L)
  expect_equal(r10$index[1], "tn10p")
  expect_equal(r90$index[1], "tn90p")
  expect_true(abs(mean(r10$value) - 10) < 5)
  expect_true(abs(mean(r90$value) - 10) < 5)
})

test_that("percentile functions error if no data in reference period", {
  dates <- seq(as.Date("2000-01-01"), as.Date("2000-12-31"), by = "day")
  tmax <- rnorm(length(dates), 15, 5)
  expect_error(
    ck_tx10p(tmax, dates, ref_start = 1961L, ref_end = 1990L),
    "reference period"
  )
})

test_that("percentile functions reject mismatched lengths", {
  expect_error(ck_tx10p(c(1, 2), as.Date("1961-01-01")), "same length")
  expect_error(ck_tn90p(c(1, 2), as.Date("1961-01-01")), "same length")
})

# Zhang 2005 in-base bootstrap (Phase D) -------------------------------------

test_that("ck_tx10p bootstrap=FALSE matches the unmodified threshold path", {
  set.seed(81)
  dates <- seq(as.Date("1961-01-01"), as.Date("1965-12-31"), by = "day")
  tmax <- 15 + 10 * sin(2 * pi * as.integer(format(dates, "%j")) / 365) +
          rnorm(length(dates))
  result <- ck_tx10p(tmax, dates, ref_start = 1961L, ref_end = 1965L,
                     bootstrap = FALSE)
  expect_s3_class(result, "data.frame")
  expect_true(all(result$value >= 0 & result$value <= 100))
})

test_that("bootstrap correction nudges in-base years toward expected 10%", {
  # Ten-year iid sample: in-sample bias suppresses exceedance under the
  # standard threshold; bootstrap should move the mean closer to 10%.
  set.seed(82)
  dates <- seq(as.Date("1961-01-01"), as.Date("1970-12-31"), by = "day")
  tmax <- 15 + 10 * sin(2 * pi * as.integer(format(dates, "%j")) / 365) +
          rnorm(length(dates))
  no_bs <- ck_tx10p(tmax, dates, ref_start = 1961L, ref_end = 1970L)
  bs    <- ck_tx10p(tmax, dates, ref_start = 1961L, ref_end = 1970L,
                    bootstrap = TRUE)
  expect_lt(abs(mean(bs$value) - 10), abs(mean(no_bs$value) - 10) + 1e-6)
})

test_that("bootstrap exceedance values stay within the unit interval after %", {
  set.seed(83)
  dates <- seq(as.Date("1961-01-01"), as.Date("1963-12-31"), by = "day")
  tmin <- 5 + 8 * sin(2 * pi * as.integer(format(dates, "%j")) / 365) +
          rnorm(length(dates))
  bs <- ck_tn90p(tmin, dates, ref_start = 1961L, ref_end = 1963L,
                 bootstrap = TRUE)
  expect_true(all(bs$value >= 0 & bs$value <= 100))
})

test_that("bootstrap argument validation", {
  set.seed(84)
  dates <- seq(as.Date("1961-01-01"), as.Date("1962-12-31"), by = "day")
  tmax <- rnorm(length(dates), 15, 5)
  expect_error(
    ck_tx10p(tmax, dates, ref_start = 1961L, ref_end = 1962L,
             bootstrap = "yes"),
    "logical"
  )
  expect_error(
    ck_tx90p(tmax, dates, ref_start = 1961L, ref_end = 1962L,
             bootstrap = c(TRUE, FALSE)),
    "logical"
  )
})

# Calendar-day-base spell duration indices (CSDI / WSDI) ---------------------

test_that("ck_wsdi returns 0 when no 6-day exceedance spell", {
  set.seed(21)
  dates <- seq(as.Date("1961-01-01"), as.Date("1962-12-31"), by = "day")
  tmax <- 15 + 10 * sin(2 * pi * as.integer(format(dates, "%j")) / 365) +
          rnorm(length(dates))
  result <- ck_wsdi(tmax, dates, ref_start = 1961L, ref_end = 1962L)
  expect_s3_class(result, "data.frame")
  expect_equal(result$index[1], "wsdi")
  expect_equal(result$unit[1], "days")
  expect_true(all(result$value >= 0))
})

test_that("ck_wsdi captures injected long warm spell", {
  set.seed(22)
  dates <- seq(as.Date("1961-01-01"), as.Date("1962-12-31"), by = "day")
  tmax <- 15 + 10 * sin(2 * pi * as.integer(format(dates, "%j")) / 365) +
          rnorm(length(dates))
  # Inject 10-day hot spell in 1962 (outside ref period 1961 only).
  hot <- which(dates >= as.Date("1962-07-15") & dates <= as.Date("1962-07-24"))
  tmax[hot] <- tmax[hot] + 20
  result <- ck_wsdi(tmax, dates, ref_start = 1961L, ref_end = 1961L)
  yr_1962 <- result$value[format(result$period, "%Y") == "1962"]
  expect_gte(yr_1962, 10)
})

test_that("ck_csdi returns 0 when no 6-day under-threshold spell", {
  set.seed(23)
  dates <- seq(as.Date("1961-01-01"), as.Date("1962-12-31"), by = "day")
  tmin <- 5 + 8 * sin(2 * pi * as.integer(format(dates, "%j")) / 365) +
          rnorm(length(dates))
  result <- ck_csdi(tmin, dates, ref_start = 1961L, ref_end = 1962L)
  expect_equal(result$index[1], "csdi")
  expect_equal(result$unit[1], "days")
  expect_true(all(result$value >= 0))
})

test_that("ck_csdi captures injected long cold spell", {
  set.seed(24)
  dates <- seq(as.Date("1961-01-01"), as.Date("1962-12-31"), by = "day")
  tmin <- 5 + 8 * sin(2 * pi * as.integer(format(dates, "%j")) / 365) +
          rnorm(length(dates))
  # Inject 10-day cold spell in 1962 (outside ref period 1961 only).
  cold <- which(dates >= as.Date("1962-02-01") & dates <= as.Date("1962-02-10"))
  tmin[cold] <- tmin[cold] - 20
  result <- ck_csdi(tmin, dates, ref_start = 1961L, ref_end = 1961L)
  yr_1962 <- result$value[format(result$period, "%Y") == "1962"]
  expect_gte(yr_1962, 10)
})

test_that("ck_wsdi / ck_csdi reject mismatched lengths and missing reference", {
  expect_error(ck_wsdi(c(1, 2), as.Date("1961-01-01")), "same length")
  expect_error(ck_csdi(c(1, 2), as.Date("1961-01-01")), "same length")
  dates <- seq(as.Date("2000-01-01"), as.Date("2000-12-31"), by = "day")
  expect_error(
    ck_wsdi(rnorm(length(dates), 15, 5), dates, ref_start = 1961L, ref_end = 1990L),
    "reference period"
  )
})