Mapping of peatlands in the forested landscape of Sweden using lidar-based terrain indices

Globally, northern peatlands are major carbon deposits with important implications for the climate system. It is therefore crucial to understand their spatial occurrence, especially in the context of peatland degradation by land cover change and climate change. This study was aimed at mapping peatlands in the forested landscape of Sweden by modelling soil data against lidar-based terrain indices. Machine learning methods were used to produce nationwide raster maps at 10 m spatial resolution indicating the presence or not of peatlands. Four different definitions of peatlands were examined: 30, 40, 50 and 100 cm thickness of the organic horizon. Depending on peatland definition, testing with a hold-out dataset indicated an accuracy of 0.89–0.91 and Matthew's correlation coefficient of 0.79–0.81. The final maps showed a national forest peatland extent of 60 292–71 996 km², estimates which are in the range of previous studies employing traditional soil maps. In conclusion, these results emphasize the possibilities of mapping boreal peatlands with lidar-based terrain indices. The final peatland maps are publicly available at https://doi.org/10.17043/rimondini-2023-peatlands-2 (Rimondini et al., 2023) and may be employed for spatial planning, estimating carbon stocks and evaluating climate change mitigation strategies.

1 Introduction

Peatlands support biodiversity, buffer hydrological cycles, and sequester and store atmospheric carbon (C). They have historically been important C sinks, and it is estimated that peatlands in the Northern Hemisphere hold 415 ± 150 Pg C (Hugelius et al., 2020). As peatlands are degraded by climate and land cover change, this long-term C sink is, however, becoming a substantial source of greenhouse gases (GHGs). Peatland draining for agricultural purposes is a large contributor to carbon dioxide emissions (e.g. Günther et al., 2020), and permafrost thaw in peatlands is projected to lead to large methane emissions (e.g. Hugelius et al., 2020). Peatlands are thus key environments for climate change mitigation strategies, and there is an urgent need to implement global peatland restoration and protection policies (Beaulne et al., 2021; Günther et al., 2020; Humpenöder et al., 2020). High-quality data on the extent and depth of peatlands are crucial to support decision-making and reach the international goals for sustainable development.

Traditionally, landscape-scale peat mapping has been conducted through aerial image interpretation combined with logistically challenging field inventories. Today, automated algorithms for digital mapping with remotely sensed data have made it possible to circumvent these limitations (Minasny et al., 2019). Optical and synthetic aperture radar (SAR) satellite imagery and lidar-based terrain indices are used for peatland delineation and depth estimation, often together with auxiliary data on climate and soil. Limitations for these modern techniques are instead set by data availability, especially lidar-based elevation data and field data for model development.

Sweden presents a good case study for digital mapping of peatlands as lidar-based elevation data are readily available, and field data on soil properties have been collected by long-standing nationwide surveys. In addition, the large extent of Swedish peatlands and their considerable share of drainage, especially in productive forest areas, lead to important GHG emissions warranting further research and mitigation policies (Lindgren and Lundblad, 2014; Skogsstyrelsen, 2021).

Swedish peatlands have been mapped with both traditional and modern methods, with estimates ranging between 15 % and 20 % of the total land area (e.g. Barthelmes et al., 2015; Joosten, 2010; Pahkakangas et al., 2016; Tanneberger et al., 2017; Ågren et al., 2022). Ågren et al. (2022) produced a series of digital Swedish peatland maps at 2 m spatial resolution representing different peat depths. These maps are reported to perform better than traditionally produced maps. Generated by fitting a cubic relationship between O-horizon thickness and a soil moisture product, these maps show greater total peatland cover than traditional maps due to the inclusion of previously unmapped smaller peatlands. In this study we used partly overlapping O-horizon data but employed multiple independent variables and machine learning modelling methods to map Swedish forest peatlands.

Study aim

We aimed at mapping the extent of forest peatlands in Sweden, by analysing O-horizon thickness measurements against lidar-based terrain indices. Machine learning methods were applied to construct models able to binary classify shallow peatlands with ≥30, 40 or 50 cm O-horizon thickness and deeper peatlands with ≥100 cm O-horizon thickness. Maps of the classifications were produced as 10 m raster grids at a national scale.

2 Materials and methods

2.1 Study area

