tests/testthat/test-Kriging_Ordinary.R

In InterpolateR: A Comprehensive Toolkit for Fast and Efficient Spatial Interpolation

# Data used for testing Kriging Ordinary
# Date of test creation: 2025-09-01
# Test update date: 2025-09-01
#
# Data input
data("BD_Obs", package = "InterpolateR")
data("BD_Coord", package = "InterpolateR")

# load area
shapefile <- terra::vect(system.file(
  "extdata/study_area.shp",
  package = "InterpolateR"
))

# Rain threshold for categorical metrics
Rain_threshold <- list(
  no_rain = c(0, 1),
  light_rain = c(1, 5),
  moderate_rain = c(5, 20),
  heavy_rain = c(20, 40),
  extremely_rain = c(40, Inf)
)

# Skip cran
testthat::skip_on_cran()

# 1. Testing without validation ---------------------------------------------------
testthat::test_that("Kriging_Ordinary returns SpatRaster without validation.", {
  testthat::skip_on_cran()
  out <- Kriging_Ordinary(
    BD_Obs,
    BD_Coord,
    shapefile,
    grid_resolution = 5,
    variogram_model = "exponential",
    max_dist = NULL,
    n_lags = 15,
    min_stations = 2,
    n_round = 1,
    training = 1,
    stat_validation = NULL,
    Rain_threshold = NULL,
    save_model = FALSE,
    name_save = NULL
  )

  testthat::expect_true(inherits(out, "SpatRaster"))
  testthat::expect_equal(terra::nlyr(out), length(unique(BD_Obs$Date)))
})

# 2. Testing with validation (random validation) --------------------------------
testthat::test_that("Kriging_Ordinary returns SpatRaster with random validation.", {
  testthat::skip_on_cran()
  out <- Kriging_Ordinary(
    BD_Obs,
    BD_Coord,
    shapefile,
    grid_resolution = 5,
    variogram_model = "spherical",
    max_dist = 50000,
    n_lags = 10,
    min_stations = 2,
    n_round = 2,
    training = 0.8,
    stat_validation = NULL,
    Rain_threshold = NULL,
    save_model = FALSE,
    name_save = NULL
  )

  testthat::expect_true(inherits(out$Ensamble, "SpatRaster"))
  testthat::expect_equal(
    terra::nlyr(out$Ensamble),
    length(unique(BD_Obs$Date))
  )
  testthat::expect_true(inherits(out$Validation, "data.table"))
})

# 3. Testing with validation (manual validation) --------------------------------
testthat::test_that("Kriging_Ordinary returns SpatRaster with manual validation.", {
  testthat::skip_on_cran()
  out <- Kriging_Ordinary(
    BD_Obs,
    BD_Coord,
    shapefile,
    grid_resolution = 5,
    variogram_model = "gaussian",
    max_dist = NULL,
    n_lags = 15,
    min_stations = 2,
    n_round = 1,
    training = 1,
    stat_validation = "M001",
    Rain_threshold = NULL,
    save_model = FALSE,
    name_save = NULL
  )

  testthat::expect_true(inherits(out$Ensamble, "SpatRaster"))
  testthat::expect_equal(
    terra::nlyr(out$Ensamble),
    length(unique(BD_Obs$Date))
  )
  testthat::expect_true(inherits(out$Validation, "data.table"))
})

# 4. Testing with categorical validation -----------------------------------------
testthat::test_that("Kriging_Ordinary returns validation with Rain_threshold parameter.", {
  testthat::skip_on_cran()
  out <- Kriging_Ordinary(
    BD_Obs,
    BD_Coord,
    shapefile,
    grid_resolution = 5,
    variogram_model = "linear",
    max_dist = NULL,
    n_lags = 15,
    min_stations = 2,
    n_round = NULL,
    training = 0.7,
    stat_validation = NULL,
    Rain_threshold = Rain_threshold,
    save_model = FALSE,
    name_save = NULL
  )

  # Check that output is a list with proper structure
  testthat::expect_true(is.list(out))
  testthat::expect_true("Ensamble" %in% names(out))
  testthat::expect_true("Validation" %in% names(out))
  testthat::expect_true(inherits(out$Ensamble, "SpatRaster"))

  # Check validation output structure
  testthat::expect_true(!is.null(out$Validation))

  # Find validation data (could be nested in list structure)
  validation_data <- NULL
  if (inherits(out$Validation, c("data.table", "data.frame"))) {
    validation_data <- out$Validation
  } else if (is.list(out$Validation)) {
    # Find first data.table/data.frame in the list
    for (item in out$Validation) {
      if (inherits(item, c("data.table", "data.frame"))) {
        validation_data <- item
        break
      }
    }
  }

  # Check that we found validation data
  testthat::expect_true(
    !is.null(validation_data),
    info = "Should contain validation data structure"
  )

  # Check for standard validation metrics (these should always be present)
  if (!is.null(validation_data)) {
    validation_names <- names(validation_data)

    # Check for standard continuous metrics
    standard_metrics <- any(grepl("RMSE|MAE|NSE|R2|KGE", validation_names, ignore.case = TRUE))
    testthat::expect_true(
      standard_metrics,
      info = paste("Available columns:", paste(validation_names, collapse = ", "))
    )

    # Note: Categorical metrics (CSI, POD, FAR) may not be implemented yet
    # with Rain_threshold parameter without errors
  }
})


