# Predicting Water Permeability of the Soil Based on Open Data

Jonne Pohjankukka, Paavo Nevalainen, Tapio Pahikkala, Eija Hyvönen, Pekka Hänninen, Raimo Sutinen, Jari Ala-Ilomäki, Jukka Heikkonen

► To cite this version:

Jonne Pohjankukka, Paavo Nevalainen, Tapio Pahikkala, Eija Hyvönen, Pekka Hänninen, et al. Predicting Water Permeability of the Soil Based on Open Data. 10th IFIP International Conference on Artificial Intelligence Applications and Innovations (AIAI), Sep 2014, Rhodes, Greece. pp.436-446, 10.1007/978-3-662-44654-6_43. hal-01391345

HAL Id: hal-01391345

https://inria.hal.science/hal-01391345

Submitted on 3 Nov 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Predicting water permeability of the soil based on open data

Jonne Pohjankukka, Paavo Nevalainen, Tapio Pahikkala, Eija Hyvönen, Pekka Hänninen, Raimo Sutinen, Jari Ala-Ilomäki, and Jukka Heikkonen

University of Turku, Computer Science Dept.
{Jonne.Pohjankukka,ptneva,Tapio.Pahikkala}@utu.fi,
{Eija.Hyvonen,Pekka.Hanninen,Raimo.Sutinen}@gtk.fi,
Jukka.Heikkonen@utu.fi, jari.ala-ilomaki@metla.fi,
http://www.utu.fi/en/units/sci/units/it/

**Abstract.** Water permeability is a key concept when estimating load bearing capacity, mobility and infrastructure potential of a terrain. Northern sub-arctic areas have rather similar dominant soil types and thus prediction methods successful at Northern Finland may generalize to other arctic areas. In this paper we have predicted water permeability using publicly available natural resource data with regression analysis. The data categories used for regression were: airborne electro-magnetic and radiation, topographic height, national forest inventory data, and peat bog thickness. Various additional features were derived from original data to enable better predictions. The regression performances indicate that the prediction capability exists up to 120 meters from the closest direct measurement points. The results were measured using leave-one-out cross-validation with a dead zone between the training and testing data sets.

**Key words:** load bearing capacity of soil, water permeability, regression, k-nearest neighbor, mobility, sub-arctic infrastructure.

## 1 Introduction

This paper is about predicting the water permeability of the soil by regression analysis using publicly available multi-source data. Water permeability (also called hydraulic conductivity) is a central soil property related to soil type and soil texture. High permeability means that soil tends to stay dry and traversable most of the year, whereas low permeability creates a risk for mobility when precipitation is high. Mobility in arctic areas is of great interest to many different parties. E.g. the mining industry is interested about the mobility estimates when placing various facilities. The forest industry is interested on the load bearing capacity of the soil, since the route solutions can be adaptive to mobility predictions.

Our input data set consists of 44 features which are publicly available. The data is in raster format with grid resolution ranging from 10 meters to 50 meters.
Water permeability of the soil has been measured in 1788 test spots at North- ern Finland provided by the Geological Survey of Finland (GTK). It is an im- portant attribute which, when combined with other features available, helps to determine the soil types. Related studies have been conducted in [1] where soil respiration rates are predicted from temperature, moisture content and soil type. Another related research was published in the paper of P. Scull, J. Franklin and O.A. Chadwick [2]. In their paper they use classification tree analysis for pre- dicting the soil type in desert landscapes. R.P.O. Schulte et al.[3] focuses on soil moisture deficit, which is a related concept but not of concern in sub-polar areas, H. Gao et al. [4] and R.A. Chapuis [5] focus on water budget model- ing, which was not yet attempted in our study. H.S. Mahmood et al. [6] uses on-site gamma-ray measurements for analysis of the farming soil. N.J. McKen- zie [7] combines gamma-ray and digital elevation model to predict the chemical composition of the farming land. Closest to our paper is [8], where several data sources (topographic and remote sensing) are combined with 85 soil samples to assess the usability of the soil within and outside the sampled area. One can try to by-pass the water permeability estimation by directly learning the dynamic coupling between the precipitation and remotely observed soil moisture. This approach must include the digital terrain model (DTM) to estimate the water catchment. An example of this approach is [9].

