Ecosystem services (ES) have progressively become an important concept in environmental planning and policy-making to bridge the science – policy interface in the management of ecosystems (Braat and de Groot 2012, Perrings et al. 2011, Groot et al. 2010). Mapping indicators of ES is one prominent approach to improve spatially explicit decision-making and land management (Burkhard et al. 2012). This ecosystem service approach can help policy-makers to target strategic areas, formulate new policies and/or evaluate impacts of previous policies (Burkhard et al. 2013, Maes et al. 2012, Martinez-Harms et al. 2015). ES maps are popular outreach and data visualisation products. They have profited from the increased availability and applications of tools such as GIS or remote sensing that helps to increase their production and distribution (Palsky 2013). Despite its popularity, mapping ES has its limitations (Eigenbrod et al. 2010van Oudenhoven et al. 2018). Specifically, spatial information used in mapping is rarely, if ever, completely accurate or verified (Heuvelink and Burrough 2002; Devendran and Lakshmanan 2014) and that spatial assessments are often based on coarse data, which decreases the confidence in the spatial products (Andrew et al. 2015; Hamel and Bryant 2017).

The importance of uncertainty in the use and analysis of spatial data has long been recognised in land-use and landcover change mapping (e.g. Dendoncker et al. 2008, Lavorel et al. 2017, Verburg et al. 2013, Verburg et al. 2011, Kuemmerle et al. 2013) and is increasingly acknowledged in the ES mapping community (Jacobs et al. 2017). The representation of uncertainties, related to missing or incomplete data, aggregation error, functional knowledge gaps or normative and value-laden indicators in these types of maps is key for the map user (Crossman et al. 2013, Jacobs et al. 2013). In principle, information on uncertainties should allow a critical look at the map and be an integral part of the decision-making process for which maps are used. As a matter of fact, there is an increasing number of scientific studies using ES maps addressing uncertainty or reliability (e.g. Bagstad et al. 2018, Brunner et al. 2017, Johnson et al. 2012, Boithias et al. 2016, Grêt-Regamey et al. 2013, Schulp et al. 2014). However, uncertainty in the valuation process still remains an exception and has not explicitly been accounted for in many of the studies.

Uncertainty is complex and there are many definitions or typologies of uncertainty in ES analyses (Hamel and Bryant 2017) and therefore there are many ways to address the issue. Most of the studies focus on the uncertainty stemming from the ecosystem complexity or on the social-ecological uncertainties from supply and demand in ES models. More recently, studies have highlighted the effects of the input data on the reliabilility of ES assessments (Kangas et al. 2018). There is, however, much less emphasis on the epistemic uncertainty in ES mapping due to the modelling process. In the field of ES mapping, there is specifically a lack in the representation of errors inherent to the extrapolation of point-based measurements to produce empirical-based ES maps. Mapping an ES indicator from point field measurements implies several causes of epistemic uncertainties, inherent in each step from measurement to extrapolation. For instance, because of the mismatch between the field sample, its GPS location information and pixel size, there are errors generated, as field and remote sensing data are linked (Réjou-Méchain et al. 2015) and there is also uncertainty generated by the application of statistical methods, i.e. related to the extrapolation process itself. A comprehensive analysis and transparent documentation of the latter are seldom provided in ES mapping approaches.

The objective of this research is to describe a simple approach to assess and spatially represent uncertainties associated with the extrapolation of measured field data using regression analysis to map an ES indicator (vegetation carbon stock). Our results express the degree of certainty that we have about our ES map and the confidence that policy-makers can have in the map. To this end, we firstly used field data measurements for vegetation carbon stock in the Brazilian Amazon forest as a response variable and various land characteristics derived from remote sensing as potential explanatory factors of the ES indicator. We implemented the mapping techniques from earlier works (Le Clec’h et al. 2017) and evaluated the variability of ES supply in space resulting from our extrapolation process. Secondly, we calculated and mapped an indicator based on the prediction interval as a measure of the uncertainty stemming from the regression. Thirdly, we represented the variation of the prediction interval within the land-use classes and within the three sites, to identify sources of uncertainty within our dataset. While such estimation and analysis of the uncertainty are common in some fields, such as geostatistics, it is still under-considered in the field of ES mapping. Therefore, in the field of ES assessment and mapping, this study can be seen as a first step in a systematic assessment of uncertainties related to ES mapping in which point-based measurements are used to construct ES maps. Systematic assessments of uncertainties should critically question the confidence and usability of ES maps. This should be critical for discussing and determining the potential use(s) of the map and/or for supporting policy and planning processes. As some studies recommend the application of the principle of parsimony in ES mapping (e.g. Jacobs et al. 2017), we favoured a simple approach in this paper.

Site description

We estimated and mapped uncertainty related to ES maps using the case of the Brazilian Amazon forest, a global hotspot of ES provision (Fearnside 2017). More specifically, we focused on the State of Pará that spans 1.25 million km², which is partly covered by tropical forest. Extensive forests, high biodiversity and rapid deforestation qualify Pará as indicative of contexts in which ES surveys are urgently required. The rapid deforestation observed in Pará is mainly caused by timber harvesting and subsequent cattle grazing, which prevents forest restoration (Fearnside 2017, Godar et al. 2012), with 1,887 km² deforested in 2014 (INPE, Instituto Nacional de Pesquisas Espaciais 2014). While deforestation has been slowing, the area deforested in Pará in 2013 still accounted for almost half of the total deforestation in Brazil.

Within Pará, we conducted field measures of ES within three local study sites that are representative of the regional variability in socio-economics, deforestation temporalities and ES change (Fig. 1, Oszwald et al. 2011). The first site, Maçaranduba belongs to the municipality of Nova Ipixuna. It covers 220 km² and has been deforested since the early 1970s. Although largely deforested, its remaining forest is relatively well preserved. The second site, named here Pacajá, belongs to the municipality of Pacajà. It covers 175 km² and is located on one of the ‘fishbone’ roads (Traverssão 338 South) of the Trans-Amazon highway. It has been deforested since the 1990s by agricultural colonisation that has resulted in the conversion of the landscape to pasture and annual and perennial crops. Despite the deforestation dynamics, forests still covered 63% of this site in 2007. The third site, Palmares II, is an assentamento (settlement) of 160 km² near the Carajas iron mine. Emblematic of the agrarian reform, this site was established following the conflict between Landless Movement workers, a fazendeiro, owner of a large farm, which was located on the site, the federal government and the Vale mining company. This conflict led to the fragmentation of the land and thus of the forest that only covered 23% of the site in 2007.

In a deforestation front context, because of the lack of law, public policies and environmental management, the study of ES is highly relevant. It allows us a better understanding of the impacts of deforestation and cultivation on the environment, in other words to evaluate the impacts of the past and current public policies (or their absence) on the ecosystems.

Data

We used two different datasets to apply predictive statistical methodologies for the ES indicator: field data related to the ES indicator (response variable) and remote sensing data (explanatory variables – Table 1).

Field data: carbon stocks

According to the definition of ES proposed by Portela and Rademacher (2001), a way to map ES consists of defining spatial indicators of the biophysical processes that provide the services (Oszwald et al. 2014). We studied here an indicator of the service of climate regulation. The chosen indicator is the vegetation carbon stock whose analysis is relevant in a pioneer front context. Even though stocks are not per se an ES, we assumed this indicator to be relevant in this context, as their main driver of change is due to human activities. It was measured for 135 sampling points that were assessed during fieldworks in 2008 before the dry season. In each of the 27 farms (nine per site), five sampling points were equally spaced (around every 200 m) along a transect corresponding to the longest diagonal of the farm or to a north-south axis.

The carbon stock was estimated using a factor of 0.5 (Markewitz et al. 2004) in aboveground dry total biomass calculated from field data shown by Costa et al. (2012). At each sampling point, the dry biomass was calculated for four strata: upper, middle, lower (same strata where the vegetation cover was inventoried) and necromass. A plot of 10 m× 50 m (500 m²) was established for the upper stratum inventory (individuals with DBH - diameter at breast height ≥ 10 cm), a subplot of 5 m× 50 m was established for the middle stratum (individuals with DBH < 10 cm and height ≥ 2.0 m) and finally 4 subplots of 1 m ×1 m were regularly distributed in the centre of the plot for the lower and necromass strata measured (individuals with 2.0 m > height ≥ 10 cm). In the upper and middle strata, dry biomass was estimated from allometric equations using diameter at breast height. In the primary forests, biomass was estimated according to Higuchi et al. (1998) using three classes of DBH: DBH > 20 cm, 5 cm < DBH < 20 cm and DBH < 5 cm. In secondary forest, biomass was calculated for Cecropia and all other trees, regardless of diameter, according to Nelson et al. (1999). In both forests, primary and secondary, biomass of lianas was estimated according to Gerwing (2002). The same equations were used in all plots to ensure a better comparison between them. Due to the low density of trees found in pastures and cultivation areas, the dry biomass of the upper and middle strata were calculated using the same formulas used in primary and secondary forests. In the lower stratum, the biomass was first calculated directly by fresh weight; thereafter, one sample was taken and also weighed; and finally, the sample was oven-dried at 70°C until a constant dry weight and the dry biomass of the sub-plot was calculated by cross-multiplication using the dry biomass of the sample. After the removal of plants from sub-plots, all remaining ground material (leaves, twigs, flowers, dead wood) was gathered to measure the necromass in a similar procedure to that used for measuring the biomass of the lower strata.

Remote sensing data

We built and applied linear regression to extrapolate and map the ES indicator, using the plot-level measures and local high-resolution satellite imagery (Table 1). Since the vegetation carbon stock is strongly related to the biophysical characteristics of the landscape (Grimaldi et al. 2014), explanatory variables were chosen to characterise different aspects of the landscape. We used remote sensing data that provides information about vegetation cover and topography. They are derived from the processing of the 1986, 1996, 2001 and 2007 Landsat TM images (30 m x 30 m spatial resolution) carried out under ENVI and the processing of the Aster Digital Elevation Model (DEM - 30 m x 30 m spatial resolution). These data are known for the three study sites and they were extracted for the 135 sampling points.

We used remote-sensing data to characterise the land-cover of the three study sites for 2007. A Landsat TM image from the dry season (30 m spatial resolution) was used to build a supervised classification by maximum likelihood, to calculate two vegetation indexes (NDVI and NDWI) and to determine a historical trajectory of land-cover (Oszwald et al. 2012). The Landsat classifications characterise six land-cover classes (forest, burned forest, juquira-capoeira [fallow lands], grasslands with trees, clean grasslands and bare soils) for 2007. Training data, used for the supervised classification, were sampled during field campaigns using a GPS. Landsat images were radiometrically and geometrically corrected prior to classification to ensure comparability between the study sites. In addition, analysis and processing of land cover classification were also extended to all Landsat TM images of the dataset (from 1986 to 2007). Using the land-use maps obtained for the four dates, five classes of land-use trajectories were determined (Oszwald et al. 2012), ranging from a homogeneous forest structure (class 1) to an agricultural dynamic of extensive breeding (class 5). Landsat TM images were also used to calculate two vegetation indices giving information about vegetation density (NDVI) and water content into plants (NDWI).

Data also provided information about the elevation (in metres) at every point. Slopes synthesised the altitude difference between two adjacent pixels and are provided as a percentage. These two variables (elevation and slope) are quantitative and are treated as continuous raw data. The "topography" variable corresponds to a synthetic characterisation of the topographic context comprising four modalities: bottom of valleys, hilltops, zones of steep slopes and zones of low slopes. Finally, we deduced the hydrographic network from the DEM which was used to determine a distance to the rivers (0 to 100 m, 100 to 200 m, 200 to 300 m, 300 to 500 m and more than 500 m).

We also used a variable, named “Site” that corresponds to the identity of the study site to which each pixel belongs. This variable aims to (1) estimate the spatial auto-correlation of the sampling points and (2) take into account the inter-sites variability, due to homogeneous biophysical conditions and socioeconomic characteristics within each location.

Statistical approach

We aimed to evaluate and map the confidence we have in our prediction, as a way to represent the uncertainty related to the ES map. We based our approach on the implementation of statistical methodologies that link field (ES indicator) and remote sensing (explanatory) data. These statistical methodologies are used to (1) extrapolate field data using remote sensing to the three study sites. They also (2) give us information about the reliability of the resulting maps, through a general and a spatialised indicator (Fig. 2).

