A New Land Cover Map of Two Watersheds under Long-Term Environmental Monitoring in the Swedish Arctic Using Random Forest Classification of Sentinel-2 Data

1. Introduction

The nature of the land cover, including vegetation covers, bare rock outcrops and surface water bodies are of major importance to understand hydrological and biogeochemical fluxes on continental surfaces [1,2,3]. It is especially true in the Arctic, where permafrost conditions may strongly control the present ecotypes and their distributions [4,5,6,7], while vegetation cover variability may in turn significantly impact thermo-hydrological conditions [8,9]. In permafrost-affected soils, strong coupling between water and heat transfers occurs, and thus the thermal buffering of the vegetation cover is a key determinants of permafrost dynamics [10,11,12,13,14]. Evapotranspiration fluxes may also be a dominant term of the water budget in permafrost regions [15,16]. For all these reasons, permafrost modeling requires detailed knowledge of up to date land cover distribution.

The vast extension and the remoteness of the Arctic regions make the establishment of field survey – based land cover maps difficult. Moreover, fine resolutions and open data maps are needed for many applications [17], including permafrost modelling. Thus there is a growing interest in airborne [18,19] and remote sensing [20,21,22] capable of producing fine resolution vegetation maps in the Arctics. These regions are experiencing intensive climate change [23]. Permafrost thawing results in methane and carbon dioxide emissions [24] which contributes to the greenhouse effect. These modifications induce changes in ecotypes [25] that are visible at the landscape level. Thus, there is a need for not only high spatial resolution maps, but also for high temporal resolution survey. In order to produce regularly updated land cover map for large areas, the use of remote sensing data from long term satellite missions combined with in-situ information is required [21].

Here we present a workflow for creating high resolution vegetation map using only open data and open source software along with dedicated field data. The workflow is applied to two watersheds in the Swedish Arctic, for parts of which previous vegetation maps at lower resolutions and/or in past climatic conditions where already available [18,26,27]. The obtained map is used for investigating a link between topography and vegetation distribution, and assessing the temporal evolution of the vegetation cover during a 14 years period (2008-2022). This information is crucially important forpermafrost modeling works of the studied sites within a permaFoam approach [16,28], and it will also provide new insights on contemporary landscape evolution in this type of environment.

2. Materials and Methods

2.1. General geographic information about the study area

We selected two watersheds close to Abisko Scientific Research Station (INTERACT Network)(Figure 1). The first one, from West to East, is Miellajokka, a sub-alpine catchment which includes the iconic mounts of Tjuonavagge. This 51.5 km² catchment presents altitudes ranging from 383 to 1731 m above sea level [29]. The most eastern watershed is Stordalen, a 16 km² catchment with a lake-rich, peat-rich Northern part, and a sub-alpine Southern part, with elevation between 350 and 770 m above sea level [30,31,32,33]. In Stordalen vegetation maps of the Northern, low elevation part has been already produced based on airborne data of 2000 [26]. Later on, another vegetation map for the whole watershed has been produced using airborne data of 2008 [27]. Both Stordalen and Miellajokka are encompassed in the area studied by Reese et al. [18], with a vegetation map established on the basis of 2010 satellite images, using also data acquired by a lidar survey.

2.2. Satellite image and digital elevation model

Obtaining images in the Arctic zone to study vegetation cover is difficult. These geographical areas are covered with snow for a large part of the year, which prevents any satellite study of the vegetation cover. In addition, frequent clouds hinder the acquisition of optical images. A single Sentinel-2 image acquired on 25 August 2022 was downloaded from https://peps.cnes.fr. Ten bands were selected for land cover classification (B02 – blue, B03 – green, B04 – red, B05 – red-edge 1, B06 – red-edge 2, B07 - red-edge 3, B08 – NIR, B08A – narrow NIR, B11 – SWIR 1, B12 – SWIR 2). The image is not corrected for atmospheric effects (Level-1C). The images are stored in the UTM34N reference coordinate system and all calculations are performed in this system to avoid altering the radiometry by re-projection.

