Observed Vegetation Greening and Its Relationships with Cropland Changes and Climate in China

Chinese croplands have changed considerably over the past decades, but their impacts on the environment remain underexplored. Meanwhile, understanding the contributions of human activities to vegetation greenness has been attracting more attention but still needs to be improved. To address both issues, this study explored vegetation greening and its relationships with Chinese cropland changes and climate. Greenness trends were first identified from the normalized difference vegetation index and leaf area index from 1982–2015 using three trend detection algorithms. Boosted regression trees were then performed to explore underlying relationships between vegetation greening and cropland and climate predictors. The results showed the widespread greening in Chinese croplands but large discrepancies in greenness trends characterized by different metrics. Annual greenness trends in most Chinese croplands were more likely nonlinearly associated with climate compared with cropland changes, while cropland percentage only predominantly contributed to vegetation greening in the Sichuan Basin and its surrounding regions with leaf area index data and, in the Northeast China Plain, with vegetation index data. Results highlight both the differences in vegetation greenness using different indicators and further impacts on the nonlinear relationships with cropland and climate, which have been largely ignored in previous studies.


Introduction
Since the 1980s, China's croplands have experienced significant changes due to policy-driven urbanization and ecological restoration (e.g., the Grain to Green Project). This has resulted in a loss of croplands and a consequent threat to food security [1,2]. Therefore, to ensure food security and sustainability for the rising population, croplands should be protected, and food production must grow substantially. However, cropland increases either in China or the other countries in the world, could conversely lead to various environmental problems such as biodiversity loss, land degradation, or pest outbreaks [3][4][5][6].
With recent advances in data availability, our ability to evaluate the impacts of cropland changes improved. Many studies used various datasets including historical survey data, reconstructed datasets, and satellite data to monitor cropland changes and inferred possible reasons driving cropland changes [7][8][9][10]. Some studies evaluated the impacts of cropland changes on production. For example, Liu et al. [11] suggested that urbanization and the Grain to Green Project (which reverted cropland to Land 2020, 9, 274 2 of 19 forests and grasslands) might lead to a decrease in the production from 1990 to 2010. Qin et al. [12] found that ecological restoration projects influenced cropland productivity in China from 2000 to 2008 but that their effects were less than those of urbanization. In contrast, the potential consequences of cropland changes on ecological aspects remain underexplored.
Currently, in the ecology-related fields, vegetation greenness is widely used to monitor vegetation status and assess the improvement or degradation of ecological environments, and thus an important indicator to evaluate the ecological impacts of cropland changes. Over the past decades, satellite remote sensing has been the only possible means of monitoring vegetation greenness on a large scale [13][14][15]. Increasing evidence from long-term satellite-derived leaf area index (LAI) or normalized difference vegetation index (NDVI) records showed that vegetation greenness increased in the northern middle and high latitudes over the past three to four decades [16][17][18]. The greening has been primarily interpreted as increased photosynthesis due to climate change. Air temperature was considered to be the main factor driving the greening trend [19,20], and atmospheric CO 2 accounted for a significant increase in the LAI of the Northern Hemisphere [21]. Despite the overall greening trend, browning trends have also been observed in some regions, probably due to the occurrence of droughts [22,23]. In addition to climatic factors, human activities (e.g., land use) have been found to contribute to changes in vegetation greenness [24][25][26][27]. However, researchers have generally attributed changes in vegetation greenness to climatic factors, assuming that little or no land use and land cover changes (LULCC) occurred in the areas studied [28,29]. This assumption is illogical, particularly in regions that have experienced significant land-cover changes [30]. In China, some studies have indicated that LULCC might be a dominant driver of vegetation dynamics but have often highlighted the significant role of China's forestry programs in explaining the spatial patterns of vegetation growth or focused on the vegetation greenness or growth in grasslands [24,[31][32][33][34], while the contribution of cropland and its changes to vegetation greening has been largely ignored.
To address both issues, this study aims to explore the relative impacts of human-induced cropland changes on greenness trends in China. To the best of our knowledge, there have been no studies that linked cropland changes with greenness trends, except for the exploration of the relationship between cropland changes in Sahel and the NDVI trends [35]. A comprehensive exploration of the linkage between cropland changes and greenness trends is still lacking for China. Therefore, the first objective of this study is to investigate the possible associations between cropland changes in China and vegetation greening or browning, compared with the impacts of climate on vegetation greening or browning.
To achieve the goal, it is essential to accurately detect greening or browning trends. Many algorithms have been developed to detect trends of greenness over large areas from remote sensing data. Most studies have adopted a monotonic linear regression approach to detect vegetation greenness trends at regional and global scales [15,36]. However, greenness trends caused by human activities or other local factors might be abrupt, and non-monotonic trend detection methods are thus needed in these scenarios. Jamali et al. [37] proposed the Detecting Breakpoints and Estimating Segments in Trend approach and applied it to detect major changes as well as the types (abrupt or non-abrupt), time, and magnitude of these changes in Iraq from time series of vegetation indices. Verbesselt et al. [38] developed the BFAST (Break detection For Additive Season and Trend) algorithm to detect disturbances in near real time. Li et al. [39] used piecewise linear regression to detect trends of vegetation greenness in the Tibetan Plateau from 1982 to 2006. Although various monotonic and non-monotonic methods were developed for trend detection analysis, few studies have compared the impacts of different algorithms on the identified greenness trends. Therefore, another objective of this study is to compare the performances of different algorithms, both linear and nonlinear, on the greenness trend results. Moreover, note that previous studies have identified trends of vegetation greenness from satellite-derived NDVI datasets using developed algorithms [15,16,29,40]. Since NDVI could be affected by environmental factors [41], some studies employed LAI to measure vegetation greenness, but did not give enough consideration to the potential discrepancies in detecting greenness Land 2020, 9, 274 3 of 19 using different metrics. [31,32]. However, in this study, both NDVI and LAI metrics were employed to detect trends of vegetation greenness (greening or browning) in China. This would enable a more robust identification of changes in the vegetation conditions and a comparison of trend results based on different metrics. Figure 1 showed the datasets, pre-processing of these datasets, trend detection techniques and statistical analysis methods for evaluating the impacts of cropland and climate on vegetation greenness. Cropland-related datasets, satellite-derived Global Land Surface Satellite (GLASS) NDVI and LAI datasets, and climate data were the main datasets used in this study. Section 2.1 described the cropland-related datasets and the pre-processing of cropland datasets to extract cropland percentage and cropland changes from 1980 to 2015. Section 2.2 introduced the satellite-derived GLASS NDVI and LAI datasets and the aggregation of these datasets from 8-day temporal resolution to the annual scale. Section 2.3 included the climate data and the pre-processing of climate data to obtain multi-year annual temperature and precipitation data, and changes in annual temperature and precipitation. Section 2.4 presented the analysis methods, including the simple linear regression, Theil-Sen slope estimator, and the polynomial fitting-based algorithm (PolyTrend) to identify greenness trends from satellite-derived NDVI and LAI data, comparison of spatial similarity of greenness trends based on GLASS NDVI and GLASS LAI using three detection techniques, and the boosted regression trees (BRT) model to evaluate the relative contributions of cropland information in Section 2.1 and climate information in Section 2.3 to vegetation greening detected from GLASS NDVI and LAI datasets.