Regression methods represent one of the possible statistical approaches to map ES indicators from field and remote sensing data (Le Clec’h et al. 2017). Several regression methods can be computed to provide ES predictions based on both qualitative and quantitative inputs (Le Clec'h et al. 2013). Obtaining prediction intervals on the predictions can be quite challenging, in particular when no analytic expression of such intervals is available. For example, non-parametric methods, such as regression trees or random forests, may require bootstrap techniques to do so. However, linear modelling is interesting, for it allows the calculation of confidence intervals on predictions (Cornillon and Matzner-Lober 2011). In the following, we present results provided by linear regression methods.

We computed a linear model on the 135 sampling points. The model linked the ES indicator (carbon stocks - field data) and remote sensing data (related to the site, land cover and topography – Fig. 2). We manually performed a selection of variables, using Mallows' Cp Mallows 2000. The selection of variables aimed to (1) identify and select the best subset amongst the remote sensing variables and to (2) avoid overfitting. Thus, this selection allows us to obtain a model, which has good predictive properties. We called it, in this paper, “final model”, the model resulting from the selection of variables and that was used for the prediction of new ES values. This final model was based on two variables: the land-cover classification and the site (split into two classes: Palmares II and the two other sites, Table 2). The model can be considered as reliable (R² = 0.75).

We used the final model to extrapolate the vegetation carbon stocks to the whole study sites. To do so, the final model, trained on the 135 sampling points, was then applied to all the study sites to predict new carbon stocks, based on the land-cover data and the typology of the study sites (for more information, see Le Clec’h et al. 2017). Thanks to these predicted values, we obtained a map of the vegetation carbon stocks.

In the second step, we used the same final model to estimate a simple indicator that gives information on the confidence we have in our prediction. In linear regression, one way to associate confidence/uncertainty to the prediction is the calculation of the prediction intervals on predictions (Cornillon and Matzner-Lober 2011). The prediction intervals of the prediction are a range that is likely to contain the mean response given specified settings of the predictors in our model. Prediction intervals assess the accuracy of the estimated ES indicator. Therefore, to characterise the confidence we have in our prediction, we evaluated and mapped an index based on the prediction intervals. This index allowed us to have spatial information about uncertainty of the estimated carbon stocks. Such information is useful to identify the areas where the prediction is very reliable or, on the contrary, where the predicted values should be considered carefully.

We proposed an index based on width of 95% prediction intervals. To do so, we calculated the difference between the upper and lower bounds of prediction confidence intervals. A high index characterised areas with high uncertainty around the prediction (up to 90 MgC/ha around the prediction). For each predicted value (pixel), we thus applied the final linear model to get (1) an estimated carbon stock and (2) the prediction interval around the predicted value. As we predicted an ES value for each pixel of the study sites, we also obtained an uncertainty index for each pixel. Thus, we can propose a spatial representation of the uncertainty. Finally, we represented the variations of the index within the two explanatory variables of the final model: the land cover and the site classifications.

Maps of carbon stocks and their global reliability

The linear model based on the land-cover classification and the site classification was used to map vegetation carbon stock. The maps of vegetation carbon stock show the influence of land-cover changes on ES supply (Fig. 3).The highest values are located in forested areas, with the lowest values in deforested areas: farms and riversides in Maçaranduba, the main road in Pacajá and the southern part of Palmares II close to the city, influenced by the railway and the road. There are differences amongst study sites, because of local specificities such as differences in terms of biophysical conditions as human pressure, political and socio-economic factors that affect vegetation conditions and thus carbon stock. In Palmares II, the highest stock is lower than the one in the two other sites because the forests in Palmares II are highly degraded.

Mapping uncertainty associated with the ES map

We mapped the uncertainty index related to the application of a statistical method to predict values of an ES indicator for the vegetation carbon stock. The resulting map gives an overview of the areas where the map is and is not reliable. High prediction intervals express low confidence and can be associated with our inability to reliably predict carbon stocks, whereas high prediction intervals express our ability to estimate the vegetation carbon stocks with much greater certainty.

The index of uncertainty takes high values in highly anthropised areas (grazed landscapes and bare soils). Grazed grasslands have very heterogeneous profiles because they can be declined from bare soil to pasture with trees. Moreover, forest areas are associated with high carbon stock and quite high variability. Forests of the three study sites are at different stages of degradation and carbon stock is therefore relatively heterogeneous within these forests. Variability in transition areas can be explained by the nature of the class itself. Transition areas consist of secondary vegetation and fallow lands. This class is relatively homogeneous in Maçaranduba (Fig. 4) and in Pacajá, whereas it is more heterogeneous in Palmares II.

Role of land cover and location

We analysed how the uncertainty index varies with the two explanatory variables of the final linear model: the land cover and the site classifications. This analysis helped us to better understand the level of uncertainty related to our assessments. To do so, we plotted the prediction intervals (1) within the six land-cover classes and (2) within the three study sites. Such analyses highlight the uncertainty related to the explanatory data that our knowledge of the study sites (socioeconomic, political and/or biophysical contexts and conditions) can explain. The range of the prediction intervals modelled, based on the final model (Table 2), varies amongst the land-cover type (Fig. 5A). Prediction intervals are low in degraded forests, even if this class is highly diversified, in terms of density, strata and type of vegetation, due to natural and human-based factors. Agricultural areas are very diverse, especially the bare soil class, which regroups two different realities that can be discriminated with the analysis of our satellite data. On the one hand, this class includes bare soils such as roads and tracks. On the other hand, it regroups as well annual crops that are actually bare soil without plant cover, during the dry season. The reliability of our prediction is also very dependent on the location, independently of the differences in land-cover between the sites (Fig. 5B). In Maçaraduba and Pacajá, our predictions are reliable with a low variability, whereas predicted values in Palmares are highly uncertain. For instance, in Maçaranduba, sharp contrasts in the landscape exist, between large areas of preserved forests and pastures. In Pacajá, despite the landscape gradient, there are also sharp contrasts in the landscape between well-conserved forests and already cultivated / pastured areas. In Palmares II, differentiation between land covers was less clear: pastures were abandoned then transformed into fallow lands and, sometimes, were transformed back into pastures. Consequently, for this site, continuity in land-cover types existed, leading to an absence of spatial contrasts that can influence our capacity to obtain reliable statistical ES predictions.

A map is a generalisation or a schematisation of reality (Palsky 2013) and, accordingly, includes uncertainties or even errors. The issue of uncertainty can become particularly important when maps aim to support policies, for example, when they are used for land planning or to allocate budget to areas considered as ES hotspots (Bagstad et al. 2018, Eigenbrod et al. 2010, Palomo et al. 2018). In ES mapping processes, there is a plurality of uncertainties, inherent in each step of the modelling, for any methods that have been implemented. Scientific literature has increasingly acknowledged the need for considering uncertainties in the ES assessments (Bagstad et al. 2018, Grêt-Regamey et al. 2013, Schulp et al. 2014). There are, however, different types of uncertainties, from epistemic to scenarios inherent in mapping ES (Lavorel et al. 2017). Quantifiable and spatially explicit information about different types of uncertainty, as we presented here, help to qualify and critically assess ES maps and reveals the limits of the presented information, irrespective of the underlying methodological approach. However, the range of uncertainties produced in an ES assessment is much broader than the uncertainties induced by the statistical prediction (Hamel and Bryant 2017) and more systematic assessments are needed.

We proposed an exploratory approach to assess and communicate epistemic uncertainties related to ES mapping and valuation and thus specifically addressed a lack in the representation of errors inherent to the extrapolation of point-based measurements to produce empirical based ES maps. It is only one of the uncertainties produced in an ES assessment (Hamel and Bryant 2017), but it is also one that can be easily identified and communicated and/or reduced in order to create reliable outputs that can support policy-making or land management. We extrapolated measured field data using regression analysis to map carbon stocks for three study sites in the Brazilian Amazon and evaluated the variability of the ES provision across space. An indicator based on prediction intervals gave an overview of the spatially explicit reliability in the statistical procedure by enabling the identification of places where the prediction is reliable and where it has to be improved and cannot currently be used to base decisions.

Overall, the regression model provided a reliable map of carbon stocks that allowed us to highlight carbon hotspots and, at the same time, to identify areas where this prediction of carbon stocks is highly uncertain. The results showed that prediction intervals for bare soils and clean pastures had higher mean values but also a wider range of values and thus emerged as highly variable classes. The map’s trustworthiness is thus partially a consequence of the uncertainty related to the explanatory data underlying these classes. From a scientific perspective, these findings have three implications. First, a better understanding of the drivers of ES provision could help to reduce the uncertainty stemming from the variability of the provision of ES per land-cover class. This aspect is very critical in the Amazon rainforest where the slowdown of deforestation rates leads to an increase in forest degradation. The consequences of forest degradation on ES provision in these areas are still largely unknown. Secondly, the land-cover classification should account for as many different classes as possible to reduce intra-class variability. This underlines the importance of studies on local scales with detailed information on land use/land cover to complement large scale maps with less accuracy in land-use classes (Foody 2015, Domac 2004). The availability of new satellite information such as PROBA-V might help to fill this gap given its relatively high spatial resolution and short processing time. Thirdly, despite this potential of new remote sensing information, data uncertainties highlight our limited ability to model some specific ES indicators. On the one hand, variations in the vegetation carbon stock, for instance, can be explained by modified properties of ecosystems such as land-cover changes. On the other hand, indicators related to soil or biodiversity are greatly influenced by inherent ecosystem properties such as biophysical processes that can rarely be inferred from remote-sensing data (Dominati et al. 2010). In addition, high uncertainty areas can vary from one ES to another and, for one specific ES, vary from one spatial scale to another (Le Clec’h et al. 2018). Thus, a systematic and transparent assessment of such uncertainties helps to reconsider and improve the methodological approach and the data used in the mapping process.

The use of prediction intervals to assess epistemic uncertainties in ES maps has also certain drawbacks. Not all methodologies allow a feedback between data and ES maps. The indicator proposed here, can only be applied in the case of certain statistical procedures. Even for some statistics, such as regression trees (CART algorithm - Breiman et al. 1984), the calculation of prediction intervals is not possible (Yohannes and Hoddinott 1999). Moreover, the indicator represents the variance in the prediction and does not contain information on the goodness of fit to the response variable. This point seems to be critical, as previous studies demonstrated that the reliability of the map in general could vary strongly from one indicator to another (Le Clec’h et al. 2017). Thus, we here used a more general indicator, such as the RMSE or the R², to complement the information on the prediction uncertainty. General information on the map´s reliability could help the cartographer to decide whether he/she should stop or proceed with the mapping process. However, it raises in turn the critical question of what an acceptable threshold of reliability is. Thus, the choice of threshold values should be transparent and be subject to sensitivity analysis. Finally, we are fully aware that the range of uncertainties in an ES assessment is much broader than the epistemic uncertainties induced by the statistical prediction (Hamel and Bryant 2017).

Uncertainty can become particularly important when maps aim to support policy development, for example, when they are used to analyse trade-offs and synergies (Bagstad et al. 2018, Eigenbrod et al. 2010, Palomo et al. 2018). Information about uncertainties in these maps can help to prevent unwanted or unintended consequences of policy measures that affect land-use decision or management (Jacobs et al. 2017). In addition, the provision of information about spatially explicit trustworthiness of ES maps helps to evaluate policies. Maps and information about their accuracy and flaws can support the spatially explicit identification of the relevant socio-economic and biophysical characteristics that lead to an (in-)effective policy implementation.

Uncertainties are inherent and unavoidable in the assessment of environmental management in general and in particular in ES mapping. They exist in different stages of the assessments and some cannot integrally be reduced. We used here a linear model to describe the method because prediction intervals can be directly calculated. This approach could be extended to other methods, especially by the use of bootstrapping. However, their identification, acknowledgement and explanation should, at least, be systematic, to raise awareness and to determine the optimal use(s) of the maps, especially in the context of environmental policy-making and governance (Jacobs et al. 2017). Information related to uncertainty, such as the map presented in this paper, might be difficult to apprehend for non-specialists. Accordingly, further research should not only provide more holistic perspectives on uncertainties in ES mapping, but also offer insights on how to communicate such information. Further work should focus on what could be the best way(s) to map such information and how to adapt these way(s) to the purpose of the map. In that sense, our approach should be seen as a first step towards a systematic consideration and acknowledgement of uncertainties that would be included as a strong and integral component of ES mapping.

This research was funded by the Institut des Amériques and by the French Agence Nationale de la Recherche through two grants: ANR AMAZ, coordinated by P. Lavelle and ANR AGES, coordinated by X. Arnauld de Sartre and by the Brazilian Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq), also through two Grants: Processes 484990/2007-1 and 490649/2006-8.