The main novelty of this paper related to the previous studies is that the prediction is based on wide-area public data on a subpolar region. The features used in this paper are basically available through-out the arctic zone.

We use regression analysis to find a mapping between the publicly available data and water permeability of the soil. In the following, we present the regression methods in Ch. 2. Then we introduce the test area, the original data sets and derived features (Ch. 3) and describe the analysis process and results of the analysis (Ch. 4). The last part is for conclusions and future approaches (Ch. 5).

## 2 Regression methods

Regularized least squares regression (RLS) is well known so we describe it mainly to introduce the variables and the notation used later in the paper. The explana- tory variables $x_1, ..., x_p$ consist of given data and dependent variable $y$ is the water permeability. We need to find a set of parameters $\mathbf{w} \in \mathbb{R}^p$ and $b \in \mathbb{R}$ such that the error function:

$$E(\mathbf{w}) = \frac{1}{n} \sum_{i=1}^{n} (y_i - \mathbf{w}^T \mathbf{x}_i - b)^2 + \frac{\lambda}{n} \mathbf{w}^T \mathbf{w} \quad (1)$$

is minimized, where $\mathbf{x}_i \in \mathbb{R}^p$ is the input vector, $y_i \in \mathbb{R}$ is the response value, $n$ is the number of observations and $\lambda$ is the regularization parameter.

The k-nearest neighbors (k-NN) approach predicts the test sample by taking the average from $k$ points nearest to it. Euclidean distance is the used metric in our analysis. Explicitly stated, if $y_1, ..., y_k$ are the response values of the k- nearest points to the test sample, then the response value for the test sample $\hat{y}_t$
is:

$$ \hat{y}_t = \frac{1}{k} \sum_{i=1}^{k} y_i . $$

## 3 Test area, data sets and features

The research area is located in the northern part of the municipality of So-dankylä, which is a part of Finnish Lapland. The size of the target area is 18432 km². The center point of the rectangular target area is at ETRS-TM35FIN coordinates 7524 kmN, 488 kmE, zone 35.

The data set consists of aerial gamma-ray spectroscopy data (referred later as gamma-ray data, AGR) combined with electromagnetic (AEM), topographical (Z), peat bog mask (PBM) and The National Forest Inventory 2011 (VMI<fnref n="1" />) data when predicting the qualities and characteristics of the soil, namely its type and water permeability (WP). Gamma-ray data is inversely related on the amount of water on the soil, which can be used to predict the type of the soil. The forest inventory data describes the profile of tree species, their maturity and foresting state. Albeit this kind of data is not directly available elsewhere in northern sub-arctic areas (e.g. Russia, Canada), several studies are underway to predict the main characteristics of the forest by remote measurement methods [10]. These methods include LiDAR and various satellite measurements.

The data providers are:

<table><caption>Table 1: Data providers, data and the grid size.</caption><thead><tr><td>Provider</td><td>Data</td><td>Grid size</td></tr></thead><tbody><tr><td>Geological Survey of Finland (GTK)</td><td>AGR, AEM<br>WP</td><td>50 m</td></tr><tr><td>Finnish Forest Research Institute (Metla)</td><td>VMI, PBM</td><td>20 m</td></tr><tr><td>National Land Survey of Finland (NLS)</td><td>Z</td><td>10 m</td></tr></tbody></table>

When considering all the derived features used in the analysis we get a total of 96 data layers.

The test site has 1788 sample points, where many mechanical and electro-chemical properties of the soil were measured, see [11]. The water permeability is a theoretical value derived from the soil particle size distribution of the soil.

We now present our data sources and donors.

### 3.1 Forest inventory data