This study was limited to the forested landscapes within the country of Sweden, located on the Scandinavian Peninsula. Sweden has a land area of 407 284 km², of which approximately 69 % is covered by forest (Statistiska Centralbyrån, 2022). The country has a pronounced south–north gradient between $\mathrm{55}{}^{\circ }{\mathrm{20}}^{\prime }$ N and $\mathrm{69}{}^{\circ }{\mathrm{3}}^{\prime }$ N, with the consequent effect on the climate ranging between temperate and subarctic. The vicinity to the Gulf Stream and continental west coast setting of Sweden enforce a warmer climate than other areas at similar latitudes. The annual rainfall generally varies between 500 and 800 mm with local extremes of 400 and 2000 mm (Sveriges Meteorologiska och Hydrologiska Institut, 2009). The topography in Sweden is mainly influenced by the Caledonian orogeny, pre-glacial deep-weathering and Quaternary glaciations. The bedrock is generally igneous and metamorphic except in the southernmost parts of Scania County, Öland, Gotland and various isolated pockets in mainland Sweden. Soils are predominately of Quaternary origins and mostly till, with coastal areas having been influenced by wave action due to isostatic depression (Fredén et al., 2009).

In this study we only produced maps of peatlands in Swedish forests and treeless wetlands, as the employed O-horizon thickness data were not available in non-forested areas (see Sect. 2.2.1). A forest mask was used to delineate areas of interest, consisting of a reclassified version of the Swedish land cover dataset from 2018, Nationella Marktäckedata 2018 (NMD2018), by the Swedish Environmental Protection Agency. See Sect. 2.2.3 for further descriptions.

2.2 Data

2.2.1 Dependent variable: O-horizon thickness

A Swedish dataset on O-horizon thickness was constructed by merging two data sources. The first source, the Swedish Forest Soil Inventory (SFSI), is a nationally covering soil survey for which O-horizon thickness has been measured in the range 0–99 cm (Olsson, 1999). The O horizon is defined as a layer of organic material with ≥20 % organic carbon. The SFSI inventories are designed to cover forested land across Sweden and exclude agricultural areas, urban land, infrastructure and mountain regions.

The second data source was the peat archive of the Geological Survey of Sweden (SGU). Peat core measurements collected by peat experts over the past 100 years are stored in this archive, with varying methods and data quality. Unlike the SFSI, no lower limit in peat depth is set.

The accuracy of the geographical positioning is 5–10 m in the SFSI (Olsson, 1999), while for the peat core data from SGU it depends on the measurement period and method. In total, the merged dataset includes 10 115 data points, of which 7175 are from the SFSI and 2940 from the SGU peat archive; their spatial distribution can be viewed in Fig. 1.

Figure 1 Locations of the O-horizon thickness sampling points included in this study. Black points are part of the Swedish Forest Soil Inventory (SFSI), whereas the red points are from the peat core archive by the Geological Survey of Sweden (SGU).

Four thresholds were used to divide the O-horizon thickness dataset into separate binary variables. The first three thresholds were set to 30, 40 and 50 cm. This is the minimum O-horizon thickness of peatlands according to a number of national and international soil classification systems (Lourenco et al., 2022). The last threshold was set to 100 cm, based on the measurement range limitations of SFSI (no data reported ≥100 cm). These binary variables (peat deeper or shallower than 30, 40, 50 or 100 cm, respectively) were used separately as dependent variables in machine learning workflows, with the final goal to produce four different geospatial products: maps over peatlands with O-horizon thickness ≥30, 40 and 50 cm and a map over deeper peatlands with O-horizon thickness ≥100 cm.

2.2.2 Features

The independent variables tested for prediction of peatland occurrence, hereby referred to as features, can be seen in Table 1. We hypothesized that these features indicate peatland occurrence at the landscape scale in Sweden and made the assumption that only lidar-based variables were sufficient for prediction. Nonetheless, peat formation and thickness are complexly controlled by climate, hydrology, topography, soil and vegetation (Limpens et al., 2008). In earlier versions of the models presented in this article, some variables on climate and vegetation were thus included, but these did not enhance the classification performance and were removed.

Table 1 Raster datasets tested as features in this study. Features denoted with an asterisk were included in the final models used for peatland mapping; the others were excluded during pre-processing or feature selection (see Sect. 2.3.2).

