Annual high-resolution grazing-intensity maps on the Qinghai–Tibet Plateau from 1990 to 2020

Grazing activities constitute the paramount challenge to grassland conservation over the Qinghai–Tibet Plateau (QTP), underscoring the urgency of obtaining detailed information regarding the extent, patterns, and trends of grazing to enable efficient grassland management and sustainable development. Here, to inform stakeholders about these issues, we provided the first annual Gridded Dataset of Grazing Intensity (GDGI), with a resolution of 100 m, from 1990 to 2020 for the QTP. The five most commonly used machine learning algorithms were leveraged to develop a livestock spatialization model, which spatially disaggregates the livestock census data at the county level into a detailed 100 m × 100 m grid based on seven key predictors from terrain, climate, vegetation, and socio-economic factors. Among these algorithms, the extreme-tree (ET) model performed the best in representing the complex nonlinear relationship between various environmental factors and livestock intensity, with an average absolute error of just 0.081 SU ha⁻² (where SU denotes sheep units), a rate outperforming the other models by 21.58 %–414.60 %. By using the ET model, we further generated the GDGI for the QTP to reveal the spatio-temporal heterogeneity and variations in grazing intensities. The GDGI indicates that grazing intensity remained high and largely stable from 1990 to 1997, followed by a sharp decline from 1997 to 2001 and fluctuations thereafter. Encouragingly, compared to other open-access datasets for grazing distribution on the QTP, the GDGI has the highest accuracy, with the determinant coefficient (R²) exceeding 0.8. Given its high resolution, recentness, and robustness, we believe that the GDGI dataset can significantly enhance our understanding of the substantial threats to grasslands emanating from overgrazing activities. Furthermore, the GDGI product holds considerable potential as a foundational source for other research, facilitating the rational utilization of grasslands, refined environmental impact assessments, and the sustainable development of animal husbandry. The GDGI product developed in this study is available at https://doi.org/10.5281/zenodo.10851119 (Zhou et al., 2024).

1 Introduction

Livestock is a crucial contributor to global food systems through the provision of essential animal proteins and fats and plays a significant role in supporting human survival and socio-economic development (Gilbert et al., 2018; Godfray et al., 2018; Humpenöder et al., 2022; Kumar et al., 2022). However, the escalating increase in human demand for meat and dairy products over recent decades has triggered a livestock boom, which, in turn, has increasingly threatened grassland ecosystems and placed a heavy burden on the environment through overgrazing and land use change (Tabassum et al., 2016; Wei et al., 2022; Minoofar et al., 2023). It is estimated that up to 300 ×10⁶ ha of land is used globally for grazing and cultivating fodder crops (Tabassum et al., 2016). Grazing activities could alter vegetation phenology and community structure (Dong et al., 2020) and could trigger deforestation (García Ruiz et al., 2020), grassland degradation (Sun et al., 2020), soil erosion (Shakoor et al., 2021), and associated direct releases of greenhouse gases that lead to climate change feedback (Godfray et al., 2018; Chang et al., 2021). Additionally, livestock are responsible for large-scale dispersion of pathogens, organic matter, and residual medications into soil and groundwater, thereby contaminating the environment (Venglovsky et al., 2009; Tabassum et al., 2016; Hu et al., 2017; Muloi et al., 2022). Consequently, more and more scholars have called for attention to provide reliable contemporary datasets to illustrate the spatio-temporal heterogeneity and variations of livestock (Petz et al., 2014; Fetzel et al., 2017; Zhang et al., 2018; Li et al., 2021).

One of the major challenges in monitoring grazing activity at regional or even larger scales is the determination of the livestock distribution pattern. Despite the importance of geographical grazing information, high spatio-temporal grazing datasets remain unavailable, posing the most critical challenge to grassland management, particularly for vulnerable grassland ecosystems in fragile regions grappling with economic- and sustainable-development contradictions (Meng et al., 2023; Pozo et al., 2021; Miao et al., 2020; He et al., 2022). In the early 2000s, the Food and Agriculture Organization of the United Nations (FAO) launched the Gridded Livestock of the World (GLW) project to facilitate a detailed evaluation of livestock production, aiming to provide pixel-scale livestock densities instead of traditional administrative unit benchmarks (Nicolas et al., 2016). Consequently, the world's inaugural dataset of livestock spatialization (GLW1) was released in 2007, providing the first globally standardized livestock density distribution map at a spatial resolution of 0.05 decimal degrees (≈ 5 km at the Equator) for 2002. It was not until 2014 that an updated GLW2 map with a 1 km resolution for 2006 was released; this updated map used a stratified regression approach, superior spatial resolution predictor variables, and more detailed livestock census data (Robinson et al., 2014). Furthermore, an evolutionary step in machine learning technology saw Gilbert et al. (2018) using random forest algorithms to forge a global livestock distribution map with a 10 km resolution for 2010 (GLW3), succeeding traditional multivariate regression methods and surpassing the precision of previous GLW1 and GLW2 maps. Beyond these global mappings, several maps with different scales have also been published, including ones on intercontinental, national, state or provincial, and local scales (Neumann et al., 2009; Prosser et al., 2011; Van Boeckel et al., 2011; Nicolas et al., 2016). However, these maps are fundamentally coarse due to constraints regarding the availability of fine-scale and contemporary census data, the grazing-spatialization method, and the identification of appropriate indicators, thereby limiting their application to local- or regional-scale studies (Nicolas et al., 2016; Gilbert et al., 2018; Robinson et al., 2014). Hence, there is an emergent demand for more refined grazing map products (Mulligan et al., 2020; Martinuzzi et al., 2021).

