GEMAS: Spatial analysis of the Ni distribution on a continental-scale using digital image processing techniques on European agricultural soil data

Highlights

• Digital image processing techniques and spatial modelling used to plot geochemical maps for the first time

• Continental-scale Ni concentration patterns are revealed in the topsoil data not seen before

• Multi-disciplinary GEMAS geochemical data reveals topsoil – bedrock interaction

• New research opened by digital terrain analysis of continental-scale geochemical bivariate surface data

Abstract

This study demonstrates the use of digital image processing for the spatial pattern recognition and characterisation of Ni concentrations in topsoil in Europe. Moving average smoothing was applied to the TIN-interpolated grid model to suppress small irregularities. Digital image processing was applied then to the grid. Several NE-SW, E-W and NW-SE oriented features were revealed at the continental scale. The dominant NE-SW linear features follow the Variscan and Alpine orogenies. The highest variability zones are in the Alps and the Balkans where mafic and ultramafic rocks outcrop. A single major E-W oriented north-facing feature runs along the last continental glaciation zone. This zone also coincides with a series of local maxima in Ni concentration along the glaciofluvial deposits. The NW-SE elongated features are located in the Pyrenees, northern Italy, Hellas and Fennoscandia. This study demonstrates the advantages of digital image processing analysis in identifying and characterising spatial geochemical patterns unseen before on conventional colour-surface maps.

Introduction

The Geochemical Mapping of Agricultural and Grazing Land Soil (GEMAS) project was established in 2008 as a joint project of the Geochemistry Expert Group of EuroGeoSurveys and Eurometaux (European Association of Metals) (Reimann et al., 2014a, Reimann et al., 2014b). The project aimed at providing harmonised geochemical background data according to the specifications of the new European Chemicals Regulation, known as REACH (Registration, Evaluation and Authorisation of Chemicals; EC, 2006). In this framework, 2108 samples of agricultural soil (Ap, 0–20 cm, regularly ploughed fields), and 2023 samples from land under permanent grass cover (Gr, 0–10 cm, grazing land soil) were collected across almost the whole European continent, at an average density of 1 sample site/2500 km², in accordance with a commonly agreed sampling protocol (EuroGeoSurveys Geochemistry Working Group, 2009). The spatial distribution of chemical elements in the GEMAS database has been investigated by several studies using interpolated concentration maps and geostatistical analysis in order to identify anomalous patterns in relation to bedrock geology and anthropogenic contamination sources (Reimann et al., 2012, Tarvainen et al., 2013, Fabian et al., 2014, Reimann et al., 2014a, Reimann et al., 2014b, Albanese et al., 2015, Birke et al., 2014, Birke et al., 2017).

In general, spatial geochemical data have long been analysed with various graphical techniques, such as proportional dots (Björklund and Gustavsson, 1987), and statistical methods, namely classification based on percentiles, boxplots (Kürzl, 1988, O'Connor and Reimann, 1993) and cumulative probability (CP) and empirical cumulative distribution function (ECDF) plots (Tennant and White, 1959, Lepeltier, 1969, Sinclair, 1974, Sinclair, 1976, Sinclair, 1983, Sinclair, 1986, Sinclair, 1991). Element concentration maps in regular grid are constructed with interpolation methods, such as polynomials functions over irregularly spaced data (Akima, 1996) or with kriging based on semi-variogram models (Reimann et al., 2012, Gosar et al., 2016, Birke et al., 2017). The advanced method of multi-fractal mapping has been widely used to decompose spatial geochemical data into anomalies and background values (Agterberg, 2001, Cheng et al., 1994, Cheng, 1999a, Cheng, 1999b, Lima et al., 2003, Zuo et al., 2015, Zuo and Wang, 2016). Additional methods are the spatially weighted singularity mapping (Xiao et al., 2017) or the spatially weighted principal component analysis for multi-element geochemical data (Cheng et al., 2011). Local neighbourhood analysis (LNA) is another widely-used method to characterise spatial geochemical patterns (Zhang et al., 2007, Cheng, 2007, Carranza, 2009, Xie et al., 2008, Zuo, 2014). Shahi et al. (2015) proposed the coupling of discrete wavelet transforms (DWT) and principal component analysis (PCA) for predicting mineralisation. Apart from a few early attempts of digital image processing of geochemical maps, such as the pioneering work by Chork and Cruikshank (1984), no systematic application of digital image processing has been applied to spatial geochemical data, despite the extensive mapping and map interpretation in applied geochemical surveys.

Geochemical maps interpolated from the original sampling points to a regular grid can be viewed as raster images and, hence, processed using digital image processing procedures to increase the apparent distinction between spatial features. There are two basic types of procedures, global and local operations. Global operations (also called point operations) act on the whole image, on all the pixels at the same time, and modify each pixel value independently from its neighbouring pixels. Local operations (also called spatial operations), on the contrary, modify each pixel value based on the neighbouring pixel values. Global operations use the image histogram and have basically two applications. Histogram or density slicing is used to classify data based on frequency distribution characteristics (Lillesand et al., 2015). Slicing of an image histogram by dividing pixel values (element concentrations in the case of geochemical maps) into specified intervals is used in the present study to display discrete categories of element concentrations, such as anomalous and background values. Contrast stretching methods (often called normalisation) enhance image contrast by expanding the range of grey levels (colours) assigned to image pixel values. These two techniques are often used in combination where grey levels of values within a class, such as the low concentration geochemical background, are streached (also called level slicing or piecewise linear contrast stretching) in order to reveal subtle concentration patterns on the geochemical map. However, contrast enhancement should be applied only after error treatment to avoid ‘over-contrasting’ by assigning different display values to pixels differing only within the error range.

Local operations, also called neighbourhood operations, change a particular pixel value based on the surrounding pixel values. The ‘surrounding’ or neighbourhood is defined by a rectangular grid (matrix; also called window or kernel) of pixels centred over the target pixel, which is slid from pixel to pixel to make the same calculation for each pixel on the image, hence the name ‘moving window’ operations. The size of the moving window is defined by the number of rows and columns of the constituting pixels, such as 3 × 3 or 9 × 9 window-size, for example. Local operations include the important class of the so called ‘low pass’ and ‘high pass’ filters, besides the many and diverse procedures. Low pass filters allow the low frequency features, such as a large area of geochemical anomalies to ‘pass’ or remain in the image, and they filter out or remove the high frequency components, such as random noise. One of the most frequently used and simplest low pass filter is the moving average filter that calculates the average value within a moving kernel and replaces the original value of the centre pixel with the calculated window average. Average filters are used to reduce noise and smooth small irregularities that can reveal broad, regional trends. High pass filters, on the contrary, allow the high frequency features, such as sudden changes in concentration values to pass or remain in the image and they remove the low frequency components, such as the spatial trend. One of the most frequently used and simplest high pass filter is the gradient filter that calculates the concentration change from pixel to pixel.

Digital image processing is only a sub-area of the broad field of mathematics called signal processing. In this study, only a few specific signal processing procedures are used, such as the digital identification of local extreme values (minima and maxima) in the Ni concentration geochemical map, local variability characterisation (relief calculation) and auto-correlation (variogram) analysis. Finally, an interesting aspect of the digital image processing of an interpolated raster geochemical map is considering the raster map as a continuous surface of a bivariate function, much like the digital terrain model of the topographical surface. Hence, the geochemical maps can be analysed for ‘slope’, i.e., the rate of concentration change (gradient) or for ‘aspect’, i.e., the direction of the maximum change to see, for example, if the regional change of concentration has a certain geographical orientation (Jordan, 2007, Jordan et al., 2005).

The objective of this study is to demonstrate the use of digital image processing techniques for reproducible spatial geochemical pattern recognition and quantitative spatial feature characterisation. The concentrations of Ni in samples of agricultural topsoil from the GEMAS project database (Reimann et al., 2014a) are used as an example to perform the detailed spatial analysis, and to relate the spatial features to possible underlying geological and geochemical processes.