Why tiling?

Why tiling in the first place? Tiling allows for massively parallelizing the use of statisticalRoughness which comes particularly handy if you have large raster or if you want to explore the evolution of the scaling over a large range of scales. In this vignette, the code is not run but walk-through as the whole vignette tasks about 40 minutes to run.


The tiling has 3 main components:

  1. Resampling the target raster
  2. Making it into polygons
  3. Using the polygons to tile the target raster


Here we use functions from statisticalRoughness to derive the values of \(L\) and \(R\). Note that Lmax corresponds to a 100 km.

Lmax <- 1E4
allLR <- get_all_R_L(Lmax, 5, only = 12, len = 32)
maxL <- max(allLR$allL)

Target DEM reading

This chunk reads a target DEM and extracts its name as fname_sans_ext. terraOptions ensures that the data are not stored on a system disk.

root.dir <- get_rootdir() # a in-house function to retrieve a path
dem.dir <- 'data/california-rivers/gis-files/NED'
terraOptions(tempdir = file.path(root.dir, dem.dir), memfrac = 0.7)
fname <- "CA_DEM.grd"
fname_sans_ext <- tools::file_path_sans_ext(fname)
DEM <- terra::rast(file.path(root.dir, dem.dir, fname))


This chunks resample the target raster to the size of maxL, the biggest value of \(L\).

maxL_raster <- aggregate(DEM, fact = maxL, fun = mean)
writeRaster(maxL_raster, file.path(root.dir, dem.dir, paste0(fname_sans_ext, "_Lmax-resampled.tif")), overwrite = TRUE)

Raster to polygons

We now transform the resampled raster into polygons.

maxL_raster <- raster::raster(file.path(root.dir, dem.dir, paste0(fname_sans_ext, "_Lmax-resampled.tif")))
maxL_polygons <- raster::rasterToPolygons(maxL_raster)
# maxL_polygons <- terra:as.polygons(maxL_raster, dissolve = FALSE, values = FALSE)
raster::shapefile(maxL_polygons, file.path(root.dir, dem.dir, paste0(fname_sans_ext, "_Lmax-polygons.shp")), overwrite = TRUE)

Tiling DEM

In this chunk, the target DEM is iteravely crop’ed into tiles.

out.dir <- file.path(root.dir, dem.dir, paste0(fname_sans_ext, "_tiles"))
if (!dir.exists(out.dir)) dir.create(out.dir)
for (n in seq(1, length(maxL_polygons))){
    fname <- file.path(out.dir, paste0(fname_sans_ext, "_tile_", n, ".tif"))
    dem <- crop(DEM, maxL_polygons[n, ])
    terra::writeRaster(dem, fname, overwrite = TRUE)

Cleaning files

Removing temporary files.

tmpFiles(old = TRUE, remove = TRUE)
raster::removeTmpFiles(h = 0)