DEM10 is a lidar-based digital elevation model (DEM) with a spatial resolution of 10 m in the geographical coordinate system SWEREF99 (EPSG: 3006). DEM10 is the base-layer of all the other features, except MarkFuktIDX.

Topographic position index (TPI; Weiss, 2001), Topographic ruggedness index (TRI; Riley et al., 1999) and Roughness (Dartnell, 2000) describe the elevation heterogeneity of a landscape. TPI is the elevation difference between a DEM pixel and the mean of a determined pixel neighbourhood surrounding it; TRI is defined as the square root of the sum of the squared difference between a DEM pixel and each pixel in its eight-cell neighbourhood; the term Roughness describes the largest elevation difference between a DEM pixel and a determined pixel neighbourhood.

The term Slope describes the inclination in a DEM neighbourhood, in our case computed using an algorithm by Horn (1981).

DTW is a topographic index developed by Murphy et al. (2007), approximating the elevation difference between a pixel in a DEM and a hydrological source. The hydrological sources are defined as pixels with an upslope area larger or equal to a set threshold, the so-called flow initiation area, generally ranging between 0.5–40 ha. Source pixels are allocated for each pixel in a slope raster, applying a least-cost-path algorithm; the accumulated slope values along the path are multiplied by the resolution of the raster to obtain the DTW index. Mathematically, DTW is defined as

$$\begin{array}{}\text{(1)}& \text{DTW}=\left[\sum {\displaystyle \frac{\mathrm{d}{z}_i}{\mathrm{d}{x}_i}}a\right]x_{\mathrm{c}},\end{array}$$

where dz_i and dx_i are the vertical and horizontal distance (in metres) between two cells; a is a constant of 1 or $\sqrt{\mathrm{2}}$ depending on whether two adjacent cells connect parallelly or diagonally; x_c is the raster cell size (in metres).

DTW has been used to identify wet areas with the advantage of being independent of spatial scale (Ågren et al., 2014). The optimal threshold for flow initiation area is, however, dependent on soil characteristics, as areas with higher hydraulic conductivity are better examined with a higher flow initiation area, and vice versa.

MarkFuktIDX is a Swedish soil moisture index product created to support delineation of wetlands in the national land cover dataset (NMD2018). It was created by merging DTW rasters with topographic wetness index (TWI) rasters, with an unbalanced weight towards DTW. The DTW raster was calculated using a 2 ha flow initiation area threshold, and the TWI raster included information on soil permeability where this was available.

2.2.3 Forest and treeless wetland mask

The forest and treeless wetland mask used to limit the extent of the analyses in this study consisted of a binary reclassified version of NMD2018. The positive class included forest classes (codes 111–128) and open wetland (code 2); the negative class (where no analyses were made) includes arable land (code 3), open land (codes 41 and 42) and artificial surfaces (codes 51–53).

2.3 Data processing

2.3.1 Raster processing

The topographic variables were calculated using a 10 m resolution DEM (DEM10; see Table 1), originally resampled with the bilinear method from the 2 m resolution Höjddata 2+ grid DEM from Lantmäteriet.

TPI, TRI, Rough and Slope were calculated using the gdaldem utility (GDAL/OGR contributors, 2022). TPI, TRI, Rough and Standard deviation of slope (SDSlope) were calculated on an eight-cell neighbourhood basis.

DTW was calculated with a Python script using tools from the Python packages RichDEM (Barnes, 2016) and Whitebox (Lindsay, 2014). Four different flow initiation area thresholds were used: 1, 5, 10 and 20 ha. These were subjectively chosen but fall within the range of thresholds used in previous studies at similar scale (e.g. Oltean et al., 2016; Ågren et al., 2014, 2021).

2.3.2 Machine learning workflow

The rasters which constituted the features (see Table 1) were sampled at the geographical position of each O-horizon data point. In total, 10 115 data points were extracted. Of these, 4230 were labelled as “peatland” using the 40 cm threshold and 3686 as “deeper peatland” (≥100 cm).

The labelled dataset was pre-processed by data imputation, feature correlation analysis, undersampling, dataset splitting and standardization. Most utilities are from the scikit-learn Python package (Pedregosa et al., 2012), hereafter referred to as Sklearn.

