A 10 m resolution land cover map of the Tibetan Plateau with detailed vegetation types

The Tibetan Plateau (TP) hosts a variety of vegetation types, ranging from broadleaved and needle-leaved forests at the lower altitudes and in mesic areas to alpine grassland at the higher altitudes and in xeric areas. Accurate and detailed mapping of the vegetation distribution on the TP is essential for an improved understanding of climate change effects on terrestrial ecosystems. Yet, existing land cover datasets for the TP are either provided at a low spatial resolution or have insufficient vegetation types to characterize certain unique TP ecosystems, such as the alpine scree. Here, we produced a 10 m resolution TP land cover map with 12 vegetation classes and 3 non-vegetation classes for the year 2022 (referred to as TP_LC10-2022) by leveraging state-of-the-art remote-sensing approaches including Sentinel-1 and Sentinel-2 imagery, environmental and topographic datasets, and four machine learning models using the Google Earth Engine platform. Our TP_LC10-2022 dataset achieved an overall classification accuracy of 86.5 % with a kappa coefficient of 0.854. Upon comparing it with four existing global land cover products, TP_LC10-2022 showed significant improvements in terms of reflecting local-scale vertical variations in the southeast TP region. Moreover, we found that alpine scree, which is ignored in existing land cover datasets, occupied 13.99 % of the TP region, and shrublands, which are characterized by distinct forms (deciduous shrublands and evergreen shrublands) that are largely determined by the topography and are missed in existing land cover datasets, occupied 4.63 % of the TP region. Our dataset provides a solid foundation for further analyses which need accurate delineation of these unique vegetation types in the TP. TP_LC10-2022 and the sample dataset are freely available at https://doi.org/10.5281/zenodo.8214981 (Huang et al., 2023a) and https://doi.org/10.5281/zenodo.8227942 (Huang et al., 2023b), respectively. Additionally, the classification map can be viewed at https://cold-classifier.users.earthengine.app/view/tplc10-2022 (last access: 6 June 2024).

1 Introduction

The Earth's surface is physically covered by various types of land cover, including forests, grasslands, croplands, lakes, wetlands, etc. The accurate mapping and classification of land cover are fundamental components of Earth observations. By understanding the distribution and characteristics of different land cover types, land cover mapping supports the assessment of carbon stocks, vegetation dynamics, and land–atmosphere interactions, contributing to the implementation of effective climate change mitigation measures (Z. Wang et al., 2022; Liu et al., 2022; Li et al., 2018).

The advent of remote-sensing technology has enabled the generation of global-scale land cover products at various resolutions. For instance, products like MCD12Q1, produced using MODIS data (Friedl et al., 2010, 2002), and the ESA CCI product derived from sensors like MERIS (Agency, 2014) have significantly contributed to the understanding of global ecosystem responses to climate change. However, their spatial resolutions are hundreds of meters, so they are unable to provide an accurate representation of the land surface conditions (Tian et al., 2021), particularly in spatially heterogeneous regions, such as the mountainous southeast Tibetan Plateau (TP) (Yang et al., 2017; Grekousis et al., 2015). In response to this limitation, several medium- to high-resolution land cover products have been created using satellite images from Landsat and Sentinel-2. Notable examples include GlobeLand30 (Chen et al., 2021, 2015), FROM_GLC30 (Gong et al., 2013), and GLC_FCS30 (X. Zhang et al., 2021) based on Landsat and FROM_GLC10 (Chen et al., 2019), Dynamic World (Brown et al., 2022), Esri Land Cover (Karra et al., 2021), and ESA WorldCover (Zanaga et al., 2022) based on Sentinel-2. However, these products use different classification systems, resulting in large divergences in certain regions (Shi et al., 2023; Hua et al., 2018), and they are often inadequate to reflect the diverse and unique land cover types in important ecosystems (C. Liu et al., 2023), such as those in the TP.

Renowned as the “Third Pole” of the world (Shukla and Sen, 2021), the TP holds a dual significance as a sensitive area and an indicator zone for global climate change (Hua et al., 2021; Li et al., 2022; Trew and Maclean, 2021; Pepin et al., 2022). It hosts a variety of vegetation types, ranging from broadleaved and needle-leaved forests at the lower altitudes and in mesic areas to alpine grassland at the higher altitudes and in xeric areas. However, many of the unique vegetation types in the TP, such as the alpine scree ecosystem in the transitional zone from alpine grasslands to bare rocks at very high altitudes and the shrubland ecosystem in the transitional zone from forests to grasslands, are not well represented in existing land cover datasets (Li et al., 2014). Furthermore, shrublands in the TP can have either evergreen leaves or deciduous leaves, depending on the local environments in which they grow, yet they are largely ignored in existing 10 m resolution land cover datasets (Venter et al., 2022). It is very important to monitor these unique ecosystems in the TP, given that the TP has experienced dramatic warming (Fu et al., 2021), increased humidity (Yang et al., 2014), rapid glacier retreat (F. Zhao et al., 2022), permafrost thawing (Gao et al., 2021), expansion of lakes (Zhang et al., 2020), and vegetation changes (Wang et al., 2020; Duan et al., 2021; Gao et al., 2014) in recent decades. Thus, detailed and accurate mapping of the diverse vegetation types in the TP is required for understanding climate change effects on the terrestrial ecosystem, but this is challenging to accomplish, given that shrublands are often confused with forests or alpine meadows and alpine grasslands are commonly misclassified as bare land in most products (Liu et al., 2021; Cai et al., 2022; Yu et al., 2014). Moreover, the extremely rough terrain in the TP results in large mountain shadows and variations in slope aspects, which complicates the accurate detection of vegetation types from satellite imagery (Pizarro et al., 2022; X. Wang et al., 2023).

