R/SOptim_RasterTiling.R

In joaofgoncalves/SegOptim: Optimization of Image Segmentation Parameters for Object-Based Image Analysis (OBIA)

#' Split a vector into n chunks
#' 
#' @param x a vector
#' @param n number of chunks
#' 
#' @author 
#'   \href{http://stackoverflow.com/users/1563634/mathheadinclouds}{mathheadinclouds}, 
#'   \href{http://stackoverflow.com/users/1737569/dis-shishkov}{Dis Shishkov}
#' @references \url{http://stackoverflow.com/questions/3318333/split-a-vector-into-chunks-in-r}
#' @export
#' @examples
#'  chunk2(1:30, 6)
#'  
chunk2 <- function(x,n) split(x, cut(seq_along(x), n, labels = FALSE)) 

#' Create a raster tiles object
#'
#' A function used to cut an input \code{SpatRaster} object into several tiles 
#' usable by \code{\link[terra]{values}} to read data by chunks.
#' 
#' @param rst A \code{SpatRaster} object.
#'
#' @param nd Number of times to slice the SpatRaster across row 
#' and column direction. The total number of tiles will be given by: 
#' \eqn{N_{tiles} = nd^{2}}.
#' 
#' @return 
#' A \code{SOptim.Tiles} object with nd^2 tiles each with the starting 
#' row and column numbers as well as the number of rows and columns to read.
#' 
#' @importFrom terra nrow
#' @importFrom terra ncol
#' 
#' @export
#' 
#' @examples
#' 
#' library(terra)
#' r <- rast(ncol = 100, nrow = 100)
#' rtiles <- createRasterTiles(r, nd = 3)
#' 

createRasterTiles <- function(rst, nd = 2){
  
  if(nd < 2){
    stop("createRasterTiles error! The number nd must be greater than two.")
  }
  
  nr <- terra::nrow(rst)
  nc <- terra::ncol(rst)
  
  blkRows <- chunk2(1:nr, nd)
  blkCols <- chunk2(1:nc, nd)
  
  tiles <- list()
  
  k <- 0
  for(i in 1:nd){
    for(j in 1:nd){
      k <- k+1
      newTile <- list()
      newTile[['rowStart']] <- min(blkRows[[i]])
      newTile[['rowLen']]   <- length(blkRows[[i]])
      newTile[['colStart']] <- min(blkCols[[j]])
      newTile[['colLen']]   <- length(blkRows[[j]])
      class(newTile) <- "SOptim.Tile"
      
      tiles[[k]] <- newTile
    }
  }
  
  class(tiles) <- "SOptim.Tiles"
  attr(tiles, "NumberSlices") <- nd
  attr(tiles, "NumberTiles") <- nd^2
  attr(tiles, "NumPixelsPerTile") <- unlist(lapply(tiles, function(x) x$rowLen * x$colLen))
  attr(tiles, "NumPixelsTotal") <- sum(attr(tiles, "NumPixelsPerTile"))
  
  return(tiles)
}

#' Read raster data for a certain tile
#' 
#' An ancilary function that uses an \code{SOptim.Tiles} object (generated by 
#' [createRasterTiles()]) and reads pixel-level data to a matrix or data.frame.
#' 
#' @param rst A \code{SpatRaster} object.
#'
#' @param tiles A \code{SOptim.Tiles} object (generated by \code{\link{createRasterTiles}}. If not 
#' provided, a \code{SOptim.Tiles} will be generated on-the-fly with nd^2 tiles.
#' 
#' @param nd Number of times to slice the SpatRaster across row 
#' and column direction. The total number of tiles will be given by: 
#' \eqn{N_{tiles} = nd^{2}}. (default: 2)
#' 
#' @param tileIndex Tile index to read (default: NULL)
#' 
#' @param as.df Return a data.frame object? (default: TRUE)
#' 
#' 
#' @importFrom terra values
#' 
#' @export
#' 
#' @return 
#' A matrix or data.frame containing pixel-level data for the whole tile. If the input is 
#' a \code{SpatRaster} then each column represents each layer.
#'    

readDataByTile <- function(rst, tiles, nd = 2, tileIndex = NULL, as.df = TRUE){
  
  if(is.null(tileIndex)){
    stop("tileIndex not defined in readDataByTile!")
  }
  
  if(missing(tiles)){
    tiles <- createRasterTiles(rst, nd)
  }else{
    if(!inherits(tiles, "SOptim.Tiles")){
      stop("Input tiles is not an object of class SOptim.Tiles!")
    }
  }
  
  if(!inherits(tiles[[tileIndex]],"SOptim.Tile")){
    stop("Object loaded is not of class SOptim.Tile (readDataByTile)")
  }
  
  #print(tiles)
  # print(tiles[[tileIndex]][["rowStart"]])
  # print(tiles[[tileIndex]][["rowLen"]])
  # print(tiles[[tileIndex]][["colStart"]])
  # print(tiles[[tileIndex]][["colLen"]])
  
  out <- terra::values(rst, 
                        row   = tiles[[tileIndex]][["rowStart"]], 
                        nrows = tiles[[tileIndex]][["rowLen"]], 
                        col   = tiles[[tileIndex]][["colStart"]], 
                        ncols = tiles[[tileIndex]][["colLen"]])
  
  if(as.df){
    return(as.data.frame(out))
  }else{
    return(out)
  }
}