*2.3. Data Processing*

#### 2.3.1. EVI Data Processing

EVI data was affected by factors such as aerosol, ice and snow, solar illumination angle, and sensor observation angle in the process of collection and processing. There may be anomalous or missing data, which caused the seasonal trend of EVI curve to be insignificant. Therefore, in order to make EVI time series data reflect seasonal variation of vegetation correctly, it was necessary to conduct a filtering process [9,37,38]. Firstly, the invalid pixels in the quality control file attached to MODIS13Q1 with a value of 65,535 were removed. With full consideration of the growth periodicity of vegetation, the 16d EVI data of IMAR were smoothed using harmonic analysis of time series (HANTS) (Figure 2). During processing, the valid range was −3000~10,000, the period was 23, and the frequency was 2 (11, 23). The time series data after reconstruction can reflect the periodic variation of the EVI curve [9,38]. Then, using the MODIS reprojection tools (MRT), the MOD13Q1 data were pre-georeferenced to the UTM zone 48 North projection WGS-84 datum resampled with the resolution of 500 m. Next, the monthly EVI data for vegetation in the growing season (May–September) were calculated using the method of maximum value composite (MVC). Subsequently, the effects of bare soil and sparsely vegetated areas were eliminated according to the following rules: (1) The annual mean value of EVI in the growing season was greater than 0.07; (2) The annual maximum value of EVI was greater than 0.10; (3) The annual maximum value of EVI should appear in July–September. Finally, the EVI pixels from May to September of the year 2000–2015, which met the requirements above, were used as the mean annual value of EVI in the growing season.

**Figure 2.** The EVI time series data before (**a**) and after (**b**) HANTS processing.

#### 2.3.2. Spatial Interpolation of Meteorological Data

Changes in meteorological sites, observation instruments and observation methods, or environmental changes in the surroundings will lead to non-uniformity of collected data. It was necessary to conduct homogeneity tests on the data of meteorological sites [39] to eliminate invalid site data. The existing research on meteorological factor interpolation adopts different methods—Kriging and IDW (Inverse Distance Weighted) [11,40] for relative humidity, Kriging and IDW [11,32] for precipitation., and TPS (Thin Plate Spline) and Kriging [40,41] for air temperature. Considering that the spatial distribution was significantly different in relative humidity, precipitation, and air temperature, and air temperature had a certain degree of altitude sensitivity [40–42], on the basis of the previous study method, after being tested by comparing meteorological site data and other methods, IDW was selected to interpolate the relative humidity, and Kriging to interpolate the precipitation. TPS was selected and DEM acted as the covariate to assist spatial interpolation of air temperature data. The spatial resolutions of meteorological element interpolation were unified to 500 m.
