Computing

Hervé Guillon

08 March 2021

Description

This script is an example of the use of the function of statisticalRoughness on a tiled DEM.

Initialization

This sets up the terraOptions and a reference coordinate system.

terra::terraOptions(tempdir = "F:/tmp/")
crs_ref <- terra::crs("+proj=aea +lat_1=34 +lat_2=40.5 +lat_0=0 +lon_0=-120 +x_0=0 +y_0=-4000000 +ellps=GRS80 + +towgs84=0,0,0,0,0,0,0 +units=m +no_defs ")

This sets up the location to look for tiles.

dem.dir <- 'data/california-rivers/gis-files/NED'
fname <- "CA_DEM.grd"
fname_sans_ext <- tools::file_path_sans_ext(fname)
out.dir <- file.path(root.dir, dem.dir, paste0(fname_sans_ext, "_tiles"))

Factoring

With this, we retrieve the scales at which the DEM has been tiled.

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

Computing

Computing can be made parallel or sequential. Here we have 80 tiles and n indicates the current tile, l indicates the current scale at which the calculation takes place. Note that the fractal dimension is calculated over the last 5 scales.

for (n in seq(80)){
    print(n)
    dem <- terra::rast(file.path(out.dir, paste0(fname_sans_ext, "_tile_", n, ".tif")))
    iters <- floor(seq(from = 1, to = 32))
    for (l in iters){
        tstart <- Sys.time()
        print(l)
        R <- allLR$allRs[[l]] %>% tail(5)
        L <- allLR$allL[[l]]
        fpath <- file.path(out.dir, paste0("Hurst_raster_", n, "_",l , ".tif"))
        if (!file.exists(fpath)){
            Hurst_rasters <- compute_Hurst_rasters(dem, L, R)
            terra::writeRaster(Hurst_rasters, fpath)
        } else {
            Hurst_rasters <- terra::rast(fpath)
        }
        # toc(log = TRUE, quiet = TRUE)
        plot(Hurst_rasters)
        print(Sys.time() - tstart)
    }
    gc()
}