An example of this need can be observed in the Qinghai–Tibet Plateau (QTP), the world's most elevated pastoral region and an important grazing area in China (Zhan et al., 2023). It is made up of abundant grassland that spans 1.5 ×10⁶ km², accounting for 50.43 % of China's total grassland area, with Yak and Tibetan sheep as the primary grazing livestock (Feng et al., 2009; Cai et al., 2014; Zhan et al., 2023). Over recent decades, the QTP has undergone escalating grassland degradation, leading to many ecological and socio-economic problems, which represent the urgent need for detailed livestock distribution datasets (M. Q. Li et al., 2022). Unfortunately, despite researchers' efforts at mapping the QTP's grazing intensity, current livestock datasets still suffer from coarse spatio-temporal resolutions and modeling accuracies. Apart from the aforementioned global grazing dataset, several other maps also cover the QTP. For instance, Liu (2021) generated annual 250 m gridded carrying-capacity maps for 2000–2019 by employing multiple linear regressions of livestock numbers, population density, net primary productivity (NPP), and topographic features. Li et al. (2021) used machine learning algorithms to produce gridded livestock distribution data at 1 km resolution for 2000–2015 in western China at a 5-year interval based on county-level livestock census data and 13 factors from land use practice, topography, climate, and socio-economic aspects, including grassland coverage, arable land coverage, forest land coverage, desert coverage, normalized difference vegetation index (NDVI), elevation, slope, daytime surface temperature, precipitation, distance to river, travel time to major cities, population density, and gross domestic product (Li et al., 2021). A contribution by Meng et al. (2023) brought forth longer annual time series grazing maps by using random forest models, integrating climate, soil, NDVI, water distance, and settlement density to decompose county-level livestock census data to a 0.083° (≈ 10 km at the Equator) grid for 1982–2015 (Meng et al., 2023). Similarly, Zhan et al. (2023) also used random forest algorithms to combine 11 influence factors to provide a winter- and summer-grazing density map at 500 m resolution for 2020 (Zhan et al., 2023).

