Geostatistics for soil mapping

DOI: doi.org/10.17027/isric-fdtz-ma12

Licence This tutorial is released under the GNU GPL v3.0 license. GNU GPL v3.0 is a strong copyleft license. This means that you may use the code and change/modify the code. If you distribute copies or modifications of the code, you are required to release these updates under the GPL v3 license.

Citation Heuvelink, G.B.M., Kempen, B., Genova, G., and Collazos Cortes, D., 2024. Geostatistics for soil mapping. Tutorial, ISRIC-World Soil Information, Wageningen. https://doi.org/10.17027/isric-fdtz-ma12

Introduction

The floodplain of the Geul river valley, located in the south of the Netherlands, is polluted by heavy metals. Historic metal mining has caused the widespread dispersal of lead, zinc and cadmium in the alluvial soil. Each time the river flooded the area, polluted sediments were deposited on the river banks. The pollutants may constrain the land use in these areas, so maps are required that delineate zones with high concentrations. In this tutorial you will create such maps from soil data collected in the Geul area. Figure 1 shows the study area and sampling locations.

Figure 1: Location of the Geul study area in The Netherlands. (a) Zoom in on the study area and (b) shows part of the Geul valley and the sampling locations. The size of the symbols is proportional to the lead concentration of the topsoil.

The tutorial consists of six chapters:

1  Analysing and modelling spatial variation Analyzing and modelling spatial variation: quantifying spatial variation with a variogram.

2  Ordinary kriging Ordinary kriging: spatial interpolation of the topsoil lead concentration using ordinary kriging.

3  Spatial stochastic simulation Spatial stochastic simulation: simulating values from the kriging probability distribution.

4  Regression kriging Regression kriging: including covariates in a geostatistical prediction model by combining linear regression with kriging.

5  Universal kriging Universal kriging: a more elegant way of regression kriging that also accounts for uncertainty in the regression model.

6  Answers Answers: answers to questions in the previous five chapters.

Note it is recommended to the run the chapters in sequence. Output objects generated in one chapter are typically used as input for the next chapter. It is possible, however, to run the chapters independently since the various output objects required by subsequent scripts are all contained in the folder with materials that is associated to this tutorial (see below for further details).

After completing this tutorial you will be able to:

• quantify and model spatial variation with a semivariogram.
• create soil maps with different forms of kriging and assess the uncertainty of the predictions.
• generate possible realities with spatial stochastic simuluation.

Materials

Download

This tutorial is supported by a set of materials that are available in a dedicated folder on ISRIC’s file service and can be accessed here. The folder contains zip file geostatistics that includes all materials. Download this file, save it in a dedicated folder on your computer and then extract the zip file. This will give you a folder named geostatistics that contains three sub-folders: code, input, output. The figure below shows the folder structure. Here one has downloaded and extracted the geostatistics zipfile in a dedicated folder named tutorials that is located on the C-drive of a Windows computer.

Figure 1.1: Folder structure.

IMPORTANT:

1. Please carefully read the README that explains the folder structure and what can be found within these folders.

2. Do not change any of the folder or file names, including the name of the ‘parent folder’ geostatistics. Changing folder and file names might cause the code stop running and produce errors.

3. Do not change any code. The entire code should execute without errors. Changing the code might cause the code to stop running and produce errors. If you wish to experiment with the code (which we encourage!), please do so in a copy of the code file. The only code line you need to modify is the setwd() function (further explanation below).

Input data

The data needed in this tutorial is stored in the input folder (all coordinates are in the Dutch grid projection system):

• geuldata.csv: CSV file with 100 lead concentration observations including their geographic coordinates;
• mask_aoi_shp: shapefile delineating the study area;
• river_line.shp: shapefile representing the river Geul;
• geul_mask.tif: raster map of the study area;
• geul_dem.tif: digital elevation model of the study area;
• geul_slope.tif: slope map of the study area derived from the digital elevation model;
• river_dist.tif: raster map with shortest distance to the Geul river.

R and Packages

Installation

The code must be executed with the R software for statistical computing. R can be downloaded here www.r-project.org.

We recommend to use a dedicated software that provides a more user-friendly Graphical User Interface (GUI) to R, such as RStudio. RStudio can be downloaded here download/rstudio-desktop.

The required R packages are:

• sf: vector (point, line and polygon) data handling
• stars: raster data handling
• gstat: spatial and spatio-Temporal Geostatistical Modelling, Prediction and Simulation

In case a package is not installed, the install.packages() function can be used to download and install it. The code example below shows how to check if a package is available on your machine and how to install a missing package.

# To check if a package is installed
system.file(package = 'gstat')

# it prints the directory where the package is installed
# otherwise, an empty directory is printed as ''

# To install a package
install.packages('stars')

Alternatively, packages can be installed using specific options in Rstudio, in case it is used to run the R code. Instructions can be found here.

Setting the working directory

To be able to execute the code you have to tell R where it can find the materials on your computer. In other words, you have to set the working directory.

This can be done with the setwd() function. Within the () you define the full file path of the folder on your computer that contains the geostatistics folder that you have downloaded.

As an example, suppose one has created a folder called tutorials on the C-drive of a Windows computer, and has saved and unzipped the downloaded geostatistics file there, then the working directory should be set as follows:

# Set the working directory pointing to the folder
# where the tutorial materials are stored
setwd("C:/tutorials/")

Note that quotation marks must be used to define the file path. Note further that the geostatistics folder must not be included in the file path if it follows the folder structure in the flowchart above.

IMPORTANT: the code line with the setwd() function is the only code line you need to modify in each of the scripts.

Alternatively, for RStudio users, the working directory can be set using specific options in the GUI. Instructions can be found here.

The working directory is set at the beginning of each chapter of this tutorial.

In case R returns an error after running the setwd() function, you probably have an error in the file path definition: check and correct.