Trends in vegetation productivity related to climate change in China’s Pearl River Delta

April and September, and cool dry winters from November to March. The natural vegetation consists of tropical and subtropical evergreen forest. Methods Data sets and pre-processing NDVI data and phenology extraction. Coarse-scale GIMMS (Global Inventory Model- ing and Mapping Studies) NDVI datasets provide the most consistent long-term data of the earth’s biosphere [22]. Its recent version (NDVI3g.v1; https://ecocast.arc.nasa.gov/data/pub/ gimms/3g.v1/) [23,24] spans more than 34 years from 1981–2015 at 8 km spatial resolution (Table 1). The datasets are based on 15-day maximum value composites that provide 24 NDVI composite images a year. For this study, we used 816 fortnightly observations for 34 years (1982 to 2015) of the NDVI data. The fortnightly Maximum Value Composites (MVC) largely removed noise, though it includes disturbances caused by cloud contamination, atmospheric variability, and bidirectional effects, which show as undesirable noise [5], therefore, the NDVI product is accompanied by a quality assessment (QA) layer which enables pixels to be excluded or weighted when reconstructing the time series. To avoid this, the Whittaker smoothing approach was implemented in R to reconstruct a continuous NDVI time-series. In evergreen ecosystems, changes in vegetation phenology and productivity are driven by the emergence of new leaves and shedding of older leaves, whereas seasonality in photosynthe- sis is explained by leaf development and demography [27]. We applied the midpoint pixel

approach to extract the Start of Season (SOS) and End of Season (EOS) time from the smoothed NDVI time series by scaling the annual cycle between 0 and 1 [28]. This is one of the most consistent methods for extraction of phenology metrics regardless of vegetation cover type, and applicable to a variety of ecosystems [7,29]. Then six additional phenology and vegetation productivity metrics (length of growing season (LOS), position of season peak (POP), position of season trough (POT), mean growing season NDVI (MGS), maximum sea- sonal NDVI (Peak), and minimum seasonal NDVI (Trough)) were determined using the Phe- nologyRaster() function in the R package ‘green-brown’. Further analysis was confined to pixels showing significant seasonality where the mean annual NDVI amplitude was greater than 0.1 [29]. Description and ecological meanings of the phenometrics are given in the sup- plementary S1 Table. Seasonal climate variables from satellite-based precipitation and temperature data. Temperature and precipitation are considered representative climate variables. Station-based point measurements cannot capture a continuous spatial distribution of precipitation and temperature [30]. Therefore, to compare the spatial and temporal variability of precipitation with the NDVI data, the CHIRPS (Climate Hazards Group Infrared Precipitation with Station) data set, developed by the United States Geological Survey (USGS) in collaboration with Earth Resource Observation and Science (EROS) centre is used [25]. CHIRPS is generated by inte-

grating satellite imagery and in-situ gauge collected observations. The daily rainfall data are distributed at 0.5˚ (~5 km) spatial resolution. In this study, we used the daily CHIRPS data Table 1. Datasets used for the study. No. Data/Product Data Type Time Range Spatial Resolution Temporal Resolution 1 NDVI GIMMS3g [23,24] Raster 1982–2015 ~8 km 15 Days 2 CHIRPS Precipitation [25] Raster 1982–2015 ~5 km Daily 3 ERA5 (mean_2m_air_temperature) [26] Raster 1982–2015 ~25 km Daily 4 MODIS LCLU product (MCD12Q1) [21] Raster 2001–2015 ~0.5 km Annual https://doi.org/10.1371/journal.pone.0245467.t001 PLOS ONE Climate and vegetation trends in the Pearl River Delta PLOS ONE | https://doi.org/10.1371/journal.pone.0245467 February 24, 2021 4 / 20