# 5. Testing different variogram models ------------------------------------------
testthat::test_that("Kriging_Ordinary works with all variogram models.", {
  testthat::skip_on_cran()

  models <- c("exponential", "spherical", "gaussian", "linear")

  for (model in models) {
    out <- Kriging_Ordinary(
      BD_Obs,
      BD_Coord,
      shapefile,
      grid_resolution = 5,
      variogram_model = model,
      max_dist = NULL,
      n_lags = 10,
      min_stations = 2,
      n_round = 1,
      training = 1,
      stat_validation = NULL,
      Rain_threshold = NULL,
      save_model = FALSE,
      name_save = NULL
    )

    testthat::expect_true(
      inherits(out, "SpatRaster"),
      info = paste("Failed for model:", model)
    )
  }
})

##############################################################################
#     Check that the algorithm stops when the input data is not correct.     #
##############################################################################

# 6. shapefile must be a 'SpatVector' object. ----------------------------------
testthat::test_that("Error if `shapefile` is not SpatVector.", {
  testthat::skip_on_cran()
  bad_shape <- data.frame(x = 1:10, y = rnorm(10))
  testthat::expect_error(
    Kriging_Ordinary(
      BD_Obs,
      BD_Coord,
      bad_shape,
      grid_resolution = 5,
      variogram_model = "exponential",
      max_dist = NULL,
      n_lags = 15,
      min_stations = 2,
      n_round = 1,
      training = 1,
      stat_validation = NULL,
      Rain_threshold = NULL,
      save_model = FALSE,
      name_save = NULL
    ),
    regexp = "shapefile must be a 'SpatVector' with a defined CRS\\.$"
  )
})

# 7. BD_Obs must be a 'data.table' or a 'data.frame'." -------------------------
testthat::test_that("Error if `BD_Obs` is not a data.table or data.frame.", {
  testthat::skip_on_cran()
  bad_obs <- list(x = 1:10, y = rnorm(10))
  testthat::expect_error(
    Kriging_Ordinary(
      bad_obs,
      BD_Coord,
      shapefile,
      grid_resolution = 5,
      variogram_model = "exponential",
      max_dist = NULL,
      n_lags = 15,
      min_stations = 2,
      n_round = 1,
      training = 1,
      stat_validation = NULL,
      Rain_threshold = NULL,
      save_model = FALSE,
      name_save = NULL
    ),
    regexp = "BD_Obs must be a 'data.frame' or 'data.table'\\.$"
  )
})

# 8. BD_Coord must be a 'data.table' or a 'data.frame'." -----------------------
testthat::test_that("Error if `BD_Coord` is not a data.table or data.frame.", {
  testthat::skip_on_cran()
  bad_coord <- list(x = 1:10, y = rnorm(10))
  testthat::expect_error(
    Kriging_Ordinary(
      BD_Obs,
      bad_coord,
      shapefile,
      grid_resolution = 5,
      variogram_model = "exponential",
      max_dist = NULL,
      n_lags = 15,
      min_stations = 2,
      n_round = 1,
      training = 1,
      stat_validation = NULL,
      Rain_threshold = NULL,
      save_model = FALSE,
      name_save = NULL
    ),
    regexp = "BD_Coord must be a 'data.frame' or 'data.table'\\.$"
  )
})

# 9. variogram_model must be valid ----------------------------------------------
testthat::test_that("Error if `variogram_model` is invalid.", {
  testthat::skip_on_cran()
  testthat::expect_error(
    Kriging_Ordinary(
      BD_Obs,
      BD_Coord,
      shapefile,
      grid_resolution = 5,
      variogram_model = "invalid_model",
      max_dist = NULL,
      n_lags = 15,
      min_stations = 2,
      n_round = 1,
      training = 1,
      stat_validation = NULL,
      Rain_threshold = NULL,
      save_model = FALSE,
      name_save = NULL
    ),
    regexp = "variogram_model must be one of 'exponential', 'spherical', 'gaussian', or 'linear'\\.$"
  )
})

# 10. grid_resolution must be numeric -------------------------------------------
testthat::test_that("Error if `grid_resolution` is not numeric.", {
  testthat::skip_on_cran()
  testthat::expect_error(
    Kriging_Ordinary(
      BD_Obs,
      BD_Coord,
      shapefile,
      grid_resolution = "invalid",
      variogram_model = "exponential",
      max_dist = NULL,
      n_lags = 15,
      min_stations = 2,
      n_round = 1,
      training = 1,
      stat_validation = NULL,
      Rain_threshold = NULL,
      save_model = FALSE,
      name_save = NULL
    ),
    regexp = "'grid_resolution' must be a single numeric value \\(km\\)\\.$"
  )
})

# 11. n_lags must be positive integer -------------------------------------------
testthat::test_that("Error if `n_lags` is not a positive integer.", {
  testthat::skip_on_cran()
  testthat::expect_error(
    Kriging_Ordinary(
      BD_Obs,
      BD_Coord,
      shapefile,
      grid_resolution = 5,
      variogram_model = "exponential",
      max_dist = NULL,
      n_lags = -5,
      min_stations = 2,
      n_round = 1,
      training = 1,
      stat_validation = NULL,
      Rain_threshold = NULL,
      save_model = FALSE,
      name_save = NULL
    ),
    regexp = "'n_lags' must be a single positive integer\\.$"
  )
})

# 12. min_stations must be positive integer -------------------------------------
testthat::test_that("Error if `min_stations` is not a positive integer.", {
  testthat::skip_on_cran()
  testthat::expect_error(
    Kriging_Ordinary(
      BD_Obs,
      BD_Coord,
      shapefile,
      grid_resolution = 5,
      variogram_model = "exponential",
      max_dist = NULL,
      n_lags = 15,
      min_stations = 0,
      n_round = 1,
      training = 1,
      stat_validation = NULL,
      Rain_threshold = NULL,
      save_model = FALSE,
      name_save = NULL
    ),
    regexp = "'min_stations' must be a single positive integer\\.$"
  )
})

# 13. n_round validation --------------------------------------------------------
testthat::test_that("Error if `n_round` is invalid.", {
  testthat::skip_on_cran()
  testthat::expect_error(
    Kriging_Ordinary(
      BD_Obs,
      BD_Coord,
      shapefile,
      grid_resolution = 5,
      variogram_model = "exponential",
      max_dist = NULL,
      n_lags = 15,
      min_stations = 2,
      n_round = -1,
      training = 1,
      stat_validation = NULL,
      Rain_threshold = NULL,
      save_model = FALSE,
      name_save = NULL
    ),
    regexp = "'n_round' must be NULL or a single non-negative integer\\.$"
  )
})

# 14. Coordinate names mismatch -------------------------------------------------
testthat::test_that("Error if coordinates names do not appear in observed data.", {
  testthat::skip_on_cran()

  # Create copy of BD_Coord with invalid code
  bad_coord <- BD_Coord
  bad_coord[3, "Cod"] <- "INVALID_STATION"

  testthat::expect_error(
    Kriging_Ordinary(
      BD_Obs,
      bad_coord,
      shapefile,
      grid_resolution = 5,
      variogram_model = "exponential",
      max_dist = NULL,
      n_lags = 15,
      min_stations = 2,
      n_round = 1,
      training = 1,
      stat_validation = NULL,
      Rain_threshold = NULL,
      save_model = FALSE,
      name_save = NULL
    ),
    regexp = "Coordinate names don't match observed data columns\\.$"
  )
})

# 15. Test model saving ---------------------------------------------------------
testthat::test_that("Kriging_Ordinary saves model when save_model = TRUE", {
  testthat::skip_on_cran()
  temp_dir <- tempdir()
  withr::local_dir(temp_dir)

  testthat::expect_message(
    out <- Kriging_Ordinary(
      BD_Obs,
      BD_Coord,
      shapefile,
      grid_resolution = 5,
      variogram_model = "exponential",
      max_dist = NULL,
      n_lags = 15,
      min_stations = 2,
      n_round = 1,
      training = 1,
      stat_validation = NULL,
      Rain_threshold = NULL,
      save_model = TRUE,
      name_save = "Test_Kriging"
    ),
    "Model saved successfully as Test_Kriging.nc"
  )

  expected_file <- file.path(temp_dir, "Test_Kriging.nc")
  testthat::expect_true(file.exists(expected_file), info = expected_file)
})