To address the aforementioned challenges, we developed a specific vegetation remote-sensing fine-classification system tailored for the TP and consisting of 12 vegetation classes and 3 non-vegetation classes. We then created a comprehensive training and validation dataset consisting of 10 242 samples through manual interpretation and field trips; based on this dataset, we performed land cover classification of the TP by integrating multiple data sources in the Google Earth Engine (GEE) platform, including satellite imagery from Sentinel-1 and Sentinel-2, topography, temperature, and precipitation. We investigated the performance of four different classification models provided in GEE and selected the highest-accuracy one to generate a 10 m resolution land cover product for the TP in 2022, referred to as TP_LC10-2022.

Figure 1 Overview of the study area colored by elevation. The black lines are the field trip routes along national roads. The photos show examples of landscape views of typical vegetation types in the Tibetan Plateau.

2 Study area and data

2.1 Study area

The TP spans from the northern foot of the West Kunlun Mountains and Qilian Mountains to the southern foot of the Himalayas and other mountain ranges and from the western edge of the Kunlun Mountains and Pamir Plateau to the eastern edge of the Hengduan Mountains (Fig. 1). It lies between latitudes 25°59^′30^′′ and 40°1^′0^′′ N and longitudes 67°40^′37^′′ and 104°40^′57^′′ E, covering a total area of 3.083 million km². Its average elevation is approximately 4320 m (Zhang et al., 2022).

Due to the combined influence of climate, topography, and human activities over time, the vegetation cover types in the TP vary significantly at different altitudes. The northwestern and central regions are characterized by extensive bare lands, alpine screes, and persistent snow cover. In the southern and eastern areas, there is a distribution of evergreen forests and mixed forests consisting of needle-leaved and broadleaved trees. The transitional zone between these regions is characterized by shrublands, alpine grasslands, and alpine meadows. We investigated the vegetation cover in a field trip carried out along national road nos. 318 and 109 in July 2023 (Fig. 1), which covered all the vegetation types in the TP.

2.2 Data

2.2.1 Satellite imagery

We used both the optical imagery from Copernicus Sentinel-2 and radar imagery from Copernicus Sentinel-1 for the classification. Sentinel-2 comprises two high-resolution multispectral imaging satellites, each equipped with a multispectral imager. It has 13 bands, with spatial resolutions of 10 m for four bands, 20 m for six bands, and 60 m for three bands. The study utilized level-2A products from the year 2022, which had undergone processing via the Sen2Cor tool at the Copernicus Scientific Data Hub (Doxani et al., 2018). Annual remote-sensing images have proven to accurately capture phenological changes in specific vegetation cover and have been successfully utilized in various large-scale land cover classification studies (Verde et al., 2020). Hence, in this study, the Sentinel-2 remote-sensing images from the entire year of 2022 were selected for band feature extraction. In this study, the initial step involved retaining the images with a cloud cover of less than 10 %. Subsequently, the quality assessment information (QA band) was utilized to exclude pixels with inadequate quality through cloud masking.

Sentinel-1 comprises two polar-orbiting satellites positioned in the same orbital plane. For this research, the ground range detected (GRD) data obtained in wide-swath (IW) mode were chosen. The GRD data consist of single-polarization (VV) and dual-polarization (VV and VH) interferometric wave modes offering 10 m resolution (Prats-Iraola et al., 2015). It enables the provision of radar images suitable for land and maritime services, regardless of weather conditions and time of day. The median-compositing method in GEE (Souza et al., 2020; Phan et al., 2020) was applied to process all bands of Sentinel-1 and Sentinel-2.

2.2.2 Topography data

The Shuttle Radar Topography Mission (SRTM) (Farr et al., 2000) was designed to generate high-quality digital elevation models (DEMs) globally using synthetic-aperture radar technology. The data collected by SRTM were used to create a global elevation model with a horizontal accuracy of 16 m and a vertical accuracy of 6 m at a spatial resolution of 30 m (Yang et al., 2011).

2.2.3 Precipitation data

The Climate Hazards Group InfraRed Precipitation with Station data (CHIRPS) (Funk et al., 2015) is a comprehensive dataset documenting global precipitation from 1981 to the present. CHIRPS integrates satellite imagery with in situ station data, allowing the generation of gridded rainfall time series suitable for trend analysis and seasonal drought monitoring at a resolution of 0.05°.

2.2.4 Temperature data

The ERA5-Land dataset (Muñoz-Sabater et al., 2021) offers a comprehensive reanalysis of land variables and presents a consistent perspective on their evolution over multiple decades at a higher resolution than ERA5. As the land component of the ECMWF ERA5 climate reanalysis, ERA5-Land combines model data and global observations to create a coherent dataset utilizing the principles of physics. Nineteen extra bands were incorporated by GEE, with each corresponding to an accumulation band, and the hourly values were calculated as the difference between two successive forecast steps (Muñoz-Sabater, 2019). For this study, hourly temperature data with a resolution of 0.1° from 2022 were used.

3 Methodology

3.1 Land cover classification