The National Forest Inventory (VMI) holds the state of Finnish forests. The data is updated once in two years. The parameters are derived from various remote sensing sources, and several spot-wise verification and calibration methods are

<footnote n="1">¹ VMI2011: http://www.metla.fi/ohjelma/vmi/vm11-info-en.html</footnote>
applied to it before publishing the data [12]. 44 numerical features include green mass, trunk dimensions and tree density per specie category. These multi-source features exhibit built-in dependencies, thus the final number of useful features is lower.

## 3.2 Aerial gamma-ray data

The aerial gamma-ray data was provided by the Geological Survey of Finland (GTK). The raster data is based on gamma-ray flux from potassium, which is the decay process of the naturally occurring chemical element potassium (K). This data indicates many significant characteristics of the soil, including the tendency to stay moist after precipitation and tendency to frost heaving. Also the soil type, especially density, porosity, grain size and humidity of the soil have an effect to gamma-ray radiation. In Fig. 1 we present the gamma-ray data from Sodankylä target area. The bright end of the gray scale is for the high gamma radiation and hence less water in the locality of the pixel.

<figure><img src="image_2.png" /><figcaption>Fig. 1: Aerial data: gamma-ray (left) and electromagnetic data (right). Air-borne electromagnetic data is sensitive to geological properties to depth of hundreds of meters, but it also indicates some features of the top soil.</figcaption></figure>

## 3.3 Electromagnetic properties of soil

The air-borne electromagnetic (AEM) data was provided by the Geological Survey of Finland (GTK). Primary AEM components, in-phase and quadrature, were transformed to apparent resistivity values by using a half-space model [13]. The apparent resistivity gives information on different kind of soil conductors. The apparent resistivity is governed by grain size distribution, water and electronic conductors content of soil and cumulative weathering.

## 3.4 Topographical height data

Topographical data provided by the National Land Survey of Finland (NLS) was included in the analysis. The data from NLS server is basically similar to aerial
laser measurements (LiDAR) except LiDAR can reach denser grid. Instead of raw height alone we used local height difference, flow accumulation area, confluence and inclination described in [14]. These four derived features are more efficient for prediction than raw height data alone.

## 3.5 Peat bog mask

Peat bog mask is created from GTK aero-radiometric data and is courtesy of NLS and METLA. The grid size is 20 m and the value 1 indicates that peat thickness is over 60 cm. Value 0 indicates thickness less than 60 cm. The limit chosen is practical for mobility prediction.

## 3.6 Derived features

The following features were derived from gamma-ray and electromagnetic data:

- Mean and variance over 3 × 3 window

- Mean and variance over Gabor filter with 8 orientations, see [15]

- Local Binary Pattern (LBP) with pixel radii $r \in \{1, 2\}$, see [16]

From topographical height we derived the following features: local height difference, ground inclination, convergence index and flow accumulation area. The definition of these features is at [17]).

There are several additional attributes possible to derive from topographical height data, and more geomorphological features will be employed in the future.

The regression methods use total of 44 original and 52 derived features, including the constant feature. The derived features are useful only if the original feature is continuous enough. E.g. the Forest Inventory data often has locally constant zones with abrupt changes and the derived features do not help much.

## 3.7 Water permeability exponent

This is the subject of prediction. Basically, the water permeability indicates the nominal vertical speed of water through the soil sample. The measurement of this quantity is indirect, based on soil particle size distribution, and the actual speed highly depends on the inhomogeneities (roots, rocks) and micro-cracks in the soil. This is why this quantity is descriptive and theoretical. In our analysis we are using a logarithmic quantity $x_{wp}$ derived from water permeability speed $v$. For purposes of this presentation it is called as the water permeability exponent and defined as:

$$x_{wp} = -\log_{10} v, \quad [v] = \frac{m}{sec}, \qquad (2)$$

This formula has $v$ as the vertical speed of water flow through the soil.
# 4 Analysis and results

