Available water-holding capacity (AWC): computation

Load packages, read data and source custom scripts

rm(list = ls())
library(stars)

#> Loading required package: abind

#> Loading required package: sf

#> Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1

path_proj <- day2day::git_path()
path_data <- file.path(path_proj, "data")
path_cleaned <- file.path(path_data, "cleaned")
path_processed <- file.path(path_data, "processed")

source(file.path(path_proj, "src", "61-vis-stars.R"))

uk <- st_read(file.path(path_cleaned, "uk_simple.gpkg"), "union")

#> Reading layer `union' from data source `/home/rstudio/Documents/Projects/land-suitability/data/cleaned/uk_simple.gpkg' using driver `GPKG'
#> Simple feature collection with 1 feature and 0 fields
#> geometry type:  MULTIPOLYGON
#> dimension:      XY
#> bbox:           xmin: -116.1923 ymin: 7054.099 xmax: 655644.8 ymax: 1218625
#> projected CRS:  OSGB 1936 / British National Grid

aws_file <- file.path(path_cleaned, "ukcp18_uk_1km_soilparams_aws_hwsd_bc.tif")
aws <- read_stars(aws_file)
aws_thick <- readRDS(sub("\\.tif$", "\\.RDS", aws_file))

Compute available water-holding capacity (AWC)

Compute available water-holding capacity based on the thickness of each layer. We use only the three main layers of available water storage. Those with thickness equals to 0.1, 0.25, 0.65. We change the units to millimeters to be comparable with precipitation and potential evapotranspiration.

aws <- aws[, , , 1:3]
thick <- aws_thick$thick[1:3]
awc <- aws[, , ,1] * thick[1] + aws[, , , 2] * thick[2] + aws[, , , 3] * thick[3]
awc <- awc / sum(thick)
awc <- adrop(awc * 1000) #  to mm
write_stars(awc, file.path(path_processed, "ukcp18_uk_1km_soilparams_awc_bc.tif"))

Visualize

awc <- read_stars(file.path(path_processed, "ukcp18_uk_1km_soilparams_awc_bc.tif"))
plot_awc(awc, uk)

Time to execute the task

Only useful when executed with Rscript.

proc.time()

#>    user  system elapsed 
#>  12.396   0.482  12.821