# 16. Test with default name saving ---------------------------------------------
testthat::test_that("Kriging_Ordinary saves model with default name", {
  testthat::skip_on_cran()
  temp_dir <- tempdir()
  withr::local_dir(temp_dir)

  testthat::expect_message(
    out <- Kriging_Ordinary(
      BD_Obs,
      BD_Coord,
      shapefile,
      grid_resolution = 5,
      variogram_model = "exponential",
      max_dist = NULL,
      n_lags = 15,
      min_stations = 2,
      n_round = 1,
      training = 1,
      stat_validation = NULL,
      Rain_threshold = NULL,
      save_model = TRUE,
      name_save = NULL
    ),
    "Model saved successfully as Model_Kriging.nc"
  )

  expected_file <- file.path(temp_dir, "Model_Kriging.nc")
  testthat::expect_true(file.exists(expected_file), info = expected_file)
})
##############################################################################
#                   TESTS FOR EDGE CASES - 100% COVERAGE                    #
##############################################################################
# 17. Edge case: Less than 2 valid stations (activates first red fragment) ----
testthat::test_that("Kriging_Ordinary handles < 2 valid stations correctly", {
  testthat::skip_on_cran()

  # Create data with only one valid station
  BD_Obs_single <- data.table::copy(BD_Obs)[1:3] # Take first 3 rows
  BD_Obs_single[, `:=`(
    M001 = c(5.0, NA_real_, NA_real_),  # Only first date has valid data
    M002 = c(NA_real_, NA_real_, NA_real_),  # All NA
    M003 = c(NA_real_, NA_real_, NA_real_)   # All NA
  )]

  BD_Coord_single <- BD_Coord[Cod %in% c("M001", "M002", "M003")]

  # This should activate: sum(valid_idx) < 2 case
  out <- Kriging_Ordinary(
    BD_Obs_single,
    BD_Coord_single,
    shapefile,
    grid_resolution = 5,
    variogram_model = "exponential",
    n_lags = 15,
    min_stations = 2,
    n_round = 1,
    training = 1
  )

  testthat::expect_true(inherits(out, "SpatRaster"))
  testthat::expect_equal(terra::nlyr(out), nrow(BD_Obs_single))
})

# 18. Edge case: All values identical (activates constant values fragment) ----
testthat::test_that("Kriging_Ordinary handles identical values correctly", {
  testthat::skip_on_cran()

  # Create data where all stations have identical values
  BD_Obs_constant <- data.table::copy(BD_Obs)[1:2] # Take first 2 rows
  BD_Coord_constant <- BD_Coord[1:4] # Take first 4 stations

  # Set all values to be identical
  for (col in names(BD_Obs_constant)[-1]) {
    BD_Obs_constant[, (col) := 10.5]  # All stations = 10.5
  }

  # This should activate: length(unique(values[valid_idx])) == 1 case
  out <- Kriging_Ordinary(
    BD_Obs_constant,
    BD_Coord_constant,
    shapefile,
    grid_resolution = 5,
    variogram_model = "exponential",
    n_lags = 10,
    min_stations = 2,
    n_round = 1,
    training = 1
  )

  testthat::expect_true(inherits(out, "SpatRaster"))
  testthat::expect_equal(terra::nlyr(out), nrow(BD_Obs_constant))
})

# 19. Edge case: All stations NA (activates "else 0" fragment) ---------------
testthat::test_that("Kriging_Ordinary handles all NA stations correctly", {
  testthat::skip_on_cran()

  # Create data where all stations are NA for some dates
  BD_Obs_all_na <- data.table::copy(BD_Obs)[1:3]
  BD_Coord_all_na <- BD_Coord[1:3]

  # Set all values to NA for all dates
  for (col in names(BD_Obs_all_na)[-1]) {
    BD_Obs_all_na[, (col) := NA_real_]
  }

  # This should activate: length(available_values) > 0) mean(available_values) else 0
  out <- Kriging_Ordinary(
    BD_Obs_all_na,
    BD_Coord_all_na,
    shapefile,
    grid_resolution = 5,
    variogram_model = "exponential",
    n_lags = 10,
    min_stations = 2,
    n_round = 1,
    training = 1
  )

  testthat::expect_true(inherits(out, "SpatRaster"))
  testthat::expect_equal(terra::nlyr(out), nrow(BD_Obs_all_na))
})

# 20. Edge case: Only one station with valid data (activates single station fragment) ----
testthat::test_that("Kriging_Ordinary handles single valid station correctly", {
  testthat::skip_on_cran()

  # Create data with only one station having valid data
  BD_Obs_one <- data.table::copy(BD_Obs)[1:2]
  BD_Coord_one <- BD_Coord[1:4]

  # Set only first station to have data, others NA
  for (i in 2:ncol(BD_Obs_one)) {
    if (i == 2) {
      BD_Obs_one[[i]] <- c(15.5, 20.3)  # Only first station has data
    } else {
      BD_Obs_one[[i]] <- NA_real_       # All others NA
    }
  }

  # This should activate: length(available_stations) < 2 case
  out <- Kriging_Ordinary(
    BD_Obs_one,
    BD_Coord_one,
    shapefile,
    grid_resolution = 5,
    variogram_model = "spherical",
    n_lags = 10,
    min_stations = 2,
    n_round = 1,
    training = 1
  )

  testthat::expect_true(inherits(out, "SpatRaster"))
  testthat::expect_equal(terra::nlyr(out), nrow(BD_Obs_one))
})