We are looking for methods which predict water permeability on areas, where there may not be direct water permeability measurements nearby. Therefore, we developed a modification of the leave-one-out cross-validation (LOOCV) for measuring the degree of spatial dependency from the nearby direct measurements, which we refer to as LOOCV with dead zone. Namely, the approach works on the measurement data just like an ordinary LOOCV in which each measurement at a time is omitted from the training set and used as a test point, except that we also remove from the training set all points that are within geographical distance *r* from the test point. This approach is illustrated in Fig. 2. By varying *r*, we can measure how far from the test area we assume the closest measurements to be at the very least. In addition, the results can be helpful in deciding how dense grid of direct measurement one should use in order to obtain a certain level of prediction performance.

We perform the regression of water permeability with the following three feature sets:

- location only

- features + location

- features only

where location refers to the geographical coordinates (e.g. latitude and longitude) and features to the ones described in Section 3. Note that one can not rely on the location information if there are no nearby direct measurements at all, and therefore we measure the prediction performance separately with these.

The prediction performances with the different feature sets as a function of the radius *r* of the dead zone are depicted in the two leftmost graphs in Fig. 3 on p. 8. The generic version based on feature data only gives weaker results, since the sample point arrangement at Sodankylä (see sample sets A and B in Fig. 2) and perhaps the phenomenon itself induce spatial dependency. No good generic regression method for this data set has been found, instead the problem is about how much additional samples are needed per target area to make the prediction useful.

The common k-NN method has one essential parameter, the number of neighbors *k*. The spatial dependency can be probed by adding the dead zone radius *r* to avoid the optimistic effect of the nearest neighbors. Fig. 2 depicts the modified leave-one-out arrangement, where *k* nearest points outside the dead zone of radius *r* are used for teaching. By varying *r* one gets a varied data set and a rough estimate on how dense it should be for it to predict well in new circumstances.

The same parameterized dead zone leave-one-out arrangement was used with regression, too.

## 4.1 Predicting water permeability

As mentioned in Sec. 3.7 before, the prediction subject is the water permeability exponent $x_{wp}$ defined by Eq. 2. The values used for regularization parameter $\lambda$ ranged from $2^{-15}$, ..., $2^{15}$. k-NN parameter had $k \in \{1, 3, 6, 12, 22\}$. Two different
<figure><img src="image_3.png" /><figcaption>Fig. 2: Left: 1788 sample points. Set A (1187 points, marked with red circles, distance to the nearest neighbor $d_{NN} \le 86m$) is tightly packed and set B is very sparse (601 points, marked with black dots, aver. $d_{NN} \approx 1.1 km$). Right: the dead zone (with radius $r$) around the leave-one-out test point (black circle). The gray circles are omitted from the training set (white circles). Both the k-NN and RLS method address the training data only, e.g. the $k$ nearest neighbours are selected from outside the circle.</figcaption></figure>

error measures were used for estimating prediction performance: mean absolute error (MAE) and concordance index (CI) [18]. Explicitly, the error measures are:

$$MAE = \frac{1}{n} \sum_{i=1}^{n} \left| \frac{y_i - \hat{y}_i}{y_i} \right| \quad CI = \frac{1}{N} \sum_{y_i < y_j} h(\hat{y}_i - \hat{y}_j) \quad (3)$$

MAE prediction baseline $\tilde{y}$ is the best possible prediction under the assumption that the prediction will be constant, thus constraining all values $\hat{y}_j = \tilde{y}, j = 1..n$ in the error minimization process. MAE baseline becomes thus:

$$MAE_b = \arg \min_{\tilde{y}} \frac{1}{n} \sum_{i=1}^{n} \left| \frac{y_i - \tilde{y}}{y_i} \right|$$

The prediction performance should be better than this to be useful. The corresponding percentage values (MAPE and $MAPE_b$) have been used in the rest of the text.

Concordance index counts the occurrences when the prediction fails to be monotonical. In equation (3) we denote $N = |\{(i, j) | y_i > y_j\}|$ as the normalization constant which equals to the number of data pairs with different label values. $h(.)$ is Heaviside step function.

