Analysing Raster Data: the rasterList Package

[This article was first published on MilanoR, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

The rasterList package has been developed to create a S4 class object that can manage complex operations and objects in a spatially gridded map, i.e. a raster map. The rasterList-Class (Cordano 2017b) object differs from a traditional raster map by the fact that each cell can contain a generic object instead of numeric values on one or more bands or layers. The raster of objects, i.e. the RasterList-class object, consists on a rasterLayer-Class (Hijmans 2016) connected to a generic list: one list element for each raster cells. It allows to write few lines of R code for complex map algebra. Furthermore, in accordance with the existing classes defined in the R environment, each cell of the raster is “occupied” by an an istanced object. This may be of utitily plenty in most applications, e.g. environmental modellig.

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 \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):

0″ />

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):

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")

Average value of Va Genuchten’s parameter per each soil type
type swc rwc alpha n m Ks_m_per_hour Ks_m_per_sec
sand 0.44 0.02 13.8 2.09 0.52153 0.15000 4.17e-05
loamy sand 0.44 0.04 11.5 1.87 0.46524 0.10000 2.78e-05
sandy loam 0.45 0.05 6.8 1.62 0.38272 0.03900 1.08e-05
loam 0.51 0.04 9.0 1.42 0.29577 0.03300 9.20e-06
silt loam 0.50 0.03 4.8 1.36 0.26471 0.00900 2.50e-06
sandy clay loam 0.40 0.07 3.6 1.56 0.35897 0.00440 1.20e-06
clay loam 0.46 0.09 3.9 1.41 0.29078 0.00480 1.30e-06
silty clay loam 0.46 0.06 3.1 1.32 0.24242 0.00130 4.00e-07
sandy clay 0.43 0.11 3.4 1.40 0.28571 0.00073 2.00e-07
silty clay 0.48 0.07 2.9 1.26 0.20635 0.00076 2.00e-07
clay 0.48 0.10 2.7 1.29 0.22481 0.00041 1.00e-07

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:


## Loading required package: raster

## Loading required package: sp

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


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):

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}{z}/{x}/{y}.png"
ref_attr <- 'Map data: © <a href="">OpenStreetMap</a>, <a href="">SRTM</a> | Map style: © <a href="">OpenTopoMap</a> (<a href="">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)


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")
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)

## 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) %>% 

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")

Points of interest: geographical coordinates (latitude andlongitude) and name
lon lat name
5.724514 50.96153 A
5.734514 50.96153 B
5.756814 50.98153 C

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:


##      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:


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

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:

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)

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)


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 <- !
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))


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


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

theta <- raster(swc_rlf(psi))


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)



The above applications illustrate how statistical analysis and other processing may be replicated at every cell of a raster map. The rasterList package has been created to answer the need to make complex operation on spatially gridded region. The shown examples go beyond their own task but they have the task of presenting the most important functions of the package and how thay can be used to solve practical problems. rasterList is able to save different type of objects within a cell of the raster map and allows that each object of a raster cell may be transformed and/or coerced to another object. Two or more RasterList-class objects can interact through operations with one another. This is of great relevance because it allows to save intermediate steps during a processing workflow (e.g. statistical analysis).


Chambers, J. M. 1992. “Linear Models.” In Statistical Models in S. Wadsworth & Brooks/Cole.

Cheng, Joe, Bhaskar Karambelkar, and Yihui Xie. 2018. Leaflet: Create Interactive Web Maps with the Javascript ’Leaflet’ Library.

Cordano, Emanuele. 2017a. LmomPi: (Precipitation) Frequency Analysis and Variability with L-Moments from ’Lmom’.

———. 2017b. RasterList: A Raster Where Cells Are Generic Objects.

Cordano, Emanuele, Daniele Andreis, and Fabio Zottele. 2017. Soilwater: Implementation of Parametric Formulas for Soil Water Retention or Conductivity Curve.

Funk, Chris, Pete Peterson, Martin Landsfeld, Diego Pedreros, James Verdin, Shraddhanand Shukla, Gregory Husak, et al. 2015. “The Climate Hazards Infrared Precipitation with Stations–a New Environmental Record for Monitoring Extremes.” Scientific Data 2 (December). The Author(s) SN -: 150066 EP.

Grolemund, Garrett, and Hadley Wickham. 2011. “Dates and Times Made Easy with lubridate.” Journal of Statistical Software 40 (3): 1–25.

Hastie, T. J., and D. Pregibon. 1992. “Generalized Linear Models.” In Statistical Models in S. Wadsworth & Brooks/Cole.

Hijmans, Robert J. 2016. Raster: Geographic Data Analysis and Modeling.

Hosking, J. R. M. 2017. L-Moments.

Hosking, J. R. M., and James R. Wallis. 1997. “L-Moments.” In Regional Frequency Analysis: An Approach Based on L-Moments, 14–43. Cambridge University Press. doi:10.1017/CBO9780511529443.004.

Kollet, Stefan, Mauro Sulis, Reed M. Maxwell, Claudio Paniconi, Mario Putti, Giacomo Bertoldi, Ethan T. Coon, et al. 2017. “The Integrated Hydrologic Model Intercomparison Project, Ih-Mip2: A Second Set of Benchmark Results to Diagnose Integrated Hydrology and Feedbacks.” Water Resources Research 53 (1): 867–90. doi:10.1002/2016WR019191.

Maidment, D. R. 1992. Handbook of Hydrology. Mc-Graw Hill.

Marsaglia, George, Wai Wan Tsang, and Jingbo Wang. 2003. “Evaluating Kolmogorov’s Distribution.” Journal of Statistical Software, Articles 8 (18): 1–4. doi:10.18637/jss.v008.i18.

Pohlert, Thorsten. 2018. Trend: Non-Parametric Trend Tests and Change-Point Detection.

VanGenuchten, M. Th. 1980. “A Closed-Form Equation for Predicting the Hydraulic Conductivity of Unsaturated Soils.” Soil Sci. Soc. Am. Jour. 44: 892–98.

The post Analysing Raster Data: the rasterList Package appeared first on MilanoR.

To leave a comment for the author, please follow the link and comment on their blog: MilanoR. offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Never miss an update!
Subscribe to R-bloggers to receive
e-mails with the latest R posts.
(You will not see this message again.)

Click here to close (This popup will not appear again)