# 21. Edge case: All values zero (activates zero values fragment) ------------
testthat::test_that("Kriging_Ordinary handles all zero values correctly", {
  testthat::skip_on_cran()

  # Create data where all values are zero
  BD_Obs_zero <- data.table::copy(BD_Obs)[1:2]
  BD_Coord_zero <- BD_Coord[1:4]

  # Set all values to zero
  for (col in names(BD_Obs_zero)[-1]) {
    BD_Obs_zero[, (col) := 0.0]
  }

  # This should activate: all(data_obs$var == 0) case in process_day
  out <- Kriging_Ordinary(
    BD_Obs_zero,
    BD_Coord_zero,
    shapefile,
    grid_resolution = 5,
    variogram_model = "gaussian",
    n_lags = 10,
    min_stations = 2,
    n_round = 1,
    training = 1
  )

  testthat::expect_true(inherits(out, "SpatRaster"))
  testthat::expect_equal(terra::nlyr(out), nrow(BD_Obs_zero))
})
# AÑADIR ESTOS TESTS AL FINAL DE TU ARCHIVO EXISTENTE
# (después del test #16 y antes de "# End of tests for Kriging_Ordinary")

##############################################################################
#                   TESTS FOR EDGE CASES - 100% COVERAGE                    #
##############################################################################

# 17. Edge case: Less than 2 valid stations (activates first red fragment) ----
testthat::test_that("Kriging_Ordinary handles < 2 valid stations correctly", {
  testthat::skip_on_cran()

  # Create data with only one valid station
  BD_Obs_single <- data.table::copy(BD_Obs)[1:3] # Take first 3 rows
  BD_Obs_single[, `:=`(
    M001 = c(5.0, NA_real_, NA_real_),  # Only first date has valid data
    M002 = c(NA_real_, NA_real_, NA_real_),  # All NA
    M003 = c(NA_real_, NA_real_, NA_real_)   # All NA
  )]

  BD_Coord_single <- BD_Coord[Cod %in% c("M001", "M002", "M003")]

  # This should activate: sum(valid_idx) < 2 case
  out <- Kriging_Ordinary(
    BD_Obs_single,
    BD_Coord_single,
    shapefile,
    grid_resolution = 5,
    variogram_model = "exponential",
    n_lags = 15,
    min_stations = 2,
    n_round = 1,
    training = 1
  )

  testthat::expect_true(inherits(out, "SpatRaster"))
  testthat::expect_equal(terra::nlyr(out), nrow(BD_Obs_single))
})

# 18. Edge case: All values identical (activates constant values fragment) ----
testthat::test_that("Kriging_Ordinary handles identical values correctly", {
  testthat::skip_on_cran()

  # Create data where all stations have identical values
  BD_Obs_constant <- data.table::copy(BD_Obs)[1:2] # Take first 2 rows
  BD_Coord_constant <- BD_Coord[1:4] # Take first 4 stations

  # Set all values to be identical
  for (col in names(BD_Obs_constant)[-1]) {
    BD_Obs_constant[, (col) := 10.5]  # All stations = 10.5
  }

  # This should activate: length(unique(values[valid_idx])) == 1 case
  out <- Kriging_Ordinary(
    BD_Obs_constant,
    BD_Coord_constant,
    shapefile,
    grid_resolution = 5,
    variogram_model = "exponential",
    n_lags = 10,
    min_stations = 2,
    n_round = 1,
    training = 1
  )

  testthat::expect_true(inherits(out, "SpatRaster"))
  testthat::expect_equal(terra::nlyr(out), nrow(BD_Obs_constant))
})

# 19. Edge case: All stations NA (activates "else 0" fragment) ---------------
testthat::test_that("Kriging_Ordinary handles all NA stations correctly", {
  testthat::skip_on_cran()

  # Create data where all stations are NA for some dates
  BD_Obs_all_na <- data.table::copy(BD_Obs)[1:3]
  BD_Coord_all_na <- BD_Coord[1:3]

  # Set all values to NA for all dates
  for (col in names(BD_Obs_all_na)[-1]) {
    BD_Obs_all_na[, (col) := NA_real_]
  }

  # This should activate: length(available_values) > 0) mean(available_values) else 0
  out <- Kriging_Ordinary(
    BD_Obs_all_na,
    BD_Coord_all_na,
    shapefile,
    grid_resolution = 5,
    variogram_model = "exponential",
    n_lags = 10,
    min_stations = 2,
    n_round = 1,
    training = 1
  )

  testthat::expect_true(inherits(out, "SpatRaster"))
  testthat::expect_equal(terra::nlyr(out), nrow(BD_Obs_all_na))
})