The values of the C-index range between 0.0 and 1.0, where 0.5 corresponds to a random predictor and 1.0 to the case where prediction is monotonically correct.

## 4.2 Results

The results for regression analysis can be seen in Fig. 3 on p. 8.
Both MAPE and C-index indicate rather good prediction performance to the distance of 120 m from the nearest soil sample point. This is seen both with k-NN and RLS methods. When MAPE is higher than the baseline, it is better to use baseline average than the prediction. MAPE baseline is the horizontal line in the lower figures in Fig. 3.

<figure><img src="image_4.png" /><figcaption>Fig. 3: Left: k-NN results with $k = 6$ and 3 different feature sets. Right: RLS results on features-only case. C-index and Pearson correlation at top and MAPE below. The prediction performance is adequate below 120 meters.</figcaption></figure>

The dead zone radius $r > 0$ simulates a situation, where the test point is at least $r$ distance away from the given training points. $r = 0$ is traditional LOO test arrangement and measures best the properties of the predicted value within the training set itself. It may be too optimistic, since we seek for generalization. A large radius $r \approx \infty$ is overly pessimistic, since it would use only tiny fragments of the training set and would completely distort the prediction.

The prediction performance near $r = 0$ seems to indicate rather good generalization ability, but the performance reduces drastically over the dead zone distance $r$. Further study, both theoretical and practical, must be done to prop-
erly address classical geoinformatics concepts such as spatial autocorrelation and
spatial semivariance together with the general prediction ability. The problem
is new, since spatial analysis in geosciences is usually applied in the sense of
interpolation and extrapolation performance, and general prediction is usually
analyzed in the terms of Machine Learning performance.

The feature selection was not attempted. There were two reasons for this:

- the number of features (96) remained modest.

- the forest inventory features are unique to Finland. They can be largely substituted by various remote measurements [12], which would extend the application scope of the method to whole sub-polar area.

# 5 Conclusions and future work

The results indicate that the chosen five data sources (forest inventory, gamma-ray, air-borne electromagnetic, topographical data and peat bog mask) can be used to estimate the water permeability to a certain range from known measurements. This range seems to be c. 120-150 m. The best results come from the k-NN method based on the location of the sample points only. This method is naturally unavailable for general prediction.

There are several possible improvements. Since the mapping from water per-
meability to soil types is not unique, see [19], a special majority rule could be
used to select the dominant soil type from neighboring grid point predictions.
Such expert rules would require additional features like sophisticated geomor-
phological categories.

The Aerial Light Detection and Ranging (LiDAR) data can substitute most
of the topographical and forest inventory data features. This would extend the
scope of the prediction to any location at the arctic zone, where only aerial and
satellite measurements are economical. LiDAR has also potential for derived
features like geological morphology [20] and soil water budget modeling [10].

The final goal is to predict the water permeability, soil types, approximate
water budget and the load bearing capacity of the terrain in relation to the given
weather forecast, while the model is based on remote measures and online learn-
ing based on measurements from the harvester fleet. The potential applications
aim to wide-area routing and location planning. In this regard, even a modest
prediction power of features-only prediction could yield a cumulative effect on
route decisions.

References

1. Azzalini, A., Diggle, P.: Prediction of soil respiration rates from temperature, moisture and soil type. Journal of the Royal Statistical Society - Series C: Applied Statistics **43** (1994) 505–526

2. Sculla, P., Franklin, J., Chadwick, O.: The application of classification tree analysis to soil type prediction in a desert landscape. *Ecological Modelling* **181** (2005) 1–15
3. Schulte, R., Diamond, J., Finkele, K., Holden, N., Brereton, A.: Predicting the soil moisture conditions of irish grasslands. Irish Journal of Agricultural and Food Research **44** (2005) 95–110

