R/cv_spatial.R

In blockCV: Spatial and Environmental Blocking for Cross-Validation

#' Use spatial blocks to separate train and test folds
#'
#' This function creates spatially separated folds based on a distance to number of row and/or column.
#' It assigns blocks to the training and testing folds \strong{randomly}, \strong{systematically} or
#' in a \strong{checkerboard pattern}. The distance (\code{size})
#' should be in \strong{metres}, regardless of the unit of the reference system of
#' the input data (for more information see the details section). By default,
#' the function creates blocks according to the extent and shape of the spatial sample data (\code{x} e.g.
#' the species occurrence), Alternatively, blocks can be created based on \code{r} assuming that the
#' user has considered the landscape for the given species and case study.
#' Blocks can also be offset so the origin is not at the outer corner of the rasters.
#' Instead of providing a distance, the blocks can also be created by specifying a number of rows and/or
#' columns and divide the study area into vertical or horizontal bins, as presented in Wenger & Olden (2012)
#' and Bahn & McGill (2012). Finally, the blocks can be specified by a user-defined spatial polygon layer.
#'
#' To maintain consistency, all functions in this package use \strong{meters} as their unit of
#' measurement. However, when the input map has a geographic coordinate system (in decimal degrees),
#' the block size is calculated by dividing the \code{size} parameter by \code{deg_to_metre} (which
#' defaults to 111325 meters, the standard distance of one degree of latitude on the Equator).
#' In reality, this value varies by a factor of the cosine of the latitude. So, an alternative sensible
#' value could be \code{cos(mean(sf::st_bbox(x)[c(2,4)]) * pi/180) * 111325}.
#'
#' The \code{offset} can be used to change the spatial position of the blocks. It can also be used to
#' assess the sensitivity of analysis results to shifting in the blocking arrangements.
#' These options are available when \code{size} is defined. By default the region is
#' located in the middle of the blocks and by setting the offsets, the blocks will shift.
#'
#' Roberts et. al. (2017) suggest that blocks should be substantially bigger than the range of spatial
#' autocorrelation (in model residual) to obtain realistic error estimates, while a buffer with the size of
#' the spatial autocorrelation range would result in a good estimation of error. This is because of the so-called
#' edge effect (O'Sullivan & Unwin, 2014), whereby points located on the edges of the blocks of opposite sets are
#' not separated spatially. Blocking with a buffering strategy overcomes this issue (see \code{\link{cv_buffer}}).
#' When using \code{\link{cv_spatial_autocor}} to inform \code{size}, note that response or covariate variograms
#' are exploratory proxies for residual autocorrelation unless model residuals have been added to \code{x} and
#' supplied through the \code{column} argument.
#'
#' @details
#' For presence-background data (\code{presence_bg = TRUE}), \code{column} holds \code{1} for presences and
#' \code{0} for \emph{background} points -- locations sampled across the study area to represent the available
#' conditions rather than confirmed absences. The fold balancing then targets only the presence records (the
#' blocks still contain all points), so the abundant background cannot dominate the split.
#'
#'
#' @param x a simple features (sf) or SpatialPoints object of spatial sample data (e.g., species data or ground truth sample for image classification).
#' @param column character (optional). Indicating the name of the column in which response variable (e.g. species data as a binary
#' response i.e. 0s and 1s) is stored. It is used to report the number of records of each class/bin in every fold and,
#' when \code{balance = TRUE}, to balance those classes across the folds. If \code{column = NULL}
#' the response variable classes will be treated the same and only training and testing records will be counted (and balanced).
#' This is used for binary (e.g. presence-absence/background) or multi-class responses (e.g. land cover classes for
#' remote sensing image classification). Continuous numeric responses are binned into quantiles using \code{num_bins}
#' before balancing.
#' @param num_bins integer; the number of quantile bins used to stratify a continuous numeric \code{column}.
#' The default is \code{4}. Set \code{num_bins = NULL} to disable binning and treat every unique value as a
#' separate class (the behaviour prior to version 3.3). If quantile breaks are tied, fewer bins may be used.
#' The raw response values are not modified; bins are only used for fold balancing and record summaries.
#' @param r a terra SpatRaster object (optional). If provided, its extent will be used to specify the blocks.
#' It also supports \emph{stars}, \emph{raster}, or path to a raster file on disk.
#' @param k integer value. The number of desired folds for cross-validation. The default is \code{k = 5}.
#' @param hexagon logical. Creates hexagonal (default) spatial blocks. If \code{FALSE}, square blocks is created.
#' @param flat_top logical. Creating hexagonal blocks with topped flat.
#' @param size numeric value of the specified range by which blocks are created and training/testing data are separated.
#' This distance should be in \strong{metres}. The range could be explored by \code{\link{cv_spatial_autocor}}
#' and \code{\link{cv_block_size}} functions. For residual-based block-size guidance, add model residuals to
#' \code{x} and pass that residual column to \code{\link{cv_spatial_autocor}}.
#' @param rows_cols integer vector. Two integers to define the blocks based on row and
#' column e.g. \code{c(10, 10)} or \code{c(5, 1)}. Hexagonal blocks uses only the first one. This
#' option is ignored when \code{size} is provided.
#' @param selection type of assignment of blocks into folds. Can be \strong{random} (default), \strong{systematic}, \strong{checkerboard}, or \strong{predefined}.
#' The checkerboard does not work with hexagonal and user-defined spatial blocks. If the \code{selection = 'predefined'}, user-defined
#' blocks and \code{folds_column} must be supplied.
#' @param iteration integer value. The number of attempts to create folds with balanced records. Only works when \code{selection = "random"}.
#' @param balance logical. When \code{TRUE} (default) and \code{selection = "random"}, the folds are chosen from \code{iteration}
#' random block assignments to balance the training/testing records (or the classes/bins of \code{column} when it is provided).
#' If \code{FALSE}, a single random assignment is returned without balancing (equivalent to \code{iteration = 1}). This argument
#' has no effect on the other selection methods.
#' @param presence_bg logical; whether to treat \code{column} as species presence-background data (0s for
#' background points and 1s for presences; see \sQuote{Details}). When \code{TRUE} (with \code{selection = "random"} and
#' \code{balance = TRUE}), the balancing search equalises only the presence (1s) records across folds so the many
#' background points cannot dominate the objective; the blocks still contain all points but the background is
#' ignored when scoring the balance. Requires a binary numeric \code{column}. The default is \code{FALSE}.
#' @param user_blocks an sf or SpatialPolygons object to be used as the blocks (optional). This can be a user defined polygon and it must cover all
#' the species (response) points. If \code{selection = 'predefined'}, this argument and \strong{folds_column} must be supplied.
#' @param folds_column character. Indicating the name of the column (in \code{user_blocks}) in which the associated folds are stored.
#' This argument is necessary if you choose the 'predefined' selection.
#' @param deg_to_metre integer. The conversion rate of metres to degree. See the details section for more information.
#' @param biomod2 logical. Creates a matrix of folds that can be directly used in the \pkg{biomod2} package as
#' a \emph{CV.user.table} for cross-validation.
#' @param offset two number between 0 and 1 to shift blocks by that proportion of block size.
#' This option only works when \code{size} is provided.
#' @param extend numeric; This parameter specifies the percentage by which the map's extent is
#' expanded to increase the size of the square spatial blocks, ensuring that all points fall
#' within a block. The value should be a numeric between 0 and 5.
#' @param seed integer; a random seed for reproducibility (although an external seed
#' should also work).
#' @param plot logical; whether to plot the final blocks with fold numbers in ggplot.
#' Defaults to \code{interactive()}.
#' You can re-create this with \code{\link{cv_plot}}.
#' @param report logical; whether to print the report of the records per fold.
#' Defaults to \code{interactive()}.
#' @param progress logical; whether to shows a progress bar for random fold selection.
#' Defaults to \code{interactive()}.
#' @param ... additional option for \code{\link{cv_plot}}.
#'
#'
#' @seealso \code{\link{cv_buffer}} and \code{\link{cv_cluster}}; \code{\link{cv_spatial_autocor}} and \code{\link{cv_block_size}} for selecting block size;
#' \code{\link{cv_plot}} to visualise, and \code{\link{cv_distance}} and \code{\link{cv_similarity}} to evaluate, the folds
#' @seealso For \emph{CV.user.table} see \code{\link[biomod2]{BIOMOD_Modeling}} in \pkg{biomod2} package
#'
#' @references Valavi, R., Elith, J., Lahoz-Monfort, J. J., & Guillera-Arroita, G. (2019). blockCV: An R package for generating spatially or environmentally separated folds for k-fold cross-validation of species distribution models. Methods in Ecology and Evolution, 10(2), 225-232. doi:10.1111/2041-210X.13107.
#'
#' Bahn, V., & McGill, B. J. (2012). Testing the predictive performance of distribution models. Oikos, 122(3), 321-331.
#'
#' O'Sullivan, D., Unwin, D.J., (2010). Geographic Information Analysis, 2nd ed. John Wiley & Sons.
#'
#' Roberts et al., (2017). Cross-validation strategies for data with temporal, spatial, hierarchical,
#' or phylogenetic structure. Ecography. 40: 913-929.
#'
#' Wenger, S.J., Olden, J.D., (2012). Assessing transferability of ecological models: an underappreciated aspect of statistical
#' validation. Methods Ecol. Evol. 3, 260-267.
#'
#' @return An object of class S3. A list of objects including:
#'    \itemize{
#'     \item{folds_list - a list containing the folds. Each fold has two vectors with the training (first) and testing (second) indices}
#'     \item{folds_ids - a vector of values indicating the number of the fold for each observation (each number corresponds to the same point in species data)}
#'     \item{biomod_table - a matrix with the folds to be used in \pkg{biomod2} package}
#'     \item{k - number of the folds}
#'     \item{size - input size, if not null}
#'     \item{block_shape - the block geometry used: \code{"hexagon"}, \code{"square"}, or \code{"user-defined"}}
#'     \item{selection - how blocks were assigned to folds: \code{"random"}, \code{"systematic"}, \code{"checkerboard"}, or \code{"predefined"}}
#'     \item{column - the name of the column if provided}
#'     \item{blocks - spatial polygon of the blocks}
#'     \item{records - a table with the number of points in each category of training and testing}
#'     }
#' @export
#'
#' @examples
#' \donttest{
#' library(blockCV)
#'
#' # import presence-absence species data
#' points <- read.csv(system.file("extdata/", "species.csv", package = "blockCV"))
#' # make an sf object from data.frame
#' pa_data <- sf::st_as_sf(points, coords = c("x", "y"), crs = 7845)
#'
#' # hexagonal spatial blocking by specified size and random assignment
#' sb1 <- cv_spatial(x = pa_data,
#'                   column = "occ",
#'                   size = 450000,
#'                   k = 5,
#'                   selection = "random",
#'                   iteration = 50)
#'
#' # spatial blocking by row/column and systematic fold assignment
#' sb2 <- cv_spatial(x = pa_data,
#'                   column = "occ",
#'                   rows_cols = c(8, 10),
#'                   k = 5,
#'                   hexagon = FALSE,
#'                   selection = "systematic")
#'
#' }
cv_spatial <- function(
        x,
        column = NULL,
        r = NULL,
        k = 5L,
        hexagon = TRUE,
        flat_top = FALSE,
        size = NULL,
        rows_cols = c(10, 10),
        selection = "random",
        iteration = 100L,
        balance = TRUE,
        presence_bg = FALSE,
        user_blocks = NULL,
        folds_column = NULL,
        deg_to_metre = 111325,
        biomod2 = TRUE,
        offset = c(0, 0),
        extend = 0,
        seed = NULL,
        num_bins = 4L,
        plot = interactive(),
        report = interactive(),
        progress = interactive(),
        ... # other arguments for cv_plot
){

    # pre-run checks ----------------------------------------------------------

    # check for availability of ggplot2
    if(plot) .check_pkgs(c("ggplot2"))
    # check for selection arg
    selection <- match.arg(selection, choices = c("random", "systematic", "checkerboard", "predefined"))
    # check x is an sf object
    x <- .check_x(x)
    # is column in x?
    column <- .check_column(column, x)
    # validate presence-background data (0/1 column) when requested
    invisible(.presence_index(x, column, presence_bg))

    # check for user_blocks format
    if(!is.null(user_blocks)){
        user_blocks <- .check_x(user_blocks, name = "user_blocks")
        # limit user defined blocks for checkerboard selection
        if(selection=="checkerboard"){
            warning("The checkerboard selection cannot be used with 'user_blocks`.\nThe random selection is used!")
            selection <- "random"
        }
    }
    # checks for pre-defined folds
    if(selection == "predefined"){
        if(is.null(folds_column) || is.null(user_blocks)){
            stop("The 'user_blocks' and 'folds_column' should be specified for 'predefined' selection")
        }
        if(!folds_column %in% colnames(user_blocks)){
            stop(sprintf("There is no column named '%s' in 'user_blocks'.\n", folds_column))
        }
        if(!is.numeric(user_blocks[,folds_column, drop = TRUE])){
            stop("The fold numbers in 'folds_column' must be integer numbers.")
        }
    }

    # change the r to terra object
    if(!is.null(r)){
        r <- .check_r(r)
        r <- r[[1]]
        # check points fall within the raster extent
        .check_within(x, r)
    }

    # if hex; selection muse random or systematic
    if(hexagon && selection %in% c("checkerboard", "predefined")){
        selection <- "random"
        message("Hexagon blocks can only be used with random or systematic selections!\nThe random selection is used.")
    }

    if(selection=="checkerboard") k <- 2

    tryCatch(
        {
            extend <- abs(extend) # this correctly identifies invalid values rather than max(0, val)
            extend <- min(5, extend) / 100
        },
        error = function(cond) {
            message("'extend' must be a numeric value between 0 and 5.")
        }
    )

    # iterations --------------------------------------------------------------

    # The iteration must be a natural number
    tryCatch(
        {
            iteration <- abs(as.integer(iteration))
            iteration <- max(1, iteration)
        },
        error = function(cond) {
            message("'iteration' must be a natural number.")
        }
    )

    # balance only applies to the random selection search
    if(selection == "random"){
        if(!balance){
            iteration <- 1L
        } else if(iteration == 1L){
            warning("'balance = TRUE' has no effect when 'iteration = 1'; increase 'iteration' to search for balanced folds.")
        }
    }

    # turn off progress if...
    if(selection != "random") progress <- FALSE
    if(iteration < 3) progress <- FALSE

    # creating blocks ---------------------------------------------------------

    if(is.null(user_blocks)){
        # create rectangular and hexagonal spatial blocks using terra and sf packages
        blocks <- .make_blocks(
            x_obj = if(is.null(r)) x else r, # select the object to make grid
            blocksize = size,
            blockcols = rows_cols[2],
            blockrows = rows_cols[1],
            hexagonal = hexagon,
            flat_top = flat_top,
            extend_perc = extend,
            degree = deg_to_metre,
            xy_offset = offset,
            checkerboard = ifelse(selection == "checkerboard", TRUE, FALSE)
        )
    } else{
        # make sure user_blocks is a data.frame/sf
        blocks <- if(methods::is(user_blocks, "sfc")) sf::st_sf(user_blocks) else user_blocks
    }

    ## subset the blocks by x and keep the intersection result for records
    blocks_intersect <- sf::st_intersects(sf::st_geometry(blocks), sf::st_geometry(x))
    blocks_has_records <- lengths(blocks_intersect) != 0
    sub_blocks <- blocks[blocks_has_records, ]
    blocks_len <- nrow(sub_blocks)

    # k must be a natural number
    tryCatch(
        {
            k <- abs(as.integer(k))
        },
        error = function(cond) {
            message("'k' must be a natural number.")
        }
    )
    # check k is not larger than len blocks
    if(k > blocks_len){
        stop("'k' is bigger than the number of spatial blocks: ", blocks_len, ".\n")
    } else if(k < 2){
        stop("'k' must be a natural number equal or higher than 2.")
    }

    # x and block intersection ------------------------------------------------

    ## use the block-to-point intersection computed above
    sub_blocks_intersect <- blocks_intersect[blocks_has_records]
    records_per_block <- lengths(sub_blocks_intersect)
    blocks_df <- data.frame(
        records = unlist(sub_blocks_intersect, use.names = FALSE),
        block_id = rep.int(seq_len(blocks_len), records_per_block)
    )
    blocks_df <- blocks_df[order(blocks_df$records, blocks_df$block_id), ]
    row.names(blocks_df) <- NULL
    # randomly remove the repeated records occurred on the edges of blocks
    if(nrow(blocks_df) > nrow(x)){
        if(!is.null(seed)){
            set.seed(seed)
        }
        blocks_df <- blocks_df[sample(nrow(blocks_df)), ]
        blocks_df <- blocks_df[!duplicated(blocks_df$records), ]
    } else if(nrow(blocks_df) < nrow(x) || anyNA(blocks_df)){
        nonoverlap <- nrow(x) - nrow(blocks_df)
        warning("At least ", nonoverlap, " of the points are not within the defined spatial blocks!\n")
        message("Consider using the 'extend' parameter to ensure points are covered by blocks e.g. extend = 0.5.")
    }

    # creating folds ----------------------------------------------------------
    # create records table
    response <- .column_response(x, column, num_bins = num_bins)

    if(selection == "random"){
        # search random block-to-fold assignments for the most balanced split
        res <- .balance_folds(
            blocks_df = blocks_df,
            blocks_len = blocks_len,
            k = k,
            iteration = iteration,
            response = response,
            seed = seed,
            biomod2 = biomod2,
            progress = progress,
            opt_cols = .balance_opt_cols(response, presence_bg)
        )
        blocks_df <- res$blocks_df
        train_test_table <- res$records
        fold_list <- res$folds_list
        fold_vect <- res$folds_ids
        biomod_table <- res$biomod_table

        # map the best folds back to the blocks
        sub_blocks$block_id <- seq_len(nrow(sub_blocks))
        blocks_df_filter <- blocks_df[, c("block_id", "folds")]
        blocks_df_filter <- blocks_df_filter[!duplicated(blocks_df_filter), ]
        sub_blocks <- merge(x = sub_blocks, y = blocks_df_filter, by = "block_id", all.x = TRUE)

    } else{
        # systematic, checkerboard and predefined selections need only one round
        train_test_table <- .records_table(k, response)
        biomod_table <- data.frame(RUN1 = rep(TRUE, nrow(blocks_df)))

        if(selection == "systematic"){
            if(hexagon){
                sub_blocks <- .fold_assign(sf::st_geometry(sub_blocks), n = k)
            } else{
                sub_blocks$block_id <- seq_len(blocks_len)
                sub_blocks$folds <- rep(1:k, length.out = blocks_len)
            }
            fold_df <- sf::st_drop_geometry(sub_blocks)

        } else if(selection == "checkerboard"){
            sub_blocks$folds <- sub_blocks$id
            sub_blocks$block_id <- seq_len(blocks_len)
            fold_df <- sf::st_drop_geometry(sub_blocks)

        } else if(selection == "predefined"){
            fold_df <- data.frame(block_id = seq_len(blocks_len), folds = sub_blocks[, folds_column, drop = TRUE])
        }

        # merge block-df once after the selection
        blocks_df <- merge(x = blocks_df, y = fold_df, by = "block_id", all.x = TRUE)

        # count the number of points in each fold
        fold_list <- list()
        fold_vect <- rep(NA, nrow(blocks_df))
        for(p in seq_len(k)){
            train_set <- blocks_df$records[which(blocks_df$folds != p)]
            test_set <- blocks_df$records[which(blocks_df$folds == p)]
            fold_vect[test_set] <- p
            fold_list[[p]] <- assign(paste0("fold", p), list(train_set, test_set))
            train_test_table <- .records_table_row(train_test_table, p, train_set, test_set, response)
            if(biomod2){ # creating a biomod2 CV.user.table for validation
                colm <- paste0("RUN", p)
                biomod_table[, colm] <- FALSE
                biomod_table[train_set, colm] <- TRUE
            }
        }
    }
    if(report){
        cat("\n")
        .print_column_bins(response)
        print(train_test_table)
    }
    # throw a warning if there are folds with zero cases
    if(any(train_test_table < 1)){
        zerofolds <- which(apply(train_test_table, 1, function(x) any(x < 1)))
        zero_text <- "class(es)"
        if(!is.null(response$bins)) zero_text <- "class/bin(es)"
        if(length(zerofolds) > 1){
            warning("Folds ", paste(zerofolds, collapse = ", "), " have ", zero_text, " with zero records")
        } else{
            warning("Fold ", zerofolds, " has ", zero_text, " with zero records")
        }
    }
    # remove the NA blocks; not for user-blocks
    if(is.null(user_blocks)){
        sub_blocks <- sub_blocks[stats::complete.cases(sub_blocks$folds), ]
    }

    # final objects for exporting
    final_objs <- list(
        folds_list = fold_list,
        folds_ids = fold_vect,
        biomod_table = if (biomod2) as.matrix(biomod_table) else NULL,
        k = k,
        size = size,
        block_shape = if(!is.null(user_blocks)) "user-defined" else if(hexagon) "hexagon" else "square",
        selection = selection,
        column = column,
        presence_bg = presence_bg,
        blocks = sub_blocks,
        records = train_test_table
    )

    class(final_objs) <- c("cv_spatial")

    # plot with the cv_plot function
    if(plot){
        plot(
            cv_plot(
                cv = final_objs,
                r = r,
                ...
            )
        )
    }

    return(final_objs)
}