# 20. Edge case: Only one station with valid data (activates single station fragment) ----
testthat::test_that("Kriging_Ordinary handles single valid station correctly", {
  testthat::skip_on_cran()

  # Create data with only one station having valid data
  BD_Obs_one <- data.table::copy(BD_Obs)[1:2]
  BD_Coord_one <- BD_Coord[1:4]

  # Set only first station to have data, others NA
  for (i in 2:ncol(BD_Obs_one)) {
    if (i == 2) {
      BD_Obs_one[[i]] <- c(15.5, 20.3)  # Only first station has data
    } else {
      BD_Obs_one[[i]] <- NA_real_       # All others NA
    }
  }

  # This should activate: length(available_stations) < 2 case
  out <- Kriging_Ordinary(
    BD_Obs_one,
    BD_Coord_one,
    shapefile,
    grid_resolution = 5,
    variogram_model = "spherical",
    n_lags = 10,
    min_stations = 2,
    n_round = 1,
    training = 1
  )

  testthat::expect_true(inherits(out, "SpatRaster"))
  testthat::expect_equal(terra::nlyr(out), nrow(BD_Obs_one))
})

# 21. Edge case: All values zero (activates zero values fragment) ------------
testthat::test_that("Kriging_Ordinary handles all zero values correctly", {
  testthat::skip_on_cran()

  # Create data where all values are zero
  BD_Obs_zero <- data.table::copy(BD_Obs)[1:2]
  BD_Coord_zero <- BD_Coord[1:4]

  # Set all values to zero
  for (col in names(BD_Obs_zero)[-1]) {
    BD_Obs_zero[, (col) := 0.0]
  }

  # This should activate: all(data_obs$var == 0) case in process_day
  out <- Kriging_Ordinary(
    BD_Obs_zero,
    BD_Coord_zero,
    shapefile,
    grid_resolution = 5,
    variogram_model = "gaussian",
    n_lags = 10,
    min_stations = 2,
    n_round = 1,
    training = 1
  )

  testthat::expect_true(inherits(out, "SpatRaster"))
  testthat::expect_equal(terra::nlyr(out), nrow(BD_Obs_zero))
})

# 22. Edge case: Very small max_dist (activates sum(valid_pairs) == 0 fragment) ----
testthat::test_that("Kriging_Ordinary handles very small max_dist correctly", {
  testthat::skip_on_cran()

  # Create data with stations that are far apart
  BD_Obs_small <- data.table(
    Date = as.Date(c("2015-01-01", "2015-01-02")),
    ST_A = c(5.0, 10.0),
    ST_B = c(8.0, 12.0),
    ST_C = c(6.5, 11.5)
  )

  # Create coordinates with large distances between stations
  BD_Coord_small <- data.table(
    Cod = c("ST_A", "ST_B", "ST_C"),
    X = c(0, 50000, 100000),     # Stations 50km apart
    Y = c(0, 50000, 100000)
  )

  # Use very small max_dist relative to station separation
  # This should activate: sum(valid_pairs) == 0 case
  testthat::expect_no_error({
    out <- Kriging_Ordinary(
      BD_Obs_small,
      BD_Coord_small,
      shapefile,
      grid_resolution = 5,
      variogram_model = "linear",
      max_dist = 100,  # Much smaller than station distances
      n_lags = 5,
      min_stations = 2,
      n_round = 1,
      training = 1
    )
  })

  testthat::expect_true(inherits(out, "SpatRaster"))
  testthat::expect_equal(terra::nlyr(out), nrow(BD_Obs_small))
})

