# Using spatial aggregation of soil multifunctionality maps to support uncertainty-aware planning decisions

Léa Courteille<fnref n="1" />, Philippe Lagacherie<fnref n="1" />, Nadia Boukhelifa<fnref n="2" />, Evelyne Lutton<fnref n="2" />, Léa Tardieu<fnref n="3" />,<fnref n="4" />

## Abstract

To ensure soil preservation, it is essential to incorporate the soil's ability to provide ecosystem services into the spatial planning process. For well-informed planning decisions, stakeholders need spatially explicit information on the state of the soils and the functions they fulfill, with sufficient spatial resolution and quantified uncertainty. It has been shown that Digital Soil Mapping (DSM) products can provide such information. However, in some cases, fine spatial resolution coupled with high levels of uncertainty may lead stakeholders to overlook the inherent uncertainties in the information. Spatial aggregation of DSM products opens up a promising avenue for obtaining maps that are more tailored to the users' scales of decision-making while facilitating uncertainty communication. In this perspective, we propose a new spatial aggregation approach relying on spatially-constrained agglomerative clustering (AC). The spatial aggregation approach is applied to a 25 m resolution soil potential multifunctionality index (SPMI) map developed for the coastal plain of the Occitanie Region. This DSM product was increasingly aggregated to obtain SPMI maps of different resolutions considering two distinct areal metrics : proportions of area above a given threshold of SPMI, and mean SPMI. Each map was evaluated through a set of indicators selected for their potential impact on user decision making: mean spatial resolution, overall predicted uncertainty, quantity of information and mean within-unit variability. The maps were compared with respect to these indicators to other maps obtained with alternative aggregation methods employed in DSM literature (maps aggregated according to some administrative units and QuadMaps). We show that all the tested aggregation methods produced a substantial decrease of the map uncertainty with moderate loss of spatial resolution. However, only AC preserved the fine spatial pattern of the initial DSM product while enabling fine tuning of the uncertainty displayed to end-users. We show that AC can simplify the identification of extensive regions characterized by low uncertainty without losing information regarding soil multifunctionality, thereby facilitating and enhancing the efficiency of planning decisions.

<footnote n="1">¹LISAH, INRAE, Institut Agro Montpellier, IRD, Univ Montpellier, Montpellier, France.</footnote>

<footnote n="2">²MIA-PS, INRAE, AgroParisTech, Université Paris-Saclay, 22 place de l'Agronomie, Palaiseau, France.</footnote>

<footnote n="3">³TETIS, INRAE, AgroParisTech, CIRAD, CNRS, Univ Montpellier, Montpellier, France.</footnote>

<footnote n="4">⁴ CIRED, École des Ponts ParisTech, AgroParisTech, université Paris-Saclay, Cirad, CNRS, EHESS, Nogent-sur-Marne, France.</footnote>
# 1. Introduction

Soils are a finite resource that ought to be conserved due to their significant contribution to the provision of ecosystem services and their primary importance for achieving many Sustainable Development Goals of the United Nations [Bouma et al. 2019]. To maintain sustainable soil use and protect the soils that provide the highest contributions to a set of targeted ecosystem services, it is necessary to acquire spatially explicit data on soil capability [Bouma et al. 2017] for each of the considered functions [Vrščaj et al. 2008]. This information must be of sufficient spatial resolution to be operationally useful at the local and regional levels, which are critical in achieving sustainability goals [Teixeira da Silva et al. 2018]. Moreover, this data should explicitly encode cartographic uncertainty allowing users to gauge the map's accuracy before usage, and thereby fostering informed decision-making and enhancing the evaluation of associated risks [Wadoux et al. 2023, Loiseau et al. 2020, Arrouays et al. 2014].

Digital Soil Mapping (DSM) [McBratney et al. 2003] has brought a major contribution to address these needs. Indeed, DSM methods allow predicting the basic soil properties that are involved in soil quality assessment, along with their associated uncertainty. These predictions are provided at the nodes of grids that span a designated area of interest. Furthermore, it has become a common practice to evaluate DSM products against independent validation datasets [Heuvelink 2014] to complete the local uncertainty assessment by statistical indicators of the overall uncertainty of the DSM products. Although such information on uncertainty is often required by the users of soil maps [Richer-de-Forges et al. 2019], there are many reasons why it remains challenging to understand and apply in practice [Voltz et al. 2020, Wadoux et al. 2021]. First, the multiplicity of metrics commonly used to quantify uncertainty (confidence interval, mean error, mean quadratic error...) might confuse map users [Loiseau et al. 2020]. Second, these metrics are not intuitive [Arrouays et al. 2020, Voltz et al. 2020] and are rarely linked to the decision-making process [Wadoux et al. 2021, Chen et al. 2022]. Third, map users may bring their own biases such as interpreting spatial information as more accurate when presented at a fine resolution compared to the same information delivered at a coarser resolution [Loiseau et al. 2020, Buzzelli 2020]. Lastly, it is unclear what visual variables (e.g. color, texture... see [Bertin 1983]) can effectively promote better understanding of uncertainty and facilitate informed decision-making [Pang 2001]. It is essential to develop methodologies not only for estimating uncertainty but also for enhancing its conceptual and visual comprehensibility. This is particularly crucial, as information that is inadequately represented tends to be disregarded by end-users [Richer-de-Forges et al. 2019, Jansen et al. 2022].