The advancement of cloud computing technology in remote sensing has revolutionized the rapid analysis and application of Earth system science on a large scale, even globally. GEE stands out among these technologies, offering online visualization, computation, and analysis capabilities for extensive Earth science data (Gorelick et al., 2017; Kumar and Mutanga, 2018). Consequently, we opted to utilize GEE for data processing and analysis. Importantly, the satellite data and auxiliary data relevant to this study can be readily accessed through GEE. Figure 2 presents our comprehensive classification system, which comprises four main steps: (1) sampling strategy, (2) data preprocessing and feature construction, (3) classification model comparison, and (4) accuracy assessment and intercomparison.

Figure 2 Flowchart of the land cover classification carried out in this study.

3.1.1 Classification system

The TP harbors the world's highest and one of the most distinctive alpine vegetation communities, posing challenges to its inclusion in both global and Chinese land cover classification systems. To address this issue, we have developed an adapted classification system specifically tailored to the alpine vegetation types found in the TP. The basis for constructing this classification system is as follows:

1. Comprehensive vegetation functional types. We have categorized the vegetation in the TP based on plant growth form (trees, shrubs, and herbs), leaf phenology (evergreen and deciduous), leaf type (broadleaved and needle-leaved), and ecosystem type. This classification system results in 12 vegetation types, including 5 types of tree cover (including evergreen needle-leaved forest (ENF), deciduous needle-leaved forest (DNF), evergreen broadleaved forest (EBF), deciduous broadleaved forest (DBF), and mixed forest (MF)); 2 types of shrub cover (including evergreen shrubland (ES) and deciduous shrubland (DS)); 2 types of herb cover (including alpine grassland (AG) and alpine meadow (AM)); 3 special vegetation cover types (including alpine scree (AS), wetland (WL), and cultivated vegetation (CV)); and 3 non-vegetation land cover types (including bare land (BL), water body (WB), and permanent ice and snow (PIS)).

2. Discriminability of different vegetation functional types in remote-sensing imagery. During the classification stage, we can effectively differentiate various land cover types, including diverse vegetation, utilizing the discriminative capabilities of the multispectral bands of Sentinel-2 (X. Liu et al., 2023). Moreover, the incorporation of high-resolution Google Earth imagery, with a spatial resolution of up to 0.3 m, enhances the distinguishability of land cover types during the sample selection phase. This ensures the feasibility of visually interpreting large-scale samples from remote-sensing imagery and obtaining reliable and up-to-date information (Gong et al., 2013).

In this study, we did not specifically select samples of built-up areas and instead categorized bare land together with built-up areas for two primary reasons. Firstly, built-up areas account for only 0.092 % of the total area in ESA WorldCover2021, highlighting their relatively small extent compared to other land cover types (Zanaga et al., 2022). Secondly, the bare land in our product exhibits spectral characteristics similar to those of built-up areas, resulting in the classification of most built-up areas as bare land (Li et al., 2017).

Figure 3 Examples of the auxiliary data used for explaining the visual interpretation, including Google Earth imagery and the Landsat monthly mean normalized-difference vegetation index (NDVI) time series for 2013–2022. The x axis represents month of the year, and the y axis represents NDVI value.

3.1.2 Sampling strategy

Supervised classification models heavily depend on a substantial number of labeled samples for effective training and validation (Foody and Mathur, 2004). While extracting samples directly from existing land cover products can save personnel, it introduces several issues: (1) extracted training samples may inherit errors from previous land cover products (Xi et al., 2022); (2) utilizing low-resolution products to extract training samples for high-resolution land cover mapping can lead to information loss and boundary effects between adjacent land parcels (X. Zhang et al., 2021; Zhang and Roy, 2017); and (3) reconciling the classification systems of different products is difficult, and global land cover products may not include specific land cover types for certain regions. Therefore, collecting samples through visual interpretation emerges as a more feasible approach (Schepaschenko et al., 2019).

Google Earth integrates high-resolution imagery from sources like QuickBird and GeoEye, providing reliable remote-sensing data sources for visual interpretation. Selecting samples in areas without Google Earth image coverage in 2022 poses a challenge. Normalized-difference vegetation index (NDVI) time series have thus been used as auxiliary data for land cover sample selection (Yang and Huang, 2021; Feng et al., 2016). To ensure the selection of stable samples, this study examines the stability of land features by reviewing the Landsat NDVI time series from 2013 to 2022. To eliminate the interference from clouds and snow in the NDVI time series, the following operations were performed on Landsat images: (1) pixels with cloud coverages of greater than 50 % were filtered out and (2) a normalized-difference snow index (NDSI) mask was applied to filter out pixels with NDSI values greater than −0.4 when selecting forest and shrub samples. To obtain a more continuous NDVI time series, the harmonic analysis of time series (HANTS) model was used for data interpolation and smoothing to remove noise and reconstruct missing data (Zhou et al., 2015). By following the steps outlined above, we detected land cover changes during 2013–2022 using Landsat NDVI time series (Fig. A2). This approach helps to avoid selecting sites where land cover change has occurred. Additionally, the monthly mean value of the NDVI time series for 2013–2022 was calculated to determine the phenological characteristics of each sample point (Chu et al., 2021). All samples were interpreted based on Google Earth images, with subsequent verification performed using NDVI time series as a supplementary measure to ensure stability and detect the phenology.

Table 1 Number of training and validation samples for the 15 land cover types.