#' @export
#' @method print cv_spatial
print.cv_spatial <- function(x, ...){
    desc <- if(!is.null(x$block_shape)){
        sprintf("spatial blocking (%s blocks)", x$block_shape)
    } else {
        "spatial blocking"
    }
    details <- list("Folds" = x$k)
    if(!is.null(x$size)) details[["Block size (m)"]] <- round(x$size)
    if(!is.null(x$selection)) details[["Fold selection"]] <- x$selection
    if(!is.null(x$column)) details[["Balancing column"]] <- x$column
    details[["Presence-background"]] <- if(isTRUE(x$presence_bg)) "yes" else "no"
    .print_cv_folds(x, desc, details)
}


#' @export
#' @method plot cv_spatial
plot.cv_spatial <- function(x, y, data = NULL, ...){
    if(!missing(y) || !is.null(data)){
        return(.plot_cv_fold_map(
            cv = x,
            y = if(!missing(y)) y else NULL,
            data = data,
            has_y = !missing(y),
            ...
        ))
    }

    if("folds" %in% names(x$blocks)){
        plot(x$blocks["folds"])
    } else{
        plot(x$blocks)
    }
    message("Supply sample data (e.g. plot(cv, samples)) to plot each fold with cv_plot.")
    invisible(x$blocks)
}