In DSM, some studies explored innovative representations of the uncertainty information focusing on map visuals [Loiseau et al. 2020]. Only a few studies [Vaysse et al. 2017, Padarian et al. 2023] have explored the use of spatial aggregation to manage uncertainty. In both of these studies, spatial aggregations of soil property
predictions were applied with the aim of decreasing the mapping uncertainty until an acceptable level that would not hamper users' decision making was found. The approaches proposed by Vaysse et al. 2017 and Padarian et al. 2023 differed with regard to (i) the aggregated units – administrative units versus super-pixels of a quadtree, (ii) the areal metrics of predictions – proportions of area above a given threshold versus means of soil property and (iii) the methods for propagating uncertainty to the areal units – analytical formula versus stochastic simulations – both accounting for the correlations of mapping errors across soil properties.

However, the downside of decreasing the uncertainty of soil predictions by spatial aggregation was a loss of information and spatial details of the final DSM products, which could also be detrimental for users' decisions. Such a loss was judged acceptable since land managers often require aggregated information over areas of interest that are much larger than the DSM grid cells [Vaysse et al. 2017]. The intentional loss of spatial details was also found useful to communicate uncertainty to users who tend to associate low spatial resolution and high uncertainty [Padarian et al. 2023]. Nevertheless, it is important to limit this loss of information by adequate spatial aggregation techniques that could preserve the initial spatial patterns of predictions.

In this perspective, we present in this paper a new spatial aggregation approach that uses spatially-constrained agglomerative clustering [Ward 1963] to generate different spatially aggregated representations of an initial DSM product. The algorithm iteratively groups neighboring pixels having the most similar predicted values of soil predictions, which allows to obtain the most homogeneous aggregated units with regard to soil predictions. The spatial aggregation approach was applied to a 25 m resolution soil potential multifunctionality index (SPMI) map developed for the coastal plain of the Occitanie Region [Angelini et al. 2023]. This DSM product was increasingly aggregated to obtain SPMI maps of different resolutions considering two distinct areal metrics : proportions of area above a given threshold of SPMI and mean SPMI. Each map was evaluated through a set of indicators selected for their impact on user decision making: mean spatial resolution, overall predicted uncertainty, quantity of information and mean within-unit variability. The maps were compared with respect to these indicators to other maps obtained with the aggregation methods mentioned above: maps aggregated according to administrative units [Vaysse et al. 2017] and QuadMaps [Padarian et al. 2023].

## 2. Study area and input data

### 2.1. Study area

The case study area is the community of municipalities of the "Hérault valley", an administrative division bringing together 28 municipalities. This kind of intermunicipal structure promotes joint urban and environmental planning between municipalities and have proved to be a major player in terms of territorial management [Cadiou 2013]. The "Hérault valley" community of municipalities covers 481 km² and is located 15 km west of Montpellier, in southwestern France. It is included in the physiographic region which was mapped in Angelini et al. 2023, defined as the "Mediterranean coast" of the Occitanie region. For the past decades, the Hérault
valley area has faced a high demographic pressure (25500 inhabitants in 1999, 41070 in 2020 [Dossier complet - Intercommunalité-Métropole de CC Vallée de l’Hérault (243400694) n.d.]), and the resulting urban sprawling is leading to a rapid loss of natural and agricultural soils (0.75% of the area have been artificialized between 2009 and 2022, while the national average is 0.50% [Cerema 2024]). This territory hosts a fair diversity of landscapes : scrublands, limestone plateaus, fertile plain of the Hérault valley planted with vineyards. As presented by Lagacherie et al. 2001, “The soil pattern of the region is extremely complex and mainly arises from geological variations. The substratum is constituted by heterogeneous Mesozoic and Cenozoic marine sediments from which have been formed either Lithic Leptosols, Calcaric Regosols or Calcaric Cambisols [Soil Map of the World at 1:5,000,000 1961] This substratum was partially overlaid during Plio-Pleistocene and Pleistocene era by successive alluvial deposits. These deposits have produced soils ranging from Calcaric to Chromic Luvisols with stone contents ranging from 20% to 70%. Calcaric Fluvisols are found over recent alluvial deposits which are located in the alluvial plains of Hérault river and of its tributaries. Recent colluvial phenomena within each main geological unit add to the soil heterogeneity”.

## 2.2. Input data

This study mainly builds upon the work of Angelini et al. 2023 who mapped a soil potential multi-functionality index for agriculture (Agri-SPMI) over the coastal plain of the Occitanie region. This Agri-SPMI had previously been developed by Rabot et al. 2022. It assesses the potential multifunctionality of a soil by evaluating whether it fulfills the four following functions, under three agricultural land uses: (i) provision of a physical habitat for plant growth, (ii) provision of chemical habitat for plant growth, (iii) retention and transfer of water and pollutants, and (iv) carbon storage. The ability of a soil to fulfill these functions was assessed using expert-defined thresholds on the values of the 14 intrinsic properties considered. The three chosen land uses were selected to be compatible with the realization of these 4 functions: annual crops, perennial crops and forest. By summing the number of functions fulfilled for the three land uses, we obtain a score for the potential multifunctionality ranging from 0 to 12.

In Angelini et al. 2023, the Agri-SPMI was mapped following a specific multi-variate DSM approach that accounted for the correlations of mapping errors across the soil properties involved in the Agri-SPMI computation. The maps at a 25 m spatial resolution for the 14 considered properties - CEC, organic carbon, silt, sand, pH, coarse fragments (all at soil depths 0-30 and 30-100), available water capacity (at soil depths 0-30, 30-60 and 60-100) and soil depth) - were predicted following a Digital Soil Mapping approach. A regression co-kriging model was developed in two steps: first the basic soil properties were individually mapped from legacy soil data and spatial soil covariates using a random forest algorithm. Next, the residuals were interpolated using cokriging and a linear model of co-regionalisation. The Agri-SPMI map was then computed from the resulting maps for the 14 considered soil properties. The mapping uncertainties of soil properties were propagated by calculating the soil quality index over 300 stochastic simulations of soil properties that account
169