On the basis of four (B03, B04, B08, B11) out of the 10 acquired channels, four derived indices were calculated: Bright, NDVI, NDWI, NDII (Table 1). The bright index is very sensitive to albedo. It distinguishes between light and dark soils. The NDWI (Normalized Difference Water Index) was used to detect water areas. The NDVI expressed the photosynthesis of the vegetation cover. The use of NDII [34,35] did not improve the results and was not retained for the final classification.

Since vegetation in mountainous areas is related to altitude, the digital terrain model is a very useful data source. ArcticDEM is an NGA-NSF public-private initiative to automatically produce a high-resolution digital surface model of the Arctic using optical stereo imagery. The majority of ArcticDEM data was generated from the panchromatic bands of the WorldView-1, WorldView-2, and WorldView-3 satellites and, for a small percentage of data, from the GeoEye-1 satellite. For this study, ArticDEM Release 7 "mosaic" format files with a spatial resolution of 2 m were downloaded at https://data.pgc.umn.edu/elev/dem/setsm/ArcticDEM/mosaic/v3.0/.

2.3. Field survey

The ground-truth survey took place from 21 July 2022 to 24 July 2022 in the Miellajokka and Stordalen watersheds, northern Sweden. We geolocalized areas of the different land cover types in field using a GPS, GLONASS, Beidou and Galileo navigation systems supported Samsung Galaxy Tab S6 Lite tablet. The Qfield software was used for data entry in the field. Its compatibility with QGIS simplifies data collection and subsequent analysis [28].

A field survey database including a color composite of Sentinel-2 image channels B08/B04/B03, the Open Street Map data and an empty vector layer was prepared in QGIS then transferred to Qfield.

Areas of observed and land cover types were highlighted as polygons on Figure 2. Each surveyed polygon serves a ground truth location. As a complement, a photo of the most characteristic observations was taken with the tablet camera.

The 270 observations conducted during the field survey only identified seven out of the 12 classes by Reese et al. [18]. “Alpine meadow” was not encountered enough to constitute an individual class. Likewise, the “Mountain birch - meadow” class was only observed in six ground truth polygons and was grouped with the “Mountain birch - moss” class to form a single “Mountain birch” class. Snowbeds were poorly represented and are not included. “Grass heaths” were not encountered. Further, “Rock” class mainly represents bed rock outcrops but may also include thin organic soil and sediment. We also added “Human infrastructure” as a new class, mainly representing the road and the railway passing through the mapped area. Shadows in the steep areas to the south of the study area hinder recognition of the landscape they cover. To avoid confusion, especially with water, a “Shadow” class was created, summarizing the total number of classes to nine (Table 2).

As far as possible, we balanced the number of survey polygon between the representative classes. Poor accessibility due to difficult terrain limited the choices of locations (Figure 2). Thus, a randomized field survey design was not possible in this natural environment, the transition from one land-use class to another is sometimes gradual, making it difficult to assign an area to a specific class. For this reason, surveys are only carried out in areas that are homogeneous in terms of the characteristics of the land cover classes.

2.4. Classification

The random forests method is widely used in image classification studies [37,38]. It uses decision trees and random draws without replacement of samples and variables to classify the Sentinel-2 image. Within each class, 30% of the surveyed polygons are randomly drawn and reserved for classification quality assessment. The classification is trained with the remaining 70%. GRASS software is used for the calculations [39]. The extension r.learn.ml2 interfaces with the Scikit-learn library written in python to perform random forests classification. Two parameters were selected in order to produce the forest trees. The number of trees was set to 100. The number of variables to be selected and tested for the best split when growing the trees was set to the square root of the number of input variables.

3. Results

3.1. Vegetation map in current climatic conditions

