Analysing Raster Data: the rasterList Package

Analysing Raster Data: Application 1

The application 1 is the estimation of the mathematical relationships between soil water content $\theta$ and soil water pressure $\psi$, i.e. soil water retention curves, in the cells of a raster map. Soil water retention curve is widely used in subsurface water/groundwater hydrological modeling because it allows to close the water balance equation. Generally some theoretical formulas are used to model this curve in function of soil properties. In this examples, given a raster map of soil properties, a map of soil water retention curves is produced. First of all, in absence of detailed soil information, as often used in hydrological modeling (Kollet et al. 2017), soil water curve is defined with the following empirical formula (VanGenuchten 1980):

$$\theta = \theta_{res}+(\theta_{sat}-\theta_{res})*(1+ |\alpha \psi|^n)^{-m} \qquad \psi>0$$

where $\theta$ is the volumetric water content (e.g. water volume per soil volume) ranging between residual water content $\theta_{res}$ and saturated water content $\theta_{sat}$. $\alpha$,$n$ and $m$ are parameters assuming to depend on soil type. An estimation of the parameter $\theta_{sat}$,$\theta_{res}$,$\alpha$,$m$,$n$ values in function of soil type is given through the following table (Maidment 1992):

library(soilwater)
library(stringr)
soilparcsv <- system.file("external/soil_data.csv",package="soilwater")
soilpar <- read.table(soilparcsv,stringsAsFactors=FALSE,header=TRUE,sep=",")
knitr::kable(soilpar,caption="Average value of Va Genuchten's parameter per each soil type")

soilpar$color <- str_sub(rainbow(nrow(soilpar)),end=7)  ## Only first 7 characters of HTML code are considered.

It it is assumed that each cell of a raster maps has its own soil water retention curve. Therefore to map the soil water retion curve and not only its parameters, rasterList package is loaded:

library(rasterList)

## Loading required package: raster

## Loading required package: sp

The study area is is the area of the meuse dataset (n.d.):

library(sp)
data(meuse)
help(meuse)
data(meuse.grid)
help(meuse.grid)

A soil map is avaialbe from soil type according to the 1:50 000 soil map of the Netherlands. In the Meuse site, the following soil type are present

• 1 = Rd10A (Calcareous weakly-developed meadow soils, light sandy clay);
• 2 = Rd90C/VII (Non-calcareous weakly-developed meadow soils, heavy sandy clay to light clay);
• 3 = Bkd26/VII (Red Brick soil, fine-sandy, silty light clay) which can be respectively assimilated as:

soiltype_id <- c(1,2,3)
soiltype_name <- c("sandy clay","clay","silty clay loam")

soilpar_s <- soilpar[soilpar$type %in% soiltype_name,]
is <- order(soilpar_s[,1],by=soiltype_name)
soilpar_s <- soilpar_s[is,]
soilpar_s$id <- soiltype_id

A geografic view of Meuse test case is available though leaflet package (Cheng, Karambelkar, and Xie 2018):

library(rasterList)
library(soilwater)
library(leaflet)
set.seed(1234)
data(meuse.grid)
data(meuse)
coordinates(meuse.grid) <- ~x+y
coordinates(meuse) <- ~x+y
gridded(meuse.grid) <- TRUE

proj4string(meuse) <- CRS("+init=epsg:28992")
proj4string(meuse.grid) <- CRS("+init=epsg:28992")
 
soilmap <- as.factor(stack(meuse.grid)[['soil']])


ref_url <- "https://{s}.tile.opentopomap.org/{z}/{x}/{y}.png"
ref_attr <- 'Map data: &copy; <a href="http://www.openstreetmap.org/copyright">OpenStreetMap</a>, <a href="http://viewfinderpanoramas.org">SRTM</a> | Map style: &copy; <a href="https://opentopomap.org">OpenTopoMap</a> (<a href="https://creativecommons.org/licenses/by-sa/3.0/">CC-BY-SA</a>)'
opacity <- 0.6
color <- colorFactor(soilpar_s$color,level=soilpar_s$id)
labFormat <- function(x,values=soilpar_s$type){values[as.integer(x)]}

leaf1 <- leaflet() %>% addTiles(urlTemplate=ref_url,attribution=ref_attr) 
leaf2 <- leaf1 %>% addRasterImage(soilmap,opacity=opacity,color=color) %>% 
addLegend(position="bottomright",colors=soilpar_s$color,labels=soilpar_s$type,title="Soil Type",opacity=opacity)

leaf2

Therefore, a map of each VanGenuchten (1980)’s formula parameter can be obtained as follow (The names of saturated and residual water contents swc and rwc mean $\theta_{sat}$ and $\theta_{res}$ according to the the arguments of swc function of soilwater package (Cordano, Andreis, and Zottele 2017))

soil_parameters_f <- function (soiltype,sp=soilpar_s) {
  o <- sp[soiltype,c("swc","rwc","alpha","n","m")]
  names(o) <- c("theta_sat","theta_res","alpha","n","m")
  return(o)
  }
soil_parameters_rl <- rasterList(soilmap,FUN=soil_parameters_f)

A RasterList-class object has been created and each cell contains a vector of soil water parameters:

soil_parameters <- stack(soil_parameters_rl)
soil_parameters