#' @export
#' @method summary cv_spatial
summary.cv_spatial <- function(object, ...){
    cat("Number of records in each training and testing fold:\n")
    print(object$records)
}


# create rectangular and hexagonal spatial blocks using terra and sf packages
.make_blocks <- function(
        x_obj,
        blocksize = NULL,
        blockcols = NULL,
        blockrows = NULL,
        hexagonal = FALSE,
        flat_top = FALSE,
        extend_perc = 0.005,
        degree = 111325,
        xy_offset = c(0, 0),
        checkerboard = FALSE){
    # xpoints and rasters inputs are checked in the parent function (cv_spatial);
    # so no need to check them here;

    # check if size/row/col are not all null
    if(all(sapply(list(blocksize, blockcols, blockrows), is.null))){
        stop("Size or the number of rows/columns should be defined for making spatial blocks.")
    }
    # extract the extent of the layer
    mapext <- terra::ext(x_obj)[1:4]
    # make sure offset is working only when size is provided
    if(is.null(blocksize)){
        xy_offset <- c(0, 0)
    } else{
        # make sure the shift match terra blocks
        if(hexagonal){
            if(!xy_offset[1] %in% 0:1) xy_offset[1] <- 1 - xy_offset[1]
            if(!xy_offset[2] %in% 0:1) xy_offset[2] <- 1 - xy_offset[2]
        }
        # ignore the integer part
        tryCatch(
            {
                xy_offset <- blocksize * (abs(xy_offset) %% 1)
            },
            error = function(cond) {
                message("Offsets should be numeric values between 0 and 1. For any higher values, only the decimal part is used.")
            }
        )
        # and make x offset for vector of length 1
        if(length(xy_offset) < 2) xy_offset[2] <- 0

        # adjust blocksize & infer the wgs crs when missing
        if(is.na(sf::st_crs(x_obj))){
            if(all(mapext >= -180) && all(mapext <= 180)){
                blocksize <- blocksize / degree
                warning("The input layer has no CRS defined. Based on the extent of the input map it is assumed to have an un-projected reference system.")
            }
        } else{
            # terra::is.lonlat is not working for sf objects currently
            if(sf::st_is_longlat(x_obj)){
                blocksize <- blocksize / degree
            }
        }
    }

    if (hexagonal) {
        # prepare offset values for hexagon
        xm <- as.numeric(sf::st_bbox(x_obj)[1])
        ym <- as.numeric(sf::st_bbox(x_obj)[2])
        xoff <- xm - xy_offset[1]
        yoff <- ym - xy_offset[2]
        # calculate the hexagon size with size or nrow parameter similar to sf package
        hexsize <- ifelse(is.null(blocksize),
                          diff(sf::st_bbox(x_obj)[c(1, 3)]) / blockrows,
                          blocksize)
        # make the hexagonal blocks
        tryCatch(
            {
                fishnet_poly <- sf::st_make_grid(
                    x_obj,
                    cellsize = hexsize,
                    offset = c(xoff, yoff),
                    square = FALSE,
                    what = "polygons",
                    flat_topped = flat_top
                )
            },
            error = function(cond) {
                message("Could not create spatial blocks! possibly because of using a very small block size.")
                message("Remember, size is in metres not the unit of the CRS.")
            }
        )
        # return sf/data.frame object
        fishnet_poly <- sf::st_sf(fishnet_poly)
        sf::st_geometry(fishnet_poly) <- "geometry"
        fishnet_poly$id <- seq_len(nrow(fishnet_poly))

    } else{
        # keep the reference extent to advise on adding extra cells
        ref_ext <- mapext
        # add 1% on both side to guarantee all points fall inside blocks
        xrange <- mapext["xmax"] - mapext["xmin"] # number of columns
        yrange <- mapext["ymax"] - mapext["ymin"] # number of rows
        mapext["xmin"] <- mapext["xmin"] - (xrange * extend_perc) # parenthesis for readability of code
        mapext["xmax"] <- mapext["xmax"] + (xrange * extend_perc)
        mapext["ymin"] <- mapext["ymin"] - (yrange * extend_perc)
        mapext["ymax"] <- mapext["ymax"] + (yrange * extend_perc)
        # make blocks based on blocksize and possible offsets
        if(!is.null(blocksize)){
            xPix <- ceiling(xrange / blocksize)
            yPix <- ceiling(yrange / blocksize)
            # calculate extra extent and divided by 2 to split on both sides
            xdif <- ((xPix * blocksize) - xrange) / 2
            ydif <- ((yPix * blocksize) - yrange) / 2
            # add both extra extent and offset to mapext
            mapext["xmin"] <- mapext["xmin"] - xdif + xy_offset[1]
            mapext["xmax"] <- mapext["xmax"] + xdif + xy_offset[1]
            mapext["ymin"] <- mapext["ymin"] - ydif + xy_offset[2]
            mapext["ymax"] <- mapext["ymax"] + ydif + xy_offset[2]
            # adding cells if needed and adjust the extent
            if(mapext["xmin"] > ref_ext["xmin"]){ # add one column by increasing the extent and number of bins
                mapext["xmin"] <- mapext["xmin"] - blocksize
                xPix <- xPix + 1
            }
            if(mapext["ymin"] > ref_ext["ymin"]){
                mapext["ymin"] <- mapext["ymin"] - blocksize
                yPix <- yPix + 1
            }
            # now update blockcols and blockrows for making the blocks
            # in this case priority is with blocksize rather than blockcols/blockrows
            blockcols <- xPix
            blockrows <- yPix
        }

        # make the raster blocks
        fishnet <- terra::rast(
            terra::ext(mapext),
            nrows = max(1, blockrows, na.rm = TRUE), # make sure it runs when one is provided
            ncols = max(1, blockcols, na.rm = TRUE),
            crs = terra::crs(x_obj)
        )

        # add default cell values
        terra::values(fishnet) <- seq_len(terra::ncell(fishnet))
        # make checkerboard folds
        if(checkerboard){
            net_rows <- seq_len(terra::nrow(fishnet))
            net_cols <- seq_len(terra::ncol(fishnet))
            net_ncol <- terra::ncol(fishnet)
            for(i in net_rows){
                row_cells <- terra::cellFromRowCol(fishnet, row = i, col = net_cols)
                if(i %% 2 == 0){
                    fishnet[row_cells] <- rep(1:2, length.out = net_ncol)
                } else{
                    fishnet[row_cells] <- rep(2:1, length.out = net_ncol)
                }
            }
        }

        fishnet_poly <- terra::as.polygons(fishnet, dissolve = FALSE)
        names(fishnet_poly) <- "id"
        fishnet_poly <- sf::st_as_sf(fishnet_poly)
    }

    return(fishnet_poly)
}


