Skip to content

Instantly share code, notes, and snippets.

@pat-s
Last active March 29, 2023 08:05
Show Gist options
  • Save pat-s/d2ab388e3c581e6db1a57aeca2ae8e4b to your computer and use it in GitHub Desktop.
Save pat-s/d2ab388e3c581e6db1a57aeca2ae8e4b to your computer and use it in GitHub Desktop.
blockCV: raster -> terra
#### 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!
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment