Detecting Vegetation Change in the Pearl River Delta Region Based on Time Series Segmentation and Residual Trend Analysis (TSS-RESTREND) and MODIS NDVI

: Time Series Segmentation and Residual Trend analysis (TSS-RESTREND) can detect an abrupt change that was undetected by Residual Trend analysis (RESTREND), but it is usually combined with the Global Inventory for Mapping and Modeling Studies (GIMMS) Normalized Di ﬀ erence Vegetation Index (NDVI), which cannot detect detailed vegetation changes in small areas. Hence, we used Time Series Segmentation and Residual Trend analysis (TSS-RESTREND) and Moderate Resolution Imaging Spectroradiometer (MODIS) NDVI (MOD-TR) to analyze the vegetation dynamic of the Pearl River Delta region (PRD) in this study. To choose the most suitable MODIS NDVI from MOD13Q1 (250 m), MOD13A1 (500 m), and MOD13A2 (1 km), whole and local comparison of results of the break year and MOD-TR were used. Meanwhile, a comparison of vegetation change at the city-scale was also implemented. Moreover, to reduce insigniﬁcant trend pixels in TSS-RESTREND, a combination method of TSS-RESTREND and RESTREND (CTSS-RESTREND) was proposed. We found that: (1) MOD13Q1 and MOD13A1 two NDVI were suitable for combination with TSS-RESTREND to detect vegetation change in PRD, but MOD13Q1 was a better choice when considering the accuracy of local detailed vegetation change; (2) CTSS-RESTREND could detect more pixels with a signiﬁcant change (i.e., signiﬁcant increase and signiﬁcant decrease) than those of TSS-RESTREND and RESTREND. Also, its e ﬀ ectiveness could be veriﬁed by Landsat data; (3) at the city-scale, the CTSS-RESTREND detected that only vegetation decreases in Shenzhen, Foshan, Dongguan, and Zhongshan were higher than vegetation increases, but, signiﬁcant vegetation changes (i.e., decreases and increases) were mainly concentrated in Huizhou, Jiangmen, Zhaoqing, and Guangzhou.


Introduction
Vegetation serves as the interaction substance between the atmosphere, water, and soil. Meanwhile, vegetation dynamics affect the atmospheric cycle [1], carbon cycle [2], and water cycle [3][4][5]. Thus, studying vegetation changes offers insight not only into the atmospheric cycle, water cycle, and carbon cycle but also the interaction between humans and the environment. In-depth knowledge of vegetation climatic raster data. Secondly, the MODIS NDVI data can be downloaded for free and many mature tools and technologies can be used to process MODIS NDVI data, for example, the MODIS Reprojection Tool (MRT). However, the detected results of MOD-TR have many insignificant pixels that increase the uncertainty of vegetation detection. Thus, we combine the significant change pixels of TSS-RESTREND and RESTREND to obtain a more accurate result.
The Pearl River Delta region (PRD) is a densely populated area which attracts a large number of people because of the job opportunities and the excellent environmental conditions. As the population has increased, the environment of the PRD has been significantly affected by humans. As such, insight into the influence of humans on the environment from a spatial perspective has become increasingly important for the Chinese government and researchers to manage the conflict between the population and environment in the PRD. Detecting the vegetation change of the PRD could provide useful information to understand the relationship between humans and the environment. Meanwhile, dense population and frequent economic activities lead to significant vegetation change (i.e., increase or decrease), which could be used to observe the differences for TSS-RESTREND results using different resolutions of NDVI data. Thus, this study focuses on the PRD.
As mentioned earlier, there is a need for using MOD-TR to detect vegetation change in small areas. Still, we were not sure which spatial resolution (i.e., 250 m, 500 m, and 1km) of MODIS NDVI was best for vegetation detection in small area; moreover, a defect we found is that MOD-TR would arise many insignificant change trend pixels which increase the uncertainty of vegetation detection. Hence, to address these two problem, the objectives of this research are: (1) to combine MODIS NDVI of 250 m, 500 m, and 1 km spatial resolutions with TSS-RESTREND, and evaluate which of them is the most suitable for PRD, which would provide a reference point for using MOD-TR to detect vegetation changes in other small areas; and (2) to use a combination of the significant pixel value of TSS-RESTREND and RESTREND to detect reliable city-scale vegetation change in the PRD.

