R/crossValidate.R

In rLakeHabitat: Interpolate Bathymetry and Quantify Physical Aquatic Habitat

#' Cross Validate Interpolated Bathymetry
#'
#' Obtain residual mean square error (RMSE) from K-fold cross validation of bathymetry interpolation.
#'
#' @param outline shapefile outline of a waterbody. Accepts a SpatVector, an sf object, or anything terra::vect() can read.
#' @param df dataframe of coordinates and depths for a given waterbody. Coordinates are assumed to be in the same CRS as 'outline'.
#' @param x character giving name of longitude column
#' @param y character giving name of latitude column
#' @param z character giving name of depth column
#' @param zeros logical describing if bounding zeros are needed (FALSE) or provided (TRUE), default = FALSE
#' @param separation number describing distance between points, in meters
#' @param k numeric value describing the number of folds to test, default = 5
#' @param res number describing desired cell resolution in meters, default = 5
#' @param seed optional numeric value used to seed the random number generator, so that fold assignment (and therefore the resulting RMSE) is reproducible across runs. Default = NULL (not seeded).
#' @details
#' Folds are assigned in two stages: points are first split into 5 depth strata (quintiles of observed depth), then within
#' each stratum, points are spatially clustered into 'k' groups (via k-means on their coordinates) and each spatial cluster
#' becomes one fold. This keeps depth ranges reasonably balanced across folds while avoiding the optimistic bias that comes
#' from randomly scattering spatially autocorrelated points across folds.
#' 'res' is required and is always in meters, regardless of the CRS 'outline' was originally supplied in.
#' @param method character describing method of interpolation, options include Inverse Distance Weighted ("IDW"), Ordinary Kriging ("OK"), or Universal Kriging ("UK"). Default = "IDW"
#' @param nmax numeric value describing number of neighbors used in interpolation, default = 20
#' @param idp numeric value describing inverse distance power value for IDW interpolation
#' @param model character describing type of model used in Ordinary Kriging, options include 'Sph', 'Exp', 'Gau', 'Sta', default = 'Sph'
#' @param psill numeric value describing the partial sill value for OK interpolation, default = NULL
#' @param range numeric describing distance beyond which there is no spatial correlation in Ordinary Kriging models, default = NULL
#' @param nugget numeric describing variance at zero distance in Ordinary/Universal Kriging models, default = 0
#' @param kappa numeric value describing model smoothness, default = NULL
#' @param trend_order numeric value (1 or 2) giving the polynomial trend order for Universal Kriging ("UK" only), default = 1
#' @param zero_threshold numeric proportion (0-1) of surface area that must interpolate to exactly 0 before the automatic zero re-interpolation pass runs - passed through to interpBathy(). Default = 0.05.
#' @details
#' For the model argument there are four different methods included here that are supported by gstat::vgm ("Sph", "Exp", "Gau", "Mat").
#' "Sph" = The default gstat::vgm method. Spherical model characterized by a curve that rises steeply to defined range then flattens, indicates no spatial correlation between points beyond that range.
#' "Exp" = Exponential model characterized by spatial correlation decaying rapidly with distance, results in a rougher surface.
#' "Gau" = Gaussian model similar to spatial model but with slower decay over distance, results in a smoother surface.
#' "Mat" = Matern model that uses kappa to define the variogram relationship. High kappa values approach a Guassian model (smooth surface), and low kappa values approach the Exponential model (kappa = 0.5 is equivalent to Exponential).
#' Three parameters (psill, range, kappa) are incorporated from a fitted variogram (default = NULL). If specified in function input, chosen values will overwrite variogram values.
#'
#' @return a named numeric value giving the mean RMSE across k folds
#' @author Tristan Blechinger, Department of Zoology & Physiology, University of Wyoming
#' @export
#' @import dplyr
#' @rawNamespace import(terra, except = c(union,intersect, animate))
#' @examples
#' #load example outline
#' outline <- terra::vect(system.file("extdata", "example_outline.shp", package = 'rLakeHabitat'))
#' #load example xyz data
#' data <- read.csv(system.file("extdata", "example_depths.csv", package = 'rLakeHabitat'))
#' #run function
#' crossValidate(outline, data, "x", "y", "z", zeros = FALSE, separation = 10, k = 5,
#' res = 50, method = "IDW", nmax = 4, idp = 1.5, seed = 123)