For instance, in Fig. 3, different color characteristics are observed for evergreen shrubs and deciduous shrubs in Google Earth imagery. Evergreen shrubs maintain their green color even during winter, while deciduous shrubs appear yellow-brown. However, during spring or summer, direct differentiation between the two from imagery is not possible. Therefore, phenological characteristics are extracted from their mean NDVI time series. Evergreen shrubs exhibit relatively stable NDVI values, whereas deciduous shrubs show a decrease in NDVI due to seasonal leaf shedding. Evergreen needle-leaved forests, evergreen broadleaved forests, and evergreen shrublands exhibit similar trends and values in NDVI time series. However, they can be discerned in Google Earth images based on their distinctive crown shapes and textures (Fig. 3).

Figure 4 Spatial distribution of the 10 242 samples for land cover classification in the Tibetan Plateau.

Google Earth imagery does not accurately determine the presence of herbaceous plant growth. Nevertheless, grasslands display a significant periodic increase in NDVI during the growing season, while bare land exhibits a relatively flat NDVI time series. This characteristic is utilized for identifying bare land. Regarding alpine grasslands and alpine meadows, judgments are based on area size, vegetation composition, moisture conditions, and terrain. Meadows typically have a smaller area compared to grasslands and better moisture conditions, and they are often accompanied by trees or shrubs in the vicinity. Grasslands have a flatter distribution area compared to meadows, as depicted in Fig. 3. Consequently, effective differentiation between alpine grasslands and alpine meadows is achieved.

Topography data (elevation, slope, and aspect) (Farr et al., 2000), a 1:1 million Chinese vegetation map (Su et al., 2020), and high-quality Google Maps photos were selected for auxiliary judgment. Ultimately, a total of 10 242 samples were collected, as illustrated in Fig. 4. Subsequently, the 10 242 samples were mixed, and the samples for each category were randomly divided into training and validation sets in an approximate 4:1 ratio, as presented in Table 1. We adjusted the ratio of training to validation samples to 4:1 instead of the commonly used 7:3 to enhance the model's fitting capability so that it could handle the complex distribution of features (Ramezan et al., 2021).

Table 2 Features used for land cover classification.

3.1.3 Feature construction for classification

The selected input bands for Sentinel-2 included B2–B8, B8A, B9, B11, and B12. Among these bands, B2–B8, B11, and B12 have been demonstrated to be effective in classifying deciduous and coniferous tree species (Immitzer et al., 2016; C. Li et al., 2021). Additionally, B8A is suitable for boreal landscape classification (Abdi, 2020), while B9 values demonstrate differences between bare soil and vegetation-covered areas (Y. Zhao et al., 2023), making them useful for classification purposes. For Sentinel-1 images, utilizing both VV and VH can enhance classification accuracy, leading to their selection as input features (Jacob et al., 2020; Steinhausen et al., 2018).

To better discern the characteristics of various land features, we calculated several indices using Sentinel-2 imagery. These included the NDVI, NDSI, normalized-difference water index (NDWI), and optimized soil-adjusted vegetation index (OSAVI). NDVI is highly sensitive to vegetation growth and is commonly used to distinguish between vegetated and non-vegetated areas (Rouse et al., 1974). NDSI effectively detects snow by utilizing the reflective properties of snow in the short infrared band, making it advantageous for studying ice and snow coverage in high mountain regions (Dozier, 1989). NDWI effectively distinguishes between water and non-water features (Xu, 2006). OSAVI improves the sensitivity and stability of vegetation indices by considering the influence of soil reflectance, providing a more accurate reflection of vegetation coverage and growth conditions, particularly in cases of bare soil or sparse vegetation (Rondeaux et al., 1996).

The topography significantly influences the vertical distribution of vegetation in high mountain areas (Zou et al., 2023). Therefore, in this study, we included elevation, slope, and aspect as input features for classification. Additionally, we incorporated annual precipitation and mean annual temperature as classification feature indicators (F. Wang et al., 2023; Shen et al., 2015). For bands with a spatial resolution different from 10 m, we employed bicubic interpolation to resample them to 10 m resolution for mapping (Liu et al., 2020). All the features and their detailed descriptions are presented in Table 2.

3.1.4 Classification model comparison

Machine learning is a technique typically employed in remote-sensing image classification. To identify the most appropriate classification model, we compared four widely-used machine learning models in GEE, including random forest (RF) (Breiman, 2001), gradient tree boosting (GTB) (Friedman, 2001), support vector machine (SVM) (Hearst et al., 1998), and minimum distance (MD) (Wacker and Landgrebe, 1972). We fine-tuned the parameters of all the classification models to achieve optimal results (Table A1). The classification model with the highest overall performance was chosen to generate the land cover map and calculate the area proportion of each land cover type.

3.1.5 Accuracy assessment and intercomparison

The accuracy of remote-sensing image classification is commonly assessed using a confusion matrix, which provides four quantitative indicators: producer's accuracy (PA), for measuring omission errors; user's accuracy (UA), for measuring commission errors; overall accuracy (OA); and the kappa coefficient.

To compare with four existing global land cover datasets, namely ESA WorldCover2021, FROM_GLC10-2017, FROM_GLC30-2015, and GLC_FCS30-2020, we merged pixels belonging to the same class (Table A2) and employed randomly sampled validation samples. Additionally, we selected three 0.1° × 0.1° grids within the TP to compare the visual classification results of TP_LC10-2022 with the existing four land cover products.

4 Results and discussion