The classified image (Figure 3) shows patterns consistent with the knowledge of the terrain. The transport infrastructures are described with precision in their continuity. The lake Torneträsk at the North is homogeneously identified.

The statistics computed from surveyed polygons reserved for classification quality assessment confirm the quality of the classification. The percentage of pixels correctly classified by random forests is 88%. The confusion matrix (Table 3) provides an analysis of the accuracy of the classification used for building our map at the class level. The largely dominant weight of the diagonal terms in this confusion matrix demonstrates the quality of the classification. The confusion between “Dry heath” and “Mesic heath” is understandable because these two formations are differentiated primarily by canopy height, a feature not accessible from the images used in our study. Likewise, the confusion between “Alpine willow” and “Wetland” is due to the difficulty of recognizing spaces occupied by a few willow plants. Besides, with such a pixel classification approach, places that are temporarily flooded at the moment of the satellite image acquisition are difficult to distinguish from true wetlands, i.e. places that are under water almost all along the active season. This could overestimate the wetland area, including the places with other vegetation types such as meadow which are temporarily flooded by ground water discharge or snowmelt water. Another important point is the detection of temporary high elevation open water bodies in several places around the Tjuonavagge lake, according to both this classification and the two indices NDII, NDVI values of the pixels. These ones may be generated by late snowmelt in the highest places of the landscape. Finally, the confusion between “Dry heath”, “Mesic heath” and “Mountain birch” are classes that can be contiguous and even associated in some places. It describes mixed spaces where several classes coexist, i.e. ecotone between these classes.

Classification accuracy was accessed by class using producer’s accuracy, user’s accuracy, class area estimates and standard error of class estimate in table 4 [40,41]. Concerning the producer's accuracy criterion, only two classes “Mesic heath” and “Alpine willow” are poorly assessed. These classes represent only 4% of the whole image (table 5). "Alpine willow" is confused with "Wetland", which is understandable given that this plant grows in very wet areas. "Mesic heath" is confused with "Dry heath". These two classes are similar and are differentiated mainly by the height of the vegetation from a remote sensing point of view. The "Dry heath" producer's accuracy criterion is poorly assessed because it is confused with "Mesic heath". This confusion is explained for the reason given above.

All the pixels classified as "Alpine willow" belong to this class. An area classified as “Alpine willow” therefore corresponds to this class. However, an area of land covered by this plant is not always assigned to this class. Thus, the area classified to “Alpine willow” is underestimated. Assuming that the distribution of surface areas follows a normal distribution, the approximate 95% confidence interval is obtained by the formula

Â_{k ∓} 1.96 S(Â_k) (1)

Â_k class area estimate

S(Â_k) standard error of class area estimate

All the class areas are estimated with low standard error. The standard error values are particularly good for classes with large areas: 17 pixels for 1331 pixels for “Mountain birch”, 19 pixels for 2977 pixels for “Water”.

The “Mountain birch” class is the most represented (Table 5). It occupies more than a quarter of the space. The “Rock” class and the wetlands (“Wetland” and “Alpine willow”) share 40% of the spaces. The other classes are much smaller.

3.2. Influence of altitude on land cover

Land cover appears to be strongly conditioned by altitude. Altitudinal zonation of land cover is encountered in various arctic contexts [42,43,44], including in the Abisko region [45,46]. In Miellajokka catchment, the seasonal variability of the hydrogeochemistry of the stream indicates a strong altitudinal control on hydrological processes, especially during the spring freshet [47], while hydrological conditions strongly interact with vegetation and carbon dioxyde fluxes [48]. Table 6, constructed by cross-referencing the land cover map with the ArticDEM, illustrates this phenomenon in the studied watersheds.

The subalpine level at altitudes below 600 m is mainly occupied by the “Mountain birch” class. The Alpine stage includes the “Dry heath”, “Mesic heath”, “Wetland” and “Alpine willow” formations. It is divided into two sub-stages: Between 600 and 800 m of altitude, the “Mountain birch” class is still very present. Above 800 m, this class gives way to the rock. The nival stage is composed only of rock, probably because it depends on harsher life conditions and more intense erosive processes at higher elevations.