crossValidate <- function(outline, df, x, y, z, zeros = FALSE, separation = NULL, k = 5, res = 5, seed = NULL, method = "IDW", nmax = 20, idp = 2, model = "Sph", psill = NULL, range = NULL, nugget = 0, kappa = NULL, trend_order = 1, zero_threshold = 0.05){

  #seed the RNG (fold assignment involves random spatial-cluster sampling
  #and k-means initialization) so results are reproducible across runs
  if(!is.null(seed)){
    if(!is.numeric(seed))
      stop("seed must be numeric")
    set.seed(seed)
  }


  #transform outline shapefile into vector
  if(!inherits(outline, "SpatVector")){
    outline <- terra::vect(outline)
  }
  else{
    outline <- outline
  }

  #checks
  if(!inherits(df, "data.frame"))
    stop("df must be a dataframe")
  if(!inherits(x, "character"))
    stop("x must be a character giving the longitude column name")
  if(!inherits(y, "character"))
    stop("y must be a character giving the latitude column name")
  if(!inherits(z, "character"))
    stop("z must be a character giving the depth column name")
  if(x %in% names(df) == FALSE)
    stop("The value of x does not appear to be a valid column name")
  if(y %in% names(df) == FALSE)
    stop("The value of y does not appear to be a valid column name")
  if(z %in% names(df) == FALSE)
    stop("The value of z does not appear to be a valid column name")
  if(!inherits(df[, x], "numeric"))
    stop("data in x column is not formatted as numeric")
  if(!inherits(df[, y], "numeric"))
    stop("data in y column is not formatted as numeric")
  if(!inherits(df[, z], "numeric"))
    stop("data in z column is not formatted as numeric")
  if(!inherits(outline, "SpatVector"))
    stop("outline is not a SpatVector or cannot be transformed")
  if(is.na(k) || is.null(k))
    stop("k must be defined as a numeric value")
  if(!is.numeric(k))
    stop("k must be defined as a numeric value")

  #rest of data checks from interpBathy function

  max_depth <- max(df[[z]])

  #assign depth strata (quintiles of observed depth)
  df <- df %>%
    dplyr::mutate(depth_stratum = dplyr::case_when(between(!!sym(z), 0, max_depth*.2) ~ 1,
                                                   dplyr::between(!!sym(z), max_depth*.2, max_depth*.4) ~ 2,
                                                   dplyr::between(!!sym(z), max_depth*.4, max_depth*.6) ~ 3,
                                                   dplyr::between(!!sym(z),  max_depth*.6, max_depth*.8) ~ 4,
                                                   dplyr::between(!!sym(z), max_depth*.8, max_depth) ~ 5))

  #within each depth stratum, spatially cluster points into k groups (rather
  #than randomly assigning fold membership) so that folds are spatially blocked
  assign_spatial_folds <- function(coords, k){
    n <- nrow(coords)
    if(n <= k){
      #not enough points in this stratum to cluster meaningfully
      return(sample(rep_len(1:k, n)))
    }
    km <- stats::kmeans(coords, centers = k)
    return(km$cluster)
  }

  df$fold <- NA_integer_
  for(s in unique(df$depth_stratum)){
    idx <- which(df$depth_stratum == s)
    df$fold[idx] <- assign_spatial_folds(as.matrix(df[idx, c(x, y)]), k)
  }

  k_values <- list()

  #conduct cross val
  for (i in 1:k) {
    testDat <- df %>% dplyr::filter(.data$fold == i)
    trainDat <- df %>% dplyr::filter(.data$fold != i) %>%
      as.data.frame()
    trainDat[[x]] <- as.numeric(trainDat[[x]])

    dem <- interpBathy(outline, trainDat, x = x, y = y, z = z, zeros = zeros, separation = separation,
                       res = res, method = method, nmax = nmax, idp = idp, model = model, psill = psill,
                       range = range, nugget = nugget, kappa = kappa, trend_order = trend_order, zero_threshold = zero_threshold)

    dem <- dem[[1]]

    preds <- terra::extract(dem, testDat[, c(x, y)], ID=F)
    testDat <- base::cbind(testDat, preds)

    testDat <- testDat %>%
      dplyr::rename(zpred = ncol(testDat))

    testDat <- testDat %>%
      dplyr::mutate(diff = (.data[[z]] - .data$zpred)^2) %>% ##
      as.data.frame() ##

    total <- base::sum(testDat$diff, na.rm=T)

    k_values[[i]] <- base::sqrt(total/nrow(testDat))
  }

  rmse_value <- mean(unlist(k_values))

  return(c(RMSE = rmse_value))
}