HiLPD-GEE: high spatial resolution land productivity dynamics calculation tool using Landsat and MODIS data

ABSTRACT Land productivity is one of the sub-indicators for measuring SDG 15.3.1. Land Productivity Dynamics (LPD) is the most popular approach for reporting this indicator at the global scale. A major limitation of existing products of LPD is the coarse spatial resolution caused by remote sensing data input, which cannot meet the requirement of fine-scale land degradation assessment. To resolve this problem, this study developed a tool (HiLPD-GEE) to calculate 30 m LPD by fusing Landsat and MODIS data based on Google Earth Engine (GEE). The tool generates high-quality fused Normalized Difference Vegetation Index (NDVI) dataset for LPD calculation through gap filling and Savitzky–Golay filtering (GF-SG) and then uses the method recommended by the European Commission Joint Research Centre (JRC) to calculate LPD. The tool can calculate 30 m LPD in any spatial range within any time window after 2013, supporting global land degradation monitoring. To demonstrate the applicability of this tool, the LPD product was produced for African Great Green Wall (GGW) countries. The analysis proves that the 30 m LPD product generated by HiLPD-GEE could reflect the land productivity change effectively and reflect more spatial details. The results also provide an important insight for the GGW initiative.


Introduction
Land degradation poses a threat to sustainable development in many countries (Cui and Li 2022)approximately 40% of global land is degraded.This degradation directly affects 50% of the global population and threatens approximately 50% of global gross domestic product (GDP) (https:// www.unccd.int/news-stories/press-releases/chronic-land-degradation-un-offers-stark-warningsand-practical).The United Nations General Assembly has adopted 'Transforming Our World: The 2030 Agenda for Sustainable Development in 2015', which specifies 17 sustainable development goals (SDGs) and 169 specific targets (Li, Lu, and Jia 2021).The target SDG 15.3 aims to combat desertification, restore degraded land and soil, including land affected by desertification, drought, and floods, and strive to achieve a land degradation-neutral world by 2030.The official indicator for reporting progress is SDG indicator 15.3.1, which is the proportion of land that is degraded over the total land area under the custodianship of the United Nations Convention to Combat Desertification (UNCCD) (UNCCD 2016;Orr et al. 2017;Cui and Li 2022).Therefore, monitoring SDG 15.3.1 is important for determining progress and implementing effective interventions against land degradation (Liniger et al. 2019).
Three key sub-indicators included in SDG 15.3.1 monitoring comprise land cover, land productivity, and soil organic carbon (SOC), and they were adopted by the eleventh session of the UNCCD Conference of the Parties (UNCCD 2013) and the Inter-Agency and Expert Group on SDG indicators (IAEG-SDGs) in 2017 (Sims et al. 2019).Land productivity is the biological capacity of land to produce and is a measure of the source of all food, fiber, and fuel that sustains humans (Running et al. 1999;Clark et al. 2001).Its value is based on the notion that the loss of vegetation in productive lands can result in land degradation and vice versa.Remote sensing is the main technology for assessing large-scale land productivity and change (Easdale et al. 2018).Among many approaches, trend analysis of time-series vegetation indicators sensitive to vegetation growth is the most robust and operational (Bai et al. 2008;Metternicht et al. 2010;Easdale et al. 2018;Cui and Li 2022).The Normalized Difference Vegetation Index (NDVI) is a useful proxy for land productivity (Bai et al. 2008;Metternicht et al. 2010;Easdale et al. 2018;Cui and Li 2022), because of a relationship identified between NDVI and NPP (Ivits and Cherlet 2013;Rotllan-Puig, Ivits, and Cherlet 2021).The present spatial resolution of the default global Land Productivity Dynamics (LPD) dataset is 1 km.It is classified changes in productivity based on the direction, magnitude, and stability of the SPOT NDVI between 1999 and 2013, developed by the European Commission Joint Research Centre (JRC) (Cherlet et al. 2018).The Food and Agriculture Organization (FAO) of the United Nations upgraded the LPD dataset between 2000 and 2021 using the same approach but with a MODIS 250 m NDVI (FAO 2022).The LPD calculation tool, LPDynR, was developed in 2021 for indicator calculation based on R (Rotllan-Puig, Ivits, and Cherlet 2021).Generally, LPD products have an extended time frame and improved resolution from 1000 to 250 m.However, LPD datasets and appropriate calculation tools for global with medium to high spatial resolution (10-30 m) are not available, although their importance has been reflected in assessing and monitoring land and soil degradation (Baskan, Dengiz, and Demirag 2017), land-use planning and management (UNCCD 2022), promoting sustainable management and equitable distribution of resources (UNCCD 2017), and assessing potential climate-derived impacts on ecosystems and human development (FAO 2022).
Increased accessibility to high spatial resolution data from satellites such as Landsat and Sentinel via Google Earth Engine (GEE) makes high spatial LPD calculations possible (Gorelick et al. 2017;Mutanga and Kumar 2019;Tamiminia et al. 2020).The major challenge is that the acquisition of consistent and comparable large-scale NDVI datasets with at least a monthly temporal resolution is hampered by low temporal resolution caused by a narrow field of view and cloud contamination (Ju and Roy 2008;Cao et al. 2020;Chen et al. 2021).Another constraint is that the ability of users to analyze data at such a temporal resolution and scale is limited by the need for either knowledge of coding and/or for downloading, storing, and processing substantial volumes of data (Gorelick et al. 2017;Tamiminia et al. 2020).
This study describes the development of a GEE-based 30 m LPD calculation tool.This tool can realize 30 m LPD calculations in any spatial range around the world and within any time window after 2013.The specific processes involved are as follows.The gap filling and Savitzky-Golay filtering (GF-SG) algorithm was used to fuse Landsat-8 and MODIS data to produce NDVI dataset with 8 days temporal, and 30 m spatial resolution, and the average annual NDVI was calculated.The LPD calculation algorithm recommended by the JRC was then applied to generate LPD product, as described by Rotllan-Puig, Ivits, and Cherlet (2021).As an example, the 2013-2020 LPD results were obtained in 11 countries involved in the African Great Green Wall (GGW) Initiative.The results were compared with 250 m LPD datasets and multi-temporal high resolution satellite images, and the LPD status was analyzed to provide insights for GGW Initiative.

HiLPD-GEE
GEE can process large amounts of data online.However, because GEE is a freely shared computing resource, user limitations are included in the operation of complex non-parallel codes (Gorelick et al. 2017).Therefore, the HiLPD-GEE tool is divided into the following sub-sections.The first part generates annual vegetation index images (https://code.earthengine.google.com/990bb85a7eaccc182c1ecd9eafc8ed7d).The second part generates the current status map (https:// code.earthengine.google.com/216be17b7b5e032dd2e4db4432aec791).The third part generates the LPD product (https://code.earthengine.google.com/915732d0891c677fce7e0d8746355a4f).The first part is required because the annual vegetation index images are the base dataset.In the second part, users can select ecoregions or land-cover data to obtain the current status map, or they can ignore this part if only long-term change map are preferred as an alternative LPD outcome.The final part is the tool used to obtain the LPD product.At the beginning of each part, parameters can be adjusted, such as position, time window, and threshold value, users can input parameters as needed.

Dataset and pre-processing
The basic vegetation index dataset was mainly generated by Landsat-8 and MODIS image datasets.Land cover or ecoregions datasets were utilized for calculating the current status map.All data preparation and pre-processing are as follows.

Landsat-8 image dataset
The Landsat-8 image collection uses the USGS Landsat 8 Level 2, Collection 2, Tier 1 dataset in the GEE platform (GEE Collection Snippet: LANDSAT/LC08/C02/T1_L2).The dataset image was started on 18 March 2013, and the product was atmospherically corrected.It contained five visible and near-infrared (NIR), one thermal infrared (TIR), and two shortwave (SWIR) bands as well as a Quality Assessment (QA) band generated by the CFMASK algorithm.The red, NIR, and QA bands were used to calculate the NDVI and remove clouds from images.

MODIS image dataset
The MODIS image collection uses the MOD09Q1.006Terra Surface Reflectance 8-Day Global 250 m product collection in the GEE platform (GEE Collection Snippet: MODIS/006/ MOD09Q1).This dataset was atmospherically corrected and synthesized with optimal values selected from eight days per pixel, including two bands of surface spectral reflectance with wavelengths ranging from 620 to 670 and 841 to 876 nm, corresponding to the red and NIR bands, respectively.These were used to calculate the MODIS 250 m spatial and 8 days temporal resolution NDVI.The dataset also included QA band of surface reflectance for cloud removal from the images.

Land cover
The land cover dataset uses the global 30 m land cover product with a fine classification system in 2020 (Zhang et al., 2021), containing 29 land cover types over a renewal period of five years.

Ecoregions
Ecoregions are region-wide ecosystems that represent different combinations of biodiversity and have natural (not political) boundaries (Dinerstein et al. 2017).The RESOLVE Ecoregions dataset (GEE Collection Snippet: RESOLVE/ECOREGIONS/2017) was updated in 2017.It provides descriptions of 846 terrestrial ecoregions worldwide and includes eight realms and 14 biomes, of which six are forests, and eight are not.
The HiLPD-GEE tool automatically pre-processes images, including cloud removal and mosaicking.Clouds are removed by the bitwise calculation of the QA band.After GF-SG fusion (described below), the generated high-quality NDVI time-series dataset with 30 m spatial and 8 days temporal resolution was automatically synthesized to obtain annually averaged NDVI dataset.Annually averaged NDVI images with different row and column numbers were combined using maximum synthesis in the regions of interest.

Gap filling and Savitzky-Golay filtering fusion
The two methods for reconstructing the 30 m NDVI time series for Landsat-8 images are interpolation and fusion with MODIS images (Chen et al. 2021).However, the number of clear and available Landsat-8 images within one year is limited, and reasonable time-series changes cannot be obtained by simple interpolation.Therefore, GF-SG was selected to produce high-quality spatiotemporally fused NDVI time-series dataset.The algorithm combined the advantages of the high temporal and spatial resolution of MODIS and Landsat-8 images.Furthermore, this method is effective and easy to implement on GEE platform (Chen et al. 2021).GF-SG has two main steps, namely gap filling and Savitzky-Golay filtering, as follows.

Gap filling
The discontinuous property of Landsat time-series data was tackled by MODIS images.GF-SG uses bicubic interpolation resampling to produce MODIS 30 m time-series data based on processed MODIS NDVI time-series data.The mixed pixel and matching problems of different source images were surmounted by GF-SG shape matching to obtain a MODIS referenced NDVI time series.This required searching for similar pixels in the neighborhood using the following formulae: (1) where w (x j ,y j ) is the weight of the jth similar pixel, M similar (x j ,y j ) is the value of the jth similar pixel, cor j (x,y) is the correlation coefficient between the jth similar pixel and the target pixel (x, y), and cor min (x,y) and cor max (x,y) are the minimum and maximum correlation coefficients among all the pixels similar to the target pixel (x, y).Although the temporal shape monitored by MODIS and Landsat images were similar, the amplitude differed between the two sensors.Thus, GF-SG corrects the reference curve (shape correction) using Landsat images as follows: where a(x, y) and a 0 (x, y) were solved using the ordinary least squares method to minimize the differences between the M_adjusted and Landsat values.At the end of the first step, GF-SG used the M_adjusted series to replace the missing values in the original Landsat time series to obtain a Landsat-MODIS synthesized time series (Chen et al. 2021).

Savitzky-Golay filtering
Thereafter, GF-SG smoothed the new NDVI time series using weighted SG filtering to reduce residual noise.The weight in SG filtering depends on adjustment factors (W i × Adjust i ), and assuming that the original Landsat value in the synthetic NDVI time series weighs the most, if the standard deviation of the NDVI of adjacent pixels within the spatial window is also small, the points filled in step 1 would have a larger weight.The adjustment factor was calculated as: where std i is the standard deviation of the neighboring pixels of (x, y) at the ith point, and std max is the maximum among std i .Finally, a smoothed Landsat-8 time series we obtained with an eight days temporal resolution (Chen et al. 2021).

LPD calculation
The current mainstream method for calculating global land productivity is dynamic monitoring of land productivity approach recommended by JRC.This method was used to calculate the LPD in the tool to ensure consistency and comparability with other global products.The method comprised a long-term change map and a current status map of land productivity (Rotllan-Puig, Ivits, and Cherlet 2021).The former characterizes long-term changes in land productivity over time.The latter represents differences in the productivity of homogeneous land during the same period and then judges the productivity status of the current land pixels.Either the former map or the combination of both can be used to calculate LPD (Rotllan-Puig, Ivits, and Cherlet 2021).
2.4.1.Long-term change map of land productivity 2.4.1.1.Steadiness index.The steadiness index was calculated by combining the slope of the linear regression with the time-series net change (NC).This method can be adapted to any type of data or time window (Ivits and Cherlet 2013), and more reasonably reflects the overall direction of ecosystem evolution to measure long-term trends in natural systems.
Slope.The linear fitting reducer ee.Reducer.LinearFit () can be used to obtain the slope in GEE, and the scale that represents the trends is reserved.The formulae are as follows: where Y is the average annual NDVI, X is the time (year of sample point), n is the number of sample points, and b 1 and b 0 represent the slope and offset, respectively, which are estimated using the least squares like Equation ( 7).Only the sign of the fitted slope in the linear regression model was retained as the slope, rather than quantity, to avoid significance tests.2.4.1.2.Baseline level.The baseline describes the NDVI level at the beginning of the monitoring time window.The average value of previous years generally serves as the baseline for high, medium, and low levels to avoid abnormal or unique values caused by extreme events in any year.In general, the division threshold is selected according to the proportion of low vegetation productivity areas within a region.Users can select various thresholds for various regions.The 2013-2020 time window was selected by default, and the NDVI average over the previous three years was selected for average value synthesis.The ee.Reducer.percentile() reducer was used in GEE, and the three-level thresholds were divided by tertiles (33%, 67%) of cumulative pixel values of three-year averaged raster images.According to the tertiles, baseline-level images can be low, medium, or high to represent different original productivity levels of intraregional ecosystems.

State change.
State change describes the difference between the starting and ending states during the monitored timeframes.Averaged values were synthesized for the first and last three years of the time series to quantify the levels of change and avoid outliers or unique values caused by extreme events in a single year.Two average images were separated into 10 levels according to the tenths of the cumulative pixel value to obtain the initial and final NDVI levels.The state change was determined by the absolute value of the difference between the two NDVI levels.Users can set the size of the difference levels according to the nature of the study area.By default, difference between initial and final levels of 0, 1-2, or >2 were assigned state change values of 0, 1, or 2, indicating, respectively, no change, small changes, or large changes.

Current status map of land productivity
The current status map of land productivity was calculated using local net production scaling (LNS) (Prince, Becker-Reshef, and Rishmawi 2009;Rotllan-Puig, Ivits, and Cherlet 2021).It represents difference between the productivity of the target pixels and the maximal productivity of homogeneous pixels during the same period, thus characterizing the potential land productivity of current target pixels.Homogeneous pixels can be clustered unsupervised by phenological or other features.The main purpose is to aggregate land units with the same type of ecosystem function and then determine the maximal productivity in the homogeneous land area.This was achieved using fine land-cover or ecoregions data in this tool.LNS was calculated as (Rotllan-Puig, Ivits, and Cherlet 2021): where AP represents the annual productivity of the target pixel (the average annual NDVI) in the monitored time window in this tool, and PP is the maximal potential productivity pixel in the homogeneous area.In general, 90% of the cumulative pixel value was selected to avoid noise and extreme values.

Comprehensive assessment of land productivity
The long-term change map of land productivity integrates the steadiness index, baseline level, and state change.The long-term change map could be used directly as LPD or combined with the current status map together (Table 1).Values 1-5 represent the five classes: declining, early signs of decline, stable but stressed, stable and not stressed, and increasing land productivity, respectively (Rotllan-Puig, Ivits, and Cherlet 2021).

Case study in Great Green Wall (GGW)
A case study of African GGW was used to demonstrate the applicability and performance of HiLPD-GEE in LPD mapping at a spatial resolution of 30 m.The GGW is an ambitious initiative launched in 2007 (Gadzama 2017).It aims to restore 100 million hectares of degraded land in 11 countries, namely Ethiopia, Chad, Niger, Mali, Burkina Faso, Eritrea, Djibouti, Sudan, Senegal, Nigeria, and Mauritania, especially in the Sahel region (Mirzabaev et al. 2022) (Figure 1), which has become a 'hotspot' for SDG 15.3.1.The annual average temperature in the region is >25°C.The elevation of this area is between −330 and 4529 m, and the plateau is concentrated mostly in Ethiopia.The annual rainfall in most areas is between 0 and 1000 mm, and the rainfall span is large and gradually increases from north to south.The climate in the northern part of the region is subtropical and tropical desert, in the southern part the climate is savanna, and a small area in the southeast has a plateau mountain climate.Being affected by climate, the land-cover types in this region are complex but regular.The northern region is mostly desert and bare land, and grasslands, croplands, and some shrubs and trees are found in the south.Across the 11 countries, the Sahel is defined as an area with 100-600 mm of rainfall, excluding Ethiopia and Djibouti (Figure 1).The Pentad precipitation dataset of the Climate Hazards Group InfraRed Precipitation with Station data (CHIRPS) was selected from the GEE platform (Collection Snippet: UCSB-CHG/CHIRPS/PENTAD) (Funk et al. 2015), which began on 1 January 1981, to obtain the Sahel boundary.The area within an isoline range of 100-600 mm precipitation was extracted and then smoothed on ArcGIS10.7 using gap filling and edge-smoothing algorithms.The Sahel is considered one of the most sensitive ecosystems in the world because of its unique geographical location and climate (Gadzama and Ayuba 2016), and the large impact from droughts causes ecological fragility (Zeng 2003).Between 2013 and 2020, 73,583 Landsat-8 images covering an area >9.34 × 10 6 km 2 were processed.The LPD product of 11 countries was analyzed to promote the monitoring and assessment of land degradation in Africa.

Annual NDVI results
The annual average NDVI results for GGW countries between 2013 and 2020 are shown in Figure 2. Because of the low or outlier NDVI time-series values in desert areas, GF-SG cannot accurately synthesize high-quality NDVI dataset.Simultaneously, the vegetation index cannot be used to evaluate land degradation in these areas because vegetation is sparse all year-round and has its own inherent climate and vegetation characteristics.Therefore, cells with null values in any year in the monitored time window were regarded as non-valued, and areas comprising such cells were considered nonvalued areas and were subsequently processed separately.Due to the influence of climatic factors such as precipitation, NDVI increases from north to south, and less vegetation grows in the northern desert areas, so the average annual NDVI is low.The 2013-2020 time series showed an overall increasing NDVI trend, particularly in the plateau of Ethiopia, and the high NDVI isoline moved northward.
The fused NDVI image eliminated the banding phenomenon in the original Landsat images.Figure 3(c,d) show the annual NDVI results calculated using only the original Landsat-8 and fusion images in specific areas.A gradient change in NDVI values is obvious in the spliced places of different scenery Landsat-8 (Figure 3(c)), resulting in banding.Sampling points numbering 203 (Figure 3 (a)) were evenly distributed, then had the 2020 original Landsat-8 NDVI time series for these points extracted and the available values counted.The results show that the original Landsat-8 images had values that were missing especially between 184 and 284 days of the year (Figure 3(e)).The high NDVI between these days directly resulted in a low annual NDVI value in the final calculation (Figure 3(b)).Therefore, the annual NDVI calculated directly using the original Landsat-8 data was less representative.The fusion image uses the time-series trends of eight days temporal resolution MODIS images to fill the gap, rendering the average annual NDVI image more consistent (Figure 3(d)).Therefore, the fusion algorithm of this tool improved the smoothness and accuracy of the annual NDVI dataset.Areas where the LNS was <50% indicate potentially lower land productivity than the maximal potential of the homogeneous area.This area was mainly concentrated in the Sahel region and eastern Ethiopia (Figure 5).The Sahel is located in the Sahelian Acacia savanna area of the ecoregions, and the north is the South Sahara Desert area.The Sahel included two baseline levels, and the regional NDVI heterogeneity was strong.Therefore, productivity for a plot of land in the north was lower than that of the homogeneous area.In contrast, fine-scale clustering is better such as using land-cover types because it is less heterogeneous across regions (Figure 5).

LPD results
Figure 6(a) shows the LPD result generated by using only the long-term change map. Figure 6(b,c) show the results combined with the current status map of ecoregions and land-cover, respectively.Areas of sparse vegetation coverage in non-valued lands were directly assigned as stable and not stressed, because LPD cannot be calculated based on NDVI.
The three results followed approximately the same trend, with no more than one grade of difference.The more obvious difference was that, after combining the current status map, most of the declining land productivity became an early sign of decline.The land productivity trends were intuitively increasing in southern Chad and most of Ethiopia, as shown in Figure 6(a).Mauritania, Senegal, Mali, Burkina Faso, Niger, and Nigeria in the west had large areas of declining land productivity in their southern parts.Some stable but stressed and declining land productivity was evident in the Sahel (Figure 6(a)).The LNS in most areas of the status map were ≥50% (Figure 5), leading to an increase in land productivity in the south (Figure 6(b,c)).The combining of the map with the ecoregions (Figure 6(b)) reveals that declining land productivity occurs in the Sahel region.

Performance analysis
Three local windows were randomly selected within the case study area to better demonstrate the LPD results (Figure 7).Window (a), located in Mali, was comprised of deciduous shrubland.The RGB images of 2013 and 2020 show that the exposed area was of greater concern in 2020 than in the 2013 image and that vegetation coverage had decreased within the seven intervening years.Land cover data in window b located in Eritrea between 2015 and 2020 shows that most of the shrubland that was evident in this area in 2015 was converted into closed deciduous broadleaved forest (fc > 0.4) by 2020. Figure 7(b) shows that the tool also monitored the increase in land productivity in window b.A comparison of RGB images between 2013 and 2020 in window c located in Chad shows that some cropland has become grassland, as well as a significant declining trend in the LPD result.A LPD product with 250 m spatial resolution was produced for comparison using the same LPD method and 2013-2020 MODIS dataset.Figure 7 shows that the 250 m LPD result had an evident block effect, and the mixed pixels significantly reduced the accuracy of the product.The 30 m LPD better described changes in land productivity at the regional scale, and the influence of mixed pixels was significantly reduced.The 30 m LPD results showed more details, allowing meticulous descriptions of the ground situation.

Statistical analysis
We counted the number of pixels in the LPD results derived from the long-term change map and the synthesis of the long-term change map and current status map with ecoregions or land-cover (Table 2).Among all classifications, most (∼50%) were in the stable and not stressed classes, followed by increasing land productivity.Because the pixel value of most of the areas in the local net scaling was ≥ 50%, compared with those in the long-term change map, the land productivity of the synthesized results was generally one level higher.Notes: A, B, and C: long-term change map, the synthesis of long-term change map and current status map merged with ecoregions, and the synthesis of long-term change map and current status map merged with land-cover, respectively.Values of 1-5, declining, early signs of decline, stable but stressed, stable and not stressed, and increasing land productivity, respectively.
The statistics of the long-term change map as the LPD result were analyzed (Figure 8).The stable and not stressed land productivity class covered an area of approximately 6.39 × 10 6 km 2 , which was the highest proportion, followed by increasing land productivity in all 11 countries.The total area of declining land productivity was 7.1 × 10 5 km 2 (7.35%).Ethiopia had the highest proportion (43.75%) of increasing land productivity, accounting for ∼502,318 km 2 whereas Djibouti had the lowest (0.42%), and 91.39% of that country's area had stable and not stressed land productivity.Burkina Faso had the highest proportion (21.31%) of declining land productivity, covering an area of 59,829 km 2 .Among all 11 countries, Eritrea had the lowest proportion (2.4%) of declining land productivity.
The highest proportions of increasing and declining land productivity in the Sahel area were in northern Nigeria, accounting for 20.95% and 18.76%, respectively.The proportions of stable and not stressed and stable but stressed land productivity throughout the Sahel were 70.56% and 11.54%, respectively.Those for declining and increasing land were 9.24% and 8.45%, respectively.In contrast, the GGW countries had increasing and declining land productivity of 16.25% and 7.35%, respectively.The proportion of increasing land productivity exceeded that of declining land productivity by 8.9% in GGW countries, however, the proportion of declining land productivity exceeded that of increasing land productivity by 0.79% in the Sahel region.This reflects the difficulty in restoring vegetation productivity in the Sahel region and the urgency and complexity of GGW initiative in Africa.

Discussion
The impact of selecting different data sources on the tool's results, differences in parameter settings, 30 m result versus existing products, and the limitations of the tool and prospects for future improvements were discussed in the following.
Parameter and data choice can lead to different results.Adding the current status map with selected data sources can reflect differences in pixels at the same level and generate more abundant results.The current status map can be calculated using other data sources, such as homogeneous pixels aggregation through phenology.Calculating with ecoregions leads to obvious regionalization and boundaries.Vegetation productivity is likely to differ within regions because ecoregions tend to divide those with biodiversity (Figure 5).In contrast, land cover divides the same type of land, and the result of the division is less heterogeneous.Although the data sources differed, the general trends were relatively consistent.Both the current status map results show that areas with LNS ≥ 50% were larger, which infers that most pixels in this area are at a higher level of productivity on homogeneous land.
Different choice of parameters and thresholds can lead to different results.For example, the baseline levels are generally selected according to the proportion of arid areas.The division of low and medium baselines is globally estimated at 40% (Rotllan-Puig, Ivits, and Cherlet 2021).In the GGW countries, 33% was selected as the low-medium baseline level division.The classification of vegetation coverage was similarly referenced to select the threshold of state change.Users can set a threshold that suits the area of investigation.
Both the JRC and FAO generated LPD products with a spatial resolution of 250-1000 m worldwide.Compared with the 2004-2020 LPD product generated by FAO using MODIS imagery (Figure 9), this study used masks for unique values in the northern arid desert regions and assigned the LPD result in these regions to the stable class, whereas most of the regions in the FAO product had decreasing land productivity.In areas where vegetation is sparse, the decreasing trend in LPD does not necessarily represent degradation but might be caused by the unique values.Except for specific desert areas in the north, the LPD results were similar to those of the FAO LPD trends in the south.Land productivity was mostly stable in the central region and declining more significantly in Nigeria and Burkina Faso.The difference between the results typically did not exceed one class.This was mainly caused by the different data sources, range of monitoring years, and selection of the thresholds.The 30 m LPD result showed increasing land productivity in the south when the trend of FAO product from 2004 to 2020 was stable, perhaps because the overall trend from 2004 to 2013 was downward.The effectiveness of the GGW project since its launch in 2007 can be visualized to some extent, such as the 11.2 million trees planted, the restoration of 25,000 hectares of degraded land in Senegal, and the restoration of 15 million hectares of degraded land in Ethiopia (Gadzama 2017).
This study also shows some limitations of the tool.The synthesis of NDVI fusion dataset is limited in specific locations with small NDVI values or large annual fluctuations (such as areas of desert or year-round snow cover).The inherent limitation of GEE computing power prevents the simultaneous generation of results over a large area.This requires sub-regional and step-by-step calculations, and users are required to create mosaics of the relevant sub-regional results.Only land cover or ecoregions were used to process the current status map, without consideration of vegetation phenology and growing seasons.The development of higher-precision algorithms for phenological extraction is anticipated to meet this demand.Future research and development will provide further opportunities and platforms for more complex algorithm execution and data processing.In the interim, the algorithm will continue to be simplified to facilitate the generation of large-scale datasets.Most importantly, improvement of the tool and its application in GGW countries and the Sahel region will be continued, and local real values will be used for accurate product verification and calibration.

Conclusions
SDG 15.3.1 progress measuring at fine spatial resolution is a prerequisite for limited resources allocation and integrated land use planning.This study describes the development of the HiLPD-GEE, the first 30 m spatial resolution LPD online calculation tool based on the GEE platform.This tool can produce the 30 m spatial resolution LPD at any global location within any time window after 2013, which could reflect the land productivity change effectively and reflect more spatial details when compared with the 250 m LPD products.The most obviously technical advances lie in the spatio-temporal inconsistency at the large scale between 30 m Landsat dataset is resolved through the fusion of Landsat and MODIS data and the feasibility of land productivity monitoring at global scale with high spatial resolution enabled by the GEE platform.With the HiLPD-GEE tool, we generated the 30 m LPD products for the GGW countries and Sahel regions on the continental scale for the first time, which shows the good performance of this tool and proves the challenges of the GGW initiative in the future through analyzing the LPD in the GGW countries and Sahel regions.
There are still many challenges for improving the HiLPD-GEE in the future.The calculation efficiency and processing capability need further improvement.Larger-scale product generation requires the simultaneous coordination of increased storage space, computing power, and quantity of multiple data-processing tasks.More field truth data on land degradation are needed to improve the accuracy and effectiveness of 30 m LPD products.More applications in other typical areas globally are necessary to validate the HiLPD-GEE performance at the global scale.Key parameters and thresholds for LPD sub-indicators calculation need to be determined and optimized from the global to region scale, which will require more in-depth and multi-angular considerations.

Disclosure statement
No potential conflict of interest was reported by the author(s).
Net change.The overall NC is calculated, and the positive (+) or negative (−) sign represents the net change in productivity.Positive and negative values indicate, respectively, increasing or decreasing dominant change trends in NDVI (Guo et al. 2008): NC = 2019 i=2013 NDVI 2020 − NDVI i (8) Steadiness index categories 1-4 represent different trends in ecosystems and are determined by the combination of positive and negative NC and Slope pixel-by-pixel values.Steadiness 1 represents negative slope, negative NC, and strong negative dynamics within the time window, thus indicating a declining shifting balance in the ecosystem.Steadiness 2 represents negative slope, positive NC, and fluctuations in the ecosystem, indicating stability within the time window.Steadiness 3 represents positive slope, negative NC, and positive trend of fluctuations in the ecosystem, but stability within the time window in general.Steadiness 4 represents positive slope, positive NC, and strong positive dynamics within the time window, indicating an increasing shifting balance in the ecosystem (Rotllan-Puig, Ivits, and Cherlet 2021).

Figure 4
shows the steadiness index, baseline level, and state change of long-term change in the GGW countries.The steadiness index (Figure 4(a)) indicates that most areas tend to move toward positive equilibrium.The baseline level (Figure 4(b)) was clearly divided into high, medium, and low from south to north, proving that the vegetation productivity base in the south was higher than that in the north.The state change (Figure 4(c)) shows few areas where the state changed by more than two classes, that is to say, the productivity change was almost two classes below or unchanged between 2013 and 2020.In the Sahel area, the steadiness index indicates that almost all dynamics were strongly positive, with few other categories of this index.The baseline level (Figure 4(b)) reveals low productivity in the Sahel area.The combining of this finding with the state change results (Figure 4(c)) indicates that the productivity of the Sahel area has remained low since 2013, mainly due to relatively low precipitation.

Figure 3 .
Figure 3.Comparison of original Landsat-8 NDVI and fusion NDVI images.(a) Result of 2020 fusion NDVI and 203 sample points.(b) Annual difference in average NDVI between fusion and original images of each sample point.(c and d) Annual average NDVI in local window of original Landsat-8 and fusion images.(e) Counts and average statistics of original Landsat-8 NDVI images.NDVI, Normalized Difference Vegetation Index; GGW, Great Green Wall.

Figure 4 .
Figure 4. Sub-indicators of long-term change map of GGW countries and Sahel.(a) Steadiness index.(b) Baseline level.(c) State change.GGW, Great Green Wall.

Figure 5 .
Figure 5. LNS of the current status map in the GGW countries and Sahel.(a) LNS with ecoregions.(b) LNS with land-cover.LNS, local net production scaling; GGW, Great Green Wall.

Figure 6 .
Figure 6.Results of LPD in the GGW countries and Sahel.(a) The long-term change map of land productivity.(b) The synthesis of long-term change map and current status map with ecoregions.(c) The synthesis of long-term change map and current status map with land-cover.LPD, land productivity dynamics; GGW, Great Green Wall.

Figure 7 .
Figure 7.Comparison of 30-and 250-m LPD results and comparison of LPD results with actual RGB images in local areas.LPD, land productivity dynamics.

Figure 8 .
Figure 8. LPD products statistics for GGW countries and Sahel.LPD, land productivity dynamics; GGW, Great Green Wall.

Table 1 .
Rules of LPD calculation with different indicators.

Table 2 .
Numbers of pixels (n) and proportions of results.