170

171

172

for their spatial correlations, since they were derived from the linear models of co-regionalisation. The final maps of both Agri-SPMI and its associated uncertainties were obtained by computing respectively the pixelwise mean and standard deviation over the 300 simulations.

173

174

175

176

177

178

179

We obtained the data for our study area using the exact same methodology as Angelini et al. 2023, except that we first enriched their input dataset (which was composed of 1229 profiles over Coastal Plain of the Occitanie region) with additional legacy soil profile data: the augmented dataset contains 9135 profiles for the same spatial extent. To carry out this study, we used the resulting 300 simulated Agri-SPMI maps at a 25 m spatial resolution and the averaged Agri-SPMI map, which, thanks to the added data, comes with a considerably lower uncertainty.

180

# 3. Methods

181

## 3.1. Spatial aggregation

182

183

184

185

186

187

188

189

The initial Agri-SPMI is a raster map of the Agri-SPMI computed as the average of 300 simulated maps. The first step of this work consists in reducing the level of spatial detail of this map by aggregating pixels that have similar averaged Agri-SPMI while optimizing a homogeneity measure of the derived regions. This problem falls under the domain of regionalization, also known under a number of other terms including zonation and spatially constrained clustering. It has been an important analysis methodology for many application domains, such as climate zoning, ecoregion analysis, and urban analysis [Guo et al. 2011].

190

191

192

193

194

195

196

197

198

199

A range of regionalization algorithms are available. We opted for the Agglomera-tive Clustering (AC) procedure based on its performance notably when dealing with large datasets. This algorithm performs bottom-up hierarchical clustering under a contiguity constraint [Murtagh 1985]. It starts with each pixel forming a separate group. At every iteration the two most similar clusters (according to the considered similarity function) are merged together into a new cluster until a desired number of clusters is reached. The contiguity constraint implies that only clusters that are neighbors (here according to 4-connectedness) can be merged. The Agglomerative Clustering algorithm is not randomized: for a given image, the order in which the pixels are merged is always the same.

200

201

202

203

204

205

206

207

208

We first constructed the 4-connected adjacency matrix of our raster image. This matrix was then fed into the scikit-learn implementation of the Agglomerative Clus-tering algorithm, with the other parameters set as default. This default implemen-tation relies on Ward's criterion Ward 1963 for recalculating the distances between clusters at each step. According to this criterion, between-cluster distance is calcu-lated as the increase in the total within-cluster variance after merging two clusters. It is defined so that within-cluster variance is minimized and inter-cluster variance maximized. This criterion therefore minimizes the information loss associated with each grouping [Guo 2009].

209

210

A map with a given level of aggregation can be obtained by stopping the AC procedure as soon as a given number of clusters *K* is reached. The smaller the value
of *K*, the coarser the spatial division of the case study area (see **Figure 1**). Starting from the 764823 pixels of the initial raster, we computed 15 maps with different levels of aggregation by taking 15 values of *K* (1, 2, 5, 10, 20, 40, 100, 200, 500, 800, 1200, 1600, 6000, 20000, 100000).

<figure><img src="image_1.png" /><figcaption>Figure 1: Examples of aggregated maps with different levels of spatial detail. The smaller the *K* value, the coarser the map. The resulting regions are irregularly shaped.</figcaption></figure>

## 3.2. Choice of the tested areal metrics

The spatial aggregation induces a loss of information concerning the distribution of SPMI values that should be mitigated as far as possible by a thorough selection of a relevant areal metric. This metric should allow users to get what they need to know while not missing important features of the data. The most straightforward spatial aggregate metric that has been tested in the literature is the mean value. It is also one aggregate users are familiar with, since block-kriging has been broadly used in geostatistics to predict the mean values of soil properties for a set of areas of interest [Burgess et al. 1980, Padarian et al. 2023, Wadoux et al. 2023]. An alternative areal metric found in the literature [Vaysse et al. 2017] was to consider a criterion (e.g. being above a threshold value) for a property of interest, and represent the proportion of area within the region where this criterion is satisfied. These two types of representations were compared in the frame of this study.

As far as the proportion metric is concerned, a prerequisite is the definition of a threshold that is relevant to the user's interest. We considered an Agri-SPMI of 10 since it emphasizes the highly multifunctional soils that should be preserved from an irrational use (e.g., soils of high environmental value are environmentally too "good" to be sealed by extensive shopping centers [Vrščaj et al. 2008]). However, this threshold value is to be seen as an example: it could instead be set by experts or end-users.

## 3.3. Map production

For each of the 15 aggregation levels we produced two separate maps – one for each areal metric. These aggregated maps were derived from the 300 equiprobable
Agri-SPMI raster maps, according to the procedure described in **Figure 2**.

Each of these 300 simulated rasters displayed the predicted SPMI value at simulation $n \in \mathcal{N}$ for every pixel $p \in \mathcal{P}$, noted $v_n^p$, where $\mathcal{P}$ is the set of pixels represented on the map, and $\mathcal{N}$ the number of simulations (here equal to 300).

For every simulation $n$ and each aggregated zone $z \in \mathcal{Z}$, we computed the two distinct metrics:

* the mean SPMI for the zone, noted $\mu_n^z$, which was computed as follows:

$$ \mu_n^z = \frac{1}{|\mathcal{P}_z|} \sum_{p \in \mathcal{P}_z} v_p^n \quad (1) $$

where $\mathcal{P}_z$ is the set of pixels contained in the zone $z$.