4.1 Comparison of classification models

Table 3 presents the evaluation results from different classification models applied to the study area using GEE. The results demonstrate that the RF model achieved the highest accuracy, with an overall accuracy (OA) of 86.5 % and a kappa coefficient of 0.854. The gradient-tree-boosting (GTB) model closely followed, with an OA of 85.6 % and a kappa coefficient of 0.844. The minimum-distance (MD) model yielded an OA of 79.7 % and a kappa of 0.781, while the support vector machine (SVM) exhibited significantly lower classification results, with an OA of 64.7 % and a kappa of 0.618.

Table 3 Comparison of classification results from random forest (RF), gradient tree boosting (GTB), minimum distance (MD), and support vector machine (SVM) at their best performance. A bold font is used to indicate the highest UA and PA for each land cover type as well as the highest OA and kappa among the four models.

BL: bare land; AS: alpine scree; AG: alpine grassland; AM: alpine meadow; ENF: evergreen needle-leaved forest; DNF: deciduous needle-leaved forest; EBF: evergreen broadleaved forest; DBF: deciduous broadleaved forest; MF: mixed forest; ES: evergreen shrubland; DS: deciduous shrubland; WB: water body; WL: wetland; CV: cultivated vegetation; PIS: permanent ice and snow.

The high accuracy achieved by RF and GTB models can be attributed to their ensemble-learning algorithms based on decision trees. These algorithms combine multiple decision trees to enhance model performance and generalization capabilities (Salditt et al., 2022). In contrast to the findings of Abdi (2020), where RF and SVM exhibited similar OAs, our SVM showed a decline of 21.8 % compared to RF (Tu et al., 2020). This discrepancy may be attributed to RF's ability to mitigate the correlation between samples and features through random sampling and feature selection, resulting in improved classification performance and robustness. Moreover, RF can effectively handle high-dimensional data and capture nonlinear relationships by integrating multiple decision trees (Tu et al., 2020; Gislason et al., 2006).

The other three classification models aside from SVM effectively distinguished water bodies from other land cover types, achieving a PA exceeding 0.99. However, all classification models performed fairly well in differentiating mixed forests. For instance, SVM achieved a low PA of only 0.269 for mixed-forest classification. Despite the integration of various machine learning models within GEE, including algorithms like RF, a distinct absence of direct support for deep learning persists. This is notable even in light of the well-established and showcased capabilities of deep learning in the fine-grained classification of land cover (Y. Wang et al., 2023). This limitation, to a certain extent, poses a hindrance to the extensive application of large-scale land cover mapping. The utilization of multisource remote-sensing data can offer a more comprehensive understanding of land cover (Xu et al., 2022; Chen et al., 2017). Given the importance of features, all features contributed to the mapping, and elevation contributed slightly more to the accuracy of the classification (Fig. A1). This is attributed to the impact of the TP's rugged terrain on the hydrothermal conditions in distinct regions, leading to notable variations in vegetation phenology (Hwang et al., 2011; Sang et al., 2024).

Figure 5 Overall comparison between TP_LC10-2022 and other land cover maps. (a) TP_LC10-2022 and the proportions of 15 land cover types. (b) An overview of four land cover products for the Tibetan Plateau, including ESA WorldCover2021, FROM_GLC10-2017, FROM_GLC30-2015, and GLC_FCS30-2020. Legend fusion rules for WorldCover2021 and FROM_GLC10-2017 are provided in Table A2, and those for FROM_GLC30-2015 and GLC_FCS30-2020 are given in Table A3.

4.2 Land cover classification map

Figure 5a provides an overview of the TP_LC10-2022 product and four global land cover products, along with the proportion of each land cover type in TP_LC10-2022. Alpine meadow and alpine grassland account for 23.76 % and 16.48 %, respectively. Alpine scree surprisingly ranks fourth, with 13.99 %, after alpine meadow, bare land, and alpine grassland. Evergreen needle-leaved forest has the largest area among the forest types, and deciduous shrubland has a larger area than evergreen shrubland, reaching 3.57 %, which surpasses all other forest types except for evergreen needle-leaved forest. Table A4 presents the statistical area results of 5 land cover products in the TP, highlighting significant discrepancies among them.

Figure 6 The confusion proportions for each of the land cover types in TP_LC10-2022. BL: bare land; AS: alpine scree; AG: alpine grassland; AM: alpine meadow; ENF: evergreen needle-leaved forest; DNF: deciduous needle-leaved forest; EBF: evergreen broadleaved forest; DBF: deciduous broadleaved forest; MF: mixed forest; ES: evergreen shrubland; DS: deciduous shrubland; WB: water body; WL: wetland; CV: cultivated vegetation; PIS: permanent ice and snow.

Table 4 Confusion matrix of the TP_LC10-2022 product extracted using the random-forest (RF) classification model. A bold font denotes correctly classified sample points.

BL: bare land; AS: alpine scree; AG: alpine grassland; AM: alpine meadow; ENF: evergreen needle-leaved forest; DNF: deciduous needle-leaved forest; EBF: evergreen broadleaved forest; DBF: deciduous broadleaved forest; MF: mixed forest; ES: evergreen shrubland; DS: deciduous shrubland; WB: water body; WL: wetland; CV: cultivated vegetation; PIS: permanent ice and snow.