Materials and Methods
Land 2020, 9, x FOR PEER REVIEW 3 of 20 NDVI and LAI metrics were employed to detect trends of vegetation greenness (greening or browning) in China. This would enable a more robust identification of changes in the vegetation conditions and a comparison of trend results based on different metrics. Figure 1 showed the datasets, pre-processing of these datasets, trend detection techniques and statistical analysis methods for evaluating the impacts of cropland and climate on vegetation greenness. Cropland-related datasets, satellite-derived Global Land Surface Satellite (GLASS) NDVI and LAI datasets, and climate data were the main datasets used in this study. Section 2.1 described the cropland-related datasets and the pre-processing of cropland datasets to extract cropland percentage and cropland changes from 1980 to 2015. Section 2.2 introduced the satellite-derived GLASS NDVI and LAI datasets and the aggregation of these datasets from 8-day temporal resolution to the annual scale. Section 2.3 included the climate data and the pre-processing of climate data to obtain multi-year annual temperature and precipitation data, and changes in annual temperature and precipitation. Section 2.4 presented the analysis methods, including the simple linear regression, Theil-Sen slope estimator, and the polynomial fitting-based algorithm (PolyTrend) to identify greenness trends from satellite-derived NDVI and LAI data, comparison of spatial similarity of greenness trends based on GLASS NDVI and GLASS LAI using three detection techniques, and the boosted regression trees (BRT) model to evaluate the relative contributions of cropland information in section 2.1 and climate information in section 2.3 to vegetation greening detected from GLASS NDVI and LAI datasets.