4. Gao, H., Tang, Q., Shi, X., Zhu, C., Bohn, T.J., Su, F., Sheffield, J., Pan, M., Lettenmaier, D.P., Wood, E.F.: Water budget record from variable infiltration capacity (vic) model. Algorithm Theoretical Basis Document for Terrestrial Water Cycle Data Records (2010)

5. Chapuis, R.: Predicting the saturated hydraulic conductivity of soils: a review. Bulletin of Engineering Geology and the Environment **71** (2012) 401–434

6. Mahmood, H., Hoogmoed, W., van Henten E.J.: Proximal gamma-ray spectroscopy to predict soil properties. Sensors **13** (2013) 16263–16280

7. McKenzie, N., Ryan, P.: Spatial prediction of soil properties using environmental correlation. Geoderma **89** (1999) 67–94

8. Emadi, M., Baghernejad, M., Pakparvar, M., Kowsar, S.: An approach for land suitability evaluation using geostatistics, remote sensing, and geographic information system in arid and semiarid ecosystems. Environ Monit Assess **164** (2010) 501–11

9. J.P., W., P.R., H.: Soil moisture estimation using remote sensing. In: Proceedings of the 27th Hydrology and Water Resources Symposium, The Institute of Engineers Australia, Melbourne (2002)

10. Leutner, B., Müller, H., Wegmann, M., Beierkuhnlein, C.: Modelling biodiversity and forest structure using hyperspectral. In: 41st Annual Meeting of the Ecological Society of Germany. (2011)

11. Hyvönen, E., Pänttäjä, M., Sutinen, M.L., Sutinen, R.: Assessing site suitability for scots pine using airborne and terrestrial gamma-ray measurements in finnish lapland. Canadian Journal of Forest Research **33-5** (2003) 796–806(11)

12. Tomppo, E., Katila, M., Mäkisara, K., Peräsaari, J.: Multi-source national forest inventory methods and applications. Volume 18 of Managing Forest Ecosystems. Springer (2008)

13. Hautaniemi, H., Kurimo, M., Multala, J., H., L., Vironmäki, J.: The three in one aerogeophysical concept ofgtk in 2004. Geological Survey of Finland, Special Paper **39** (2005) 21–74

14. Schwanghart, W., Kuhn, N.: Topotoolbox: a set of matlab functions for topographic analysis. Environmental Modelling & Software **25** (2010) 770–781

15. Weldon, T.P., Higgins, W.E., Dunn, D.F.: Efficient gabor filter design for texture segmentation. Pattern Recognition (1996) 2005–2015 This could be removed, a more arcaic Gabor reference instead!

16. Ojala, T., Pietikäinen, M., Mäenpää, T.: Multiresolution gray-scale and rotation invariant texture. IEEE Transactions on Pattern Analysis and Machine Intel- **24** (2002) 971987

17. Kiss, R.: Determination of drainage network in digital elevation models, utilities and limitations. Journal of Hungarian Geomathematics **2** (2004) 16–29

18. Gönen, M., Heller, G.: Concordance probability and discriminatory power in pro- portional. Biometrika **92** (2005) 965–970

19. Pohjankukka, J., Nevalainen, P., Pahikkala, T., Hyvönen, E., Sutinen, R., Hänninen, P., Heikkonen, J.: Arctic soil hydraulic conductivity and soil type recog- nition based on aerial gamma-ray spectroscopy and topographical data. Proceed- ings of the 22nd International Conference on (ICPR 2014) (2014) to appear.

20. Sutinen, R., Hyvönen, E., Middleton, M., Ruskeeniemi, T.: Airborne lidar detection of postglacial faults and pulju moraine. Global and Planetary Change (2014) 24–32
ACKNOWLEDGEMENTS

This work is done as a part of ULJATH project, which is funded by the *Finnish Funding Agency for Technology and Innovation* (TEKES). ULJATH stands for *New Computational Methods For The Efficient Utilization of Public Data*.

The authors would like to thank the anonymous reviewers for valuable com-
ments and suggestions.