According to Fig. 5b, the ESA WorldCover2021, FROM_GLC10-2017, and FROM_GLC30-2015 products overestimate the area of bare land in the TP, which is similar to the issues observed in the FROM_GLC-agg and ESA CCI land cover products (Liu et al., 2021; Yu et al., 2014). This may be due to the misclassification of alpine grassland as bare land because these products captured less spectral information during the growing season of alpine grasslands. GLC_FCS30-2020 exhibits the highest consistency with TP_LC10-2022 regarding bare land (Table A4 and Fig. 5), and it classified more grasslands while failing to differentiate between grasslands and meadows. Additionally, GLC_FCS30-2020 assigned 61.44 % of the total TP area to grassland, indicating an overestimation of grassland extent (Table A4).

The TP exhibits significant variations in annual rainfall and land surface temperature across its diverse regions, resulting in distinct hot and cold spots (Rao et al., 2019; Wu et al., 2019). Consequently, leveraging climate data can prove beneficial when categorizing alpine meadows in the southeastern TP and alpine grasslands in the northwestern TP at regional climatic scales, given their high sensitivity to changes in annual precipitation and land surface temperature (Su et al., 2020; Wang et al., 2021). Our study also found that incorporating resampled coarse-resolution climate data can help improve the classification accuracy of finer-resolution data (Jia et al., 2014). However, this may cause potential loss of spatial information (Xu et al., 2020), which has not been observed in the TP_LC10-2022 dataset. Table 4 illustrates the confusion matrix of TP_LC10-2022, with an overall accuracy of 86.5 % and a kappa coefficient of 0.854. Water body achieved a PA of 100 %, while mixed forest only reached 46.2 %. Figure 6 shows that most mixed forests are challenging to differentiate from other forest types, with over 25 % of mixed forests misclassified as evergreen needle-leaved forests and 11.5 % misclassified as evergreen broadleaved forests or deciduous broadleaved forests. The classification accuracy of wetlands is also unsatisfactory, with a PA of only 61.1 %. Over 16 % of wetlands were classified as bare land, and over 15 % were incorrectly classified as alpine meadows. The UA for the water body reached 98.3 %, while evergreen needle-leaved forests had the lowest UA at 76.4 %.

In addition, the spectral variations within urban areas have also resulted in substantial uncertainties. Our approach of categorizing built-up areas and bare land may lead to the misclassification of urban pixels. To minimize the uncertainties in urban areas on our final map, we applied the ESRI land cover map for 2022 to mask off urban pixels (Karra et al., 2021).

Although we employed the Sentinel-2 median composition method in this study, we acknowledge the potential enhancement that time-series analysis could bring to our research. In comparison to median composition, time-series analysis has the potential to more comprehensively capture phenological information on vegetation, thereby yielding more accurate land cover classification results (Xie et al., 2019; Nguyen et al., 2020). However, time-series methods also have their limitations, such as the requirement for a greater number of valid observations (Hemmerling et al., 2021). For example, for the summer of 2022 (June–August), when we set the “CLOUDY_PIXEL_PERCENTAGE” parameter to 10 %, 20 %, 30 %, and 40 % and applied QA band masking, we lost 13.59 %, 5.81 %, 2.44 %, and 1.32 % of the Sentinel-2 image area in the TP. The removed pixels were concentrated mainly in the cloudy southeastern TP (shown for a 10 % threshold in Fig. A3) (Tang et al., 2022). This constraint can preclude the attainment of desired outcomes in regions where cloud-free image availability is low (Chu et al., 2021; Coluzzi et al., 2018).

The blue, red-edge, and shortwave-infrared (SWIR) bands of mono-temporal median Sentinel-2 imagery have proven effective for vegetation classification, distinguishing between crop types and tree species (Immitzer et al., 2016). As shown in Fig. A4, both evergreen and deciduous vegetation exhibit similar trends in Sentinel-2 multispectral bands, but they display significant differences in spectral reflectance values. This indicates that median-composited bands of Sentinel-2, along with constructed spectral indices, can be used to distinguish between evergreen and deciduous vegetation. Median composites are affected by the number of available images, so we ensured that there was a minimum of three high-quality observations across the entire TP when preprocessing the annual Sentinel-2 images. The composites from at least three Sentinel images make it possible to achieve the seamless effect shown in Fig. A4 in various locations over large areas of the TP. The integration of multiple satellite images over time helps capture the phenology of different vegetation types while mitigating the influence of outliers (Carrasco et al., 2019; Pizarro et al., 2022; Tu et al., 2020; Verde et al., 2020; Xie et al., 2019).

However, relying solely on median-composited bands of Sentinel-2 and constructed spectral indices may not suffice to achieve a high classification accuracy, emphasizing the importance of multisource data. Notably, elevation emerges as the most important feature among all the ancillary ones (Fig. A1), reflecting the natural distribution of vegetation types, which is predominantly shaped by latitudinal zonation in the mountainous TP (Sherman et al., 2008) (Figs. A5 and 7). Conversely, in flat areas where the vegetation distribution is minimally influenced by the topography or in urban areas where the vegetation distribution is affected by anthropogenic activity, topographic information may exhibit limitations in land cover classification (Zeng et al., 2019). Thus, leveraging features derived from multisource data allows us to amplify and capture differences between evergreen and deciduous vegetation as well as between shrubs and woodlands, ultimately leading to a high classification accuracy (Xu et al., 2018; Yan et al., 2023).