## class       : RasterStack 
## dimensions  : 104, 78, 8112, 5  (nrow, ncol, ncell, nlayers)
## resolution  : 40, 40  (x, y)
## extent      : 178440, 181560, 329600, 333760  (xmin, xmax, ymin, ymax)
## coord. ref. : +init=epsg:28992 +proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel +towgs84=565.4171,50.3319,465.5524,-0.398957,0.343988,-1.87740,4.0725 +units=m +no_defs 
## names       : theta_sat, theta_res,   alpha,       n,       m 
## min values  :   0.43000,   0.06000, 2.70000, 1.29000, 0.22481 
## max values  :   0.48000,   0.11000, 3.40000, 1.40000, 0.28571

A map of saturated water content (porosity) can be visualized as follows:

theta_sat <- soil_parameters[["theta_sat"]]
color <- colorNumeric("Greens",domain=theta_sat[])

leaf3 <- leaf1 %>% addRasterImage(theta_sat,color=color,opacity=opacity) %>% 
addLegend(position="bottomright",pal=color,values=theta_sat[],opacity=opacity,title="Porosity")

Let consider 3 points of interest in the Meuse located (latitude and logitude) as follows:

lat <- 50.961532+c(0,0,0.02)
lon <-  5.724514+c(0,0.01,0.0323)
name <- c("A","B","C")

points <- data.frame(lon=lon,lat=lat,name=name)
knitr::kable(points,caption="Points of interest: geographical coordinates (latitude andlongitude) and name")

coordinates(points) <- ~lon+lat
proj4string(points) <- CRS("+proj=longlat  +ellps=WGS84") 
leaf3 %>% addMarkers(lng=points$lon,lat=points$lat,label=points$name)

<

They are indexed as follows in the raster map:

points$icell <- cellFromXY(soil_parameters,spTransform(points,projection(soil_parameters)))

where icell is a vector of integer numbers corresponding to the raster cells containing the points A,B and C. Information on these points can be extracted as follows:

soil_parameters[points$icell]

##      theta_sat theta_res alpha    n       m
## [1,]      0.48      0.10   2.7 1.29 0.22481
## [2,]      0.48      0.10   2.7 1.29 0.22481
## [3,]      0.43      0.11   3.4 1.40 0.28571

Using swc , a function that given the soil parameters returns a soil water function is defined:

library(soilwater)

swc_func <- function(x,...) {
         o <- function(psi,y=x,...) {
         
              args <- c(list(psi,...),as.list(y))
              oo <- do.call(args=args,what=get("swc"))
         }     
        return(o)
}

And in case that it is applied to point A (i.e. setting the soil parameters of point A), it is:

soilparA <- soil_parameters[points$icell[1]][1,]
swcA <- swc_func(soilparA)

where scwA is a function corresponding to the Soil Water Retention Curve in point A that can be plotted as follows:

library(ggplot2)
psi <- seq(from=-5,to=1,by=0.25)
title <- "Soil Water Retention Curve at Point A"
xlab <- "Soil Water Pressure Head [m]"
ylab <- "Soil Water Content (ranging 0 to 1) "
gswA <- ggplot()+geom_line(aes(x=psi,y=swcA(psi)))+theme_bw()
gswA <- gswA+ggtitle(title)+xlab(xlab)+ylab(ylab)
gswA

The procedure is replicated for all cells of the raster map, saving the soil water retention curves as an R object for all the cells, then making use of rasterList function:

swc_rl <- rasterList(soil_parameters,FUN=swc_func)

The RasterList-class object contains each soil water retention curve per each cell. Since the list element referred to each raster cell are function class, swc_rl is coerced to a single function defined all raster extension through rasterListFun function:

swc_rlf <- rasterListFun(swc_rl)

In such a way, the function swc_rlf can be used to estimate to estimate soil water function given a spatially unimorm value of soil water pressure head $\psi$ on all Meuse region.

psi <- c(0,-0.5,-1,-2)
names(psi) <- sprintf("psi= %1.1f m",psi)
soil_water_content <- stack(swc_rlf(psi))
names(soil_water_content) <- names(psi)
plot(soil_water_content,main=names(psi))

Finally, assuming that $\psi$ may vary with space ranging between -0.5 m and -1.5 m with a spatial exponetial dacay profile from point A, a map for $\psi$ is calculated as follows:

region <- raster(swc_rlf(0))
mask <- !is.na(region)
region[] <- NA
region[points$icell[1]] <- 1
dist <- distance(region)
dist[mask==0] <- NA

mdist <- max(dist[],na.rm=TRUE)

psi_max <- 0 
psi_min <- -2
psi <- psi_max+(psi_min-psi_max)*(1-exp(-dist/mdist))
plot(psi)

color <- colorNumeric("Blues",domain=psi[])

leaf_psi <- leaf1 %>% addRasterImage(psi,color=color,opacity=opacity) %>% 
addLegend(position="bottomright",pal=color,values=psi[],opacity=opacity,title="Psi") %>% addMarkers(lng=points$lon,lat=points$lat,label=points$name)

## Warning in colors(.): Some values were outside the color scale and will be
## treated as NA

leaf_psi

The obtained map of $\theta$ is then visualized through leaflet:

theta <- raster(swc_rlf(psi))
plot(theta)

color <- colorNumeric("Blues",domain=theta[])

leaf_psi <- leaf1 %>% addRasterImage(theta,color=color,opacity=opacity) %>% 
addLegend(position="bottomright",pal=color,values=theta[],opacity=opacity,title="Psi") %>% addMarkers(lng=points$lon,lat=points$lat,label=points$name)

leaf_psi