A new global land productivity dynamic product based on the consistency of various vegetation biophysical indicators

ABSTRACT Changes in land productivity have been endorsed by the Inter Agency Expert Group on Sustainable Development Goals (IAEG-SDGs) as key indicators for monitoring SDG 15.3.1. Multiple vegetation parameters from optical remote sensing techniques have been widely utilized across different land productivity decline processes and scales. However, there is no consensus on indicator selection and their effectiveness at representing land productivity declining at different scales. This study proposes a fusion framework that incorporates the trends and consistencies within the four commonly used remote sensing-based vegetation indicators. We analyzed the differences among the four vegetation parameters in different land cover and climate zones, finally producing a new global land productivity dynamics (LPD) product with confidence level degrees. The LPD classes indicated by the four vegetation indicators(VIs) showed that all three levels (low, medium, and high confidence) of increasing area account for 23.99% of the global vegetated area and declining area account for 7.00%. The Increase high-confidence(HC) area accounted for 2.77% of the total area, and the Decline-HC accounted for 0.35% of the total area. This study demonstrates the accuracy of the high-confidence (HC) area for the evaluation of land productivity decline and increase. The “forest” landcover type and “humid” climate zone had the largest increasing and declining area but had the lowest high-confidence proportion. The data product provides an important and optional reference for the assessment of SDG 15.3.1 at global and regional scales according to the specific application target. The “Global Land Productivity Dynamic dataset” is available in the Science Data Bank at http://www.doi.org/10.11922/sciencedb.j00076.00084.


