Analysing Raster Data: the rasterList Package

Analysing Raster Data: an introduction

Nowadays the availability of spatially gridded coverages , often obtained by distributed biophysical/environmental models or remote sensed observations, e. g. CHIRPS rainfall (Funk et al. 2015), lead to the fact that normal operations on time series analysis and statistics must be replicated in each pixel or cell of the spatially gridded domain (raster). Based on raster package (Hijmans 2016), a S4 class has been created such that results of complex operations or speficfic R objects (e.g, S3 or S4) can be executed on each cells of a raster map. The raster of objects contains the traditional raster map with the addition of a list of generic objects: one object for each raster cells. In such a way, different kinds of objects, e.g. the probability distrubution of the mapped random variable or an htest, lm ,glm or function class objects (Chambers 1992)(Hastie and Pregibon 1992), can be saved for each cell, in order that the user can save intermediate results without loosing any geographic information. This new oblect class is implemented in a new package rasterList (Cordano 2017b). Two applications are presented:

• Creation of a 2D/3D modelled map of soil water retention curve based on existing soil properties map and state-of-art empirical formulation (VanGenuchten 1980);
• Processing spatial gridded datasets of daily or momthly precipitation (Funk et al. 2015) for a trend or frequancy analysis (estimation of the return period of critical events).

Analysing Raster Data: Application 1

The application 1 is the estimation of the mathematical relationships between soil water content and soil water pressure , 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):

0″ />

where is the volumetric water content (e.g. water volume per soil volume) ranging between residual water content and saturated water content . , and are parameters assuming to depend on soil type. An estimation of the parameter ,,,, 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 and 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 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 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 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 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