Figure 7 Detailed comparison of TP_LC10-2022 with four other global land cover products. Panels (a), (b), and (c) present 0.1° × 0.1° grids used for detailed comparisons. Legend fusion rules for WorldCover2021 and FROM_GLC10-2017 are provided in Table A2, and those for FROM_GLC30-2015 and GLC_FCS30-2020 are given in Table A3.

4.3 Intercomparison with other products

The land cover samples selected remained stable across the years from 2013 to 2022 for the other four land cover products, thus making them comparable to our TP_LC10-2022 map. Therefore, we validated the aggregation of samples into eight categories and assessed the performance of TP_LC10-2022 and the four other land cover products in the TP region, as depicted in Table 5.

For shrubland, the classification performance of the four global land cover products is remarkably low. Notably, ESA WorldCover2021 achieves a PA and UA of 0 for shrubland classification. Among these land cover products, FROM_GLC30-2015 exhibits the highest UA for shrubland classification, albeit a mere 59.3 %. This suggests substantial shortcomings in the precise classification of shrubland in the TP region by the current land cover products.

Table 5 Comparison of mapping accuracy based on validation samples merged into eight land cover types. The fusion rules for the five land cover products are provided in Table A2. A bold font indicates the highest UA and PA for each land cover type as well as the highest OA and kappa among the five products.

BL: bare land; GL: grassland; FST: forest; SHR: shrubland; WB: water body; WL: wetland; CV: cultivated vegetation; PIS: permanent ice and snow.

A simultaneous visual comparison of the five products was conducted. In Fig. 7a, TP_LC10-2022 and FROM_GLC30-2015 exhibit superior performance, revealing more intricate forest details compared to the other products. Notably, other products largely disregard vast areas of high-elevation alpine shrublands above the timberline, while TP_LC10-2022 delineates them (shown in brown) and exhibits distinct vertical zonation. In Fig. 7b, the other four products, particularly FROM_GLC30-2015 and GLC_FCS30-2020, tend to misclassify shrublands as forests, whereas TP_LC10-2022 accurately differentiates between forests and shrubs. In Fig. 7c, both FROM_GLC10-2017 and FROM_GLC30-2015 depict scattered shrublands but lack continuity. The four products other than TP_LC10-2022 overestimate grasslands and underestimate the extent of shrubland areas. This discrepancy may stem from the similar phenological characteristics of deciduous shrublands and meadows, which poses difficulties when distinguishing them based on spectral features (X. Li et al., 2021). However, TP_LC10-2022 integrates topographic and climatic factors as classification features, facilitating precise differentiation between shrublands and grasslands.

Lakes and glaciers are the sentinels of global climate change and constitute the foundation of the TP as a crucial water source for surrounding regions (Zhang et al., 2017; Zhang and Duan, 2021). Precisely extracting the boundaries of lakes and glaciers is imperative for quantitatively monitoring lake expansion and glacier melting as well as understanding the dynamic relationship between them and precipitation (R. Zhao et al., 2022; Tong et al., 2016; J. Zhang et al., 2021). Our land cover data, samples, and mapping methodology can serve as baseline support for these endeavors (Yan et al., 2020; Korzeniowska and Korup, 2017), facilitating the effective utilization of available water resources and promoting the sustainable development of the economy and society in the Greater Tibetan Plateau area and downstream regions of rivers originating from the TP (Ding et al., 2019).

Alpine forests play a crucial role in carbon storage and sequestration, thereby enhancing ecosystem services in the TP (Lin et al., 2023; Z. Wang et al., 2022; H. Zhao et al., 2023). Our study revealed that TP_LC10-2022 yielded the smallest forested area (8.60 %), while GLC_FCS30-2020 and FROM_GLC30-2015 yielded the largest and second-largest areas of alpine forest, respectively (12.86 % and 11.89 %) (Table A4). Conversely, the area of shrubland exhibited nearly the opposite trend (Table A4). Confusion also arises between alpine grassland and bare land, potentially leading to variations in carbon storage estimation for each vegetation type. These discrepancies could impact efforts related to forest resource protection and grassland management for animal husbandry (Li et al., 2020; Yu et al., 2022).

Alpine screes are extensively distributed across the TP, but they are frequently disregarded from other products. Our product presents the initial description of alpine scree vegetation locations, which will contribute to environmental monitoring and biodiversity research in the periglacial zone of the TP (Li et al., 2014). Shrublands play a vital role as carbon sinks in ecosystems and have substantial implications for biomass estimation and global carbon cycling (Ma et al., 2021; Nie et al., 2018). TP_LC10-2022 accurately predicts the spatial distribution of shrublands, which is of considerable importance when forecasting the impact of future changes in the biomass and carbon cycle on global-scale ecosystems (Chang et al., 2022).

High-resolution and accurate land cover data encompassing diverse vegetation types are crucial for monitoring large-scale alpine vegetation dynamics (F. Wang et al., 2023; Z. Wang et al., 2022; Wang et al., 2020). For instance, relying on land cover maps such as ESA WorldCover as the foundation to examine tree lines and vegetation lines in the TP may lead to the underestimation of tree lines due to misclassifications of grasslands and shrublands (Fig. 7) (Zou et al., 2023). Additionally, the vegetation line may also be underestimated because of the absence of alpine scree (Fig. 7). In our future work, we aim to leverage Sentinel-2, Sentinel-1, and other multisource data to annually generate TP_LC10 products. This approach will facilitate alpine vegetation monitoring and change detection, thereby enriching our comprehension of the dynamic TP amidst intensifying global climate change (Y. Wang et al., 2022).

