Predicting Abiotic Soil Characteristics Using Sentinel-2 at Nature-Management-Relevant Spatial Scales and Extent

1. Introduction

Knowledge of local plant community characteristics is imperative for practical nature planning and management, and for understanding local-scale plant diversity and distribution drivers. Today, such data of high enough quality for these purposes can only be derived through targeted field campaigns, which are costly both in time and money. Here we use multispectral high-to-medium resolution (10–60 m) satellite data [1] to predict the local composition of plant species over a large extent (approx. 30.000 km²) in terms of their preference for soil moisture, soil fertility, and pH which are among the most important factors for local plant community composition [2,3].

In recent years, massive vegetation plot databases have provided large amounts of data enabling studies of local plant diversity at national and continental extents [4,5,6]. Combining these with indicator values such as Ellenberg opens for novel insights into how local plant communities are distributed at broad spatial scales [7,8]. In Denmark, the monitoring of protected natural and semi-natural habitats (protected under the EU Habitats Directive Annex I, i.e., the nature deemed the most important to the EU) generates a significant amount of vegetation plot data each year [9,10,11].

At the same time remote sensing techniques are slowly maturing and are approaching a stage where they can be used to predict local plant diversity over large areas [12,13]. A large number of studies have attempted to use remote sensing techniques to predict ecological factors like plant species richness, plant phenology, plant traits and habitat characteristics [4,12,14,15]. However, most of these studies have done this only at local scale or involving only one or a few species, traits or habitat types, preventing broad scale generalization. At the same time no study has succeeded in providing sufficiently detailed insight into a given area based solely on remotely sensed information to use remote sensing as a stand-alone tool for monitoring and management of nature. This is likely related to the fact that retrieving some of the most important abiotic factors for plant species distribution – namely soil chemical factors such as soil fertility, soil moisture and pH [16] – with remote sensing is notoriously hard as most remote sensing signals do not penetrate into the soil and hence do not gather data directly from the soil.

Predicting soil moisture, fertility and pH using remotely sensed data in nature areas has been attempted several times, but often at too coarse a spatial resolution for being relevant for practical nature management [17] or across fairly small spatial extents preventing generalization of results [18]. The hypothesis behind using multispectral satellite imagery to generate such predictions is that soil chemistry is thought to affect the plants structure and leaf chemistry causing differences in reflectance of plants depending on the underlying soil chemistry and moisture [19,20] and hence the plants can act as indicators of soil chemical properties. For example, Schmidtlein (2005) and Möckel et al. (2016) [21,22] used imaging spectroscopy to map these factors at high spatial resolution, but the study areas were relatively small. On the other hand, studies combining relatively fine resolution and large extents are starting to improve our ability to predict these important soil factors at relevant scales and extents for nature monitoring and management [23].

Here we take the next major step to bridge this gap by using plant community data from > 50,000 small (5 m circular) vegetation plots with indicator values for soil moisture, soil fertility and pH combined with current multispectral satellite data to show how this can be used to predict some of the most important abiotic soil factors for local plant diversity and distribution – soil moisture, soil fertility and pH – at sufficient detail to be relevant for nature monitoring and management. More specifically, we investigate (1) if 9 bands of Sentinel-2 data are linked to plants’ preference for soil moisture, -fertility and -pH, and if so (2) how well can we predict these plant-indicated soil factors using current state-of-the-art remote sensing techniques.

2. Materials and Methods

Vegetation Data

In this study, we used vegetation data from the Danish national program for monitoring of nature (NOVANA, freely available from naturdata.miljoeportal.dk [9,10,11]). On 16 September 2016 we extracted higher plant species lists for 97,334 vegetation plots (5-m diameter radius). The plots are randomly placed within sites belonging to 32 different habitat types, all of them listed in the EU Habitats Directive’s Annex I [24]. These natural and semi-natural habitats range from grasslands, over heathlands, dunes, and shrublands to wetlands, but in this analysis, we excluded forests (Figure S1). The plot-data were collected by trained field botanists from 2004–2015 and cover all parts of Denmark (Figure 1, but here only the plots used for this study are shown). For analysis, we used 57,871 plots of which 31,001 were unique, meaning that some plots were revisited several times over the years (See Figure S2). Though, since the vegetation plots were georeferenced in the field using traditional GPS having a vertical uncertainty of 5-10 m, revisits were not exactly at the same place at every revisit. The exact GPS models used differed depending on inventory year and region.

To represent the three main drivers of local plant diversity; soil moisture, -fertility, and -pH, in each plot, we calculated – based on the plant species list for each plot – the mean Ellenberg indicator values (EIVs [20]) for soil moisture (EIV_F), soil nutrients (EIV_N) and soil reaction (EIV_R). These indicator values represent the plants’ preferences for soil moisture, the fertility of the soil and its pH respectively. Earlier studies have shown that these EIVs are typically rather tightly correlated with corresponding measured abiotic variables [8]. Soil nutrients and pH is often highly correlated, so to represent the nutrient availability independently from pH, were also calculated the EIV_N/EIV_R ratio [25], referred to as the N/R ratio in the following. Not all species had an Ellenberg indicator value – in that case they were not included in the mean calculation.

Each plot also holds information on which habitat type was present. This information is based on a preliminary assignment in the field and a subsequent verification process involving the species found in the plot [11].

Remote Sensing Data

In this study, we used multispectral imagery from the European Space Agency Sentinel-2 mission [26]. To match the plants’ growing season, we retrieved all images overlapping the study area from 1 June to 31 August 2016 through the Copernicus Open Access Hub (https://scihub.copernicus.eu) during the weeks 19 and 20 in 2018. The total dataset consisted of 38 image products (see Appendix A for specific product names). To obtain surface reflectance, all images were atmospherically corrected using the MAJA ver. 1.0 software tool [27] with default settings. Unlike most other software tools for satellite imagery processing, MAJA takes advantage of using imagery time series to estimate cloud mask and aerosol optical thickness. Thanks to the use of this multi-temporal information, MAJA cloud masks are highly reliable [28] and the accuracy of the surface reflectance estimates computed is among the best when applied to Sentinel-2 [29,30]. This processing was conducted in a Docker installation of CentOS version 7.4.

Preprocessing

To provide a basis for testing how well Sentinel-2 imagery can predict local soil moisture, fertility and pH levels, we extracted the raw values for the band number 2, 3, 4, 5, 6, 8, 8a, 11 and 12 (after atmospheric correction) corresponding to the position of each vegetation plot. Due to missing coverage or clouds that overlapped the vegetation plots, we decided to use the satellite data from one single cloud free day (16 August 2016) as that resulted in the largest possible dataset with 58,071 vegetation plots having band-data from the Sentinel-2 images and the above-mentioned average Ellenberg indicator values. For that particular day, satellite data was not available for the whole country and therefore our study is only based on plots from approximately 2/3 of the country (Jutland, Figure 1). However, data does cover all main soil- and habitat types despite of this (with the exception of rock habitats that constitutes a tiny part of the Danish landscape), so probably the model is quite broadly applicable both in the rest of Denmark and in other parts of the North-European temperate region.

Statistical Analysis

To test if the Sentinel-2 data is linked to the abiotic environment in the plots (the EIVs were used as response variables), we used supervised learning AI-algorithms (see below) on a training data set consisting of 30,000 randomly selected vegetation plots.

The nine intensity values from Sentinel-2 were treated as a vector with 9 numeric elements. Using the raw intensity values instead of the traditional indices like NDVI (see e.g., reference [31]) has previously been shown to work well for this kind of modelling [22]. To investigate if the links between the imagery and the EIVs were influenced by habitat type, this information was also included as a predictor (categorical variable).

To make sure to get the best model, we let the function “Predict” in Wolfram Mathematica version 14 automatically select which machine learning algorithm that provided the best predictions for each EIV and for satellite data only, habitats only and these two data sources combined respectively (Table 1). The “random forest” algorithm was the most frequently selected method, but also the algorithms “decision trees” and “nearest neighbor” were found to provide the best predictions in some cases (Table 1).

To investigate the within-plot among-year variation of the calculated Ellenberg indicator values, the indicator values were analyzed in in a mixed linear model with habitat type as a fixed factor and plot as random factors assuming that the residual variation was normally distributed. The R (ver. 4.3.2) procedure lmer (ver. 1.1-35.1) in the package lme4 [32] were used for this part of the analysis.

3. Results

Overall, we found that the Sentinel-2 imagery explained a significant part of the variation in the plot-level abiotic environment (Figure 2, Table 1), with > 26, 59 and 54% of the variation explained in soil moisture, fertility and pH levels indicated by the plant species composition, respectively. This means, that satellite data is linked to and can indeed be used to predict plant indicated local scale soil moisture, nutrient levels and pH across numerous habitat types with some uncertainty. We also found that combining satellite data with information about the plots’ habitat type considerably improved this predictive power (over 70% of variation explained in all cases mentioned above), but also that habitat – when used alone as predictor – consistently had no or lower predictive power than using satellite data alone (Table 1). We found significant within-plot among-year variation in the calculated Ellenberg indicator values (Table S1). See result details in Figure 2 and Table 1. For examples showing the actual vs. the predicted values in selected sites see Figure 3.

4. Discussion

Soil Fertility and pH

Earlier studies have been able to explain about 35–40% of the variation in plant indicated soil fertility using remotely sensed spectral data [31,33]. However, we were able to explain almost 60% which is in good agreement with Möckel et al. 2016 [22], who reached a similar result. Both Möckel et al. (2016) [22] and we tried to utilize the information in all available spectral bands, whereas previous studies have typically used NDVI – which is only based on two bands – as the only spectral predictor. This suggests that there is indeed important information for predicting soil fertility not only in the red and infrared bands but in other bands as well. We therefore recommend that this full suite of information (i.e., all spectral bands) is used in future studies and applications of multi- and hyperspectral remote sensing data when applied to issues related to soil fertility indicated by plants. In our study, the nutrient ratio was predicted less well. This parameter integrates both nutrient and soil reaction and hence indicates the actual soil fertility independently from soil reaction and this is often important in ecological questions as it mirrors the actual nutrients available for plant growth [25]. Consequently, this result means that it can be quite complicated to use remotely sensed spectral data to gain insight into plant-relevant soil fertility dynamics at local scale. That said, our relatively good ability to predict soil nutrient status (i.e., only EIV_N) is indeed encouraging in that respect, and we do believe there is a potential for further development of methods, that could improve our results (see “Perspectives”). When adding habitat information to our models, the predictive power for both soil fertility and pH rose to 70% or more, indicating that knowledge of habitats further refines our models and should be used whenever available.

Soil Moisture

In contrast to our soil fertility and pH (EIV_N and EIV_R) results, we achieved poorer results when it comes to predicting soil moisture than earlier studies (R² = 29%). For example, both Möckel et al. (2016) [22], Weber et al. (2018) [33] and Löfgren et al. (2018) [31] were able to explain >50% of the variation in soil moisture indicated by the plants. However, their studies were all conducted solely in dry grasslands, and that could explain this difference: In dry grasslands the soil moisture gradient is quite short compared to when all the different habitat types are considered from bogs and fens, over fresh- and saltwater meadows to moist and dry grasslands like we do here. We believe that having this longer moisture gradient causes the spectral signals to vary much more independently of the soil moisture. For example, it seems intuitive that in dry grasslands, the plants in moister areas are greener than those in drier areas. In other habitats there may not be this difference in the appearance of the plants, which could explain our lower explanatory power. Hence, for endeavors to remotely sense soil moisture it may be necessary to include other techniques, such as lidar-based topographical wetness indices [34]. Also, our modelling highlighted the importance of having access to habitat type data to make good predictions of soil moisture from spectral data.

Limitations and Uncertainties

Following from the predictive power of our models being lower than what would be acceptable for everyday nature management and monitoring, our approach does have limitations when it comes to practical applications. Our models may not be precise enough for such uses. On the other hand, plant indicator values themselves are also uncertain and not necessarily mirroring abiotic soil conditions one to one [35,36,37], and clearly, predicting uncertain factors is bound to come with some uncertainty to the results. That said, plant indicators do integrate abiotic conditions over time and space like no other data source [8] and therefore our ability to predict them with remotely sensed data is core for developing future digital nature mapping and monitoring solutions.

Our models were built on satellite imagery from only one day. This means that first of all, we are missing information on how the spectral signal changes over the year, and secondly, that the spectral data does not necessarily match the recording year of the floristic data from the vegetation plot. Therefore, future works should test these models on data from multiple days over some years to make sure they are not overfitted to data from a certain day, and hence are sufficiently general for applications under day-to-day and year-to-year varying phenology and long-term weather conditions (e.g., moist and dry summers).

The within-plot among-year variation in the EIVs is probably partly due to the GPS-uncertainty rendering it impossible to find the same position when revisiting vegetation plots. It could also be due to changes in plant community composition over the years. Considered together with the fact that our satellite data only captures a glimpse of conditions in time (see above), this means that we cannot expect perfect accuracy of the local predictions of our models.

For soil moisture and pH (EIV_F and EIV_R), including habitat type in our models strongly improved predictive power while models using solely habitat type had relatively poor predictive power. This probably means that the abiotic and biotic information represented by the habitat type is important to make good models of soil abiotic conditions, and at least for some habitat types getting this information from remote sensing is notoriously hard [38].

Perspectives

For soil nutrients (EIV_N) and pH (EIV_R), we achieved considerable predictive power (R²: 0.5 – 0.6) using satellite data alone. This means that our method could have potential for further development, e.g., by using finer resolution data (both spectrally and spatially), for example from drones or airplanes [39,40]. Several of our models were machine learning based (random forests or gradient boosted regression trees). Another possible future improvement could be to try implementing deep-learning-based methods as they may be better at capturing patterns in the data [41].

Including habitat type improved our models significantly so being able to remotely sense this information and feed it into models like ours would be clearly desirable. Currently, researchers and technicians from several different parts of the World are or have been working on such remote-sensing-based classification of habitat types over large areas [42,43,44]. While these models are still far from perfect, they are slowly improving and that could possibly enhance our ability to improve models like the ones we developed here.

While satellite data often suffers from poor resolution compared to drone or lidar data, they have a clear strength in that they offer time series data (i.e., repeated data recording of the same area). This opens up for developing models that can predict changes in soil characteristics over time. Such change-detection models will possibly experience the same kind of error as we showed here, but we see clear perspectives in such models and urge researchers to test this and pursue this research direction in the future.