Introduction
Land degradation is one of the most serious ecological and environmental problems affecting most countries on Earth. Sustainable development goals (SDGs) were implemented in January 2016 (UNCCD, 2016). SDG target 15.3 is that we need to "by 2030, combat desertification, restore degraded land and soil, including land affected by desertification, drought, and floods, and strive to achieve a land degradation-neutral world." An increased number of countries have begun to monitor indicator 15.3.1, i.e. the "proportion of land that is degraded over the total land area," to achieve SDG target 15.3 in 2030(UNCCD, 2016. Thus,evaluating SDG 15.3.1 at the global scale has become an important recent focus. Land productivity reflects the efficiency of solar energy conversion and land fixation, driven by plant photosynthesis that forms terrestrial vegetation cover. It can be quantified in terms of the net primary production (NPP). All other organisms (e.g. humans, other species of animals, bacteria, and fungi) depend directly and indirectly on this primary production for their health and well-being. Globally, humans appropriate a constantly increasing proportion of this NPP, affecting the structure and functioning of ecosystems; in many cases, this proportion exceeds the natural variability and dynamics of the ecosystems. Therefore, land productivity is an essential variable for detecting and monitoring active land transformations typically associated with land degradation processes. Trends in land productivity have thus been adopted by the United Nations Convention to Combat Desertification (UNCCD) as one of three biophysical progress indicators (UNCCD, 2013). Trends in land productivity have been proposed as a sub-indicator for the global indicator to monitor progress toward achieving SDG target 15.3 for land degradation neutrality (LDN) (UNCCD, 2016).
Global land productivity monitoring currently relies on multi-temporal evaluations of long-term remotely sensed vegetation index time series computed from continuous spectral measurements of photosynthetic activity. The normalized difference vegetation index (NDVI) is the most commonly used vegetation index; other indices have been proposed and used for different scales, such as the enhanced vegetation index (EVI), leaf area index (LAI), and net primary productivity (NPP) (Chen, Park, et al., 2019;de Jong, Verbesselt, Zeileis, & Schaepman, 2013;Fensholt et al., 2012;Pan et al., 2018;Watts & Laffan, 2014;Wu et al., 2015;Zhao & Running, 2010). These indices reflect different vegetation characteristics and have been reported to perform better than other indices under specific vegetation conditions (Baret & Guyot, 1991); however, there is no consensus on the selection of appropriate indicators.
Retrieving meaningful metrics based on time series vegetation indices is another challenge for land degradation assessments owing to observations of highly variable land productivity between different years or vegetation growth cycles, even under stable conditions. The Joint Research Centre (JRC) developed a land productivity dynamics (LPD) dataset for land degradation and improvement assessments, which was adopted by the third edition of the World Atlas of Desertification (WAD3) (Cherlet et al., 2018) and the UNCCD LDN Target Setting Program (Gonzalez-Roglich et al., 2019). This LDP dataset defines five classes by combining three metrics, i.e. trends, states, and performance, based on the SPOT NDVI from 1999 to 2013. Here, trends reflect the significance of the magnitude and persistence of changes in plant productivity over time, states compare the current productivity levels with the historical productivity levels for the same region, and performance compares the local productivity levels with the productivity levels for similar land units at a national scale. Although the state and performance metrics provide new important insights into land degradation, their effective determination faces many technical challenges. Therefore, the LPD class mainly depends on the objective trend metric in a real assessment (Sims et al., 2019).
Studies have used different metrics to explore the time series potential in LPD calculations. However, relatively few studies have focused on metrics based on the difference between various remote sensing-based land productivity indicators. This study aimed to produce a new global LPD dataset by combining the trend metric of multiple stable time series vegetation indicators (i.e. the NPP, LAI, EVI, and NDVI), with a special focus on confidence classification based on consistencies among different indicators. Additionally, we used ESA-CCI landcover in 2015 and arid index data to analyze the LPD differences between indicators and in different land cover and climate regions to facilitate the selection of LPD results, which finally contributed to an improved assessment of the SDG 15.3.1 target.

NDVI and EVI data
The MODIS MOD13A1 Collection6 products provide two vegetation indices (VIs): NDVI and EVI (GEE product id: MODIS/006/MOD13A1). The algorithm for this product selects the best available pixel value from all the acquisitions over a 16-day period. The NDVI and EVI were calculated from the red, near-infrared, and blue bands after atmospheric correction. The NDVI is significantly affected by soil brightness in areas with low vegetation coverage. The EVI can minimize changes in the canopy soil and sensitivity under dense vegetation conditions. The MODIS Collection 6 products are now widely used in ecosystem, climate, and resource management research.

NPP data
The MODIS Terra NPP product (GEE product id: MOD17A3HGF.006) was provided by the NASA LP DAAC at the USGS EROS Center. The NPP was calculated from the sum of the eight-day net photosynthetic products using the gross primary productivity (GPP) and maintenance respiration (MR).

LAI data
The MODIS MOD15A2H dataset provides the LAI index (GEE product id: MODIS/006/ MOD15A2H), which is an eight-day composite dataset with a pixel size of 500 m. The LAI is defined as the one-sided green leaf area per unit ground area across broadleaf canopies, as well as one-half of the total needle surface area per unit ground area across coniferous canopies. It is widely used to calculate photosynthesis, surface evapotranspiration, NPP, the C and water cycle processes, and the vegetation biogeochemical cycle.

Landcover data
The European Space Agency (ESA) Climate Change Initiative Land Cover (CCI-LC) maps were released as the first annual global land cover time series maps at 300 m, spanning 28 years from 1992 to 2019, which is recommended by the Good Practice Guidance (GPG) SDG Indicator 15.3.1 (Sims et al., 2019). These maps were constructed by merging multiple Earth observation products based on the GlobCover products produced by the ESA, with an overall accuracy of approximately 71.1% (Liu et al., 2018). The Land Cover Classification System (LCCS) contains 22 categories at a global scale and detailed classifications at a regional scale. We ignored the dynamics of the vegetation productivity in the "settlement" and "water" classes. In general, these two types of landcover are considered stable areas. According to Trend.Earth (https://trends.earth) SDG 15.3.1 documentation and GPG (Sims et al., 2019), "Sparse vegetation" and "Shrubland" were reclassified into grassland categories; therefore, we used "bare land" instead of the "other" category. The 2015 landcover was used and resampled to 500 m spatial resolution by using the nearest neighbor method (Bilal et al., 2018). Figure 1 shows the distribution of the reclassified LCCS-LC data.

Aridity geospatial data
The Consultative Group for International Agriculture Research (CGIAR) Global-Aridity Geospatial Dataset was obtained from Trabucco & Zomer (2009). The Global-Aridity dataset uses the standard ESRI grid geospatial format, with a spatial resolution of 1 km. The Global-Aridity dataset was used and resampled to 500 m spatial resolution by using the nearest neighbor method (Bilal et al., 2018). These data are widely used to support sustainable development, biodiversity, and environmental protection. Figure 2 shows a global aridity map.

Data preprocessing
The LAI and NPP indicators were produced at eight-day intervals, whereas the EVI and NDVI indicators were produced at 16-day intervals. To integrate different indicators, the time resolution for the four indicators was unified via preprocessing. For the three indicators NDVI, EVI, and LAI, we calculated the annual average value of each indicator to be 1.1 to 12.31 from 2000 to 2015. The annual average values of the four indicators were calculated using the Google Earth Engine (Gorelick et al., 2017) to analyze interannual trends, where each indicator had 16 images for their annual average values from 2000 to 2015 at the global scale. For the NPP indicator, we used the annual average product. Identical resolution datasets were used for the spatial resolution; the spatial resolution of the four indicators was unified to 500 m. The spatial distributions of the four indicators were inconsistent. For example, some pixels have NDVI indicator values, but the NPP indicator is masked. The 16-year averages of the four indicators were used to mask each index, and the intersection range was used as the study area. Table 1 shows the time period and processing method of the four indicators.

Trend analysis
Trend analysis was performed using the Mann-Kendall trend test (Hamed & Rao, 1998;Mann, 1945), which is a non-parametric test used to detect a monotonic trend in time series data. Compared with parametric tests, non-parametric tests do not require samples that follow a specific distribution and are not disturbed by outliers. This method is more suitable for the detection of changes in all types of distribution (Hamed, 2009).
The Mann-Kendall correlation for the two sets of observations was formulated to calculate S as follows: where Figure 2. Global aridity climate classification. AI denotes aridity index value that increases for more humid conditions and decreases with more arid conditions.
(2) and the S statistics represent trends in the time series. A positive trend occurs if S > 0. A negative trend occurs if S < 0. Mann (1945) provided proof of the asymptotic normality of the S statistic. The standardized test statistic Z was calculated as follows: The significance of the trends was tested by comparing the standardized test statistic Z with the standard normal variate at the desired significance level. Trends characterized by P < 0.05 represented a significant level, indicative of increasing or declining areas. Pixels without significant trends were considered stable.

Land productivity confidence classification
The Mann-Kendall trend detection method ignores the autocorrelation and cross correlation for every pixel, which could lead to an overestimation of trend estimation and regional land degradation statistics. Therefore, this study used the number trends as the upper limiting conditions to avoid the overestimation in singal indicators. In this study, we assume that the single indicator can reflect land productivity decline but lacks accuracy. We introduce a method to accurately determine land productivity trends. Identical pixels may exhibit different trends based on the various indicators. For each pixel, more indicators showing the same trend yield increased confidence. According to the count of indicators with the same trends, the confidence for increasing and declining areas can be divided into three levels, which are the combination of the tendency and the confidence level. For the trend in single pixels, the pixel that contained two trends of increase and stability was considered as the Increase class. The pixel that only contained stable trends was considered as the Stable class. The pixel that only contained trends of increase and stability trends was considered as the Increase class. The pixel that contains both trends of increase and decline was considered as the Conflict class. For the confidence level in single pixels, the pixel that contained four increase or decline indicators was considered at the high confidence (HC) level. The pixel that contained two or three increase or decline indicators and one or two stable indicators was considered at the medium confidence (MC) level. The pixel that contained one increase or decline and three stable indicators was considered at the low confidence (LC) level. Pixels that have both positive and negative trends with different indicators are considered conflict areas. Table 2 lists the classifications for productivity confidence. In Table 2, i is the number of increased pixels for each indicator, j is the number of stable pixels for each indicator, and k is the number of declined pixels for each indicator. Table 2 also lists the productivity confidence classification conditions. Figure 3 shows the land productivity declining and increasing based on the NPP, EVI, NDVI, and LAI indices. Significant differences were observed between the increasing and declining areas on a global scale. For the statistics of each indicator, the results for the NPP indicator had the highest proportion of increasing area (13.03%), followed by the EVI (12.47%), NDVI (11.32%), and LAI (10.78%). The NPP also had the highest proportion of declining area (3.85%), followed by the NDVI (3.41%), EVI (2.79%), and LAI (2.76%).   Figure 4 shows the spatial distribution of declining and increasing areas based on the NPP, EVI, NDVI, and LAI indicators, indicating that the spatial distribution of the various indicator results is different. The declining and increasing areas based on the EVI, LAI, and NDVI indices are consistent; the results show that increasing areas mainly occur in southcentral China, India, and central South America, whereas declining areas mainly occur in Central Asia. The NPP-based LPD results indicated differences, mainly located in North America, South America, and Central Africa. Compared with the NDVI, EVI, and LAI, no declining areas were observed in North America. Increasing areas were identified in Central Africa and central South America. A significantly declining area was found in the northern region of South America. Table 3 lists the proportion and areas of each and productivity dynamic (LPD) class. The LPD classes indicated by the four vegetation indicators(VIs) showed that all three levels (low, medium, and high confidence) of increasing area account for 23.99% of the global vegetated area and declining area account for 7.00%. The Increase-HC accounted for 2.77% of the total area, and Decline-HC accounted for 0.35% of the total area. Only 1.13% of the total area had the opposite indicator.   Figure 5 shows the global patterns of the LPD classes. Areas identified as Increase-LC were distributed mainly in central South America, Central Africa, and Northern Europe. The most conspicuous region in terms of Increase-HC was distributed mainly on the Loess Plateau of China, West India, and Eastern Europe. The most conspicuous region in terms of Decline-HC was distributed mainly in Central Asia. South America had the highest concentration of conflict areas. There was an increase in confidence and a decrease in the proportion of the declining and increasing areas.

LPD classes for different landcover and climate zone
For the LPD classes in different land cover types, the total declining area (i.e. Decline-LC, Decline-MC, and Decline-HC) proportions for the forest, grassland, cropland, wetland, and bare land accounted for 47.53%, 24.14%, 25.32%, 0.97%, and 2.04%, respectively. The increasing (i.e. Increase-LC, Increase-MC, and Increase-HC) proportions accounted for 43.49%, 25.17%, 28.88%, 1.21%, and 1.25%, respectively. Forest land cover had the highest increasing and declining proportions and areas. Table 4 lists the statistics for the different land cover types.
For the confidence in the different land cover types, bare land had the highest proportion of Decline-HC and the highest consistency in declining areas. This was followed by grassland, cropland, wetland, and forest. Cropland had the highest proportion of Increase-HC and the highest consistency in increasing areas, followed by grassland, Figure 5. Global land productivity increasing (green) and declining (red) area and confidence level (HC, MC, and LC). A higher green degree indicates a high confidence level in productivity increasing areas. A higher red degree indicates a high confidence level in a productivity declining area. Darker green, medium green, and light green indicate Increase-HC, Increase-MC, and Increase-LC, respectively. The darker red, Orange, and yellow indicate Decline-HC, Decline-MC, and Decline-LC, respectively. Purple indicates at least two indicators conflict; gray indicates all four indicators are stable; white indicates at least one indicator has no data. LC, low confidence; MC, medium confidence; HC, high confidence. bare land, wetland, and forest. Forest cover had the largest increasing and declining areas, but the lowest high-confidence proportion. Figure 6 shows the proportions of each climate zone with different confidence levels.
For the LPD classes in different climate zones, the majority of the hyper-arid areas were masked during the data preprocessing process. The total proportion of declining area (i.e. Decline-LC, Decline-MC, and Decline-HC) in the humid, semi-arid, dry sub-humid, and arid climate zones accounted for 60.11%, 10.58%, 19.89%, and 9.42%, respectively. The total proportion of increaseing area (i.e. Increase-LC, Increase-MC, and Increase-HC) in these zones accounted for 59.90%, 13.99%, 21.40%, and 4.72% respectively. Humid climates had the highest proportion of increasing and declining areas. Table 5 lists the statistics for the different climate zone.
For the confidence in the different climate zones, semi-arid regions had the highest proportion of Decline-HC and the highest consistency in the declining areas, followed by arid, dry sub-humid, and humid. Semi-arid zones had the highest proportion of Increase-HC and the highest consistency in increasing areas, followed by dry sub-humid, arid and  humid. Humid regions had the largest increasing area, but the lowest HC proportion. Figure 7 shows the proportions of each climate zone as a function of the different confidence levels.

Characteristic of different indicators for representing land productivity
This study assumes that a single indicator has the potential to assess land productivity changes but lacks accuracy. In particular, various aspects of vegetation productivity dynamics and phenology can lead to land degradation (Ivits, Horion, Erhard, Fensholt, & Penuelas, 2016), which cannot be reflected by a single indicator. Therefore, we selected four indicators representing greenness (NDVI, EVI), cover (LAI), and productivity (NPP). NDVI has been widely used for productivity change mapping and monitoring. EVI, which is a revised NDVI, is sensitive to high vegetation areas and topography (Jiang, Huete, Didan, & Miura, 2008). Vegetation structure  parameters such as LAI indicate solar radiation conditions and canopy structures (Kang et al., 2016). NPP indicates the improvement in the photosynthesis of vegetation (Zhao, Heinsch, Nemani, & Running, 2005) and was selected to reflect the productivity change. The NPP trend is consistent with previous studies in South America (Tum, Zeidler, Günther, & Esch, 2016). The inconsistency between NPP and other indicators has resulted in conflict zones in northern South America, and Increase-LC in southern South America and Central Africa. Each of the four indicator trends can provide insights into land productivity. However, the different indicators may provide varied results for the land productivity declining assessment, which may be related to indicator characteristics and data quality. Based on the characteristics of each indicator, this study provides insights into land productivity dynamics for specific landcover and climate zones. Forest cover had the largest increasing and declining area but the lowest high-confidence proportion. Fang, Baret, Plummer, and Schaepman-Strub (2019) and Shi et al. (2017) indicated that EVI and LAI are more sensitive than NDVI for detecting changes in areas of high vegetation cover. Therefore, the HC classification indicated the most significant changes, and the MC and LC classifications explained the fluctuation in regions with high forest cover. Semi-arid areas had the highest proportion of HC areas and the highest consistency. In arid climate zones, vegetation cover is greatly affected by factors such as precipitation, solar radiation, and phenology. Therefore, the HC statistics indicated the significant trends of the climate and phenological factors, whereas, the LC statistics indicated the different topography, the algorithm for each indicator, and the data quality. Through statistics on the dynamic distribution of land productivity at different levels of confidence, we can understand the consistency of productivity indicators in different regions, and provide a better reference for selecting suitable land productivity indicators in different regions and interpreting land degradation.  (Zhang et al., 2017), (c) global land productivity dynamics from 1999 to 2013 (Cherlet et al., 2018), (d) NPP trend from 2000 to 2015 (Ding et al., 2020), (e) LAI trend from 2000 to 2017 (Chen, Park, et al., 2019), and (f) (Cherlet et al., 2018) global productivity increase and declining and confidence level.

Consistency analysis with other related datasets
A comparison with other land productivity dynamic datasets is shown in Figure 8, including global NDVI and NPP trends from 2000 to 2015 (Ding, Peng, Qiu, & Zhao, 2020), EVI trends from 2001 to 2015 (Zhang, Song, Band, Sun, & Li, 2017), LAI trends from 2000 to 2017 (Chen, Park, et al., 2019), and JRC's LPD dataset from 1999 to 2013 (Cherlet et al., 2018). The time window may affect the trend estimation since the time window for comparing datasets needs to be as consistent as possible. Compared with existing LDN assessment results, the time windows of these datasets are relatively close. We comprehensively considered the consistency of these previous studies on global products with the LPD dataset and comment on the convergence among them. The comparative analysis indicates that some inconsistencies do exist between different single indicators, while an evident consistency could be observed. Eastern China, Central Africa, Western India, and Europe have a high consistency of productivity increase. The degradation results with high consistency are mainly concentrated in Central Asia and southern South America. In the LDP map, the western regions of India and China were where the Increase-HC areas were primarily distributed, whereas the Decline-HC areas mainly occurred in Central Asia, central Australia, and southern Africa. Central Africa, Western Europe, and Central America showed an Increase-MC trend. The consistency of each dataset is basically the same as the LPD high confidence area, and the comparison results indicate the accuracy of the LPD dataset and verify the correctness of the data.
The distribution of the increasing and declining areas is consistent with that reported by other local studies (de Jong et al., 2013;Fang et al., 2018;Wu et al., 2015). The Loess Plateau in China had the main concentration of Increase-HC areas. China has implemented a series of measures to address land desertification (Su & Shangguan, 2019), soil erosion (Zhao, Mu, Wen, Wang, & Gao, 2013), and other environmental issues. The implementation of natural rehabilitation, terracing, tree planting, and warp land dam strategies (Ping, Anbang, & Xinbao et al., 2013) in 1950 on the Loess Plateau may the reason for the significant Increase-HC area in China. Additionally, cultivated fields in China and India show large Increase-HC areas. The use of chemical fertilizers and multiple crops has led to a substantial increase in food production in both India and China (Chen, Park, et al., 2019). In Central Asia, Zhang et al. (2018) showed that an increasing area of grasslands was converted to sparsely vegetated lands from 2000 to 2014. Persistent droughts may be the main factor causing grassland productivity declining and desertification in Central Asia. Jiang et al. (2019) also showed that convergence from 2008 to 2015 in Central Asia had a negative trend owing to decreased precipitation. These findings support the results of the LPD dataset.

Comparision of our LPD and JRC LPD products for selected areas
The JRC developed another LPD dataset (Cherlet et al., 2018), representing the period from 1999 to 2013. Compared with the JRC LPD dataset (Figure 8c), the areas in this study characterized by inconsistent productivity statuses were mainly distributed in the following regions. In Southern Africa, the JRC LPD shows a declining state, but our LPD shows a stable state. In Australia, the JRC LPD shows a declining status whereas the LPD shows a commonly stable state, with relatively few declines at low and medium confidence.
Our LPD results for Australia are consistent with the Australian vegetation assets, states, and transitions (VAST) classification (Metcalfe & Bui, 2016), which indicates that the state of vegetation in most regions has remained stable with the simultaneous acceleration of deforestation. The land degradation economics in Africa (ELD Initiative & UNEP, 2015) show that Southern Africa is in a stable state with respect to biomass, which is consistent with the LPD results. The ELD Initiative & UNEP (2015) report also showed that there is no consensus on the exact quantity of productivity losses owing to land degradation in eastern Africa. The mainland cover is grassland in southern Africa and central Australia, where the productivity assessment fluctuates significantly based on remote sensing analyses. Grassland productivity is highly inflexible based on climate factors (Moore et al., 2018) and radiation (Moore, Beringer, Evans, Hutley, & Tapper, 2017). Therefore, changes in the actual situation are difficult to interpret using the true ground degradation state. These areas can be characterized as stable and declining. LC is based on the results of this study, which is consistent with the actual situation. Figure 9(a) and (b) show the land degradation in Africa and the VAST classification for Australian vegetation, respectively.
The LPD dataset in this study differs from the JRC LPD dataset in terms of the spatial distribution of the productivity dynamics. This inconsistency is mainly due to different methods and LPD concerns. The JRC LPD uses the multi-temporal image differencing (MTID) method to analyze the trends and changes in the time series, and adopt the Local Net Scaling method to retrieve performance, then integrates the three indicators (trend, change, and performance) to determine the increases in productivity and degradation levels. The JRC LPD is divided into five levels (persistent severe decline, persistent moderate decline, stable but stressed, stable, and persistent increase), indicating the levels of decline and increase with respect to productivity. The LPD dataset considered different parameters representing greenness, cover, and productivity and generated an LPD class with different confidence levels based on the convergence, which provides a new prospect for land productivity dynamic assessment.
The Mann-Kendall algorithm used in the present study is recommended for analyzing the productive trend in the GPG framework. In the land degradation assessment framework SDG 15.3.1, land productivity is a combination of performance, trend, and state. The LPD dataset only uses the land productivity trend without considering the two indicators of performance and state, which is a limitation of the land degradation assessment described in the present study. Even so, both datasets can provide a reference for the evaluation of SDG 15.3.1 and serve to promote target land degradation neutrality.

Conclusions
Evaluations of different productivity declining processes and scales have widely used four indicators (i.e. the NDVI, LAI, EVI, and NPP). The NDVI, EVI, and LAI resulted in identical increasing and declining areas, but the NPP results were different from the three other results in South America and Africa. By combining the four indicators, the results can provide information on both the productivity and confidence levels, which can reflect the actual land productivity declining and increasing situation. Therefore, this study developed a new LPD dataset with a confidence level by combining the time series trends from the NDVI, NPP, EVI, and LAI based on the consistencies between them; the dataset has significant potential for SDG 15.3 monitoring and reporting.