5 Data availability

The TP_LC-2022 product generated in this paper is available at https://doi.org/10.5281/zenodo.8214981 (Huang et al., 2023a). The TP_LC-2022 product covering the entire Tibetan Plateau is grouped into 54 3° × 3° tiles in the GeoTIFF format (EPSG: 4326), which are named “TP_LC10-2022_E**N**.tif”, where “E**N**” corresponds to the longitude and latitude information for the upper left corner of each regional land cover map. The multisource data used in this study, including the Sentinel-2 data, can be directly accessed from Google Earth Engine.

The corresponding sample dataset, produced by manual interpretation and field trips, is available at https://doi.org/10.5281/zenodo.8227942 (Huang et al., 2023b). The classification map can be viewed at https://cold-classifier.users.earthengine.app/view/tplc10-2022 (Huang et al., 2024).

6 Conclusions

We present a detailed land cover map of the Tibetan Plateau that includes 12 vegetation types and 3 non-vegetation types at 10 m spatial resolution for the year 2022 (TP_LC10-2022). This was obtained by integrating multisource data, including Sentinel-1, Sentinel-2, SRTM, CHIRPS, and ERA5-Land data, and applying and comparing four classification models via GEE. TP_LC10-2022 achieved an overall accuracy of 86.5 % and a kappa coefficient of 0.854 % when using the RF model, which outperformed other classification models, including GTB, MD, and SVM. Comparisons between TP_LC10-2022 and four widely used land cover products (GLC_FCS30-2020, FROM_GLC30-2015, FROM_GLC10-2017, and WorldCover2021) demonstrated that TP_LC10-2022 has a higher overall accuracy and reflects the local-scale variations of vegetation types with latitude. In particular, TP_LC10-2022 incorporated unique land cover types like alpine scree, alpine grassland, and alpine meadow, which accounted for 54.23 % of the total coverage. Moreover, it accurately depicted the distribution of shrubland, which occupied 4.63 % of the TP and was underestimated in the other products. The proposed vegetation classification system for the TP can serve as a foundation for land cover mapping in this region and a reference approach for mapping shrubland globally. The developed TP_LC10-2022 product can facilitate the monitoring of vegetation changes and the study of the response to climate change in the TP.

Appendix A

Table A1 Optimal parameters for random forest (RF), gradient tree boosting (GTB), minimum distance (MD), and support vector machine (SVM) in this study.

Table A2 Table for cross-walking between different land cover products.

This table includes only land cover types present within the study area.

Table A3 Table for cross-walking between FROM_GLC30-2015 and GLC_FCS30-2020.

This table includes only land cover types present within the study area.
The “cloud” class in FROM_GLC30-2015 and “shrubland” class in GLC_FCS30-2020 have been omitted from the table due to their small areas.

Table A4 Area statistics from land cover products for the Tibetan Plateau.

BL: bare land; AG: alpine grassland; AM: alpine meadow; ENF: evergreen needle-leaved forest; DNF: deciduous needle-leaved forest; EBF: evergreen broadleaved forest; DBF: deciduous broadleaved forest; MF: mixed forest; ES: evergreen shrubland; DS: deciduous shrubland; WB: water body; WL: wetland; CV: cultivated vegetation; PIS: permanent ice and snow; AS: alpine scree.
The unit of area is 10⁴ km², and the unit of proportion is %.
Please refer to Table A3 for the rules for comparing land cover between FROM_GLC30-2015 and GLC_FCS30-2020.
The “cloud” class in the FROM_GLC30-2015 product and the “shrubland” class in the GLC_FCS30-2020 product have been omitted from the table due to their small areas.
All built-up pixels have been merged with bare land.

Figure A1 Statistical chart of the importance of different features in the random forest classification model. AP: annual precipitation; AMT: annual mean temperature.

Figure A2 Landsat NDVI time series and HANTS-filtered NDVI time series for stability verification. Panel (a) depicts a deciduous needle-leaved forest, while panel (b) shows a transition from forest to farmland at the edge of the deciduous broadleaved forest in 2015; this area was annually cultivated following deforestation.

Figure A3 Number of available observations with cloud cover <10 % from Sentinel-2 optical data for the Tibetan Plateau during summer 2022 (1 June 2022 to 31 August 2022).

Figure A4 Sentinel-2 spectral curves for forest and shrubland types. The spectral curve for each type was derived by calculating the average and standard deviation of the surface reflectance across all samples for the processed cloud-free Sentinel-2 median composite for 2022 in the Tibetan Plateau. ENF: evergreen needle-leaved forest; DNF: deciduous needle-leaved forest; EBF: evergreen broadleaved forest; DBF: deciduous broadleaved forest; MF: mixed forest; ES: evergreen shrubland; DS: deciduous shrubland.

Figure A5 Box plot derived from the SRTM for the distribution of sample elevation across different land cover types in the Tibetan Plateau. BL: bare land; AS: alpine scree; AG: alpine grassland; AM: alpine meadow; ENF: evergreen needle-leaved forest; DNF: deciduous needle-leaved forest; EBF: evergreen broadleaved forest; DBF: deciduous broadleaved forest; MF: mixed forest; ES: evergreen shrubland; DS: deciduous shrubland; WB: water body; WL: wetland; CV: cultivated vegetation; PIS: permanent ice and snow.