# 23. Edge case: Numerical extremes (activates error fallback fragment) -----
testthat::test_that("Kriging_Ordinary handles numerical extremes correctly", {
  testthat::skip_on_cran()

  # Create data with extreme values that might cause numerical issues
  BD_Obs_extreme <- data.table(
    Date = as.Date(c("2015-01-01", "2015-01-02")),
    ST_A = c(5.0, 10.0),
    ST_B = c(8.0, 12.0),
    ST_C = c(6.5, 11.5)
  )
  BD_Coord_extreme <- data.table(
    Cod = c("ST_A", "ST_B", "ST_C"),
    X = c(0, 50000, 100000),     # Stations 50km apart
    Y = c(0, 50000, 100000)
  )

  # Set extreme observation values
  for (i in 2:ncol(BD_Obs_extreme)) {
    BD_Obs_extreme[[i]] <- c(1e8 * (i-1), -1e8 * (i-1))  # Very large values
  }

  # This might activate the tryCatch error handler and fallback
  testthat::expect_no_error({
    out <- Kriging_Ordinary(
      BD_Obs_extreme,
      BD_Coord_extreme,
      shapefile,
      grid_resolution = 5,
      variogram_model = "exponential",
      max_dist = 1000,  # Small distance to potentially cause issues
      n_lags = 5,
      min_stations = 2,
      n_round = 1,
      training = 1
    )
  })
})

# 24. Edge case: Test with validation to ensure coverage in validation paths --
# Test corregido para casos extremos con validación habilitada
testthat::test_that("Edge cases work correctly with validation enabled", {
  testthat::skip_on_cran()

  # Test edge case with validation to ensure validation paths are also covered
  BD_Obs_edge_val <- data.table::copy(BD_Obs)[1:5]  # More rows for validation
  BD_Coord_edge_val <- BD_Coord[1:6]  # More stations for validation

  # Obtener los nombres de las columnas de estaciones (excluyendo Date)
  station_cols <- names(BD_Obs_edge_val)[-1]

  # Verificar que tenemos exactamente 6 estaciones como esperamos
  if (length(station_cols) != 6) {
    # Si no tenemos exactamente 6, ajustar los datos
    # Tomar solo las primeras 6 estaciones
    station_cols <- station_cols[1:6]
    BD_Obs_edge_val <- BD_Obs_edge_val[, c("Date", station_cols), with = FALSE]
    BD_Coord_edge_val <- BD_Coord_edge_val[1:6]
  }

  # Create mixed scenario: some NA, some constant, some varying
  # Usando una forma más segura de asignar valores
  BD_Obs_edge_val[1, (station_cols) := list(5.0, 5.0, 5.0, NA_real_, NA_real_, 7.0)]
  BD_Obs_edge_val[2, (station_cols) := list(NA_real_, 5.0, 5.0, 8.0, NA_real_, 9.0)]
  BD_Obs_edge_val[3, (station_cols) := list(6.0, NA_real_, 5.0, 8.5, 10.0, NA_real_)]
  BD_Obs_edge_val[4, (station_cols) := list(7.0, 6.0, NA_real_, NA_real_, 11.0, 12.0)]
  BD_Obs_edge_val[5, (station_cols) := list(NA_real_, NA_real_, NA_real_, 9.0, 12.5, 13.0)]

  # Test with validation to ensure edge cases work in validation context too
  out <- Kriging_Ordinary(
    BD_Obs_edge_val,
    BD_Coord_edge_val,
    shapefile,
    grid_resolution = 5,
    variogram_model = "spherical",
    n_lags = 8,
    min_stations = 2,
    n_round = 1,
    training = 0.7,  # Enable validation
    stat_validation = NULL,
    Rain_threshold = NULL
  )

  testthat::expect_true(is.list(out))
  testthat::expect_true("Ensamble" %in% names(out))
  testthat::expect_true("Validation" %in% names(out))
  testthat::expect_true(inherits(out$Ensamble, "SpatRaster"))
})

##############################################################################
#                        END OF COVERAGE TESTS                              #
##############################################################################
# Tests adicionales

# Test 25: Caso con datos insuficientes (< 2 estaciones válidas)
test_that("Kriging handles insufficient data cases correctly", {
  # Crear datos con solo 1 estación válida
  BD_Obs_insufficient <- data.table::data.table(
    Date = as.Date(c("2015-01-01", "2015-01-02")),
    ST001 = c(5.0, 3.0),
    ST002 = c(NA, NA),  # Todas las demás estaciones son NA
    ST003 = c(NA, NA)
  )

  BD_Coord_insufficient <- data.table::data.table(
    Cod = c("ST001", "ST002", "ST003"),
    X = c(500000, 501000, 502000),
    Y = c(9500000, 9501000, 9502000)
  )

  # Ejecutar Kriging con datos insuficientes
  result <- Kriging_Ordinary(
    BD_Obs = BD_Obs_insufficient,
    BD_Coord = BD_Coord_insufficient,
    shapefile = shapefile,
    grid_resolution = 10,
    variogram_model = "exponential"
  )

  # Verificar que devuelve un resultado válido
  expect_s4_class(result, "SpatRaster")
  expect_true(terra::nlyr(result) == 2)
})