3.3. Comparison between past and present vegetation maps

Three maps of parts of our study area have been produced by different authors. A map was constructed from aerial images of 8 August 1970 and 29 July 2000 [26], but the area covered is too small to allow comparison with our data. Another map of the Stordalen watershed by Lundin et al. [31] was obtained from images obtained from a helicopter flight on 1 August 2008. The most recent map of the Miellajokka watershed was produced by Reese et al. [18,49] from SPOT5 images of 28 July 2011 and laser data acquired under leaf-on conditions from two scanning dates (20 August 2010 and 9 September 2010). As the semantics of these maps are not identical to ours, an analysis of the variable typologies is necessary prior to the study of the landscape evolution. In order to overlay the maps and then calculate statistics, the Lundin and Reese maps extracted from the publications were georeferenced from control points identified in the landscape.

Comparison with the map of Reese et al. (based on data acquired in 2010)

Reese et al. [18,49] produced a land cover map with a larger number of classes. In order to compare the Reese map with our data, some classes of the Reese map are merged. “Snow ice” and “Snow bed” classes are grouped together. “Dry heath”, “Extremely dry heath”, “Grass heath” are also grouped together. The “Human infrastructure” class is not considered because it does not exist in Reese's work. "Alpine meadow” and “Tall Alpine meadow" were not confirmed by ground observations during our field trip.

The spatial distribution of the classes is slightly different (Table 7). The “Rock” class accounts for 36% for our classification and only 14% for the Reese map.

The change matrix shown in Table 8 encompasses the differences in semantics of classes and the changes in the landscape between the dates of the two maps (i.e.: 2014 and 2022), and may be also discrepancies due to the use of different methodologies. Table 8 and Figure 3 show that some “Rock” areas in the middle our map are covered by “Grass heath”, “Dry heath” and “Alpine willow” on the Reese map. The forest is also growing slightly to the south in sparse patches.

On the other hand, “Alpine willow” is also more represented on Reese map without any conclusion being drawn because this class is misclassified by Reese [18]: the confusion matrix indicates 20 of 44 pixels misclassified. In addition, “Alpine meadow” which we have not taken into account, is identified as “Wetland” on our map (Figure 4), maybe because most of “Alpine meadow” places were temporarily flooded at the time of observation(see also section 3.1.). This class “Wetland” exists on the map published by Borgelt [49] but it is not present in the confusion matrix published by Reese [18]. The other classes do not show significant difference in in their proportion and spatial distribution.

Despite the differences between the two maps, figure 3 shows a great consistency in the distribution of space. Three elements can explain these differences. 1) The landscape is natural. There are no parcels to structure it. Between areas occupied by two vegetation classes there is often a transition zone which is difficult to assign to a class. 2) The class definition of Reese [18] takes into account the height of the stratum using metrics derived from laser acquisitions, a technology we did not employ. 3) The landscape has evolved between 2010 and 2022. For example, some areas of rock have become vegetated.

Comparison with the map of Lundin et al. (based on data acquired in 2008)

The definition of classes in Lundin et al. [31] is significantly different from the one presented in this work, with a smaller number of classes in the map of Lundin. So the classes we used have to be modified to achieve a consistency between two maps. Grouping “Dry heath”, “Mesic heath”, “Alpine willow” classes of our classification allows them to be compared to the “Alpine tundra” class of the map of Lundin. Similarly, our “Wetland” class is compared to the “Peatland” class of the map of Lundin. The classes “Human infrastructure” and “Non vegetated” correspond to road, railway and building. The latter class is much more represented on the map of Lundin (Table 9). The differences are related to a larger road and railway footprint, which does not affect the landscape dynamics. This observation shows the satisfactory overlay of the two maps.