# generate fold number in systematic selection for hexagonal blocks
.fold_assign <- function(blocks, n, checkerboard=FALSE){
    # solve problems of digits in R
    old <- options("digits") # save away original options
    options(digits = 22) # change the option
    on.exit(options(old))
    # compute centroids
    cent <- sf::st_centroid(blocks)
    xy <- as.data.frame(sf::st_coordinates(cent))
    # to avoid problem of digits precision
    xy$X <- as.factor(xy$X)
    xy$Y <- as.factor(xy$Y)
    # get the dimension of blocks
    # xlev <- levels(xy$X)
    ylev <- levels(xy$Y)
    ny <- length(ylev)

    len <- nrow(xy)
    xy$ids <- seq_len(len)
    xy <- xy[order(xy$Y), ]
    xy$z <- 0

    if(checkerboard){
        for(i in rev(seq_len(ny))){
            wyi <- which(xy$Y == ylev[i])
            nx <- length(wyi)
            if(i %% 2){
                xy$z[wyi] <- rep(1:2, length.out = nx)
            } else{
                xy$z[wyi] <- rep(2:1, length.out = nx)
            }
        }
    } else{
        xy$z <- rep(1:n, length.out = len)
    }

    blocks <- sf::st_sf(blocks)
    blocks$block_id <- 1:nrow(blocks)
    xy <- xy[order(xy$ids), ]
    blocks$folds <- xy$z

    return(blocks)
}