blockCV: raster -> terra

test.R

#### blockCV old

library(raster)
x <- raster(system.file("external/rlogo.grd", package="raster")) 
resolution = 10000
ybin = 3
xbin = 3
checkerboard = FALSE

ext <- raster::extent(x)
extRef <- raster::extent(x)
if(is.na(raster::projection(x))){
  mapext <- raster::extent(x)[1:4]
  if(all(mapext >= -180) && all(mapext <= 180)){
    resolution <- resolution / 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 {
    resolution <- resolution
    warning("The input layer has no CRS defined. Based on the extent of the input map it is assumed to have a projected reference system")
  }
} else{
  if(raster::isLonLat(x)){
    resolution <- resolution / degree
  } else{
    resolution <- resolution
  }
}
if(!is.null(xbin) && is.null(ybin)){
  rasterNet <- raster::raster(ext, nrow=1, ncol=xbin, crs=raster::projection(x))
} else if(is.null(xbin) && !is.null(ybin)){
  rasterNet <- raster::raster(ext, nrow=ybin, ncol=1, crs=raster::projection(x))
} else if(!is.null(xbin) && !is.null(ybin)){
  rasterNet <- raster::raster(ext, nrow=ybin, ncol=xbin, crs=raster::projection(x))
} else if(is.null(xbin) && is.null(ybin) && !is.null(resolution)){
  xrange <- raster::xmax(x) - raster::xmin(x) # number of columns
  yrange <- raster::ymax(x) - raster::ymin(x) # number of rows
  xPix <- ceiling(xrange / resolution)
  yPix <- ceiling(yrange / resolution)
  xdif <- ((xPix * resolution) - xrange) / 2 # the difference of extent divided by 2 to split on both sides
  ydif <- ((yPix * resolution) - yrange) / 2
  ext@xmin <- raster::xmin(x) - xdif
  ext@xmax <- raster::xmax(x) + xdif
  ext@ymin <- raster::ymin(x) - ydif
  ext@ymax <- raster::ymax(x) + ydif
  if(!is.null(xOffset)){
    if(xOffset > 1 || xOffset < 0){stop("xOffset should be between 0 and 1")}
    ext@xmin <- ext@xmin + (resolution * xOffset)
    ext@xmax <- ext@xmax + (resolution * xOffset)
  }
  if(!is.null(yOffset)){
    if(yOffset > 1 || yOffset < 0){stop("yOffset should be between 0 and 1")}
    ext@ymin <- ext@ymin + (resolution * yOffset)
    ext@ymax <- ext@ymax + (resolution * yOffset)
  }
  # adding cells if needed
  if(ext@xmin > extRef@xmin){ # add one column by increasing the extent and number of bins
    ext@xmin <- ext@xmin - resolution
    xPix <- xPix + 1
  }
  if(ext@ymin > extRef@ymin){
    ext@ymin <- ext@ymin - resolution
    yPix <- yPix + 1
  }
  rasterNet <- raster::raster(ext, nrow=yPix, ncol=xPix, crs=raster::projection(x))
} else stop("A value should be specified for the block size")

if(checkerboard){
  raster::values(rasterNet) <- seq_len(raster::ncell(rasterNet))
  m <- raster::as.matrix(rasterNet)
  for(i in seq_len(raster::ncol(rasterNet))){
    if(i %% 2 == 0){
      m[,i] <- rep(1:2, nrow(m))[seq_len(nrow(m))]
    } else{
      m[,i] <- rep(2:1, nrow(m))[seq_len(nrow(m))]
    }
  }
  raster::values(rasterNet) <- m # rasterNet[] <- m
} else{
  raster::values(rasterNet) <- seq_len(raster::ncell(rasterNet))
}
rasterNet_pol <- raster::rasterToPolygons(rasterNet)

### Terra approach -------------------------------------------------------------

ras <- raster(system.file("external/rlogo.grd", package="raster")) 
spatras = terra::rast(ras)

ext_terra <- terra::ext(spatras)
extRef_terra <- terra::ext(spatras)
if(is.na(terra::crs(spatras))){
  mapext <- terra::ext(spatras)[1:4]
  if(all(mapext >= -180) && all(mapext <= 180)){
    resolution <- resolution / 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 {
    resolution <- resolution
    warning("The input layer has no CRS defined. Based on the extent of the input map it is assumed to have a projected reference system")
  }
} else{
  if(terra::is.lonlat(spatras)){
    resolution <- resolution / degree
  } else{
    resolution <- resolution
  }
}
if(!is.null(xbin) && is.null(ybin)){
  terraNet <- terra::rast(ext_terra, nrow=1, ncol=xbin, crs=terra::crs(spatras))
} else if(is.null(xbin) && !is.null(ybin)){
  terraNet <- terra::rast(ext_terra, nrow=ybin, ncol=1, crs=terra::crs(spatras))
} else if(!is.null(xbin) && !is.null(ybin)){
  terraNet <- terra::rast(ext_terra, nrow=ybin, ncol=xbin, crs=terra::crs(spatras))
} else if(is.null(xbin) && is.null(ybin) && !is.null(resolution)){
  xrange <- terra::xmax(spatras) - terra::xmin(spatras) # number of columns
  yrange <- terra::ymax(spatras) - terra::ymin(spatras) # number of rows
  xPix <- ceiling(xrange / resolution)
  yPix <- ceiling(yrange / resolution)
  xdif <- ((xPix * resolution) - xrange) / 2 # the difference of extent divided by 2 to split on both sides
  ydif <- ((yPix * resolution) - yrange) / 2
  ext_terra@xmin <- terra::xmin(spatras) - xdif
  ext_terra@xmax <- terra::xmax(spatras) + xdif
  ext_terra@ymin <- terra::ymin(spatras) - ydif
  ext_terra@ymax <- terra::ymax(spatras) + ydif
  if(!is.null(xOffset)){
    if(xOffset > 1 || xOffset < 0){stop("xOffset should be between 0 and 1")}
    ext_terra@xmin <- ext_terra@xmin + (resolution * xOffset)
    ext_terra@xmax <- ext_terra@xmax + (resolution * xOffset)
  }
  if(!is.null(yOffset)){
    if(yOffset > 1 || yOffset < 0){stop("yOffset should be between 0 and 1")}
    ext_terra@ymin <- ext_terra@ymin + (resolution * yOffset)
    ext_terra@ymax <- ext_terra@ymax + (resolution * yOffset)
  }
  # adding cells if needed
  if(ext_terra@xmin > ext_terra@xmin){ # add one column by increasing the extent and number of bins
    ext_terra@xmin <- ext_terra@xmin - resolution
    xPix <- xPix + 1
  }
  if(ext_terra@ymin > ext_terra@ymin){
    ext_terra@ymin <- ext_terra@ymin - resolution
    yPix <- yPix + 1
  }
  terraNet <- terra::ras(ext_terra, nrow=yPix, ncol=xPix, crs=terra::crs(spatras))
} else stop("A value should be specified for the block size")

if(checkerboard){
  terra::values(terraNet) <- seq_len(terra::ncell(terraNet))
  m_terra <- terra::as.matrix(terraNet)
  for(i in seq_len(terra::ncol(terraNet))){
    if(i %% 2 == 0){
      m_terra[,i] <- rep(1:2, nrow(m_terra))[seq_len(nrow(m_terra))]
    } else{
      mm_terra <- rep(2:1, nrow(m_terra))[seq_len(nrow(m_terra))]
    }
  }
  terra::values(terraNet) <- m 
} else{
  terra::values(terraNet) <- seq_len(terra::ncell(terraNet))
}
terraNet_pol <- terra::as.polygons(terraNet)

# OPTIONAL: convert terra to sf object
terraNet_sf = sf::st_as_sf(terraNet)

# OPTIONAL: compare raster and terra
rasterNet_terra = terra::rast(rasterNet)
all.equal(terraNet, rasterNet_terra) # besides the layer name, everything matches!