Sklearn SimpleImputer was employed for imputation, replacing missing values in the features with their respective mean. This process was applied to avoid data loss, and because some machine learning algorithms cannot handle missing values.

Features with a Pearson correlation coefficient ≥0.8 between them were then removed from the labelled dataset, as highly correlated features can negatively affect the performance of machine learning algorithms, and redundant data enhance processing times and computational costs.

RandomUnderSampler from the imbalanced-learn Python package (Lemaitre et al., 2016) was applied to avoid detrimental class imbalance caused by overemphasizing of the most populous class. This algorithm removes randomly selected samples from the over-represented class to reach complete class balance.

The Sklearn train_test_split algorithm was used to randomly split the labelled dataset into a training and a testing dataset, with a $\mathrm{70}/\mathrm{30}$ split ratio. A standardization algorithm, Sklearn StandardScaler, was lastly fitted to the features of the training dataset and used to transform the features of all datasets by mean removal and data scaling.

In machine learning, feature selection is performed to find a balance between computational cost, model interpretability and model performance (Brank et al., 2011). In this study, the feature selection algorithm RFECV from Sklearn was run following pre-processing and prior to model training. RFECV is a wrapper method testing recursively smaller subsets of features by cross-validation and eliminating the weakest feature at each iteration.

Four machine learning algorithms from the Sklearn library were tested: support vector classifier (SVC), logistic regression (LR), random forest classifier (RF) and multi-layer perceptron classifier (MLP). Sklearn RandomizedSearchCV was used for hyperparameter tuning, which randomly tests subsets of parameters by cross-validation and selects the best-performing configuration.

Model configurations were evaluated using accuracy, precision, recall and Matthew's correlation coefficient (MCC) as performance metrics, all computed on the test dataset. These performance metrics were used to identify the most suited model configurations for final prediction mapping, referred to as Peat30 (O-layer thickness ≥30 cm), Peat40 (O-layer thickness ≥40 cm), Peat50 (O-layer thickness ≥50 cm) and Peat100 (O-layer thickness ≥100 cm). We regarded high precision as a positive model trait, as for peatland management we believe it is relatively more important to minimize false positives rather than false negatives. Therefore, if two models had similar accuracy and MCC, we prioritized the one with higher precision.

Permutation feature importance was calculated for each final model, which is a measure of the impact of each feature on predictive power. Permutation feature importance is calculated as performance lost, in our case accuracy loss, as the feature is permuted (i.e. shuffled) in the test dataset. Hereafter it is referred to as feature importance.

2.3.3 Final prediction mapping

Final predictions were executed with Peat30, Peat40, Peat50 and Peat100 on a raster grid of 10 m resolution, corresponding to the raster variables used as features (see Table 1). Four maps were thus produced: three for peatlands (O-layer thickness ≥30, 40 or 50 cm) and one for deeper peatlands (O-layer thickness ≥100 cm). Some areas in the mountainous regions are missing in the final products, as lidar data are not available there.

Gotland and Öland were also excluded in the final maps. These islands were scarcely sampled with 64 sampling points in Gotland and 16 points in Öland, of which only 1 point per island was classified as peatland using O-layer thickness ≥40 cm as definition. The signal from the sampling on mainland Sweden could therefore have adverse influence on the quality of prediction on these islands, which have relatively distinct geological characteristics.

Non-forested land was masked in the final maps using the mask described in Sect. 2.2.3. The maps were compressed using the gdal_translate Lempel–Ziv–Welch compression algorithm, and overview images were built internally using gdaladdo (GDAL/OGR contributors, 2022).

Map visualizations were created using ESRI ArcGIS Pro 2.9.2.

3 Results

Prior to model training, the same features were selected by the selection algorithm RFECV in all classification tasks: DEM10, DTW1ha, DTW10ha, Rough, SDSlope, TPI, MarkFuktIDX.

The test results are shown in Table 2, demonstrating that RF and MLP have marginally higher accuracy and MCC than LR and SVC in most cases.

Table 2 Performance metrics computed on the test set for the models operating on the various binary O-layer thickness thresholds: (a) 30, (b) 40, (c) 50 and (d) 100 cm. Asterisk and bold numbers denote the models selected for final prediction mapping, i.e. Peat30, Peat40, Peat50 and Peat100. LR = logistic regression, SVC = support vector classifier, RF = random forest and MLP = multi-layer perceptron.

In the 30 cm classification task, the RF model scored the highest MCC and accuracy values with balanced precision and recall values, thus being selected as Peat30 (Table 2a). In the 40 and 100 cm classification tasks, the RF and MLP models scored the same performance metrics (Table 2b and d); however, for the sake of consistency between the final maps, the RF models were selected as Peat40 and Peat100. In the 50 cm classification task, the SVC and RF models scored the same accuracy and MCC values; however, the RF model scored higher precision, thus being selected as Peat50 (Table 2c).

Permutation feature importance for all final models can be viewed in Fig. 2. The two most prominent features were the same in all models, MarkFuktIDX and Rough. In Peat100, MarkFuktIDX and Rough showed highly similar results compared to the other models showing larger gaps between them (Fig. 2). It is worth noting that a low feature importance does not necessarily mean that the feature is insignificant for a given classification task. Features with low feature importance in one model can be of greater importance in a differently configured model trained on the same data.

Figure 2 Boxplots showing the permutation feature importance of (a) Peat30, (b) Peat40, (c) Peat50 and (d) Peat100. Permutation feature importance is calculated as the loss of accuracy when a feature is permuted, in this case iterated 10 times. The boxes represent the lower and upper quartiles with an orange line at the median; whiskers indicate data ranges, and isolated point are outliers.

The final prediction maps are available at https://doi.org/10.17043/rimondini-2023-peatlands-2 (Rimondini et al., 2023) and excerpts of Peat40 can be viewed in Figs. 3 and 4 (Rimondini et al., 2023). Peat areas clearly follow topographic basins, and shallow peat often surrounds areas of deeper peat or is found in smaller patches and furrows. The estimated national extent of peatlands according to Peat30, Peat40, Peat50 and Peat100 can be viewed in Table 3.

Figure 3 Excerpt from the final Peat40 map. Orange areas are classified as having an O horizon ≥40 cm. The map covers an area with extensive aapa mire complexes in Udtja, Norrbotten County, Sweden. Base map: Ortofoto, 0.5 RGB Tiles 735_67_50_2014 and 736_67_00_2014 © Lantmäteriet.

Figure 4 Excerpt from the final Peat40 map. Orange areas are classified as having an O horizon ≥40 cm. In the centre of the map is an area with an extensive fen complex in Kråketorp, Kronoberg County, Sweden. Base map: Ortofoto, 0.5 RGB Tile 634_48_05_2011 © Lantmäteriet.

Table 3 Estimated national extent of peatlands in forested and treeless wetland areas according to the models developed in this study. %_Tot is the percentage of extent compared to the total Swedish land area (407 284 km²; Statistiska Centralbyrån, 2022); %_Forest is the percentage of extent compared to the Swedish forested and treeless wetlands area (NMD2018 Forest and treeless wetland mask; see Sect. 2.2.3).

4 Discussions

This a product-driven study, with the aim to map the extent of peatlands in the forested landscape of Sweden. The performance metrics (Table 3) indicate that the models we developed have good predictive power and that the resulting maps are of adequate quality.

Some errors in the features may have been caused by the spatial inaccuracies of the soil data, especially in the case of the SGU peat archive which is a collection of historical data with varying quality. However, we believe that this has had a minor effect on our results, as the data from the SGU archive mostly include points from the inner parts of peatlands with an O-layer depth >100 cm and a relatively narrow value range in the geodata used as features. Only a minority of the SGU archive data points around Stockholm part from this pattern, meaning that a small portion of the points may lead to any important errors in the features.

It would have been beneficial to have extensive O-layer thickness measurements from all land cover types, in which case a land cover-based mask would not have been needed, and the resulting maps could include all land cover types. We also note that a mask of forest areas based on SFSI would have been more consistent with the utilized soil data than a reclassified version of NMD2018, but it was unavailable at the time of this study.

We hypothesize that the greater importance attributed to Rough in Peat100 (Fig. 2) is due to deeper peat being generally more concentrated in the inner and flat parts of peatlands. The overarching importance of MarkFuktIDX (Fig. 2), we believe, is caused by it being a high-quality descriptor of soil moisture, an essential edaphic factor for peat formation.

Our estimates of national forest peatland extent are within the range of studies employing traditional peatland maps for data merging and upscaling (Barthelmes et al., 2015; Pahkakangas et al., 2016; Tanneberger et al., 2017). These studies were based on maps of Quaternary deposits by SGU and estimated that 62 073–69 155 km² of Sweden is covered by peat; the range is caused by the utilization of different techniques to fill the gaps in the SGU maps and by using (or not) the SGU map with coarsest resolution. These estimates are, however, not entirely comparable to ours, as SGU does not exclude non-forested land as we did. Also, the SGU maps often exclude smaller peatlands by instead reporting the dominant soil type in a given area (Karlsson et al., 2021), probably leading to an overall underestimation of peatland extent.

Assuming that agricultural land on peat is around 2300–2500 km² (Berglund and Berglund, 2010; Lindahl and Lundblad, 2021) and that the area of wetlands in the mountainous areas with no lidar scans is approximately 650 km² (and on Gotland and Öland 150 km² together) (Wetland class in NMD2018, Naturvårdsverket), our national peatland extent estimates are similar to the higher end in the range of the traditional maps (e.g. Barthelmes et al., 2015).

On the other hand, Ågren et al. (2022) produced 2 m resolution peatland maps using modern automated mapping methods. They estimated that 70 000–94 000 km² of Sweden is covered by peat, depending on peatland definition, of which 68 000–88 000 km² is in forests. If defining peat as having an O-layer thickness ≥40 cm, the study estimated 79 000 km² total peatland extent and 76 000 km² forest peatland extent, considerably higher estimates than ours. The authors concluded that the difference in estimates between their study and previous ones was caused by the inclusion of smaller peatlands associated with streams and narrow pits not included in the SGU maps. The difference between our estimates of total peatland extent and those of Ågren et al. (2022) is mainly because they do not exclude agricultural land from their final products, as we did.

On our test dataset, the 40 cm peatland map by Ågren et al. (2022) achieved accuracy, precision, recall and MCC values of 0.93, 0.93, 0.93 and 0.86, respectively. Comparatively, Ågren et al. (2022) employed a simpler model algorithm as us but obtained better performance metrics. We believe that these results can be deduced to a high-quality soil moisture map being the foundation of their model, supported by our results indicating that soil moisture datasets are the most important DEM-derived products for peatland delineation. Furthermore, the considerably higher performance metrics reported in this paper compared to Ågren et al. (2022) are probably due to the inclusion of the SGU peat core archive as input dataset, which is skewed towards easily classifiable areas of deep peatland.

5 Data availability

The peatland maps Peat30, Peat40, Peat50 and Peat100 are available as GeoTIFF files in the Bolin Centre Database (Rimondini et al., 2023; https://doi.org/10.17043/rimondini-2023-peatlands-2).

6 Conclusions

Peatland area, as mapped in this study, is a crucial variable for estimating C stock (Yu, 2012). In turn, such estimates are important for understanding global C cycles and calculating GHG emissions from peatlands (Beaulne et al., 2021; Hugelius et al., 2020). Our peatland maps may, therefore, be used for estimating the C stock of Swedish peatlands, e.g. by coupling them with information on peat thickness and C density. This is especially important in the context of land cover and climate change, which may further contribute to peatland degradation and oxidation (Hugelius et al., 2020; Humpenöder et al., 2020).

In conclusion, our results highlight the possibilities of producing boreal peatland maps using lidar-based terrain indices and spatially and numerically extensive field data. We have demonstrated that an approach such as ours, readily replicable and using only open-source software, can yield accurate and uniform products surmounting the limitations of traditional mapping techniques. Ågren et al. (2022) and Karlson et al. (2019) also concluded that terrain indices are highly valuable for peatland mapping in Sweden; we therefore emphasize that the availability of high-quality lidar and field data is a cornerstone to improve mapping of boreal peatlands. Such advancements could in turn reduce the uncertainties in C stock estimations attributed to poor data on peatland extent, and thus improve estimates of GHG fluxes. Information on the location and depth of peatlands is also crucial for protection and restoration initiatives, making high-quality geospatial products essential for climate change mitigation strategies and decision-making.