The change matrix (Table 10) shows that the "Rock" class is also more represented on the map of Lundin, which covers twice as much space. The areas of this class that are not classified as “Rock” for our map are located in the south of the map (Figure 5). They are contiguous to the “Rock” areas of our map. It is also possible that there has been a forest expansion between 2008 and 2002 in this area. Indeed, the forest has grown between 2008 and 2022. It has gained some space in all land cover categories. Furthermore, some areas of “Peatland” appear to be transformed into “Wetland” but Lundin et al. [31] indicate "Peatlands were subdivided into wet areas (fen) and dry areas (bog) proportionally to what was found by Malmer et al. (2005)". It is therefore not possible to draw a conclusion from this observation. In summary, it is possible to compare our map with the map of Lundin after a semantic analysis of the categories. The main differences between the two maps concerns the south of the Stordalen watershed where Alpine toundra and forest are intermixed. The forest seems to have taken over areas previously occupied by tundra.

4. Discussion

The methodology implemented in this study requires little sampling effort to quickly obtain a land cover map. This operation is feasible every year to monitor the evolution of land cover. These results are particularly important in a region subject to climate change which is currently undergoing major upheaval.

The results of the Sentinel-2 image classification show a structuring of the landscape as a function of elevation, with different vegetation level along with altitude [50]. In the Miellajokka and Stordalen catchments, three levels are present. The subalpine level (< 600 m) is mainly occupied by birch forest. The alpine level [600 m , 1100 m] is characterized by heath and willow. This level could be split into two sub-levels according to the respective abundance of birch forest at lower altitude and outcrops at upper altitude. The upper level, the nival level (> 1100 m), is only composed of rock. Some temporary open water bodies were localised at high altitude, which is a surprising feature, may be linked to late melt of high elevation snow bodies.

Comparison between past and present vegetation maps is not straightforward due to a lack of a common typology field survey protocols. Nevertheless, it was possible to identify a change in the landscape between 2008 and 2022. Analysis comparison of the maps of Lundin [31] and Reese [18] with the one presented in this work demonstrated an extension of the forest on the tundra towards the south (i.e., toward higher elevations) during the 2008-2022 period. This finding is in agreement with Rundqvist's work [51] which shows an upward movement of species observed over a study period between 1976 and 2010. Nevertheless one should be careful with this possible interpretation because of the statistical uncertainty of the different classifications. This trend could be a consequence of the on-going climate warming, demonstrated across the Arctic [52,53].

5. Conclusion

In this work we provided a new land cover map for two watersheds located nearby the Abisko Scientific Research Station, to be used in future permafrost modelling work in the framework of the HiPerBorea project (hiperborea.omp.eu) based on permaFoam approach [16,28]. This new map also provides some insights on recent land cover changes in the studied area, by comparison with previous maps based on data acquired in 2008 [31] and 2010 [18]. High elevation temporary water bodies have been detected, which requires further investigation in the future.

The proximity of the Abisko observation station makes the Miellajokka and Stordalen watersheds a privileged study area for the evolution of landscapes in the Arctic zone, in particular, where thawing of the permafrost at high altitudes is attested. This monitoring requires annual surveys according to a unified protocol in terms of sampling, definition of classes and method of recording in order to monitor the evolution of this region under on-going climate change. The present study presents a protocol that would be suitable for such purpose.

At the same time, the study of land cover in the Arctic zone poses a number of difficulties. The areas are covered by snow for a large part of the year, which limits observations to a few months of summer. Acquisition by passive optical satellites can only be made during the six months of perpetual daylight, from 21 June to 21 December. Fortunately, this period includes the summer months when snow cover is at its minimum extension. Cloud cover frequently hampers optical acquisitions. For the year 2022, only one Sentinel-2 image could be used, which illustrates the problems to rely on optical images only in such environments. Future work could involve radar images whose acquisition is not affected by clouds and polar night.