* the proportion of highly multifunctional soils, noted $y_n^z$, which was computed as the ratio of the number of pixels with $v_n^p > 10$ over the total number of pixels contained in the spatial unit. Defining $\mathcal{H}_z$ as $\{p \in \mathcal{P}_z | v_p^n > 10\}$, we get: $y_n^z = \frac{|\mathcal{H}_z|}{|\mathcal{P}_z|}$. In order to compute this zonal statistic we had to introduce a thresholding step, where only pixels with $v_p^n > 10$ were kept on each of the simulated maps (see **Figure 2**).

<figure><img src="image_2.png" /><figcaption>Figure 2: Procedure for the spatial aggregation of Agri-SPMI maps. For the sake of readability only 3 simulated maps (instead of 300) are represented.</figcaption></figure>

We generated two sets of 300 aggregated maps – one for each zonal metric. We computed for every zone the average values of these metrics over the 300 simulations, which we noted $\mu_z$ and $y_z$.
$$ \mu_z = \frac{1}{\mathcal{N}} \sum_{n=1}^{\mathcal{N}} \mu_z^n \qquad (2) $$

$$ y_z = \frac{1}{\mathcal{N}} \sum_{n=1}^{\mathcal{N}} y_z^n \qquad (3) $$

## 3.4. Map evaluation indicators

To evaluate and compare the different aggregated map products, a set of quantitative indicators of map performances were computed.

### 3.4.1. Mean spatial resolution of polygons

In order to characterize the level of spatial detail of any given aggregated map, we first computed the mean area of polygons, which is the average size area of polygons on the ground (in squared meters) [Rossiter 2003, cited by Hengl 2006]. We then took the square root of these values, so that the resolution of our aggregated maps can be directly compared to the resolution of the initial raster map, expressed in meters.

It is defined as follows:

$$ r = \sqrt{\frac{1}{|Z|} \sum_{z \in Z} A_z} \qquad (4) $$

where $A_z$ is the area of zone z.

As expected, the mean spatial resolution consistently increases with higher levels of aggregation.

### 3.4.2. Mean uncertainty of prediction

The mean uncertainty of prediction was computed to provide an indicator of the expected accuracy of the maps as calculated by the regionalization algorithm. The uncertainty of prediction associated with the displayed $\mu_z$ and $y_z$ was estimated by computing the variance of these metrics across the 300 simulations (see the two maps at bottom right on Figure 2).

$$ Var(\mu_z) = \frac{1}{\mathcal{N}} \sum_{n=1}^{\mathcal{N}} (\mu_z^n - \mu_z)^2 \qquad (5) $$

$$ Var(y_z) = \frac{1}{\mathcal{N}} \sum_{n=1}^{\mathcal{N}} (y_z^n - y_z)^2 \qquad (6) $$

In general, this variance decreases as the size of the spatial support increases: since the mean is computed over a larger sample of pixels, it tends to vary less across simulations. This variance is therefore minimal in the extreme case where the support of prediction is the whole study area.

For each level of spatial aggregation, the mean uncertainty of prediction over the whole map was computed as the area-weighted sum of prediction uncertainties:

$$ V_{\mu} = \frac{1}{\mathcal{A}} \sum_{z \in Z} A_z Var(\mu_z) \qquad (7) $$
$$V_y = \frac{1}{\mathcal{A}} \sum_{z \in \mathcal{Z}} \mathcal{A}_z Var(y_z) \qquad (8)$$

### 3.4.3. Mean information content of the predictions

Spatial aggregation induces a loss of spatial detail over the study area, yet this level of detail can be seen as a determining factor of the map quality. We quantified how this information content varies with spatial aggregation using a distinct metric. To elaborate a formula for this metric, we took note of the fact that the more the data is aggregated, the more smoothed and similar the displayed SPMI values are, which decreases the information content of the map. In the extreme case where the SPMI is averaged over the whole study area the information content is minimal. On the contrary we assume that a map where the displayed SPMI values are highly variable is more reflective of the true spatial distribution of SPMI values and therefore has a high information content. The information content was therefore defined as the variance of the displayed values. For each zone, we computed how much its zonal metric $\mu_z$ or $y_z$ differed from the mean metric value for the whole map $\bar{\mu}$ or $\bar{y}$, defined as:

$$\bar{\mu} = \frac{1}{\mathcal{A}} \sum_{z \in \mathcal{Z}} \mathcal{A}_z \mu_z \qquad (9)$$

$$\bar{y} = \frac{1}{\mathcal{A}} \sum_{z \in \mathcal{Z}} \mathcal{A}_z y_z \qquad (10)$$

where $\mathcal{A}$ is the area of the whole mapped region.

The information content of each map was then computed as the area-weighted sum of the squared deviation from $\bar{\mu}$ or $\bar{y}$ of all $\mu_z$ and $y_z$.

### 3.4.4. Mean within-unit variances of predictions

As the degree of spatial aggregations increases, the generated areal units become increasingly large and heterogeneous with regard to the SPMI values. Consequently the mean values of SPMI selected as the areal metric become progressively less relevant for estimating the real SPMI values of each location within the unit. This is to be taken into account as a distinct error term relative to each displayed zonal mean SPMI value, which adds up to the above-presented uncertainty of prediction. Borrowing from Leenhardt et al. 1994, we quantified this uncertainty source by computing the within-polygon variance of the predicted SPMI.

From there we derived the area-weighted metric $\mathcal{H}_\mu$, which measures the total heterogeneity over the whole map.

As opposed to the “mean-value” zonal metric, the relevance of the “proportion of pixels with SPMI>10” metric to summarize the spatial unit information is insensitive to the degree of heterogeneity of this spatial unit.
## 3.5. Comparison to existing aggregation methods

In order to compare our aggregation method to the ones that were already introduced in the DSM literature [Padarian et al. 2023, Vaysse et al. 2017], we produced two additional sets of maps:

* 6 maps regionalized using the quadtree-based approach described in Padarian et al. 2023. Using their code, we computed 6 “QuadMaps” of different spatial resolutions. We then considered each formed super-pixel as a separate spatial cluster. It should be highlighted that this approach aggregates pixels into super-pixels, for prediction uncertainty to come below a target threshold: this differs from the optimizing criteria of the AC procedure that prioritizes the within-cluster homogeneity of SPMI values.

* 3 maps whose regions are administrative units: Three levels of administrative units were chosen: plot-level, land registry-level and municipality-level.

The resulting polygon maps served as a basis for carrying out the procedure described in Figure 2. The selected evaluation metrics were then computed for the 15 maps regionalized using Agglomerative Clustering, the 6 QuadMaps and the 3 maps of administrative units.

# 4. Results

## 4.1. Initial Digital Soil Mapping product

The initial Agri-SPMI information was represented through a set of two maps, as usually practiced in the DSM community: a map showing the predicted Agri-SPMI values and another showing the associated widths of their 90% confidence intervals (Figure 3). The latter metric allows end-users to directly perceive the extent of the variation in the model's predictions at a given point [Loiseau et al. 2020]. The first map depicts the spatial distribution of soil multifunctionality, which shows some fine- and coarse-scale patterns, notably related to landform elements (valleys, hills...). It is worth noting that the lowest Agri-SPMI predictions are associated with the largest errors, as already observed by Angelini et al. 2023.

## 4.2. Evaluation of the aggregation method

Figure 4 shows that the total prediction uncertainties exhibit the same decrease regardless of the chosen spatial aggregation units and the areal metrics. It is interesting to notice that most of the decrease occurs in the first steps of aggregation ($r < 500$ m) for the two considered spatial aggregations.

Concomitantly, as the initial map is regionalized into fewer and fewer zones, its information content also decreases. In the least aggregated maps, only the short-range variability is filtered out. Conversely, when the data is aggregated into few zones, only a small fraction of the spatial variability displayed in the initial map remains. It is worth noting that this information loss is significantly faster for QuadMaps and maps of administrative units than for spatial aggregations using AC.
<figure><img src="image_3.png" /><figcaption>Figure 3: Initial raster maps of SPMI values (left) and their associated prediction errors (right) for the case study area.</figcaption></figure>

As foreseen, the uncertainty stemming from zonal heterogeneity, consistently rises with spatial aggregation. We observe that this increase is significantly slower when the data is aggregated using AC. The heterogeneity of Administrative Units maps rises just as fast as the heterogeneity of QuadMaps.

The examples of SPMI choropleth maps of similar spatial resolution (r ≈ 1,500m) obtained by the three considered spatial aggregation techniques (Figure 5) revealed very different spatial patterns both for SPMI values (upper line) and the prediction uncertainties (bottom line). Aggregated maps obtained from AC are by far those that less deform the genuine spatial pattern of predicted SPMI and uncertainties (Figure 3). The formed regions are very irregularly shaped and of variable sizes. The spatial resolution is unevenly distributed across the mapped area: some local structures are captured by small spatial clusters and delivered with a high level of detail, whereas much larger homogeneous structures (plains, plateaus...) tend to be described by a single large cluster. The aggregated maps using administrative units or quadtree reveal a great impact of the geometries of the a priori chosen aggregation units and many losses of interesting variation details shown in the genuine map. It is also interesting to notice that the uncertainty map generates aggregated units with more contrasted uncertainty values than the two others.

# 5. Discussion

## 5.1. Decreasing map prediction uncertainties

The top graphs of Figure 4 show the evolution of the total map prediction uncertainty when spatial resolution becomes coarser. The fact that spatial aggregation leads to a decrease in prediction uncertainty was expected in light of previous experiments (Burgess et al. 1980, Vaysse et al. 2017). We note that in this study, the decrease in prediction uncertainty is very fast at early aggregation stages. This is linked to the spatial structure of the prediction errors computed through the 300 stochastic simulations. This structure may be explained by the fact that most soil
<figure><img src="image_4.png" /><figcaption>Figure 4: Evolution of the three chosen map evaluation metrics with spatial aggregation for maps aggregated using Agglomerative Clustering, Administrative Units or Quadtree algorithm, for the two considered areal metrics: mean SPMI (left column) and proportion of highly multifunctional soils (right column).</figcaption></figure>

properties involved in the SPMI assessment show high short-range spatial variability that cannot be captured by the explanatory environmental covariates. This constitutes a favorable setting for our study since even slightly aggregated maps with a fine spatial resolution display a considerably lower prediction error than the initial raster map. This present case is representative of many areas in the world, as evidenced by the compilation of variograms gathered by McBratney et al. 1999. Decreasing prediction error may help decision-makers since the uncertainty values associated with predictions at pixel support may be so high that the displayed information is hard to interpret, and as a result tends to be ignored [Richer-de-Forges et al. 2019].

380

381

382

383

384

385

386

387

388
<figure><img src="image_5.png" /><figcaption>Figure 5: Examples of aggregated maps with similar spatial resolutions (r ≈ 1,500m) showing the zonal proportions of highly multifunctional soils, obtained with the 3 zonation methods considered. The AC map was obtained by specifying K=200, the Administrative Units map was obtained by aggregating SPMI information at the land registry level and the QuadMap was obtained by specifying an uncertainty threshold of 0.01.</figcaption></figure>

## 5.2. Selecting the aggregation technique

