## 2 Study Region, Methods and Datasets
### Datasets  
Data from multiple sources are utilized in this study, including ecosystem types, meteorological data, elevation data and soil data. Th ecosystem type data are interpreted from remote-sensing image data, mainly Landsat TM/ETM, GF-1 CCD, GF-2 CCD images from 2000, 2010 and 2020. The interpretation basis is image feature information such as the color tone, texture, geometric shape, etc. Ecosystems are classified into seven land types, including forests, shrubs, grasslands, wetlands, farmland, towns, and other land types. According to the classification of land space, forests, shrubs, wetlands, and others are merged into ecological space; farmland ecosystems such as paddy fields and dry lands are merged into agricultural space; and urban ecosystems such as cities, rural residential areas, and industrial and mining land are merged into urban space.  
The ecological parameter data used for the ecosystem quality are obtained from the website of the National Aeronautics and Space Administration (NASA) ([[http://modis.gsfc.nasa.gov/data/](http://modis.gsfc.nasa.gov/data/)]([http://modis.gsfc.nasa.gov/data/](http://modis.gsfc.nasa.gov/data/))) (accessed on 1 January 2022). The MODIS product is used in this research. The Fractional Vegetation Cover (FVC) is based on the pixel binary model inversion [48]. The Leaf Area Index (LAI) is derived from the inversion of the relationship model from preprocessed data, MOD09A1, to the LAI value, generated from an artificial neural network [49; 50]. The Gross Primary Productivity (GPP) is estimated based on the Bayesian multi-algorithm integration method, combining eight mainstream GPP algorithms. By checking the simulation accuracy of various light energy utilization models and using the observation data of 155 vortex-related carbon flux stations worldwide, the simulation accuracy of the eight light energy utilization models involved in the algorithm is calculated [51; 52; 53; 54]. The data utilized in this research have a resolution of 500 m, a time resolution of one year, and cover the period from 2000 to 2020.  
In addition, among other data required for the ecosystem services assessment, meteorological data are obtained from the China Meteorological science data-sharing service  
\begin{table}
\begin{tabular}{c c c c} \hline \hline \multirow{2}{*}{**Evaluation Contents**} & \multirow{2}{*}{**Ecosystem Quality**} & \multicolumn{2}{c}{**Ecosystem Services**} \\ \cline{3-4}  & & **Water Conservation** & **Soil Conservation** \\ \hline \multirow{3}{*}{Methods} & \multirow{2}{*}{EQI\({}_{\text{ij}}\) = \(\frac{\text{LAI}_{\text{ij}}+\text{FVC}_{\text{ij}}+\text{GPP}_{\text{i}}}{3} \times 100\%\)} & \multirow{2}{*}{Q\({}_{\text{str}}\) = \(\sum\limits_{i=1}^{n}\text{A}_{i}\times(\text{P}_{i}-\text{R}_{i}-\text{ET}_{i}) \times 10^{-3}\)} & \multirow{2}{*}{Q\({}_{\text{str}}\) = \(\sumnetwork; elevation data are obtained from the geospatial data cloud; and soil data are obtained from the 1:1,000,000 soil types of the Chinese Academy of Sciences [41]. The boundary of important ecological protection areas is obtained from the satellite environmental application center of the Ministry of Ecological Environment. Finally, the data are preprocessed by splicing and fusion, resampling, and unified coordinates.