Chinese Cropland and Its Changes
China is a large country with nine major agricultural zones (http://www.resdc.cn/data.aspx? DATAID=275, Figure 2b). The cropland area was estimated to be, in total, 141.1 million hectares in China in 2000, with uneven spatial distributions [42]. Large cropland areas are generally located in provinces with high population density. Major agricultural zones include the Northeast Plain, the Huang-Huai-Hai Plain, Sichuan Basin, and the Yangtze River Basin, whereas cropland covers in the plateaus and hills such as the Yunnan-Guizhou, Tibet, and Loess Plateaus are relatively lower [43,44]. Additionally, during the past decades, China has experienced significant climate changes that also affect the cropping system [45]. The largest increase in temperature occurred in northeast China, and regional shifts in precipitation were also observed, with a decrease in summer-autumn precipitation in northeastern areas and increases in summer and winter precipitation in southeastern areas [46].
Several datasets have been published to describe the spatial distribution of croplands in China; however, comparison studies have suggested that their consistency is poor [43,[47][48][49]. In this study, the Chinese land-use/land-cover datasets (CLUDs) were selected to characterize Chinese cropland changes due to their high accuracy and long-time series [50]. The CLUDs have been available every 5 years since 1980. To build each CLUD, over 500 Landsat Thematic Mapper images were georeferenced, orthorectified, and classified into 25 land-use/land-cover categories at a spatial scale of 1:100,000, then further converted to the 1 km gridded data by calculating area percentages for each land use category within every cell [42,51,52]. The 1-km gridded CLUDs for 1980,1990,1995,2000,2005,2010, and 2015 were used in this study. They were downloaded from the Data Center for Resources and Environmental Sciences, Chinese Academy of Sciences (http://www.resdc.cn/). The CLUDs adopt a hierarchical classification system. The 25 land-cover categories were grouped into six aggregated classes: cropland, woodland, grassland, water bodies, unused land, and built-up areas [42]. Since we focus on cropland changes, each CLUD was first reprocessed into 2 broad categories, cropland and non-cropland, labelled 1 and 0, respectively. Cropland changes during two consecutive periods were quantified by subtracting the cropland dataset for the former period from that for the latter period. Using this method, cropland changes for 1980-1990, 1990-1995, 1995-2000, 2000-2005, 2005-2010, and 2010-2015 were quantified.
Although the accuracy of six classes of land use with cropland included was above 94.3% when validated with field survey materials [50], the potential misclassification of cropland as non-cropland, or non-cropland as cropland, in the CLUDs could lead to spurious land-use differences and inaccurately identified cropland changes. To improve this situation, we assumed that local pixels or areas could not experience conversions to other land-use types and then be immediately transformed to their original land-use types; thus, pixels labelled as 1-0-1 and 0-1-0 in 3 consecutive periods suggested their misclassification as non-cropland and cropland, respectively, in the second period. The potentially misclassified pixels in CLUDs for 1990CLUDs for , 1995CLUDs for , 2000CLUDs for , 2005, and 2010 would be identified and masked out in the analysis, which could improve the accuracy of cropland changes during the periods 1980-1990, 1990-2000, 2000-2010, and 2010-2015. Moreover, the accuracy would be improved for cropland changes for the entire 1980-2015 period quantified by the sum of cropland changes for 1980-1990, 1990-2000, 2000-2010, and 2010-2015. To match the satellite-derived greenness datasets described in Section 2.2, cropland changes from 1980 to 2015 were then averaged to a 0.05 • spatial resolution.
We also created the cropland masked dataset for the period 1980-2015. Pixels classified as cropland at least twice in the land-use datasets in 1980,1990,1995,2000,2005,2010, and 2015 were taken as cropland pixels and labelled as 1; otherwise, they were labelled as 0. The cropland masked data were then aggregated to 0.05 • spatial resolution. Pixels with a cropland percentage of no less than 4% at the 0.05 • resolution were classified as cropland, whereas other pixels were classified as non-cropland and masked in the subsequent analysis.

Vegetation Greenness
Greenness trends were identified from both satellite derived NDVI and LAI products in this study. Datasets used to detect trends of vegetation greenness (greening or browning) in Chinese croplands include the GLASS Advanced Very High Resolution Radiometer (AVHRR) NDVI [53] and the GLASS LAI dataset [54,55]. The GLASS products have longer time series and are thus highly suitable for monitoring environmental changes [56]. They can be downloaded from the National Earth System Science Data Sharing Infrastructure (http://www.geodata.cn) as well as from the University of Maryland Global Land Cover Facility (http://glcf.umd.edu).
The GLASS NDVI and LAI data from 1982 to 2015, with a temporal resolution of 8 days and a spatial resolution of 0.05 • , were used in this study. The 8-day datasets were first aggregated to a monthly scale and then aggregated to an annual scale by averaging the monthly NDVI values no less than 0.1 to minimize the influence of bare or sparsely vegetated pixels on greenness trends [57]. Since one of the objectives is to compare the performances of NDVI and LAI in monitoring greenness, the annual LAI values within a year were obtained by averaging the monthly LAI values with the corresponding monthly NDVI values of no less than 0.1 to ensure that annual NDVI and LAI described the cropland greenness covering the same period.

Climate Data
To fully investigate the linkage between cropland changes and vegetation greening, climatic factors were also included in the following contribution analysis of factors driving cropland greenness. Annual mean temperature and precipitation data were used in this study. The annual mean temperature and precipitation data were downloaded from the Data Center for Resources and Environmental Sciences, Chinese Academy of Sciences (http://www.resdc.cn). They were originally interpolated from meteorological data using ANUSPLIN software and were available on a 1 × 1 km grid. We transformed the gridded dataset from 1 km to a resolution of 0.05 • by re-projection and then averaged the multiyear average annual temperature and precipitation data from 1982 to 2015. Changes in annual temperature and precipitation from 1982 to 2015 were derived using a simple linear regression method.

Analysis Methods
To investigate the relationship between cropland changes and vegetation greenness trends, we first identified trends of annual mean GLASS NDVI and LAI from 1982 to 2015 in Chinese croplands by estimating the slope coefficients of simple linear regression models, with a significance level of 0.05 [58]. We also included the non-parametric Theil-Sen slope estimator and a nonlinear PolyTrend algorithm to identify vegetation greening or browning in China. The Theil-Sen slope is a non-parametric trend estimator that is insensitive to noise and outliers in time series data. Similar to linear regression, the Theil-Sen median slope estimator is limited in detecting abrupt or non-monotonic changes in vegetation greenness trends caused by land use changes [39]. The PolyTrend could identify the possible non-monotonic trends of vegetation greenness and classify the identified greenness trends into linear, quadratic, cubic, concealed, and no-trend types [59].
We compared the similarities and differences in the magnitude and spatial distribution of greenness trends detected using the above three techniques. The percentage of cropland pixels that experienced greening and browning as detected from GLASS NDVI and LAI by simple linear regression, the Theil-Sen slope estimator, and the PolyTrend algorithm were first statistically analyzed in the major agricultural zones of China. Both the linear and nonlinear trends detected by PolyTrend were considered. The corresponding greenness trends using PolyTrend algorithm comprised the positive linear, quadratic, and cubic trends, and the browning trends comprised the negative linear, quadratic, and cubic trends. Then the similarity of spatial distributions of greenness trends detected using the Land 2020, 9, 274 6 of 19 above 3 trend detection techniques was assessed using fuzzy numerical (FN) maps [60]. FN maps represent the numerical similarity of detected trends at the pixel level, and were obtained as follows: where the cell values (a and b) represent the trends of greenness detected from two different datasets or with two different methods, and s represents the numerical similarity between a and b. The numerical similarity ranges from 0 (fully distinct trends) to 1 (fully identical trends). Using Equation (1), we calculated the FN maps that measured the similarity of LAI trends and NDVI trends detected using simple linear regression, the similarity of LAI trends and NDVI trends using Theil-Sen slope estimator, the similarity of NDVI trends detected by simple linear regression and NDVI trends detected by Theil-Sen slope estimator, and the similarity of LAI trends detected by simple linear regression and LAI trends detected by Theil-Sen slope estimator, in terms of spatial distribution. The former two FN maps represented the similarity of greenness trends detected from different datasets (LAI and NDVI), and the latter two represented the similarity of greenness trends detected using different techniques (simple linear regression and Theil-Sen slope estimator). The distribution of FN maps in nine major agricultural zones were graphically depicted using boxplots that were based on a five-number summary including the minimum, the maximum, the median, the first and third quantiles. On each box, the central mark indicates the median, and the bottom and top edges of the box indicate the 25th and 75th percentiles of FN values corresponding to an agricultural zone, respectively.
Finally, the BRT model was used to evaluate the relative importance of climate and cropland changes on vegetation greenness. The BRT model combines the strengths of regression trees and boosting technique, and can provide high prediction accuracies and good explanation of the underlying relationships between response and predictors, compared with other data mining methods [61][62][63][64]. In this study, we developed the BRF model for each major agricultural zone separately using the gbm library in R, to explore the relationships of cropland and climate with vegetation greening. Predictor variables of a BRF model included the cropland percentage, cropland changes, mean annual temperature and precipitation from 1982 to 2015, and annual changes in temperature and precipitation from 1982 to 2015. The fitting of a BRF model was controlled by the learning rate which determined the contribution of each tree to the growing model, the number of trees to be grown in the sequence, and tree complexity [65]. Therefore, we first estimated the optimal learning rate, number of trees, and tree complexity with cross-validation, and then established the BRF model with optimal meta-parameters.
The relative contributions of cropland change and climate to vegetation greening in a BRT model were assessed based on the number of times that a variable was selected for splitting the BRT, weighted by the squared improvement to the model as a result of each split, and averaged over all of the trees [65,66]. The most important variable per BRT model was selected and the corresponding partial dependency plot was used to describe the relationships between vegetation greening and the most important predictor variable while controlling for the average effect of all other predictors [63,66]. Figure 2 shows the spatial distribution of average cropland percentage from 1980 to 2015 and its changes during this period based on CLUDs at the 0.05 • resolution. Cropland increases occurred in most regions of northern China, whereas decreases in cropland were mainly located in Northwest and South China. Summarizing the cropland changes in the nine major agricultural zones, we found that northern arid and semi-arid regions experienced the highest cropland losses during 1980-2015, followed by the middle-lower Yangtze Plain region. More than 8% of croplands in both agricultural zones have been lost. Similar phenomena were also found in Southern China. As suggested by previous studies, the decreases in cropland areas in the middle-lower Yangtze Plain region and in Southern China were mainly due to the conversion of croplands to urban areas, whereas the losses of croplands in Northwest China were related to the "Grain for Green" policy [50]. The significant cropland increases in northern China could be largely attributed to the reclamation of grassland for cropland and additionally by climate change [10,50,67].

Vegetation Greenness Trends in Chinese Croplands
Annual changes in GLASS NDVI and LAI from 1982 to 2015 detected by simple linear regression, the Theil-Sen slope estimator, and the PolyTrend algorithm all revealed widespread cropland greening, particularly in central and southern China (Figures 3 and 4). Identified browning trends were mainly located in urbanized regions, such as the Yangtze River Basin, the Pearl River Delta region, and the Sichuan Basin (Figures 3 and 4).

Vegetation Greenness Trends in Chinese Croplands
Annual changes in GLASS NDVI and LAI from 1982 to 2015 detected by simple linear regression, the Theil-Sen slope estimator, and the PolyTrend algorithm all revealed widespread cropland greening, particularly in central and southern China (Figures 3 and 4). Identified browning trends were mainly located in urbanized regions, such as the Yangtze River Basin, the Pearl River Delta region, and the Sichuan Basin (Figures 3 and 4).  Greenness trends detected by the Theil-Sen slope estimator were basically the same as those detected by simple linear regression algorithms and thus were not shown in this study. The PolyTrend algorithm identified both linear and nonlinear trends of time series LAI and NDVI data (Figure 4). The monotonic linear trends accounted for at least 50% of the significant trends detected. The three trend detection methods produced similar results in terms of the proportions of cropland pixels that experienced greening or browning in each major agricultural zone ( Figure 5) and consistent spatial patterns of greenness trends suggested by higher FN indices of the majority of pixels (Figure 6b). All these results showed that greening or browning trends derived from satellite data time series were basically independent of the trend detection algorithm. Greenness trends detected by the Theil-Sen slope estimator were basically the same as those detected by simple linear regression algorithms and thus were not shown in this study. The PolyTrend algorithm identified both linear and nonlinear trends of time series LAI and NDVI data ( Figure 4). The monotonic linear trends accounted for at least 50% of the significant trends detected. The three trend detection methods produced similar results in terms of the proportions of cropland pixels that experienced greening or browning in each major agricultural zone ( Figure 5) and consistent spatial patterns of greenness trends suggested by higher FN indices of the majority of pixels ( Figure 6b). All these results showed that greening or browning trends derived from satellite data time series were basically independent of the trend detection algorithm.    However, we found that greenness trends in Chinese croplands were greatly affected by the selected proxies for greenness. More pixels in cropland areas were detected to exhibit greening trends characterized by NDVI, compared with those detected by LAI (Figure 7). At least 54% of cropland areas in all agricultural zones were found to have experienced increasing NDVI trends, and the proportion reached above 90% in the Loess Plateau and the Yunnan-Guizhou Plateau. The proportion of croplands that exhibited increasing LAI trends in nine major agricultural zones ranged from 16.93% in the Northeast China Plain to 71.15% in the Loess Plateau. In Southern China croplands, the proportion of browning pixels was 13.1% for GLASS LAI and 7.76% for GLASS NDVI, respectively, slightly higher than in other agricultural zones. Vegetation browning was also found in the middlelower Yangtze Plain croplands, and the corresponding browning proportions indicated by GLASS LAI and NDVI were 7.45 and 6.88%, respectively. In addition to the large discrepancies in the However, we found that greenness trends in Chinese croplands were greatly affected by the selected proxies for greenness. More pixels in cropland areas were detected to exhibit greening trends characterized by NDVI, compared with those detected by LAI (Figure 7). At least 54% of cropland areas in all agricultural zones were found to have experienced increasing NDVI trends, and the proportion reached above 90% in the Loess Plateau and the Yunnan-Guizhou Plateau. The proportion of croplands that exhibited increasing LAI trends in nine major agricultural zones ranged from 16.93% in the Northeast China Plain to 71.15% in the Loess Plateau. In Southern China croplands, the proportion of browning pixels was 13.1% for GLASS LAI and 7.76% for GLASS NDVI, respectively, slightly higher than in other agricultural zones. Vegetation browning was also found in the middle-lower Yangtze Plain croplands, and the corresponding browning proportions indicated by GLASS LAI and NDVI were 7.45 and 6.88%, respectively. In addition to the large discrepancies in the proportions of observed greening and browning pixels derived from GLASS LAI and GLASS NDVI datasets, the spatial distributions of greenness change also differed (Figures 3 and 6). Cropland greenness calculated from GLASS LAI and NDVI data even changed in opposite directions in the Northeast China Plain, where croplands probably experienced nonlinear trends (e.g., quadratic, and cubic, or even concealed) as suggested by the polynomial fitting-based algorithm results (Figure 4). proportions of observed greening and browning pixels derived from GLASS LAI and GLASS NDVI datasets, the spatial distributions of greenness change also differed (Figures 3 and 6). Cropland greenness calculated from GLASS LAI and NDVI data even changed in opposite directions in the Northeast China Plain, where croplands probably experienced nonlinear trends (e.g., quadratic, and cubic, or even concealed) as suggested by the polynomial fitting-based algorithm results (Figure 4).

Relationship between Cropland Change and Vegetation Greenness Trends
As was shown in Section 3.2, greenness trends were less sensitive to the trend detection algorithms. For simplicity, the following revealed relationships between cropland changes and greenness trends were all based on trends calculated using the simple linear regression method.
The results of BRT models showed the diverse contributions of cropland changes and climate to greenness trends in major Chinese agricultural zones and, in addition, the dominant factors driving greenness changes differed when different metrics were used to describe vegetation greening. In the Loess Plateau and Qinghai Tibet Plateau, greenness changes were predominantly affected by temperature with LAI data, but by precipitation with NDVI data (Figures 8 and 9). However, with either LAI or NDVI data, greenness changes were predominantly affected by temperature as opposed to cropland changes in the Huang-Huai-Hai Plain, northern arid and semi-arid regions, and Southern

Relationship between Cropland Change and Vegetation Greenness Trends
As was shown in Section 3.2, greenness trends were less sensitive to the trend detection algorithms. For simplicity, the following revealed relationships between cropland changes and greenness trends were all based on trends calculated using the simple linear regression method.
The results of BRT models showed the diverse contributions of cropland changes and climate to greenness trends in major Chinese agricultural zones and, in addition, the dominant factors driving greenness changes differed when different metrics were used to describe vegetation greening. In the Loess Plateau and Qinghai Tibet Plateau, greenness changes were predominantly affected Land 2020, 9, 274 11 of 19 by temperature with LAI data, but by precipitation with NDVI data (Figures 8 and 9). However, with either LAI or NDVI data, greenness changes were predominantly affected by temperature as opposed to cropland changes in the Huang-Huai-Hai Plain, northern arid and semi-arid regions, and Southern China. The relative contribution of temperature in these regions ranged from 40.47% in Huang-Huai-Hai Plain to 48.59% in northern arid and semi-arid regions with LAI data, and ranged from 40.30% in Southern China to 62.35% in Huang-Huai-Hai Plain with NDVI data. While in the middle-lower Yangtze Plain, temperature change played a dominant role in vegetation greening, contributing to 41.95% of greenness changes characterized by LAI data and 42.66% with NDVI data. Precipitation contributed the most to vegetation greening characterized by either LAI or NDVI data in the Yunnan-Guizhou Plateau. However, using LAI to describe the greenness changes, cropland percentage had a large relative contribution of 26.66%, slightly lower than the 27.76% contributed by precipitation. Land 2020, 9, x FOR PEER REVIEW 12 of 20   Compared with climate information, cropland could be more closely linked to vegetation greening in the Sichuan Basin and its surrounding regions with LAI data and in the Northeast China Plain with NDVI data. In the Sichuan Basin and its surrounding regions, cropland percentage accounted for 32.73% of vegetation greening characterized by LAI data, followed by precipitation with a relative contribution of 27.08% and temperature with a relative contribution of 21.13%. While in the Northeast China Plain, cropland percentage accounted for 28.91% of vegetation greening with NDVI data, and relative contributions of precipitation and temperature to vegetation greening were 26.60% and 22.85%, respectively.
In all the major agricultural regions, the impacts of precipitation and cropland change were not strong, and their relative contributions were no larger 15.71% and 9.79%, respectively.   The partial dependence plots showed the nonlinear relationships between vegetation greening and its most important predictor in Chinese agricultural zones (Figures 10 and 11). Vegetation greenness increased almost linearly with the increase of temperature when temperature was within 0 to 1 • C in the northern arid and semi-arid regions. Significant non-monotonic relationships between vegetation greening and temperature were observed in the Huang-Huai-Hai Plain and Southern China. Vegetation greening increased very rapidly as temperature increased, kept stable when it reached 0.5 • C, and changed dramatically when it exceeded 1.0 • C in the Huang-Huai-Hai Plain. While in Southern China, vegetation greening was basically constant when temperature was less than 1.5 • C, and then slowly decreased when temperature was within 1.5 to 2.0 • C. Precipitation was negatively associated with vegetation greening in the Yunnan-Guizhou Plateau, particularly when precipitation values were between 100 and 200 mm. While in Qinghai Tibet Plateau, vegetation greenness characterized by either LAI or NDVI data, basically kept stable when precipitation was within about 40-80 mm. In the Loess Plateau, with the increase in temperature, greenness trends from LAI data were overall increased, and with the increase in precipitation, NDVI trends were increased. Cropland was more likely to promote vegetation greening in Sichuan Basin and its surrounding regions ( Figure 10). However, in Northeast China Plain, greenness decreased with an increase in cropland percentage (Figure 11), which suggested that cropland to other land cover types (e.g., forests) more likely led to vegetation greening [68].
values were between 100 and 200 mm. While in Qinghai Tibet Plateau, vegetation greenness characterized by either LAI or NDVI data, basically kept stable when precipitation was within about 40-80mm. In the Loess Plateau, with the increase in temperature, greenness trends from LAI data were overall increased, and with the increase in precipitation, NDVI trends were increased. Cropland was more likely to promote vegetation greening in Sichuan Basin and its surrounding regions ( Figure  10). However, in Northeast China Plain, greenness decreased with an increase in cropland percentage (Figure 11), which suggested that cropland to other land cover types (e.g., forests) more likely led to vegetation greening [68].

Uncertainties in Cropland Changes
Previous studies have explored cropland changes in China but with inconsistent results [9,69]. One primary reason for the uncertainties in estimates of cropland areas and cropland changes was that data sources used to quantify cropland related information were multiple and inconsistent. In

Uncertainties in Cropland Changes
Previous studies have explored cropland changes in China but with inconsistent results [9,69]. One primary reason for the uncertainties in estimates of cropland areas and cropland changes was that data sources used to quantify cropland related information were multiple and inconsistent. In this study, satellite data-derived CLUDs, which have been widely validated and identified as reliable datasets to depict land cover status, were used to extract cropland and its changes in China. However, the 1-km grid land cover dataset was inferior to continuous land cover datasets that could report the percentage of each grid cell that was cropland [30].
Since the 1-km grid land cover datasets were publicly available from 1980 to 2015 while the 1-km continuous versions were publicly available only from 1990 to 2005, we used the discrete rather than continuous land cover datasets in this study. When aggregating both types of datasets to the 0.05 • resolution, we found that there were no significant differences in cropland changes extracted from the two types of land cover datasets from 1990 to 2005. Therefore, the 1-km grid land cover datasets used to characterize cropland changes could be reliable.

Metrics to Depict Vegetation Greenness
The NDVI is a measure of chlorophyll abundance in leaves and is widely used as an indicator of vegetation coverage and greenness [70]. Existing studies have used NDVI as a proxy for several key biophysical properties, such as LAI, fraction of absorbed photosynthetically active radiation (FAPR), biomass, and net primary production (NPP) [16,71,72], and have monitored trends in vegetation activities using NDVI. Moreover, previous studies have often emphasized the strong correlations among these biophysical variables and thus selected one of them as the greenness indicator to evaluate vegetation growth [73][74][75]. The inherent discrepancies in using different biophysical variables to monitor vegetation greenness have been largely ignored.
Results of this study showed that greenness trends as characterized by LAI and NDVI were inconsistent, or even opposite in some regions where cropland changes occurred, although widespread increasing trends were observed by both metrics in major agricultural zones. This inconsistency indicated the discrepancies of various biophysical variables in measuring vegetation greenness, rather than the strong correlations among them reported by previous studies. Moreover, this study revealed that the choice of metrics and the corresponding datasets predominantly affected the magnitude and spatial distribution of greenness trends, which were basically insensitive to the trend detection methods employed. The importance of biophysical variables was highlighted in this study, and more biophysical variables, such as biomass, NPP, fractional vegetation cover, FPAR, or even phenological metrics derived from vegetation indices, could be incorporated to comprehensively assess vegetation greening from different perspectives in future studies [76]. Moreover, since diverse indicators provide vegetation greening from different perspectives, complex indices combining these variables could be developed to improve the evaluation of vegetation greenness and the impacts of anthropogenic activities on greenness in the future.

Implications of the Relationship between Vegetation Greenness and Cropland Changes as well as Climate
Cropland changes in China and vegetation greening are two topics that have been receiving much attention from scientists and policymakers. Given the above ambiguity about the environmental consequences of Chinese cropland changes as well as the impacts of a specific LULCC on vegetation greening or browning, we addressed both issues by exploring the possible impacts of human-induced cropland change as well as climate on greenness trends in China.
We found that the relative contribution of cropland to greenness trends from 1982 to 2015 characterized by LAI was 32.73% in the Sichuan Basin and its surrounding regions, and 28.91% in the Northeast China Plain with NDVI data, both of which were slightly higher than the relative contributions of temperature and precipitation. However, in most Chinese agricultural zones, greenness changes were predominantly affected by temperature as opposed to cropland percentage and cropland changes. Previous studies also suggested that in most parts of China, vegetation greenness is more sensitive to temperature rather than precipitation [77,78], which was consistent with our findings. Hua et al. [24] found that in southern China, the relationship between NDVI and temperature was positive, suggesting that warmer climate could enhance vegetation activity. While in our study, the results showed that temperature played a dominant role in vegetation greening, but the relationships between NDVI and temperature in southern China were nonlinear, and thus warmer climate might not lead to the increase in vegetation activity. Liu et al. [63] found the nonlinear relationships of vegetation greening in southwest China with climate and human factors, in agreement with our results. Li et al. [79] reported that both climate variability and human activities contributed to vegetation change in Loess Plateau, but human activities account for 55% of vegetation changes over the 2000-2015 period. However, results of this study showed that the relative contribution of precipitation to vegetation greening accounted for more than 60% for the period 1982-2015.
Additionally, the results of this study provided observational evidence and improved our understanding of the relative contributions of cropland changes and climate to vegetation greening in China. Since croplands in China have generally decreased since 2000 and may continually decrease under the influence of constructing an ecological security barrier, civilization, and urbanization [11], the results of this study are significant for assessing the environmental impacts or issues associated with cropland related programs in China.

Conclusions
For the first time, this study comprehensively explored the relationships between Chinese cropland changes and vegetation greening or browning. We initially detected greenness trends in Chinese cropland areas from satellite derived NDVI and LAI datasets with three trend detection techniques and then used boosted regression trees to explore underlying relationships of cropland changes and climate with observed vegetation greening or browning from satellite data. The results showed that the greenness trends detected from different datasets could be quite different, and that the nonlinear relationships between cropland changes and climate and vegetation greening varied with the metrics used. However, using either satellite-derived NDVI or LAI datasets, climate predominantly contributed to vegetation greening in most Chinese agricultural zones, and the relative importance of temperature accounted for more than 40.30% in the Huang-Huai-Hai Plain, northern arid and semi-arid regions, and Southern China. Cropland percentage was the most important predictor of vegetation greening in the Sichuan Basin and its surrounding regions with LAI data and in the Northeast China Plain with NDVI data, and their relative contributions were 32.73% and 28.91%, respectively, while the impacts of cropland changes on vegetation greening in Chinese agricultural zones were not as strong as climate factors.
For future studies, our results highlight the necessity of carefully examining the applicability of biophysical variables and associated datasets for monitoring environmental change and exploring its underlying driving factors. However, it should be noted that this investigation was an exploration of the direct impacts of cropland changes and percentage on greening or browning based on satellite observations, without considering indirect effects such as pest outbreaks or agricultural management practices (e.g., irrigation). Both direct and indirect factors that drive vegetation greening or browning should be incorporated in future research.