The value of the considered evaluation metrics and the visual examination of the map products converge to indicate the major advantage that AC has over other methods to optimize spatial aggregation with regard to the underlying soil variability. To the best of our knowledge, it is the first time that this method has been employed in DSM. The two other techniques start from spatial entities that were selected a priori without explicit regard to the underlying soil variability. Furthermore, the quadtree technique did not aim to generate homogeneous clusters but to aggregate neighbouring pixels until they reach an user-defined level of uncertainty. They yield maps with less information content that distort the initial soil variation patterns. Admittedly, choosing units of interest (e.g. plots, municipalities) as starting points allows the soil map to fit the map of land management units – as recommended in Vaysse et al. 2017. However, these spatial units might be highly heterogeneous when it comes to soil, and the resulting choropleth maps can be affected by the Modifiable Areal Unit Problem (MAUP) [Openshaw 1984]. This issue refers to the sensitivity of analytical results (statistical hypothesis tests, regression coefficients...) to the arbitrary choice of units at which data is aggregated [Röthlisberger et al. 2017].

The AC procedure is also more flexible: the spatial clusters can be computed
from one or several soil descriptors (soil property, soil function fulfillment...), and
the resulting boundaries will coincide with the underlying patterns of the data. This,
coupled with the fact that our partition of data minimizes intra-cluster variance,
makes our maps robust to the MAUP [Fotheringham et al. 1991].
407
408
409
410

Decisions about soils typically relate to spatial objects that are much larger than
a pixel (plots, water catchments...) and often require the delineation of areas with
sharp boundaries (e.g. zoning decisions, taxation decisions). Spatial aggregation of
DSM products helps relieve decision-makers from the cognitive or technical burden
of aggregating pixel values and may simplify the identification of distinctive areas.
A major advantage of AC is that it guarantees that the resulting regions reflect
the underlying soil spatial patterns. This helps identify homogeneous regions for
conservation or management activities, while still allowing for carrying out spatial
quantitative analyses. It should be noted however that the resulting maps could
exhibit very variable levels of uncertainty across clusters, which differed from the
results obtained by Padarian et al. 2023 that optimized the uncertainty homogeneity
within the map.
411
412
413
414
415
416
417
418
419
420
421
422

### 5.3. Selecting the aggregation metric

We tried out two aggregation metrics introduced in the literature. Their uncer-
tainty decreases in a similar way as spatial aggregation is carried out. In the case
where we summarize pixel values by taking their mean, this gain in prediction error
is canceled out by the rise of an error linked to the heterogeneity of the aggregated
zones [Leenhardt et al. 1994], which has a much larger magnitude. Furthermore, it
has been shown that the average value tends to be a poor describer of individual
observations and hence misleading [Wong 2004] notably in the case where distribu-
tions of pixel values within regions are skewed or contain outliers [Andrienko et al.
2006]. Summarizing the initial spatial information using metrics that are insensitive
to such issues and still deliver valuable information to end-users might therefore be
the best approach. Some commonly-used metrics could have been computed, like the
intra-zonal median or quantiles of the SPMI, but we chose to represent the proportion
of soils satisfying a given criterion (SPMI > 10). Such a metric has the advantage of
being easily interpretable for end-users. Although our heavily aggregated map prod-
ucts still show some sharp contrast (see **Figure 5**), we can expect that generally for
advanced aggregation stages, the displayed proportion value will converge towards an
averaged proportion which does not carry much information.
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440

## 5.4. Tuning spatial resolution for decision making

The great flexibility offered by AC provides numerous representations of a single
SPMI spatial distribution, among which the best one is hard to identify. The fact that
the information content of aggregated maps decreases more slowly than prediction
errors suggests that there could be an optimal representation with an almost minimal
error and a still high information content. But the decision-maker's preferences and
the nature of the decision to be made should play a major role in the selection of
the most suitable spatial resolution. Moreover, the computed map quality indicators
only inform about the overall picture. The AC map of uncertainty in **Figure 5**
441
442
443
444
445
446
447
448
449
450 clearly shows that these indicators do not reflect the uneven distribution of prediction
451 uncertainty across the study area. In some circumstances, decision-makers might
452 notice that they can have it all: their management area may be described with fine
453 spatial resolution and low uncertainty. This points towards the need to present the
454 aggregated maps in a dynamic interface, allowing the user to interactively explore the
455 trade-off between resolution and prediction error: within some regions, increasing the
456 level of spatial detail amounts to having to deal with highly uncertain information,
457 which makes decision-making a much more challenging and risky task. The user could
458 therefore choose the most adequate map product with respect to both their needed
459 level of detail and their tolerance to uncertainty.

460
5.5. Selecting the optimization criterion for spatial aggregation

The initial selection and the testing of our aggregation algorithm (AC) were based
461 on the underlying hypothesis that the users would prefer aggregated maps that best
462 represent the known soil variability. However, Padarian et al. 2023 built their ap-
463 proach from a very different hypothesis, stating that the users would relate the un-
464 certainty level with the map resolution. This motivated the aim to produce QuadMaps
465 with limited and homogeneous uncertainties. Although the AC procedure provided
466 a similar decrease of the overall uncertainty (Figure 4, top graphs), it did not fully
467 fulfill this alternate objective since the resulting maps exhibited very variable levels
468 of uncertainty across clusters (Figure 5, bottom left). Although we believe that pre-
469 serving the visualization of the soil variability is more important than obtaining maps
470 with homogeneous levels of uncertainty, the final selection among these two alternate
471 targets belong to the map users.
472

6. Conclusion and perspectives