Study Area
The PRD is located in central of Guangdong Province, with a coordinating range of 111 • 21 -115 • 25 E, 21 • 34 -24 • 23 N (Figure 1). There are four main land-cover types in the PRD: cropland, tree cover, urban land, and water bodies. Area percentages of these four land-covers in PRD are as follows: 37.42% of PRD is for Cropland, 46.1% for Tree cover, 11.3% for urban land and 5.18% for Water bodies. There are nine cities: Guangzhou, Shenzhen, Foshan, Dongguan, Zhuhai, Zhongshan, Jiangmen, Huizhou, and Zhaoqing. According to the Guangdong Statistical Yearbook for 2018 (http://stats.gd.gov.cn/gdtjnj/ content/post_1424903.html), the nine cities in the PRD have a total area of 54.770.21 km 2 , accounting for less than one-third of the land area of Guangdong Province. However, they account for 55.53% of the population of the Guangdong Province and 80.23% of the total economic output.

NDVI Data and Its Preprocessing
In this study, NDVI data includes GIMMS NDVI and MODIS NDVI. The reason we chose GIMMS NDVI is that it is the classic data used in TSS-RESTREND to detect vegetation change. The classic data is also a good comparison with MODIS data to show that MODIS NDVI can detect more breakpoints and better reflect the detailed vegetation change in small areas. The GIMMS NDVI data (i.e., GIMMS NDVI3g v1) are available for the period from 1982 to 2015, and MODIS NDVI data are available for the period from 2000 to 2020; therefore, we selected the data for the period 2000-2015.
The GIMMS NDVI data were downloaded from the Ecological Forecasting Lab at NASA Ames Research Center (https://ecocast.arc.nasa.gov/data/pub/gimms). The NDVI and flag data in the nc4 files were extracted using the "rasterizeGimms" function in R-package "gimms" [31] downloaded from Remote Sens. 2020, 12, 4049 4 of 25 The Comprehensive R Archive Network (CRAN; https://CRAN.R-project.org); their projection was performed using the "projectRaster" function in R-package "Raster" [32] from CRAN. Meanwhile, all data were resampled with a resolution of 8 km; the GIMMS NDVI was reconstructed using the weighted SG filter [33] in the R-package "phenofit" [34] from CRAN.

NDVI Data and Its Preprocessing
In this study, NDVI data includes GIMMS NDVI and MODIS NDVI. The reason we chose GIMMS NDVI is that it is the classic data used in TSS-RESTREND to detect vegetation change. The classic data is also a good comparison with MODIS data to show that MODIS NDVI can detect more breakpoints and better reflect the detailed vegetation change in small areas. The GIMMS NDVI data (i.e., GIMMS NDVI3g v1) are available for the period from 1982 to 2015, and MODIS NDVI data are available for the period from 2000 to 2020; therefore, we selected the data for the period 2000-2015.
The GIMMS NDVI data were downloaded from the Ecological Forecasting Lab at NASA Ames Research Center (https://ecocast.arc.nasa.gov/data/pub/gimms). The NDVI and flag data in the nc4 files were extracted using the "rasterizeGimms" function in R-package "gimms" [31] downloaded from The Comprehensive R Archive Network (CRAN; https://CRAN.R-project.org); their projection was performed using the "projectRaster" function in R-package "Raster" [32] from CRAN. Meanwhile, all data were resampled with a resolution of 8 km; the GIMMS NDVI was reconstructed using the weighted SG filter [33] in the R-package "phenofit" [34] from CRAN.
To choose the suitable MODIS NDVI to detect vegetation change in small areas, the MODIS NDVI (i.e., MOD13Q1, MOD13A1, and MOD13A2) data for the period 2000-2015 were downloaded from the Level-1 and Atmosphere Archive & Distribution System (LAADS) Distributed Active Archive Center (DAAC) (https://ladsweb.modaps.eosdis.nasa.gov/); their NDVI data and pixel To choose the suitable MODIS NDVI to detect vegetation change in small areas, the MODIS NDVI (i.e., MOD13Q1, MOD13A1, and MOD13A2) data for the period 2000-2015 were downloaded from the Level-1 and Atmosphere Archive & Distribution System (LAADS) Distributed Active Archive Center (DAAC) (https://ladsweb.modaps.eosdis.nasa.gov/); their NDVI data and pixel quality layers were extracted using the MRT; their re-projection was performed using ArcGIS and all data were resampled with a resolution of 250 m, 500 m, and 1 km. The MODIS NDVI were reconstructed using the weighted SG filter in the R-package "phenofit". The reconstructed GIMMS NDVI and MODIS NDVI data were all converted to an NDVI monthly value by MVC (maximum value composite) [35] and were clipped by the mask of the PRD to obtain the NDVI monthly value in the study area.

Climatic Data and its Preprocessing
Climatic data used in this study included temperature and precipitation. The temperature and precipitation (i.e., 0.5 • × 0.5 • grid data set of surface monthly precipitation and monthly temperature in China) from 2000 to 2015 were downloaded from the China Meteorological Data Sharing Service System (http://data.cma.cn). The data which was from 30 stations ( Figure 2) were interpolated into Remote Sens. 2020, 12, 4049 5 of 25 raster data of 250 m and 500 m using Ordinary kriging interpolation, respectively, and then were clipped by the mask of the PRD.
Climatic data used in this study included temperature and precipitation. The temperature and recipitation (i.e., 0.5° × 0.5° grid data set of surface monthly precipitation and monthly temperature n China) from 2000 to 2015 were downloaded from the China Meteorological Data Sharing Service ystem (http://data.cma.cn). The data which was from 30 stations ( Figure 2) were interpolated into aster data of 250 m and 500 m using Ordinary kriging interpolation, respectively, and then were lipped by the mask of the PRD. Climate Change Initiative (CCI) program website https://climate.esa.int/en/esa-climate/esa-cci/). All three data sources were re-projected, resampled o 250 m, and clipped using the PRD mask.  in China) from 2000 to 2015 were downloaded from the China Meteorological Data Sharing Service System (http://data.cma.cn). The data which was from 30 stations ( Figure 2) were interpolated into raster data of 250 m and 500 m using Ordinary kriging interpolation, respectively, and then were clipped by the mask of the PRD. Climate Change Initiative (CCI) program website (https://climate.esa.int/en/esa-climate/esa-cci/). All three data sources were re-projected, resampled to 250 m, and clipped using the PRD mask.   To facilitate the follow-up study, we combined the land types of ESA CCI LC in 2015 into four categories, namely: cropland, tree cover, urban land, and water bodies. Moreover, to analyze the effect of urban land change, cropland change, and tree cover change on vegetation change in the PRD, the urban land, cropland, and tree cover were extracted from the ESA CCI LC in 2015, and the urban land factor layer, cropland factor layer, and tree cover factor layer were constructed. The urban land factor layer included urban land and non-urban land, the cropland factor layer included cropland and non-cropland, and the tree cover factor layer included tree cover and non-tree cover. The images were downloaded from Google Earth Engine (https://code.earthengine.google.com/), and were reprojected and clipped using the PRD mask. The auxiliary data was mainly used for verification and comparison in this study.

RESTREND
RESTREND is a common method used to detect vegetation change. It calculates the trend of residuals to reflect the vegetation change. The residuals result from the linear regression between the NDVI and climatic factors (i.e., temperature and precipitation), and its Equations are as follows: where y is the NDVI value, y is NDVI value predicted by Equation (1), ε is the residual of Equation (1), a 1 , a 2 , and a 3 are regression coefficients, in which a 3 is the trend (slope), b 1 and b 2 are the intercepts, x 1 and x 2 are climatic factor values, and t is the time range of NDVI. We only select the residual and slope, whose significance level reaches 99%. To facilitate the subsequent writing, the combination of MODIS NDVI and RESTREND is abbreviated as MOD-R, the combination of MOD13Q1 NDVI and RESTREND are abbreviated as Q1-R, and the combination of MOD13A1 NDVI and RESTREND is abbreviated as A1-R.

TSS-RESTREND and its Validation
TSS-RESTREND is an extended algorithm of RESTREND. It combines the breakpoint with linear regression to study the segmented vegetation climate relationship and construct segmented RESTREND (seg-RESTREND). Its Equation are as follows: y s = a 1 x 1 + a 2 x 2 + a 3 z + a 4 x 1 z + a 5 x 2 z + b 3 (4) Equation (4) is the segmented linear regression equation that is used to fit the relationship between vegetation and climatic factors, in which y s is the NDVI value, x 1 and x 2 denote temperature and precipitation, a 1~a5 are regression coefficients, and z is a dummy variable (0 or 1). The value of z before the breakpoint is 0 and after the breakpoint is 1. The breakpoint is calculated by BFAST [20,21] and its significance is calculated using the Chow test [36]. Finally, b 3 is the intercept.
Remote Sens. 2020, 12, 4049 where ε s is the residual of segmented linear regression and y s is the NDVI value predicted by Equation (4). ε s = a 6 t + a 7 z + a 8 tz + b 4 (6) where a 6~a8 are regression coefficients, in which a 7 is the slope required, t is the time of NDVI, z is a dummy variable (0 or 1), and b 4 is intercept. In this study, we only select the residuals, breakpoints, and slopes whose significance level reaches 99%. To facilitate the subsequent writing, the combination of MODIS NDVI and TSS-RESTREND is abbreviated as MOD-TR, the combination of MOD13Q1 NDVI and TSS-RESTREND are abbreviated as Q1-TR, and the combination of MOD13A1 NDVI and TSS-RESTREND is abbreviated as A1-TR. Validation of MOD-TR is very vital for scientific rigor and its application. In our work, a local comparison between MOD-TR, MOD-R, and Landsat was used to verify the effectiveness of MOD-TR.

CTSS-RESTREND and its Validation
According to the Equation of TSS-RESTREND to detect the long-term vegetation change, there are many pixels with an insignificant change. In theory, the pixels with an insignificant change in TSS-RESTREND represent an insignificant abrupt vegetation change caused by breakpoints or denote that the vegetation time series curve of those pixels do not have breakpoints. Thus, they are considered to have a linear change trend that can neglect the abrupt vegetation change. This is the theoretical basis of combining TSS-RESTREND and RESTREND.
According to this theoretical basis, the pixels with a significant change in RESTREND can be applied to overcome the issues caused by the pixels with an insignificant change in TSS-RESTREND. However, it should meet a requirement in that the pixels with an insignificant change in TSS-RESTREND should have a significant change trend value in RESTREND. The combination is easy to implement using ArcGIS; the conditional function (con) and math logical function (IsNull) of the Raster Calculator was used to combine the significant pixels in TSS-RESTREND and RESTREND. The IsNull function was used to select the pixels with an insignificant change in TSS-RESTREND, the con function was used to select the pixels with a significant change in RESTREND to overcome the pixels with an insignificant change in TSS-RESTREND. Equation (7) shows the related Equation in the Raster Calculator of ArcGIS and the operation process is shown in Figure 4. For the convenience of the following description, we refer to the combination of TSS-RESTREND and RESTREND as CTSS-RESTREND. The combination of MODIS NDVI and CTSS-RESTREND is abbreviated as MOD-CTR, the combination of MOD13Q1 NDVI and CTSS-RESTREND is abbreviated as Q1-CTR, and the combination of MOD13A1 NDVI and CTSS-RESTREND is abbreviated as A1-CTR.
Validation of MOD-CTR is very important for scientific rigor. In our study, a local comparison between MOD-CTR, MOD-TR, MOD-R, and Landsat image was used to verify the effectiveness of MOD-CTR. Result = con(IsNull(Raster 1 ), Raster 2 , Raster 1 ), .
where Result indicates the combination result of TSS-RESTREND and RESTREND, Raster 1 indicates the raster of TSS-RESTREND result, Raster 2 indicates the raster of the RESTREND result. Raster 1 and Raster 2 are the significant results of TSS-RESTREND and RESTREND, respectively.
RESTREND to overcome the pixels with an insignificant change in TSS-RESTREND. Equation (7) shows the related Equation in the Raster Calculator of ArcGIS and the operation process is shown in Figure 4. For the convenience of the following description, we refer to the combination of TSS-RESTREND and RESTREND as CTSS-RESTREND. The combination of MODIS NDVI and CTSS-RESTREND is abbreviated as MOD-CTR, the combination of MOD13Q1 NDVI and CTSS-RESTREND is abbreviated as Q1-CTR, and the combination of MOD13A1 NDVI and CTSS-RESTREND is abbreviated as A1-CTR.

Geographical Detector
A geographical detector was proposed by Wang et al. [37] to detect the health risk factors and then was applied to quantify influences of physiographic factors on vegetation [38]. The geographical detector (GD) is composed of four detectors, namely: Risk detector, Factor detector, Ecological detector, and Interaction detector. We only used the Factor detector to quantify the influences of human factors on vegetation change in the PRD. We applied the geographical detector's R-package "GD" [39] to optimize the classification of GDP and population and implemented the Factor detector to analyze the impact of human factors on vegetation change in the PRD. Factor detector is an algorithm used to detect the power of determination (PD) of factor X to dependent variable Y, which is measured by q value. Its Equation is as follows: where N and Nm are the total number of samples and the number of samples in sub-area m, respectively; σ 2 and σ 2 m are the variance of Y over the study area and the variance of Y in sub-area m, respectively; L is the number of sub-area (categories) of factor X. The q value will be between 0 and 1; the closer the q value is to 1, the bigger the influence of this factor is.
where Y j is the jth sample value and Y is the mean of Y over the study area.
where Y m,i and Y m are the value of the ith sample and the mean of Y in sub-area m, respectively.

Net Change
Exclusively using significant vegetation change pixel statistics is insufficient to reflect the vegetation change as it neglects the actual pixel values. Therefore, there is a need to consider both the number and value of pixels to evaluate vegetation change. The net change is an indicator that involves pixel numbers and pixels value. Its Equation is as follows: where Net is the net change, v i is the ith pixel value, n i is the number of the ith pixel value. When the pixels calculated using the Equation (11) show a significant increase in vegetation pixels, the net change is positive; the larger the value of positive net change (PNC), the more the vegetation has grown. When the pixels calculated using the Equation (11) show a significant decrease in vegetation pixels, the net change is negative. The smaller the value of the negative net change (NNC), the more the vegetation has decreased. The sum of PNC and NNC in a study area (i.e., country, province, or city) is called the total net change (TNC); TNC included the negative total net change (NTNC) and positive total net change (PTNC). If the TNC is negative, the study area is dominated by vegetation decrease. If the TNC is positive, the study area is dominated by vegetation increase.

The Process for Choosing Suitable MODIS NDVI
Choosing the best suitable MODIS NDVI is one of the purposes of our study. The selection methods and steps are as follows: Firstly, detailed breakpoints information are very important to TSS-RESTREND to detect vegetation change. Hence, the whole and local breakpoints information of NDVI (i.e., MOD13Q1, MOD13A1, MOD13A2, and GIMMS NDVI) are compared to choose the suitable NDVI and exclude the coarse resolution of NDVI. In the local comparison, some typical significant vegetation changes would be chosen reference to compare detailed vegetation change (i.e., an airport). Secondly, the whole and the local result of the combination of remaining MODIS NDVI with TSS-RESTREND would be compared to show their ability to detail vegetation change. Thirdly, pixel number and net change in every city were calculated to show the difference of MOD-CTR, MOD-TR, and MOD-R at the city-scale, including a comparison of the pixel number of significant decrease (PNSD), pixel number of significant increase (PNSI), and total pixel number (i.e., the value of PNSI minus PNSD) (TPN); comparing of PNC, NNC, NTNC and PTNC; comparison of the order of PNSD, PNSI, TPN, PNC, NNC and TNC in cities. Figure 5, MODIS NDVI (i.e., MOD13Q1, MOD13A1 and MOD13A2) detected more breakpoints than GIMMS NDVI. The proportion of pixels with breakpoints was 62.18% for GIMMS NDVI, 83.60% for MOD13Q1 NDVI, 85.67% for MOD13A1 NDVI and 85.49% for MOD13A2 NDVI, which meant the MOD-TR detected more vegetation changes than GIM-TR. From the Figure 5a-c, we noticed that distribution of break year for MOD13Q1, MOD13A1 and MOD13A2 were highly similar, their area of break year did not vary much across different spatial resolution (i.e., resolution changed from 500 m to 1 km) (Figure 6a), but the pixel numbers of break year varied much across different spatial resolution (Figure 6b), which indicated detail information of vegetation change would be less with the spatial resolution become coarser. When the spatial resolution changed to 8 km, the area of break year for GIMMS NDVI showed an obvious difference with those of MODIS NDVI (Figure 6a,c), and its pixel numbers of break year dropped dramatically (Figure 6d). To show the advantage of MODIS to GIMMS and the difference between MODIS NDVI, we selected Guangzhou Baiyun International Airport as the comparison site. The airport was cropland in 2000 ( Figure 7a) and was transformed to the airport in 2004 (Figure 7b); the airport boundary could be detected by the break year pixels of MODIS NDVI (Figure 7c-e), in which MOD13Q1 detected the most accurate location of the airport (Figure 7c), followed by MOD13A1 (Figure 7d) and MOD13A2 (Figure 7e). Only part of the airport was detected by break year pixels of GIMMS NDVI (Figure 7f). In terms of pixel number, pixel area and accuracy of feature location, GIMMS NDVI, and MOD13A2 NDVI were too coarse to detail the vegetation change in PRD. Hence, GIMMS NDVI and MOD13A2 NDVI were not suitable for detailed vegetation detection in small areas like PRD. Therefore, we only showed the result of vegetation detection using MOD13Q1and MOD13A1 throughout the rest of this article.

A1-TR and Q1-TR Comparison in the PRD
As shown in Figure 8, the spatial distribution of vegetation change detected by Q1-TR and A1-TR were highly similar. The amount of vegetation change detected by the two is also similar. The vegetation change in Q1-TR detected that 24.55% of the pixels showed a significant decrease (Figure 8a), which is approximately 2.54 times more than those detected by Q1-R (9.64%) (Figure 8b). A1-TR detected that 22.46% of the pixels showed a significant decrease (Figure 8c), which is approximately 2.59 times more than those detected by A1-R (8.65%) (Figure 8d). This indicated that MOD-TR detected a more significant decrease in vegetation pixels than MOD-R. It also indicated that TSS-RESTREND detected vegetation decreases that were not detected by RESTREND. Although their spatial distribution and amount of vegetation change are highly similar, compared with A1-TR, the Q1-TR shows more detailed vegetation change. The land of Guangzhou Baiyun International Airport was mostly farmland in 2000 ( Figure 9a); in 2015, the airport had already been built for 11 years (Figure 9b). The Q1-TR detected a significant decrease in pixels caused by the airport construction, and boundaries of the airport were all detected (Figure 9c). However, A1-TR did not detect the change caused by the airport construction in detail; this method just identifies a part of the significant change for the airport area ( Figure 9d). Moreover, we found that TSS-RESTREND (i.e., Q1-TR [30.06%] and A1-TR [28.04%]) classified more pixels as 'no significant change' than RESTREND (i.e., Q1-R [20.81%] and A1-R [18.85%]), which increases the uncertainty detecting a change.  airport construction, and boundaries of the airport were all detected (Figure 9c). However, A1-TR did not detect the change caused by the airport construction in detail; this method just identifies a part of the significant change for the airport area ( Figure 9d). Moreover, we found that TSS-RESTREND (i.e., Q1-TR [30.06%] and A1-TR [28.04%]) classified more pixels as 'no significant change' than RESTREND (i.e., Q1-R [20.81%] and A1-R [18.85%]), which increases the uncertainty detecting a change.

A1-CTR and Q1-CTR Comparison in the PRD
As mentioned in Section 3.2, some no significant change pixels would lead to a greater uncertainty of vegetation detection. To overcome this drawback and obtain a more accurate indication of vegetation change, we combined the significant change pixels of TSS-RESTREND and RESTREND (i.e., CTSS-RESTREND). As shown in Figure 10, the spatial distribution of vegetation change detected by A1-CTR and Q1-CTR were quite similar, and the amount of no significant change pixels decreased significantly after they were combined. As shown in Figures 8 and 10, for the Q1-CTR, the proportion of no significant change pixels decreased from 30.06% to 7.43%; for the A1-CTR, the proportion of no significant change pixels decreased from 28.04% to 6.34%.

A1-CTR and Q1-CTR Comparison in the PRD
As mentioned in Section 3.2, some no significant change pixels would lead to a greater uncertainty of vegetation detection. To overcome this drawback and obtain a more accurate indication of vegetation change, we combined the significant change pixels of TSS-RESTREND and RESTREND (i.e., CTSS-RESTREND). As shown in Figure 10, the spatial distribution of vegetation change detected by A1-CTR and Q1-CTR were quite similar, and the amount of no significant change pixels decreased significantly after they were combined. As shown in Figure 8 and Figure 10, for the Q1-CTR, the proportion of no significant change pixels decreased from 30.06% to 7.43%; for the A1-CTR, the proportion of no significant change pixels decreased from 28.04% to 6.34%. With the decrease in no significant change pixels, both "significant decrease pixels" and "significant increase pixels" increased. For the Q1-CTR, the proportion of "significant decrease pixels" increased from 24.55% to 27.40%, and the proportion of "significant increase pixels" increased from 45.39% to 65.16%. For the A1-CTR, the proportion of "significant decrease pixels" increased from 22.46% to 25.18%, and the proportion of "significant increase pixels" increased from 49.5% to 68.48%. Thus, the amount of vegetation change detected by Q1-CTR and A1-CTR were actually very similar. With the decrease in no significant change pixels, both "significant decrease pixels" and "significant increase pixels" increased. For the Q1-CTR, the proportion of "significant decrease pixels" increased from 24.55% to 27.40%, and the proportion of "significant increase pixels" increased from 45.39% to 65.16%. For the A1-CTR, the proportion of "significant decrease pixels" increased from 22.46% to 25.18%, and the proportion of "significant increase pixels" increased from 49.5% to 68.48%. Thus, the amount of vegetation change detected by Q1-CTR and A1-CTR were actually very similar.
The results of the city-scale vegetation change are shown in Figure 11. The difference between TSS-RESTREND and RESTREND was that city-scale vegetation change of RESTREND was dominated by vegetation increase (Figure 11e, f). Namely, the city-scale vegetation significant decrease in RESTREND was far less than the decrease in TSS-RESTREND. For Q1-R and A1-R, the PNSD in all nine cities was < 20,000 and < 3200 (Figure 11e, f). However, for Q1-TR and A1-TR, there were five cities where the PNSD was > 20,000 ( Figure 9c) and six cities where the PNSD was > 3200 (Figure 9d). Actually, for Q1-TR and A1-TR, their PNSD in the remaining cities were all more than those detected by Q1-R and A1-R (Figure 11c, d, e, f), which indicates the RESTREND underestimated the vegetation decrease. However, the distribution of significant vegetation increase for TSS-RESTREND and RESTREND was similar. Their significant vegetation increase was mainly concentrated in Zhaoqing, Huizhou, Jiangmen, Guangzhou, and Foshan city.
The results of the city-scale vegetation change are shown in Figure 11. The difference between TSS-RESTREND and RESTREND was that city-scale vegetation change of RESTREND was dominated by vegetation increase (Figure 11e,f). Namely, the city-scale vegetation significant decrease in RESTREND was far less than the decrease in TSS-RESTREND. For Q1-R and A1-R, the PNSD in all nine cities was < 20,000 and < 3200 (Figure 11e,f). However, for Q1-TR and A1-TR, there were five cities where the PNSD was > 20,000 ( Figure 9c) and six cities where the PNSD was > 3200 (Figure 9d). Actually, for Q1-TR and A1-TR, their PNSD in the remaining cities were all more than those detected by Q1-R and A1-R (Figure 11c-f), which indicates the RESTREND underestimated the vegetation decrease. However, the distribution of significant vegetation increase for TSS-RESTREND and RESTREND was similar. Their significant vegetation increase was mainly concentrated in Zhaoqing, Huizhou, Jiangmen, Guangzhou, and Foshan city.
The difference between using CTSS-RESTREND and TSS-RESTREND was that the PNSI in the nine cities in CTSS-RESTREND were more than those in TSS-RESTREND. But their results were remarkably similar. Q1-TR and A1-TR identified three cities in which a decrease in vegetation was greater than the increase in vegetation, namely Shenzhen, Dongguan, and Zhongshan (Figure 11c,d). MOD-CTR (i.e., Q1-CTR and A1-CTR) also identified two cities (Dongguan and Zhongshan) as having more vegetation decrease than vegetation increase (Figure 11a,b). Besides, although the spatial resolution of MOD13Q1 and MOD13A1 were different, the order of TPN, PNSI, and PNSD detected by Q1-CTR and A1-CTR; Q1-TR and A1-TR; Q1-R and A1-R were almost the same. It indicated that MOD13Q1 and MOD13A1 two NDVI were suitable for vegetation detection in small areas like PRD. The difference between using CTSS-RESTREND and TSS-RESTREND was that the PNSI in the nine cities in CTSS-RESTREND were more than those in TSS-RESTREND. But their results were remarkably similar. Q1-TR and A1-TR identified three cities in which a decrease in vegetation was greater than the increase in vegetation, namely Shenzhen, Dongguan, and Zhongshan (Figure 11c, d). MOD-CTR (i.e., Q1-CTR and A1-CTR) also identified two cities (Dongguan and Zhongshan) as having more vegetation decrease than vegetation increase (Figure 11a, b). Besides, although the spatial resolution of MOD13Q1 and MOD13A1 were different, the order of TPN, PNSI, and PNSD detected by Q1-CTR and A1-CTR; Q1-TR and A1-TR; Q1-R and A1-R were almost the same. It indicated that MOD13Q1 and MOD13A1 two NDVI were suitable for vegetation detection in small areas like PRD.

Comparison of NDVI Net Change between CTSS-RESTREND, TSS-RESTREND, and RESTREND
There were obvious differences in the city-scale net change between the use of RESTREND, TSS-RESTREND (i.e., Q1-TR and A1-TR), and CTSS-RESTREND (i.e., Q1-CTR and A1-CTR) ( Figure 12). Compared with RESTREND (i.e., Q1-R and A1-R), four cities (i.e., Shenzhen, Foshan, Dongguan, and Zhongshan) had a NTNC when using Q1-TR and Q1-CTR, three cities (i.e., Shenzhen, Dongguan, and Zhongshan) had a NTNC when using A1-TR and A1-CTR. Compared with the other two methods, the TNC detected using RESTREND were all positive (Figure 12e,f). Using the Q1-R and A1-R, they two all detected a PNC in all nine cities (Figure 12e,f). Moreover, the PNC was mainly seen in Zhaoqing, Huizhou, Jiangmen, and Guangzhou. However, RESTREND also detected a small amount of NNC (i.e., vegetation decrease). When using Q1-R, a NNC which was detected in nine cities ranged from −7.99 to −2.19 (Figure 12e). When using A1-R, the NNC, which was detected in nine cities ranged from −1.53 to −0.39 (Figure 12f). two methods, the TNC detected using RESTREND were all positive (Figure 12e, f). Using the Q1-R and A1-R, they two all detected a PNC in all nine cities (Figure 12e, f). Moreover, the PNC was mainly seen in Zhaoqing, Huizhou, Jiangmen, and Guangzhou. However, RESTREND also detected a small amount of NNC (i.e., vegetation decrease). When using Q1-R, a NNC which was detected in nine cities ranged from −7.99 to −2.19 (Figure 12e). When using A1-R, the NNC, which was detected in nine cities ranged from −1.53 to −0.39 (Figure 12f). For TSS-RESTREND and CTSS-RESTREND, their detected city-scale NNC were all less than when using RESTREND. This means that TSS-RESTREND and CTSS-RESTREND detected a higher city-scale vegetation decrease than when using RESTREND. In terms of Q1-CTR and A1-CTR, their detected city-scale NNC varied from −39.96 to −7.48, and from −7.58 to −1.45, respectively ( Figure  12a, b). In terms of Q1-TR and A1-TR, their detected NNC varied from −38.14 to −6.63, and from −7.24 to −1.25, respectively (Figure 12c, d). The NNC detected by Q1-CTR, A1-CTR, Q1-TR and A1-TR at the city-scale were all smaller than those detected by Q1-R and A1-R in every city ( Figure  12).
Besides, as shown in Figure 12, we noticed that Q1-CTR, Q1-TR, A1-CTR, and A1-TR had the same city's order of detected TNC. Q1-CTR and A1-CTR had the same order of detected PNC for For TSS-RESTREND and CTSS-RESTREND, their detected city-scale NNC were all less than when using RESTREND. This means that TSS-RESTREND and CTSS-RESTREND detected a higher city-scale vegetation decrease than when using RESTREND. In terms of Q1-CTR and A1-CTR, their detected city-scale NNC varied from −39.96 to −7.48, and from −7.58 to −1.45, respectively (Figure 12a,b). In terms of Q1-TR and A1-TR, their detected NNC varied from −38.14 to −6.63, and from −7.24 to −1.25, respectively (Figure 12c,d). The NNC detected by Q1-CTR, A1-CTR, Q1-TR and A1-TR at the city-scale were all smaller than those detected by Q1-R and A1-R in every city ( Figure 12).
Besides, as shown in Figure 12, we noticed that Q1-CTR, Q1-TR, A1-CTR, and A1-TR had the same city's order of detected TNC. Q1-CTR and A1-CTR had the same order of detected PNC for nine cities, Q1-TR and A1-TR had the same order of detected PNC for nine cities, their difference focus on the order of detected NNC. Overall, in terms of net change, characteristics of vegetation change detailed by Q1-CTR and A1-CTR, Q1-TR and A1-TR were almost same.

Human-Driven Factor Analysis of Vegetation Change in the PRD
We applied the GD method to analyze the influence of human-driven factors (i.e., population, GDP, and land-cover change) on vegetation change in the PRD. The result of MOD-CTR was used to reflect the vegetation change. The grid spatial data of population, GDP and Land Cover raster were chosen as the human-driven factors.
The q value of factor detection can measure the influence of human factors. The q value of human factors on vegetation in the PRD were ranked as follows: tree cover change (0.001457) < cropland change (0.006932) < urban land change (0.034328) < land change (0.035532) < GDP (0.041599) < population (0.041646) ( Table 1). In terms of the rank of the q value, Population and GDP were the dominant factors affecting the vegetation change, followed by land change, and then urban land change, a form of land change. This means urban land change is an important driver causing the vegetation to decrease. However, most land changes were caused by economic development and population growth [40]; thus, the influence of land-cover change was determined to be an indirect impact of economic development and population growth to vegetation change.

Validation of the MOD-TR in PRD
Compared with similar studies [19,23,25], here, we used MOD-R to study vegetation changes in small areas. The results showed that MOD-TR had detected more pixels with a significant decrease than when using MOD-R in small areas. But the validity of the MOD-TR needed to be tested. To confirm the result of MOD-TR, part of the result of MOD-TR was chosen to compare with MOD-R and Landsat. As shown in Figure 13, a part of Dongguan was chosen to prove the validity of the MOD-TR. In a location which were circled by blue line, Q1-R could not detect pixels with a significant decrease (Figure 13a), but some pixels with a significant decrease could be detected by Q1-TR ( Figure 13b). As shown in Figure 13c

Validation of the CTSS-RESTREND in PRD
CTSS-RESTREND is a combination of the significant pixels in TSS-RESTREND and RESTREND. Although CTSS-RESTREND could detect more pixels with a significant change (i.e., decrease and increase) trend than those of TSS-RESTREND and RESTREND, it needs to demonstrate its validity. Hence, we chose a location in Zhaoqing to show its validity ( Figure 14). As showed in Figure 14, Q1-CTR detected more pixels with a significant decrease and increase trend than when using Q1-TR and Q1-R.

Validation of the CTSS-RESTREND in PRD
CTSS-RESTREND is a combination of the significant pixels in TSS-RESTREND and RESTREND. Although CTSS-RESTREND could detect more pixels with a significant change (i.e., decrease and increase) trend than those of TSS-RESTREND and RESTREND, it needs to demonstrate its validity. Hence, we chose a location in Zhaoqing to show its validity ( Figure 14). As showed in Figure 14, Q1-CTR detected more pixels with a significant decrease and increase trend than when using Q1-TR and Q1-R. The pixels which were circled by the blue line in Q1-CTR demonstrated a significant decrease trend ( Figure 14a). Meanwhile, they were some pixels with no significant change in Q1-TR ( Figure  13b), but they were pixels with a significant decrease in Q1-R (Figure 14c). The Land-cover circled by the blue line most was vegetation in 2000 ( Figure 14d); this vegetation changed to wasteland or decreased in 2015 (Figure 14e), which proved the rationality of Q1-CTR. Hence, it was reasonable to fill these pixels with no significant change in Q1-TR with the pixels with a significant decrease in Q1-R.
In the same way, the vegetation, which was circled by the purple line in Q1-CTR demonstrated significant increase (Figure 14a). Meanwhile, they were some pixels with no significant change in Q1-TR (Figure 14b), but they were pixels with a significant increase in Q1-R (Figure 14c). The Land-cover circled by purple line most were vegetation in 2000 ( Figure 14d); this vegetation increased in 2015 (Figure 14e), which indicated Q1-CTR was effective. Therefore, it was rational to fill these pixels with no significant change in Q1-TR with the pixels with a significant increase in Q1-R. Above all, the comparison in Figure 14 showed the CTSS-RESTREND was rational; meanwhile, it indicated that CTSS-RESTREND could detect more accurate vegetation significant change than when using TSS-RESTREND and RESTREND.

Limitations of MOD-TR Applied to Small Areas
However, MOD-TR also has some limitations when applied to small areas. Firstly, there are significantly more pixels in the study area as the image resolution increased from 8 km (i.e., GIMMS NDVI) to 250 m (i.e., MODIS NDVI). This increases the detection time of the breakpoints of The pixels which were circled by the blue line in Q1-CTR demonstrated a significant decrease trend ( Figure 14a). Meanwhile, they were some pixels with no significant change in Q1-TR (Figure 13b), but they were pixels with a significant decrease in Q1-R (Figure 14c). The Land-cover circled by the blue line most was vegetation in 2000 ( Figure 14d); this vegetation changed to wasteland or decreased in 2015 (Figure 14e), which proved the rationality of Q1-CTR. Hence, it was reasonable to fill these pixels with no significant change in Q1-TR with the pixels with a significant decrease in Q1-R.
In the same way, the vegetation, which was circled by the purple line in Q1-CTR demonstrated significant increase (Figure 14a). Meanwhile, they were some pixels with no significant change in Q1-TR (Figure 14b), but they were pixels with a significant increase in Q1-R (Figure 14c). The Land-cover circled by purple line most were vegetation in 2000 ( Figure 14d); this vegetation increased in 2015 (Figure 14e), which indicated Q1-CTR was effective. Therefore, it was rational to fill these pixels with no significant change in Q1-TR with the pixels with a significant increase in Q1-R. Above all, the comparison in Figure 14 showed the CTSS-RESTREND was rational; meanwhile, it indicated that CTSS-RESTREND could detect more accurate vegetation significant change than when using TSS-RESTREND and RESTREND.

Limitations of MOD-TR Applied to Small Areas
However, MOD-TR also has some limitations when applied to small areas. Firstly, there are significantly more pixels in the study area as the image resolution increased from 8 km (i.e., GIMMS NDVI) to 250 m (i.e., MODIS NDVI). This increases the detection time of the breakpoints of vegetation, which is a drawback for a regional study. Thus, faster algorithms to detect the breakpoints should be applied or developed in future studies. Secondly, the time range (i.e., 2000-2020) of MODIS NDVI data is shorter than in GIMMS NDVI , which is not conducive to the longer-term study of vegetation change. An intra-ensemble of MODIS NDVI and GIMMS NDVI should be conducted in future studies.

Advantages and Limitations of CTSS-RESTREND
CTSS-RESTREND is a combination of TSS-RESTREND and RESTREND. Its advantage was as following: CTSS-RESTREND has fewer "no significant pixels" than TSS-RESTREND and RESTREND. Meanwhile, CTSS-RESTREND was able to detect some significant vegetation changes (i.e., decrease and increase) that were not detected by TSS-RESTREND and RESTREND. As shown in Figures 8 and 10, the significant increase of pixels was greater than the increase in a significant decrease in pixels. This shows that TSS-RESTREND may underestimate the number of a significant increase in pixels. Interestingly, the proportion of significant decrease in pixels also increased, indicating that CTSS-RESTREND and RESTREND also detected some vegetation decrease pixels, which were undetected by TSS-RESTREND.
However, CTSS-RESTREND also has its limitations. Namely, there were also some no significant pixels exits in MOD-CTR. This is because those pixels were no significant pixels in TSS-RESTREND and RESTREND. If there were to be a high number of these pixels, which are insignificant in TSS-RESTREND and RESTREND, the CTSS-RESTREND results might not be reliable.

Human Factors Have a Double Effect on Vegetation Change in the PRD
We observed that there were both vegetation decreases and vegetation increases in the urban land in the PRD (Figure 15). For example, 9246 pixels showed a significant increase, and 9108 pixels showed a significant decrease in Guangzhou's urban land (Table 2). Guangzhou's PNC was 6.76 and its NNC was −10.91 (Table 2). Moreover, as shown in Table 2, the remaining eight cities also had a significant vegetation increase and significant vegetation decrease. This means that urban land change leads to both vegetation decreases and vegetation increases. The vegetation increase in urban areas can lead to a better environment for humans, and economic development and human activities can promote vegetation increase, as found in a previous study [41]. showed a significant decrease in Guangzhou's urban land (Table 2). Guangzhou's PNC was 6.76 and its NNC was −10.91 (Table 2). Moreover, as shown in Table 2, the remaining eight cities also had a significant vegetation increase and significant vegetation decrease. This means that urban land change leads to both vegetation decreases and vegetation increases. The vegetation increase in urban areas can lead to a better environment for humans, and economic development and human activities can promote vegetation increase, as found in a previous study [41].  In general, the dominant forces affecting vegetation change in the PRD are population and GDP. These caused the land-use change (the indirect impact of population and GDP), land-use changes enabled vegetation change, which includes vegetation decrease and vegetation increase, and thus, population and GDP have both a positive effect and negative effect. Vegetation increase in urban land may be caused by the construction of urban greenbelts; vegetation decrease in urban land was mainly caused by urban expansion and traffic road construction, but it is impossible for cities and traffic roads to develop without restriction, so vegetation growth may exceed vegetation reduction in urban areas.

Conclusions
In this study, we rewrote the TSS-RESTREND code to be adapted to monthly MODIS NDVI data with a resolution of 250 m, 500 m, and 1 km. The new TSS-RESTREND code could use MODIS NDVI to detect more accurate vegetation change data. However, there were many pixels with no significant change. Thus, CTSS-RESTREND was proposed to overcome this limitation and used the GD method to analyze the human-driven factors affecting vegetation change. We draw the following main conclusions:

5.
Although MOD-TR could detect more detail, it takes more time to detect vegetation change; meanwhile, MODIS NDVI is available for a short period. Thus, faster methods for breakpoint detection should be developed for MOD-TR to detect detailed vegetation change; an intraensemble of MODIS NDVI and GIMMS NDVI should be conducted in future studies to do a longer-term study.