# Test 26: Caso con valores constantes (todas las estaciones tienen el mismo valor)
test_that("Kriging handles constant values correctly", {
  # Crear datos donde todas las estaciones tienen el mismo valor
  BD_Obs_constant <- data.table::data.table(
    Date = as.Date(c("2015-01-01", "2015-01-02")),
    ST001 = c(10.0, 15.0),
    ST002 = c(10.0, 15.0),  # Mismo valor que ST001
    ST003 = c(10.0, 15.0)   # Mismo valor que ST001
  )

  BD_Coord_constant <- data.table::data.table(
    Cod = c("ST001", "ST002", "ST003"),
    X = c(500000, 501000, 502000),
    Y = c(9500000, 9501000, 9502000)
  )

  # Ejecutar Kriging con valores constantes
  result <- Kriging_Ordinary(
    BD_Obs = BD_Obs_constant,
    BD_Coord = BD_Coord_constant,
    shapefile = shapefile,
    grid_resolution = 10,
    variogram_model = "exponential"
  )

  # Verificar que devuelve un resultado válido
  expect_s4_class(result, "SpatRaster")
  expect_true(terra::nlyr(result) == 2)

  # Los valores del raster deberían ser constantes (iguales al valor de entrada)
  values_layer1 <- terra::values(result[[1]], na.rm = TRUE)
  expect_true(all(abs(values_layer1 - 10.0) < 0.01, na.rm = TRUE))
})

# Test 27: Caso que fuerza el cálculo del variograma empírico con datos límite
test_that("Kriging handles edge cases in empirical variogram calculation", {
  # Crear datos con muy poca variabilidad espacial
  BD_Obs_edge <- data.table::data.table(
    Date = as.Date("2015-01-01"),
    ST001 = 0.001,  # Valores muy pequeños para forzar variograma plano
    ST002 = 0.002,
    ST003 = 0.001
  )

  BD_Coord_edge <- data.table::data.table(
    Cod = c("ST001", "ST002", "ST003"),
    X = c(500000, 500100, 500200),  # Estaciones muy cercanas
    Y = c(9500000, 9500100, 9500200)
  )

  # Ejecutar con n_lags específico para cubrir la secuencia en el código
  result <- Kriging_Ordinary(
    BD_Obs = BD_Obs_edge,
    BD_Coord = BD_Coord_edge,
    shapefile = shapefile,
    grid_resolution = 5,
    variogram_model = "exponential",
    n_lags = 10  # Específicamente para cubrir la línea del lag_distances
  )

  expect_s4_class(result, "SpatRaster")
})

# Test 28: Caso con una sola estación disponible (available_values de longitud 1)
test_that("Kriging handles single available station", {
  # Crear datos donde solo una estación tiene datos válidos por fecha
  BD_Obs_single <- data.table::data.table(
    Date = as.Date(c("2015-01-01", "2015-01-02")),
    ST001 = c(5.0, NA),    # Solo válida el primer día
    ST002 = c(NA, 8.0),    # Solo válida el segundo día
    ST003 = c(NA, NA)      # Nunca válida
  )

  BD_Coord_single <- data.table::data.table(
    Cod = c("ST001", "ST002", "ST003"),
    X = c(500000, 501000, 502000),
    Y = c(9500000, 9501000, 9502000)
  )

  result <- Kriging_Ordinary(
    BD_Obs = BD_Obs_single,
    BD_Coord = BD_Coord_single,
    shapefile = shapefile,
    grid_resolution = 8,
    variogram_model = "spherical"
  )

  expect_s4_class(result, "SpatRaster")
  expect_true(terra::nlyr(result) == 2)
})

# Test 29: Caso específico para cubrir la condición when available_values is empty
test_that("Kriging handles completely empty available values", {
  # Crear datos donde todas las estaciones son NA para una fecha
  BD_Obs_empty <- data.table::data.table(
    Date = as.Date(c("2015-01-01", "2015-01-02")),
    ST001 = c(NA, 5.0),
    ST002 = c(NA, 3.0),
    ST003 = c(NA, 4.0)
  )

  BD_Coord_empty <- data.table::data.table(
    Cod = c("ST001", "ST002", "ST003"),
    X = c(500000, 501000, 502000),
    Y = c(9500000, 9501000, 9502000)
  )

  result <- Kriging_Ordinary(
    BD_Obs = BD_Obs_empty,
    BD_Coord = BD_Coord_empty,
    shapefile = shapefile,
    grid_resolution = 10,
    variogram_model = "gaussian"
  )

  expect_s4_class(result, "SpatRaster")
  # El primer layer debería tener valor 0 (caso else en el código)
  values_layer1 <- terra::values(result[[1]], na.rm = TRUE)
  expect_true(all(values_layer1 == 0, na.rm = TRUE))
})