We successfully applied the Agglomerative Clustering procedure to aggregate a
474 soil multifunctionality raster map into a choropleth map with irregularly-shaped zones
475 that reflect the underlying spatial distribution of the data. By varying the number
476 of clusters, the information can be rendered at any desired level of detail. The aggre-
477 gation metric we chose to summarize the soil information within the zones was their
478 proportion of highly multifunctional soils. We computed the uncertainties associated
479 with each metric value and displayed them on a separate choropleth map. We hope
480 that our method will pave the way for alternative visualizations of soil information
481 that are more tailored to the needs of decision-makers than conventional display of
482 DSM products. Given the widespread need for data abstraction and uncertainty vi-
483 sualization for better decision-making, it could be useful in similar fields (ecosystem
484 services modeling, species distribution modeling...).
485

Finally, users must be allowed to choose the level of detail which is the most
486 relevant to their decision tasks, ideally through an interactive interface. But making
487 the map fully adaptable to the needs of different actors and stakeholders requires an
488 in-depth examination of their preferences and information needs when they have to
489 make decisions under uncertainty.
490
7. Acknowledgements

The authors would like to thank Dr Cássio Fraga Dantas and Dr Dino Ienco for the insightful discussions on the mechanics of the Agglomerative Clustering algorithm, and Dr Gerard Heuvelink for his useful comments at the early stages of this work.

8. Conflict of interest statement

The authors declare no conflicts of interest.

9. Data availability statement

The data that support the findings of this study are available from the corresponding author upon request.

References

Andrienko, N. and G. Andrienko (2006). *Exploratory analysis of spatial and temporal data: a systematic approach*. Springer Science & Business Media. DOI: 10.1007/3-540-31190-4.

Angelini, M.E., G.B.M. Heuvelink, and P. Lagacherie (2023). "A multivariate approach for mapping a soil quality index and its uncertainty in southern France". In: *European Journal of Soil Science* 74.2, e13345. DOI: 10.1111/ejss.13345.

Arrouays, D., L. Poggio, Osvaldo A Salazar Guerrero, and V.L. Mulder (2020). "Digital soil mapping and GlobalSoilMap. Main advances and ways forward". In: *Geo-derma Regional* 21, e00265.

Arrouays, D. et al. (2014). "GlobalSoilMap: Toward a fine-resolution global grid of soil properties". In: *Advances in agronomy* 125, pp. 93-134. DOI: 10.1016/B978-0-12-800137-0.00003-0.

Bertin, Jacques (1983). *Semiology of graphics*. University of Wisconsin press. DOI: 10.5555/1095597.

Bouma, J., L. Montanarella, and G. Evanylo (2019). "The challenge for the soil science community to contribute to the implementation of the UN Sustainable Development Goals". In: *Soil Use and Management* 35.4, pp. 538-546. DOI: 10.1111/sum.12518.

Bouma, J., M.K. Van Ittersum, J.J. Stoorvogel, N.H. Batjes, P. Droogers, and M.M. Pulleman (2017). "Soil capability: exploring the functional potentials of soils". In: *Global soil security*, pp. 27-44. DOI: 10.1007/978-3-319-43394-3_3.

Burgess, T.M. and R. Webster (1980). "Optimal interpolation and isarithmic mapping of soil properties: II. Block kriging". In: *Journal of Soil Science* 31.2, pp. 333-341.

Buzzelli, M. (2020). "Modifiable areal unit problem". In: *International encyclopedia of human geography*, p. 169. DOI: 10.1016/B978-0-08-102295-5.10406-8.

Cadiou, S. (2013). "Intercommunality, or the broken promises of local democracy". Trans. by O. Waine. In: *Métropolitiques*. URL: http://www.metropolitiques.eu/Intercommunality-or-the-broken.html.
Chen, S. et al. (2022). "Digital mapping of GlobalSoilMap soil properties at a broad scale: A review". In: *Geoderma* 409, p. 115567. DOI: 10.1016/j.geoderma.2021.115567.

Soil Map of the World at 1:5,000,000 (1961). URL: https://www.fao.org/soils-portal/data-hub/soil-maps-and-databases/faounesco-soil-map-of-the-world/en/.

Fotheringham, A.S. and D.W.S. Wong (1991). "The modifiable areal unit problem in multivariate statistical analysis". In: *Environment and planning A* 23.7, pp. 1025-1044. DOI: 10.1068/a231025.

Guo, D. (2009). "Greedy optimization for contiguity-constrained hierarchical clustering". In: *2009 IEEE International Conference on Data Mining Workshops*. IEEE, pp. 591-596. DOI: 10.1109/ICDMW.2009.75.

Guo, D. and H. Wang (2011). "Automatic region building for spatial analysis". In: *Transactions in GIS* 15, pp. 29-45. DOI: 10.1111/j.1467-9671.2011.01269.x.

Hengl, T. (2006). "Finding the right pixel size". In: *Computers & geosciences* 32.9, pp. 1283-1298. DOI: 10.1016/j.cageo.2005.11.008.

Heuvelink, G.B.M. (2014). "Uncertainty quantification of GlobalSoilMap products". In: *GlobalSoilMap; CRC Press: Boca Raton, FL, USA*, pp. 335-340.

Dossier complet - Intercommunalité-Métropole de CC Vallée de l'Hérault (243400694) (n.d.). https://www.insee.fr/fr/statistiques/2011101?geo=EPCI-243400694. Accessed: 2023-08-02.

Jansen, Jan, Skipton NC Woolley, Piers K Dunstan, Scott D Foster, Nicole A Hill, Marcus Haward, and Craig R Johnson (2022). "Stop ignoring map uncertainty in biodiversity science and conservation policy". In: *Nature Ecology & Evolution* 6.7, pp. 828-829. DOI: 10.1038/s41559-022-01778-z.

Lagacherie, P., J.-M. Robbez-Masson, N. Nguyen-The, and J.P. Barthès (2001). "Mapping of reference area representativity using a mathematical soilscape distance". In: *Geoderma* 101.3-4, pp. 105-118. DOI: 10.1016/S0016-7061(00)00101-4.