However, although these maps have provided good help in understanding grazing conditions on the QTP, there are currently still no maps that can satisfy the need for fine-scale grassland management with a long time span. In addition, the available livestock distribution maps of the QTP still need improvement in terms of modeling techniques and factor selection to obtain high-precision livestock spatialization data. For example, traditional methods like multiple linear regression, while having been proven to be fundamental and widely applicable for livestock spatialization (Robinson et al., 2014; Ma et al., 2022), are being challenged by the development of computational science in recent years. Among them, machine learning technology is providing new opportunities for more accurate predictions of livestock distribution (García et al., 2020). Random forest regression, for instance, is currently widely used to construct global, national, and regional livestock spatialization datasets and has been proved to have much better accuracy than traditional mapping techniques (Rokach, 2016; Nicolas et al., 2016; Gilbert et al., 2018; Dara et al., 2020; Chen et al., 2019; Li et al., 2021). Nevertheless, other more advanced machine learning methods with superior feature learning and more robust generalization capabilities remain largely untapped for modeling geographic data (Ahmad et al., 2018; Heddam et al., 2020; Long et al., 2022). Thus, exploring the potential application of new advanced machine learning technologies in livestock spatialization remains a critical task. Furthermore, selecting suitable factors that influence livestock grazing preferences is also the other critical challenge for enhancing the precision of grazing-distribution datasets (Meng et al., 2023). Livestock grazing activities are often affected by abiotic and biotic resources, including climatic and environmental factors (Waha et al., 2018), herd foraging and grazing behaviors (Garrett et al., 2018; Miao et al., 2020), and conservation-oriented policies (Li et al., 2021). For instance, regions exceeding elevations of 5600 m or with slopes greater than 40 % are customarily unsuitable for grazing (Luo et al., 2013; Mack et al., 2013; Robinson et al., 2014; Chen et al., 2019). Moreover, the livestock generally prefer areas abundant in water and pasture resources for foraging (Li et al., 2021). Besides, ecological conservation policies also exert a substantial influence, significantly affecting grazing distribution relative to the level of conservation priority. In addition, the health status of the grassland is an important factor influencing whether or not livestock choose to feed (Li et al., 2021). Consequently, indicators related to the above aspects are often employed to gauge the spatial heterogeneity of livestock distribution (Allred et al., 2013; Sun et al., 2021; Meng et al., 2023). Nonetheless, some of the most commonly used indicators like NPP or NDVI can result in misconceptions as they may not fully characterize the grazing intensity. For example, grasslands with high NPP or NDVI are often preferred by livestock, but this does not necessarily correlate with grazing intensity in nature reserves due to strict policy restrictions (Veldhuis et al., 2019; O'Neill and Abson, 2009; Y. Zhang et al., 2021). Conversely, areas with sparse grassland cover may support considerable livestock numbers despite evidence of degradation (B. H. Zhang et al., 2021; Guo et al., 2015). Accordingly, further investigation into novel indicators is imperative to enhance the correlation between grassland and grazing intensity, thereby optimizing the integration of such influencing factors into grazing-spatialization models.

In summary, the QTP is in pressing need of a high-spatio-temporal-resolution grazing dataset to address urgent and realistic challenges. However, the existing livestock datasets specific to the QTP are fraught with several insufficient and predominantly concerning rough resolutions, relatively backward census data, and conventional methods in livestock spatialization. Moreover, the discrepancies in terms of predictive indicators and modeling approaches within these datasets discourage their application in time series analysis. Consequently, the generation of high-resolution and high-quality grazing map products has emerged as the most pressing challenge for the QTP. Here, we aim to (1) establish a methodological framework by using more rational models and indicators than traditional studies to achieve fine-scale livestock spatialization; (2) select the grazing spatialization model with good performance by incorporating multi-source data with advanced machine learning techniques; and (3) ultimately, provide an annual grazing-intensity dataset with 100 m resolution spanning from 1990 to 2020. These maps can not only provide fundamental datasets with finer spatio-temporal resolutions to address the limitations of existing grazing-intensity maps but also enable a better understanding of sustainable management practices, as well as of other grassland-related issues, across the QTP.

Figure 1 The geographic zoning map of the Qinghai–Tibet Plateau (QTP) superposed with grassland vegetation.

2 Data and methods

2.1 Study area

Known as Asia's water tower and the world's third pole, the QTP is geographically situated between 26°00^′–39°47^′ N and 73°19^′–104°47^′ E, with a total area of about 2.61 ×10⁶  km² (Fig. 1). Its jurisdiction encompasses 182 counties within six provincial regions of China, including the Tibet Autonomous Region, Qinghai Province, Xinjiang Uygur Autonomous Region, Gansu Province, Sichuan Province, and Yunnan Province (Meng et al., 2023). The elevation of the QTP predominantly ranges between 3000 and 5000 m, with an average altitude exceeding 4000 m. With grasslands constituting over half of its land cover, the QTP emerges as one of the most important pastoral areas in China. Alpine steppe, alpine meadow, and temperate steppe characterize the main grassland types on the QTP (Han et al., 2019; Zhai et al., 2022; Y. Y. Zhu et al., 2023). The complex geographical and climatic conditions of the QTP contribute to the markedly heterogeneous grassland distribution, which correspondingly leads to the high heterogeneity in livestock distribution. Moreover, social and economic developments, coupled with policy initiatives directed towards grassland restoration, have noticeably impacted the livestock numbers on the QTP over recent decades (Li et al., 2021; Li et al., 2016).

Table 1 Summary of the livestock data used in this study.

2.2 Data source

2.2.1 Census livestock data

The county-level census livestock data for the period between 1990 and 2020 were obtained from the Bureau of Statistics of each county across the QTP (Table 1). The data include the number of cattle, sheep, horse, and mule, with the exception of counties in Yunnan Province, which lack data for the years from 1990 to 2007; counties in Ganzi Prefecture in Sichuan Province, which lacks data for the years from 1990 to 1999; and Muli County in Sichuan Province, which lacks data for the years from 1990 to 2007. For these counties belonging to the same prefecture, including counties in Ganzi and Aba prefectures in Sichuan Province, we used the livestock census data at the prefecture level to carry out spatialization. For these counties in Yunnan Province, since they belong to different municipalities, it is not reasonable to replace them with municipal-level data. For these counties without livestock census data for some years, we supplemented the missing data by linear interpolation with grazing-density data in available years. In total, livestock data were available for 182 counties, and 4998 independent records were finally generated. Furthermore, the respective quantities of different livestock types are converted to standard sheep units (SU), in compliance with the Chinese national regulations (Meng et al., 2023).

Due to the difficulty of collecting township-level census livestock data, the validation data at the township scale collected in this study only involved the townships of Baching County (2010–2018) and Gaize County (2018–2020) in Tibet and Hongyuan County in Sichuan Province (2008). The township-level census livestock data cumulatively involve 18 townships with a total of 112 records and were only used for auxiliary validation of the simulation results.

The validation data at the pixel scale also encompass a total of 112 records from 68 sites, which were collected from the literature, questionnaires, and field surveys. Specifically, 93 records at 49 sites spanning the 1990–2021 period were obtained from 17 literature sources, and 19 records from 19 sites were obtained from the questionnaires and the field survey in 2021. The detailed information for these records can be found in Fig. S3 and Table S3 in the Supplement.

Table 2 Summary of factors affecting grazing activities on the QTP.

2.2.2 Factors affecting grazing activities

Livestock grazing activities are often affected by abiotic and biotic resources, including climatic and environmental factors (Waha et al., 2018), as well as herd foraging and grazing behaviors (Garrett et al., 2018; Miao et al., 2020). For instance, high-altitude and steep hillsides are unsuitable for grazing due to terrain constraints, and the distribution of herders directly affects the grazing areas (Luo et al., 2013; Mack et al., 2013; Robinson et al., 2014; Chen et al., 2019). Moreover, the livestock generally prefer areas abundant in water and pasture resources for foraging (Li et al., 2021). Therefore, in this study, topography, climatic, environmental, and socio-economic impacts were considered to be influential factors in grazing activities (Li et al., 2021; Meng et al., 2023).

We utilized correlation analysis and the random forest importance-ranking tool to eliminate redundant environmental factors and to determine the contribution of each factor. Ultimately, altitude, slope, distance to water source, population density, air temperature, precipitation, and human-induced impacts on NPP (HNPP) were selected as indicators (Table 2). Specifically, elevation is derived from the digital elevation model (DEM) dataset accessible via the Resource and Environmental Data Cloud Platform of the Chinese Academy of Sciences (https://www.gscloud.cn, last access: 12 July 2022), which also facilitated slope calculation. Data on rivers and lakes were obtained from the National Tibetan Plateau Data Center (https://data.tpdc.ac.cn, last access: 21 October 2022), and the nearest Euclidean distance from each pixel to rivers or lakes is calculated accordingly. Meteorological elements such as daily air temperature and precipitation were downloaded from the China Meteorological Data Service Center (http://data.cma.cn, last access: 18 April 2022). For the grid dataset that is not conditionally available, including population density, temperature, precipitation, and HNPP, we detailed the creation process in the Supplement. All datasets utilized in this study were harmonized in relation to consistent coordinate systems and resolutions (WGS 1984 Albers, 100 m).

Figure 2 Flowchart of creating grazing-intensity maps using different methods and source products.

2.3 Methodological framework

We adopted a comprehensive methodological framework for mapping high-resolution grazing intensity on the QTP. This study applied FAO's assumption that the relationship between environmental factors and livestock intensity is identical at both the administrative and pixel level. Three major steps are included to predict the distribution pattern of grazing intensity: (1) identifying factors affecting grazing activities and extracting theoretically suitable areas for livestock grazing, (2) building a grazing-spatialization model, and (3) filtering the model and correcting the grazing map. An exhaustive explanation of each step is provided in Fig. 2.

2.3.1 Identifying factors and theoretically suitable areas for grazing

In this study, we assumed that grazing activities are confined solely to grassland. Consequently, the potential grazing areas for each year were identified on the basis of grassland boundaries, which were extracted from the 30 m annual land cover dataset (CLCD) (Yang and Huang, 2021). Furthermore, grasslands with slopes over 40 % and elevations higher than 5600 m were considered to be unsuitable for grazing and were therefore excluded from the potential grazing area in the subsequent simulations (Robinson et al., 2014). In addition, the grasslands with a population density greater than 50 inhabitants per square kilometer were also excluded (Li et al., 2018). The remaining isolated grasslands were thus categorized as theoretically feasible grazing regions.

The spatial patterns of abiotic and biotic resources, incorporating food availability, environmental stress, and herder preference, critically affect grazing activities (Meng et al., 2023). In light of this, seven influencing factors covered by four categories were selected for grazing-intensity mapping (Fig. 2.1).

2.3.2 Building a grazing-spatialization model

By assessing regional statistics, the annual average values for each grazing influence factor were extracted from the theoretically suitable grazing areas at the county scale and were further used as independent variables in the model construction. The dependent variable for the model was acquired by determining the livestock density within each county, followed by a logarithmic transformation of the values to normalize the distribution of the dependent variable. Consequently, a total of 4998 samples were derived from the aforementioned independent and dependent variables. Of these samples, 70 % were allocated to model training, while the remaining 30 % comprised the test sets, serving to validate the model's performance. Subsequently, we built grazing-spatialization models using five machine learning algorithms at the county scale, including support vector (SV) regression (Cortes and Vapnik, 1995; Lin et al., 2022), K-nearest neighbors (KNNs) (Cover and Hart, 1967), gradient-boosting (GB) regression (Friedman, 2001; Pan et al., 2019), random forests (RFs) (Breiman, 2001), and extra-tree (ET) regression (Geurts et al., 2006; Ahmad et al., 2018) (see the Supplement for details). Lastly, to assess the accuracy of the spatialized livestock map, the predicted livestock intensity values were juxtaposed with the livestock statistical data from each respective county.

2.3.3 Correcting the grazing map

We further used the optimal model to predict the geographical distribution of grazing density across the QTP. To maintain better consistency between the predicted livestock number and the census data, the estimated results were adjusted using the census livestock numbers at the county scale as a control according to Eq. (1). Consequently, the corrected and refined map is presented as the final grazing-intensity map in this study.

$$\begin{array}{}\text{(1)}& L_{\text{correction}}={\displaystyle \frac{L_{\text{CCensus}}}{L_{\text{Cgrid}}}}\times L_{\text{grid}}\end{array}$$

In the above, L_correction is the predicted pixel-scale livestock number after adjustment, L_Cgrid represents the estimated livestock number for each county, L_CCensus is the census livestock number for each county, and L_grid refers to the predicted livestock number at the pixel scale.

2.4 Accuracy evaluation

We used three accuracy validation indices, namely coefficients of determination (R²), mean absolute error (MAE), and root mean square error (RMSE), to evaluate the performance of five machine learning algorithms through a comparison of the predicted value with the census data. The definitions of three metrics are presented in Eqs. (2) to (4).

$$\begin{array}{}\text{(2)}& {\displaystyle }& {\displaystyle }{R}^{\mathrm{2}}=\mathrm{1}-{\displaystyle \frac{\sum _{i=\mathrm{1}}^n(C_i-P_i{)}^{\mathrm{2}}}{\sum _{i=\mathrm{1}}^n(C_i-\stackrel{\mathrm{‾}}{C}{)}^{\mathrm{2}}}}\text{(3)}& {\displaystyle }& {\displaystyle }\text{MAE}={\displaystyle \frac{\mathrm{1}}{n}}\sum _{i=\mathrm{1}}^n\mathrm{|}{C}_i-P_i\mathrm{|}\text{(4)}& {\displaystyle }& {\displaystyle }\text{RMSE}=\sqrt{{\displaystyle \frac{\mathrm{1}}{n}}\sum _{i=\mathrm{1}}^n(C_i-P_i{)}^{\mathrm{2}}}\end{array}$$

In the above, C_i and P_i are the census livestock data and the predicted value for county i, respectively; $\stackrel{\mathrm{‾}}{C}$ represents the mean census value for all counties; and n gives the total number of counties.

2.5 Uncertainty evaluation

Uncertainty in our grazing intensity maps can stem from multiple sources, such as the constraints of cross-scale modeling and the intrinsic inaccuracies of the input data. To quantify these uncertainties, we utilized the Monte Carlo (MC) method, conducting 100 iterations of simulation. Subsequently, we evaluated uncertainty through the mean relative error (MRE) and assessed the model's robustness using the standard deviation (SD) following established methodologies (Yang et al., 2020; Alexander et al., 2017; McMillan et al., 2018). The definitions for these metrics are delineated in Eqs. (5) to (7).

$$\begin{array}{}\text{(5)}& {\displaystyle }& {\displaystyle }\text{MC}={\displaystyle \frac{\mathrm{1}}{n}}\sum _{i=\mathrm{1}}^{n}f\left(x_i\right)\text{(6)}& {\displaystyle }& {\displaystyle }\text{MRE}={\displaystyle \frac{\mathrm{1}}{n}}\sum _{i=\mathrm{1}}^n\left|{\displaystyle \frac{x_i-\stackrel{\mathrm{‾}}{x}}{\stackrel{\mathrm{‾}}{x}}}\right|\text{(7)}& {\displaystyle }& {\displaystyle }\text{SD}={\displaystyle \frac{\mathrm{1}}{n}}\sum _{i=\mathrm{1}}^{n}f\left(x_i\right)\sqrt{{\displaystyle \frac{\mathrm{1}}{n}}\sum _i^n(x_i-\stackrel{\mathrm{‾}}{x})}\end{array}$$

In the above, x_i are random samples, f(x_i) is the function evaluated at x_i, and n is the number of simulations. $\stackrel{\mathrm{‾}}{x}$ represents the mean value for all simulation maps.

Figure 3 Scatterplots of model-predicted livestock numbers and census grazing data (scaled by logarithm) at the county level. The solid red line and the solid black line are the fitting line and the 1:1 diagonal line, respectively.

3 Results

3.1 Performances of models

Table 3 summarizes the efficiency of the five machine learning models used, with a consideration of all three accuracy evaluators of R², MAE, and RMSE. It can be seen that the ET model performs the best, with its R² exceeding 0.955 and with its MAE (0.081 SU ha⁻²) and RMSE (0.164 SU ha⁻²) being significantly lower than the values of RF, GB, KNN, and support vector machine (SVM) models. Figure 3 illustrates the correlation between the census livestock data and the livestock numbers predicted by the model for each county from 1990 to 2020. It is demonstrated that the ET-predicted data displayed a distribution pattern consistent with that of other models, but the scatter points of the ET model were more convergent in relation to the 1:1 diagonal line, indicating a superior fit compared to the other models. These comparisons suggest that the ET model possesses superior robustness and can, therefore, provide stable estimations of livestock intensity on the QTP.

Table 3 Comparison of mapping accuracy for five machine learning models based on the same validation datasets.

Note that the MAE and RMSE have been calculated using inverse logarithmic transformation, representing the actual values.

Figure 4 Accuracy of the ET-predicted grazing-intensity results at a spatial resolution of 100 m from 1990 to 2020. (a) Absolute error (AE) between the predicted and the census data at the county scale from 1990 to 2020; (b) comparison of the predicted and census data of the whole QTP from 1990 to 2020; (c) spatial distribution of the mean absolute error (MAE) during 1990 to 2020 for each county.

Using the ET model, we projected the spatio-temporal distribution of grazing intensity across the QTP from 1990 to 2020 at a 100 m × 100 m resolution. To validate the accuracy of these predictive maps, we upscaled the pixel-level predictions to the county level and compared them against livestock census data (Fig. 4a and b). The results clearly show a high degree of consistency between the predicted livestock intensity and the county-level census data, especially in areas with lower grazing intensity (Fig. 4a and b). Specifically, while the mean census data indicated 2.983 SU ha⁻² for livestock intensity, our county-level predictions yielded an average of 3.106 SU ha⁻², with an MAE of 0.123 SU ha⁻², an RMSE of 0.580 SU ha⁻², and an R² value of 0.669. Additionally, 76.31 % of the counties (n = 3814) exhibited data discrepancies of no more than 0.6 SU ha⁻², and 91.74 % (n = 4585) had discrepancies under 1.0 SU ha⁻². Regarding the spatial distribution, areas with data discrepancies of less than 0.3 SU ha⁻² were predominantly located in the northwestern and southeastern regions of the QTP. In certain counties of the northeast and southwest, the variations were even below 1.0 SU ha⁻² (Fig. 4c).

Figure 5 Uncertainty analysis of grazing-intensity maps based on ET and Monte Carlo methods. (a) MRE of grazing-intensity maps from 1990 to 2020, (b) spatial distribution of MRE, and (c) spatial distribution of SD.

3.2 Evaluation of uncertainties

We have chosen the mean relative error (MRE) as a key metric for evaluating the simulation accuracy of grazing intensity within the QTP. Employing Monte Carlo simulations spanning the period from 1990 to 2020, our research findings demonstrate that the average MRE for grazing intensity across the QTP ranged between 6.84 % and 9.08 % (Fig. 5a). The spatial distribution of MRE indicates that the majority of the plateau exhibits low error margins. For example, in 2020, areas with an MRE of less than 5 % accounted for 35.86 % of the total grassland area, while those with an MRE below 10 % constituted 75.84 %. Only 3.38 % of the grasslands had an MRE exceeding 20 %, with these regions being primarily located in the southwestern portion of the QTP (Fig. 5b). Moreover, the robustness analysis suggests that the majority of regions within the QTP display relatively stable grazing-intensity trends. For instance, the overall standard deviation (SD) in 2020 was 0.059 SU ha⁻², with the northwestern region demonstrating remarkable stability, reflected in an SD of less than 0.005 SU ha⁻². Although some areas within the Yarlung Zangbo River basin and in the eastern part of Qinghai Province experienced higher variability, their SD was still maintained below 0.3 SU ha⁻² (Fig. 5c).

Figure 6 Validation of the GDGI dataset using 112 field grazing-intensity records at the pixel scale: (a) linear fitting results, (b) absolute error (AE) distribution.

3.3 Validation of the GDGI dataset

After employing the county-level livestock census as a benchmark for quality control, we obtained the annual Gridded Dataset of Grazing Intensity (GDGI) across the QTP, spanning 31 years from 1990 to 2020. We firstly confirmed the accuracy of the GDGI dataset based on 112 field grazing-intensity records at 68 sites (see Table S3 in the Supplement for details), which ranged from 0 to 5.61 sheep units per hectare (SU ha⁻²) and covered three main grasslands on the QTP: the alpine steppe (N = 62), alpine meadow (N = 46), and alpine desert steppe (N = 4). The GDGI dataset was assessed by undertaking a comparative accuracy assessment with the field grazing-intensity data (Fig. 6a). It can be seen that, in general, our dataset was highly consistent with the reference ground-truth validation data, with R² = 0.804, MAE = 0.572 SU ha⁻², and RMSE = 0.953 SU ha⁻². Moreover, the absolute errors between the GDGI data and the field grazing-intensity data were relatively small, with more than half of the records having an error below 0.3 SU ha⁻², 78.57 % having an error below 1.0 SU ha⁻², and 89.29 % having an error below 1.5 SU ha⁻² (Fig. 6b).

Figure 7 Validation of the GDGI dataset using census livestock data at the township level: (a) linear fit of predicted number and census data, (b, c) logistic fit of grazing-intensity data and grassland area.

We further validated the precision of the GDGI dataset using the township-level livestock census data. Encouragingly, the evaluation results showed that the GDGI dataset demonstrates excellent performance at the township scale (Fig. 7a), with R² of 0.867, MAE of 0.208 SU ha⁻², and RMSE of 0.276 SU ha⁻². In addition, similarly to the census data, the GDGI dataset indicated that some townships with few grasslands areas are still under high grazing pressure (Fig. 7b and c).

Figure 8 Validation of the GDGI maps using the census grazing data from 1990 to 2020: (a) violin plot of the census data and the predicted value, (b, c) spatial distribution (in SU per pixel), (d) temporal change (in SU per year) (only including 124 counties with livestock census data), (e–f) spatial distribution of SU changes tested by Sen's slope and Mann–Kendall tests. Note that ESI stands for extremely significant increase (slope > 0 and p < 0.01), SI stands for significant increase (slope > 0 and p < 0.05), NSI stands for non-significant increase (slope > 0 and p > 0.05), ESD stands for extremely significant decrease (slope < 0 and p < 0.01), SigDec stands for significant decrease (slope < 0 and p < 0.05), and NSD stands for non-significant decrease (slope < 0 and p > 0.05).

3.4 Spatio-temporal variations in grazing intensity

In terms of the temporal trends of grazing intensity, overall, the GDGI dataset exhibited consistent trends with the livestock census data (Fig. 8d–f). Specifically, the census data indicated that the livestock numbers remained high and largely stable from 1990 to 1997, followed by a sharp decline from 1997 to 2001 and then a period of fluctuation post-2001, which was successfully captured by the GDGI dataset. Moreover, the spatial heterogeneity of grazing intensity within the counties over the QTP was also effectively reflected by the GDGI dataset, a characteristic not illustrated by the census dataset. For example, areas of high grazing intensity were concentrated in the northeastern and south-central regions of the plateau, mainly including the eastern part of Qinghai Province, the southwestern part of Gansu Province, the northwestern part of Sichuan Province, and the eastern region of the Tibet Autonomous Region (Fig. 8e and f).

Over the past 31 years, 63.95 % of the plateau's grassland showed a decreasing trend in grazing intensity, with 49.80 % showing significant decreases, primarily located in the eastern Sanjiangyuan area and the southwestern region of the QTP (Fig. 8e and f). Meanwhile, grazing intensity saw an increase in 36.05 % of the grassland, but most of these areas (60.16 %) did not reach the level of significance and were mainly distributed in the northeastern plateau (Fig. 8e and f).

Table 4 Summary of map-derived parameters for this study and the seven other public gridded livestock datasets covering the QTP.

The abbreviations used in the table are as follows: LRA – linear regression analysis; MLR – multivariate linear regression; DNN – deep neural network.

Figure 9 Comparisons of different grazing datasets for the years 2001 and 2015: (a) spatial patterns, (b) predicted livestock number and census data at county scale, and (c) accuracy evaluation between predicted livestock number and census data at township scale.

4 Discussion

4.1 Comparison with other grazing-intensity maps

To further assess the effectiveness and reliability of the developed GDGI dataset, the mapping results were juxtaposed with seven publicly available grazing-intensity maps covering the QTP (Table 4). It can be seen that, despite their public availability, these maps were lacking in both spatial and temporal resolutions when juxtaposed with the GDGI maps. Our analysis was extended to four openly accessible gridded livestock datasets, including GI-Sun (Sun et al., 2021), ALCC (Liu, 2021), GI-Meng (Meng et al., 2023) and GLW (Gilbert et al., 2018). A commonality among all five maps was the consistency of the spatial patterns of grazing intensity, with prevalent high and low intensities in the northeastern and northwestern regions, respectively (Fig. 9). However, these maps differed significantly in terms of accuracy. As the grazing intensity maps of GLW and ALCC were produced based on the livestock census data in 2001 and 2015, an accuracy comparison for the corresponding years was conducted among the five datasets at both the county and township scale. Observations at the county scale indicate that all four datasets, with the exception of GI-Sun, are largely in alignment with the county census data (Fig. 9b). When examined at the township scale, GI-Sun and GLW demonstrate the most significant discrepancies, with MRE values surpassing 68 %. ALCC and GI-Meng follow, recording MREs of 30.69 % and 38.80 %, respectively. Additionally, the GDGI shows the highest degree of accuracy in relation to the township census data, as indicated by the lowest MAE and RMSE values (Fig. 9c). Moreover, the GDGI dataset spanning 31 years (1990–2020) was earmarked as a more suitable choice for long-term studies in comparison to the other four datasets. Regarding spatial distribution, the overall patterns of these grazing maps are largely consistent, exhibiting higher density patterns in the southeast and lower in the northwest. However, notable discrepancies are still apparent in the finer details. In general, in terms of visually representing the spatial distribution of livestock, the GDGI maps exhibit the best performance.

Several potential factors may contribute to the improved accuracy of the GDGI. First, the livestock census data used in GDGI are more detailed, aiding in enhancing the accuracy of the estimation results. Specifically, GI-sun, ALCC, GI-Meng, and GDGI all use county-level livestock statistics to map grazing intensity, whereas GLW3 and GLW4 are based on provincial-level census data, which results in their accuracy lagging significantly behind the other four datasets (Nicolas et al., 2016; Sun et al., 2021). Second, after ruling out areas unsuitable for grazing, we estimate grazing intensity by dividing the livestock numbers from statistical data by the area of suitable grazing regions. However, these maps differ in terms of their definitions of suitable grazing areas. In this study, as with the GI-sun and GI-Meng maps, we considered grazing to occur only on grasslands and further excluded unsuitable areas such as high elevations and steep slopes. This kind of definition is clearly more reasonable than the GLW series, which removed only waterbodies, urban core areas, and protected areas with relatively tight regulations regarding human activity (McSherry and Ritchie, 2013; He et al., 2022). However, the GI-Meng dataset considers the core areas of protected areas to be grazing-free regions; this does not match the actual situation on the QTP (Jiang et al., 2023; T. Li et al., 2022; Zhao et al., 2020). Those different thresholds for the definition of suitable grazing areas account for the fact that each map has different theoretical grazing regions. Third, the selection of models and environmental factors may also be a significant contributing factor, leading to variations in predictive accuracy. For instance, different algorithms were employed, including linear regression and machine learning methods (Nicolas et al., 2016; Li et al., 2021). Additionally, the environmental factors considered varied; specifically, the GDGI utilized the human-induced net primary productivity (HNPP) to represent grasslands, whereas other maps relied on net primary productivity (NPP) and the normalized difference vegetation index (NDVI) (Allred et al., 2013; Sun et al., 2021; Meng et al., 2023).

4.2 Spatial heterogeneity of grazing intensities

In general, the multiyear average grazing intensity on the QTP increased from west to east during 1990 to 2020, with broad spatial heterogeneity (Fig. 8). Highest grazing intensity was found mainly in the northeastern and south-central regions of the Plateau (mostly higher than 5.0 SU ha⁻²), while they were lowest in the northwest (mostly less than 1.0 SU ha⁻²). Over the past 31 years, the average grazing intensity decreased across most of the Plateau, but 36.05 % of the entire QTP grassland still encountered a continuous grazing-intensity increase, especially in the northeastern regions (Fig. 8).

The spatial heterogeneity of grazing intensities on the QTP may be attributed to the following reasons. First, complex geographic and climatic conditions on the QTP determine the heterogeneity of grasslands, which in turn affects livestock distribution (Wang et al., 2018; Wei et al., 2022). In general, the grazing-intensity patterns shown in the GDGI maps are basically consistent with the stocking-rate threshold patterns in the QTP grasslands, which both decreased from east to west (Q. Zhu et al., 2023). This phenomenon partially reflects the heterogeneity of the grasslands as the alpine meadows and the steppes are mainly distributed in the east and the west, respectively. Second, the dynamics of socio-economic development are obviously another important factor determining grazing intensity. In areas falling behind in terms of the socio-economic indicators, herders prefer to increase livestock in an effort to improve household incomes, leading to greater pressure on grasslands in these regions (Fang and Wu, 2022). In addition, the perceived increases in human population also resulted in the considerably increased need for more livestock (Wei et al., 2022).

The grazing-intensity dynamics across the QTP are partly reflective of the impacts of various management policies that have been implemented over distinct periods. For example, a significant increase in grazing intensity on the QTP was observed in the early 1990s, potentially a direct result of the introduction of the household contract responsibility system. Moreover, the grazing intensity experienced a pronounced decline from 1997 to 2001, as illustrated in Fig. 8d, indicative of the adverse effects of natural disasters. Notably, the severe snowstorms that struck Naqu in the central QTP during 1997–1998 are documented to have caused the mortality of over 820 000 livestock (Ye et al., 2020). Figure 8d further delineates a considerable upsurge in grazing intensity on the QTP between 2000 and 2010, aligning with the trends reported by Sun et al. (2021) and Li et al. (2021). This observed increase may be attributed to a rebound in grazing activity following the aforementioned natural disasters. In addition, Fig. 8d indicates a sustained decrease in grazing intensity post-2010 across the plateau, which can be predominantly ascribed to the implementation of extensive ecological conservation projects.

4.3 Implications for grazing management

Nearly half of the grasslands on the QTP have been reported to be degraded over the past 4 decades (Wang et al., 2018; Dong et al., 2020), with some reports even indicating that the degraded grassland has reached 90 % (Wang et al., 2021). It is widely recognized that overgrazing is the predominant and most pervasive unsustainable human activity continuing to drive grassland degradation on the QTP (Wang et al., 2018; Chen et al., 2019). Generally, these degraded grasslands on the QTP can be effectively restored by adaptive management (Wang et al., 2022). However, better management of grasslands requires a deeper understanding of the anthropogenic activities, which still remains an important challenge and can be effectively addressed by the GDGI dataset.

According to the GDGI maps generated in this study, high-intensity grazing activities are mainly concentrated in the northeastern and south-central parts of the QTP, with the grazing intensity in some areas even being nearly more than 10 times the average value of the entire plateau (Fig. 6b), and have exceeded the stocking-rate threshold of these grasslands (Q. Zhu et al., 2023). Population growth and the related increasing livelihood demands constitute one of the main reasons for this increase. To meet daily needs and to enhance household incomes, the herders have endeavored to increase livestock, thereby intensifying grazing pressures on the grasslands over the QTP (Fang and Wu, 2022; Abu Hammad and Tumeizi, 2012). Although the current average grazing intensity in the northwestern QTP (around 1.0 SU ha⁻²) is below the average stocking-rate threshold (around 1.5 SU ha⁻²) (Q. Zhu et al., 2023), the grassland management should still be given adequate attention. The northwestern QTP, the driest area on the plateau with low stocking-rate thresholds, has seen an increase in grazing intensity in recent years. Nevertheless, it must be noted that the stocking-rate threshold may exceed the carrying capacity because it is predicted to lead to extreme grassland degradation (Q. Zhu et al., 2023). The GDGI dataset also showed a similar pattern between the grazing-intensity data and the WorldPop data near the built-up areas, indicating higher grazing intensity around settlements compared to other regions on the QTP. In addition, the GDGI dataset also indicates that, from 1990 to 2020, although the grazing intensity of the Plateau has generally decreased, the hotspot areas for grazing activities have remained almost unchanged. This implies that these regions should be the focus of adaptive grassland management to effectively prevent grassland degradation, mainly based on the grass–livestock balance which varies with time and space.

Encouragingly, the GDGI dataset shows that the grazing intensity for two-thirds of the entire QTP grassland has decreased over the past 31 years, which is also consistent with the findings of other studies (Sun et al., 2021; Li et al., 2021). Recent decades of biodiversity protection; active restoration projects; and management measures, such as nature reserves, grazing exclusion, partial grazing bans, and fencing enclosure, are believed to have driven these decreases (Deng et al., 2017; Li and Bennett, 2019). In addition, most grassland in eastern Sanjiangyuan, mid-eastern Changtang, and the northern foothills of the Himalayas showed a significant decrease in grazing intensity (Fig. 6e), indicating the importance of protected areas in preventing overstocking and grassland degradation. Meanwhile, the GDGI maps also show that the grazing density varies greatly among protected areas, possibly owing to the difference in policy implementation. For instance, it can be seen from the GDGI maps that grazing intensity is increasing in some protected areas, especially in several wetland nature reserves on the Zoigê Plateau (Fig. 6e). Moreover, the average grazing intensity in all nature reserves on the QTP has increased overall from 1990 to 2020, although their increase rate is much lower than that of the non-protected areas (0.0125 SU ha⁻² ⋅ 10a vs. 0.0304 SU ha⁻² ⋅ 10a), which implies that grassland management in protected areas still needs to be strengthened on the QTP.

The grazing initiatives in alignment with the Sustainable Development Goals (SDGs) on the QTP can benefit from the GDGI dataset. Firstly, determining a reasonable stocking rate is vital to prevent overstocking of the pastures, which will possibly induce extreme grassland degradation (Q. Zhu et al., 2023). Stocking-rate determination can be optimized by using our grazing-intensity maps and the stocking-rate threshold maps of the QTP. Secondly, the GDGI maps can contribute to strategic placement of fence, which is a common practice adopted to prevent grassland degradation on the QTP. Building fences in areas with high grazing intensity and exceeding the carrying capacity can improve the effectiveness of fence construction (Zhou et al., 2023; Zhang et al., 2023). Thirdly, the GDGI dataset can provide a solid support for promoting effective nature reserve management, with nature reserves covering, in total, nearly one-third of the entire QTP. For example, the GDGI maps showed that grazing activities still exist in most nature reserves on the Plateau, although most of them have significantly lower grazing intensities compared with their adjacent non-protected areas. By using the GDGI maps, the conflict between ecological protection and grazing activities in nature reserves can be alleviated. Finally, our grazing-intensity maps can act as a basic dataset to support other grassland-related policies. Currently, these policies on the QTP often adopt a one-size-fits-all approach to determine the carrying capacity and to carry out ecological compensation, which may lead to overstocking or an unfair financial distribution (Wang et al., 2022). The grassland management strategies balancing carrying capacity and stocking rates are more likely to result in optimal management choices for policymakers and stakeholders, and our GDGI maps can contribute to this decision-making processes.

4.4 Uncertainties and limitations

Although this study has collected datasets that are as reliable as possible, users of the GDGI products should be cognizant of the inherent uncertainties and limitations within these datasets. Notably, the mean relative error of the GDGI dataset spanning 1990 to 2020 was recorded at 4.2 % (Fig. 4a), calculated from the average errors across 182 counties within the QTP that had accessible livestock census data. Furthermore, approximately 8.26 % of grassland areas exhibited a relative error exceeding 1.0 SU ha⁻² (Fig. 4b). Such discrepancies arise from several limitations that were subsequently propagated to the final grazing-intensity maps, thereby contributing to the dataset's overall uncertainties.

Firstly, the estimations of grazing intensities were fundamentally conservative, primarily due to the lack of comprehensive input data. Livestock numbers, derived from year-end data at the county level, inadvertently led to underestimations of grazing intensity by not accounting for livestock off-take rates. Likewise, the evaluation focused solely on livestock grazing intensity, excluding wild herbivores and forage-dependent livestock, which potentially results in an underestimation of actual grazing pressures on the QTP. Additionally, despite identifying seven main factors influencing livestock distribution, the study did not encompass all potential factors, such as fencing, forage availability, road proximity, and seasonal transformation in grazing practices. Moreover, to align with county-scale livestock census data, we averaged the environmental factors at the county scale. Although this approach has been widely used based on the hypothesis that a consistent causal relationship between livestock intensity and environmental factors persists across various scales (Robinson et al., 2014; Nicolas et al., 2016; Li et al., 2021; Meng et al., 2023), this might oversimplify the intricate dynamics between grazing intensity and lead to a certain degree of estimation inaccuracies. In addition, the reliance on linear extrapolation to supplementary missing gridded 100 m population density data from 1990–1999 introduced further uncertainties due to the limited resolution (1 km) and interval (5 year) of the ChinaPop dataset.

Secondly, the modeling process for mapping grazing intensity also suffered from several challenges. Specifically, this study adopted the FAO's assumption that the relationship between environmental factors and livestock intensity is uniform across both administrative and pixel levels. However, it is unlikely that these relationships are entirely consistent across scales, and the county-level model's approach inevitably smooths spatial details, potentially reducing the precision of the data. Furthermore, the ET model was trained with a limited sample size of 4998 and was applied to a vast area consisting of 150 million pixels, which could compromise the model's accuracy. In addition, despite the ET model's design being intended to reduce overfitting risks by using randomly selected features and partitioning decisions, the potential for overfitting effects still remains, particularly when faced with a high number of output classes or insufficient sample sizes (Geurts et al., 2006; Galelli and Castelletti, 2013). In fact, this limitation was evident in this study as the generalization capability of the ET model was restricted by the disparity between the number of training samples and the total number of pixels, leading to predictions that often exceed actual livestock census (Fig. 4a).

Thirdly, our methodological framework for high-resolution gridded grazing dataset mapping was developed based on the assumption that all grasslands were accessible to livestock. However, in reality, the amount of available grassland was less due to fencing and grazing bans on the QTP (Zhan et al., 2023). Moreover, transhumant herders generally follow a seasonal calendar for summer pastures and winter pastures on the QTP. However, we did not consider these seasonal movements due to data limitations, which further restrict the analysis of seasonal livestock distribution patterns (Kolluru et al., 2023). Additionally, the model's reliance on the human population as a proxy for livestock locations overlooked the possibility of high grazing intensity in areas with low human populations on the QTP, particularly in regions designated for summer pastures.

Finally, it is important to note that gathering livestock census data in the Qinghai–Tibet Plateau presents significant challenges, leading to a scarcity of livestock validation data in this study, particularly at the township and pixel scales. This limitation may, to some extent, impact the reliability of the grazing-intensity data we have presented.

In summary, all these limitations associated with input data, the modeling process, and the methodological framework collectively contribute to the uncertainties and reduce the accuracy of the GDGI maps. We henceforth recommend that future research should aim to incorporate more detailed data, consider additional influential factors, enhance key datasets' time series consistency, and refine the methodological framework to improve the accuracy of grazing-intensity mapping.

5 Data availability

The annual gridded grazing-intensity maps of the QTP spanning from 1990 to 2020 are accessible at the following link: https://doi.org/10.5281/zenodo.10851119 (Zhou et al., 2024). Each map is catalogued by year and recorded in GeoTIFF format, with values represented in $\mathrm{SU}{\mathrm{ha}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$. These datasets, with a spatial resolution of 100 m and annual temporal resolution, utilize the WGS 1984 Albers geographic coordinate system. To streamline data transfer and download processes, the comprehensive 31-year dataset has been compressed into a ZIP file, readily available for download and compatible with geographic information system (GIS) software for viewing.

6 Conclusions

In this study, we introduce a framework utilizing ET machine learning algorithms to achieve fine-scale livestock spatialization, subsequently generating the GDGI dataset across the QTP. The GDGI has a spatial resolution of 100 m and spans 31 years, from 1990 to 2020. It is consistent with county livestock census data of the QTP, and the accuracy evaluations at both the pixel level and township level underscore the outstanding reliability and applicability of the GDGI dataset, which can successfully capture the spatial heterogeneity and variations in grazing intensities in greater detail. Moreover, comparisons between the GDGI dataset and other existing grazing map products further proved the robustness and efficiency of our dataset and demonstrated the validity of the proposed framework in the research of livestock spatialization. Nonetheless, it is imperative for data users to recognize that the GDGI may still contain inherent uncertainties. Our Monte Carlo simulations have estimated the average MRE for grazing intensity across the QTP to vary between 6.84 % and 9.08 %. The GDGI dataset, as presented in this study, can enhance the understanding of grazing activities on the QTP. This, in turn, can aid in the rational utilization of grasslands and facilitate the implementation of informed and sustainable management practices.