Leenhardt, D., M. Voltz, M. Bornand, and R. Webster (1994). "Evaluating soil maps for prediction of soil water properties". In: *European Journal of Soil Science* 45.3, pp. 293-301. DOI: 10.1111/j.1365-2389.1994.tb00512.x.

Loiseau, T., A.C. Richer-de-Forges, P. Roudier, C. Ducommun, S. Chen, P. Lagacherie, and D. Arrouays (2020). "Essais de représentation cartographique de l'incertitude pour les utilisateurs de cartographies des sols par modélisation statistique". In: *Étude et Gestion des Sols* 27.1, pp. 257-275.

McBratney, A.B. and M.J. Pringle (1999). "Estimating average and proportional variograms of soil properties and their potential use in precision agriculture". In: *Precision Agriculture* 1, pp. 125-152. DOI: 10.1023/A:1009995404447.

McBratney, A.B., M.L.M. Santos, and B. Minasny (2003). "On digital soil mapping". In: *Geoderma* 117.1-2, pp. 3-52. DOI: 10.1016/S0016-7061(03)00223-4.

Murtagh, F. (1985). "A survey of algorithms for contiguity-constrained clustering and related problems". In: *The computer journal* 28.1, pp. 82-88. DOI: 10.1093/comjnl/28.1.82.

Openshaw, S. (1984). "The modifiable areal unit problem". In: *Concepts and techniques in modern geography*.
Padarian, J. and A.B. McBratney (2023). "QuadMap: Variable resolution maps to better represent spatial uncertainty". In: *Computers & Geosciences* 181, p. 105480.
DOI: 10.1016/j.cageo.2023.105480.

Pang, A.T. (2001). "Visualizing uncertainty in geospatial data". In: *Proceedings of the workshop on the intersections between geospatial information and information technology*. Vol. 10. 1.20, p. 3823.

Cerema (2024). *Indicateurs consommation d'espaces produites dans le cadre du Pnb Action7 à partir des Fichiers Fonciers*. Accessed: 2024-04-05. URL: %5Curl%7Bhttps://cerema.app.box.com/v/pnb-action7-indicateurs-ff%7D.

Rabot, E., M. Guiresse, Y. Pittatore, M.E. Angelini, C. Keller, and P. Lagacherie (2022). "Development and spatialization of a soil potential multifunctionality index for agriculture (Agri-SPMI) at the regional scale. Case study in the Occitanie region (France)". In: *Soil Security* 6, p. 100034. DOI: 10.1016/j.soisec.2022.100034.

Richer-de-Forges, A., D. Arrouays, M. Bardy, A. Bispo, P. Lagacherie, B. Laroche, B. Lemercier, J. Sauter, and M. Voltz (2019). "Mapping of soils and land-related environmental attributes in France: analysis of end-users' needs". In: *Sustainability* 11.10, p. 2940. DOI: 10.3390/su11102940.

Rossiter, D.G. (2003). *Methodology for Soil Resource Inventories*.

Röthlisberger, V., A.P. Zischg, and M. Keiler (2017). "Identifying spatial clusters of flood exposure to support decision making in risk management". In: *Science of the total environment* 598, pp. 593-603. DOI: 10.1016/j.scitotenv.2017.03.216.

Teixeira da Silva, R., L. Fleskens, H. van Delden, and M. van der Ploeg (2018). "Incorporating soil ecosystem services into urban planning: status, challenges and opportunities". In: *Landscape Ecology* 33, pp. 1087-1102. DOI: 10.1007/s10980-018-0652-x.

Vaysse, K., G.B.M. Heuvelink, and P. Lagacherie (2017). "Spatial aggregation of soil property predictions in support of local land management". In: *Soil Use and Management* 33.2, pp. 299-310. DOI: 10.1111/sum.12350.

Voltz, M., D. Arrouays, A. Bispo, P. Lagacherie, B. Laroche, B. Lemercier, A.C. Richer-de-Forges, J. Sauter, and N. Schnebelen (2020). "Possible futures of soil-mapping in France". In: *Geoderma Regional* 23, e00334. DOI: 10.1016/j.geodrs.2020.e00334.

Vrščaj, B., L. Poggio, and F.A. Marsan (2008). "A method for soil environmental quality evaluation for management and planning in urban areas". In: *Landscape and urban planning* 88.2-4, pp. 81-94. DOI: 10.1016/j.landurbplan.2008.08.005.

Wadoux, A. M.J.-C. and G.B.M. Heuvelink (2023). "Uncertainty of spatial averages and totals of natural resource maps". In: *Methods in Ecology and Evolution* 14, pp. 1320-1332. DOI: 10.1111/2041-210X.14106.

Wadoux, A. M.J.-C., G.B.M. Heuvelink, R.M. Lark, P. Lagacherie, J. Bouma, V.L. Mulder, Z. Libohova, L. Yang, and A.B. McBratney (2021). "Ten challenges for the future of pedometrics". In: *Geoderma* 401, p. 115155. DOI: 10.1016/j.geoderma.2021.115155.
618 Ward, J.H. Jr (1963). "Hierarchical Grouping to Optimize an Objective Function".

619 In: *Journal of the American Statistical Association* 58.301, pp. 236–244. DOI: 10.

620 1080/01621459.1963.10500845.

621 Wong, D.W.S. (2004). “The modifiable areal unit problem (MAUP)”. In: *World-

622 *Minds: geographical perspectives on 100 problems: commemorating the 100th an-
623 anniversary of the association of American geographers 1904–2004*. Springer, pp. 571–

624 575. DOI: 10.1007/978-1-4020-2352-1_93.