An Open Data Approach for Estimating Vegetation Gross Primary Production at Fine Spatial Resolution

Abstract

Accurate simulations of the spatial and temporal changes in vegetation gross primary production (GPP) play an important role in ecological studies. Previous studies highlighted large uncertainties in GPP datasets based on satellite data with coarse spatial resolutions (>500 m), and implied the need to produce high-spatial-resolution datasets. However, estimating fine spatial resolution GPP is time-consuming and requires an enormous amount of computing storage space. In this study, based on the Eddy Covariance-Light Use Efficiency (EC-LUE) model, we used Google Earth Engine (GEE) to develop a web application (EC-LUE APP) to generate 30-m-spatial-resolution GPP estimates within a region of interest. We examined the accuracy of the GPP estimates produced by the APP and compared them with observed GPP at 193 global eddy covariance sites. The results showed the good performance of the EC-LUE APP in reproducing the spatial and temporal variations in the GPP. The fine-spatial-resolution GPP product (GPP_L) explained 64% of the GPP variations and had fewer uncertainties (root mean square error = 2.34 g C m⁻² d⁻¹) and bias (−0.09 g C m⁻² d⁻¹) than the coarse-spatial-resolution GPP products. In particular, the GPP_L significantly improved the GPP estimations for cropland and dryland ecosystems. With this APP, users can easily obtain 30-m-spatial-resolution GPP at any given location and for any given year since 1984.

1. Introduction

Vegetation gross primary production (GPP), as the principal component of terrestrial ecosystems, plays an important role in supporting global sustenance and energy cycles, regulating the carbon balance and alleviating the rise in atmospheric CO₂ concentrations and global climate change. However, the modeling of GPP is one of the largest sources of uncertainty in global carbon cycle [1,2]. The current satellite-derived GPP products are based on coarse-spatial-resolution data (500 m~8 km), which cannot reproduce the spatial heterogeneity of GPP and leads to uncertainties in GPP estimates [3]. Previous studies showed that the spatial heterogeneity in farmland led to a 30% underestimation in GPP from the agroecosystem [4], while Chen et al. (2010) showed a 5% overestimation in forest ecosystems [5]. Therefore, the improved capturing of the spatial heterogeneity of landscapes can reduce GPP simulation uncertainty.

With fine-spatial-resolution remote sensing data as the input, light use efficiency (LUE)-based GPP models, which account for both physiological process and vegetation structure change, can precisely characterize the GPP in heterogeneous regions [6,7]. Studies showed that the GPP derived from fine-spatial-resolution satellite-based datasets has a better performance than coarse-resolution GPP products because fine-resolution satellite data can capture the heterogeneity of vegetation cover [6]. For example, a LUE model driven by fine-spatial-resolution remote sensing data can reduce GPP simulation uncertainties by 20% [4]. Gitelson et al. (2012) showed that fine-spatial-resolution GPP driven by a LUE model captured the spatial and temporal changes in GPP in a typical agroecosystem [8]. Past researches focused on how fine-spatial-resolution data improved GPP estimation [9,10], but the lack of multi-site studies prevented the assessment of how spatial heterogeneity influenced GPP modeling globally.

More importantly, the collection of fine-scale remotely sensed data requires a significant amount of storage space, and computing fine-spatial-resolution GPP is time-consuming. For example, the Landsat 8 project acquired around 500,000 images at a 30-meter spatial resolution per year, requiring an enormous data storage capacity (~600 TB) [11]. Before users can execute a study, a large amount of time is spent accessing the data and completing the necessary preprocessing steps. Thus, it is difficult to download such a significant amount of data to derive fine-spatial-resolution GPP for various sites. Using the Google Earth Engine (GEE) platform to generate such fine-spatial-resolution GPP products on demand can save significant time and data storage. Furthermore, the revisit periods of fine-spatial-resolution-satellite remote sensing data (for example, ~16 days for Landsat, ~5 days for Sentinel-2), are much longer than those of coarse spatial resolution, which can provide almost two observations each day in the case of moderate-resolution-imaging spectro-radiometer (MODIS) and advanced very-high-resolution radiometer (AVHRR) series. Therefore, the availability of fine spatial remote sensing data may fail to represent the growing condition of vegetation and lead to high uncertainties in GPP modeling. However, few studies have focused on this topic.

In this study, we used the Eddy Covariance-Light Use Efficiency (EC-LUE) model [12] and long-time-series Landsat data as the input to derive fine-spatial-resolution GPP and compared it to the GPP driven by coarse-spatial-resolution remote sensing data. We asked the following questions: (1) How and where did the fine-spatial-resolution remote sensing data improve the GPP estimation compared to coarse-spatial-resolution GPP? (2) How did the data missing from Landsat affect GPP modeling accuracy? To address these questions, we built a web application (APP) on GEE to generate GPP with Landsat 30-meter-spatial-resolution data and compared the fine spatial resolution GPP with the existing coarse-spatial-resolution GPP products. With this APP, researchers can easily obtain fine-spatial-resolution GPP for their purposes and use it for further research.

2. Materials and Data

2.1. Revised EC-LUE Model

The EC-LUE model framework has the following form [12]:

$$\left.\mathrm{GPP}=\mathrm{FPAR}\times \mathrm{PAR}\times {\mathsf{\epsilon }}_{\max }\times {\mathrm{C}}_{\mathrm{s}}\times \min {\left(\mathrm{T}\right.}_{\mathrm{s}},{\mathrm{W}}_{\mathrm{s}}\right)$$

(1)

where FPAR is the fraction of absorbed photosynthetically active radiation calculated by fine-spatial-resolution normalized difference vegetation index (NDVI), calculated as:

$$\mathrm{FPAR}=1.24\times \mathrm{NDVI}+0.168$$

(2)

The NDVI in Equation (2) is calculated from red and near-infrared band reflectance from different-spatial-resolution remote sensing data (e.g., Landsat, MODIS, AVHRR). The PAR in Equation (1) is the photosynthetically active radiation (MJ m⁻²); ε_max (LUEmax, g C MJ⁻¹) is the maximal light use efficiency, which is biome-specific; Cs, Ts, and Ws represent the downward-regulation scalars for the respective effects of atmospheric CO₂ concentration (ppm), temperature (°C), and atmospheric water demand (vapor pressure deficit, VPD, in kPa) on LUE, respectively; min denotes the minimum value of Ts and Ws. The Ts is calculated as:

$$\mathrm{Ts}=\frac{\left(\mathrm{Ta}-\mathrm{Tmin}\right)\times \left(\mathrm{Ta}-\mathrm{Tmax}\right)}{\left(\mathrm{Ta}-\mathrm{Tmin}\right)\times \left(\mathrm{Ta}-\mathrm{Tmax}\right)-{\left(\mathrm{Ta}-\mathrm{Topt}\right)}^2}$$

(3)

where Ta is the air temperature. The Tmin, Tmax, and Topt are set as 0, 40, and 20.33 °C in the model, respectively. The Ws is calculated as:

$$\mathrm{Ws}=\frac{\mathrm{VPD}0}{\mathrm{VPD}0+\mathrm{VPD}}$$

(4)

where VPD0 is biome-specific half-saturation coefficient of the VPD (in k Pa). The Cs is calculated as:

$$\mathrm{Cs}=\frac{\mathrm{Ci}-\vartheta }{\mathrm{Ci}+2\vartheta }$$

(5)

with Ci as leaf internal CO₂ concentration [12]. The $\vartheta$ in Equation (5) is the CO₂ compensation point in the absence of dark respiration, which is also biome-specific [3]. Since the NDVI data sources are different from those used by Yuan et al. (2019) [12], we needed to optimize the model coefficients based on the time series of continuous Landsat NDVI. The fine-spatial-resolution GPP model framework is shown in Figure 1.

To evaluate the merit of fine-spatial-resolution GPP based on reconstructed Landsat data (GPP_L), we compared GPP_L with the GPP measured at FLUXNET towers (GPP_FLUX). Next, we compared GPP_L with two other coarse-spatial-resolution GPP products based on the revised EC-LUE model with MODIS and AVHRR data as input, named GPP_M and GPP_A [3], respectively. The GPP_M is driven by MODIS-based 500-m LAI product [13]. The GPP_A is based on AVHRR data at 0.05 degrees (~8km) of spatial resolution.

2.2. Reconstructed Landsat 30-m NDVI Data

Around 35% of fine-spatial-resolution Landsat data are missing because of cloud cover and long revisit periods [14]. Before we optimized the model coefficients, we used the following steps to reconstruct spatially and temporally continuous Landsat NDVI data (Figure 2).

We chose the surface reflectance data of Landsat Collection 2 Tier 1. We set a data selection rule to minimize the scan line corrector (SLC)-off effect in Landsat 7 Enhanced Thematic Mapper Plus (ETM+): (i) from 1984 to 2002, we selected all the available Landsat data from Landsat 5 and Landsat 7; (ii) from 2003 to 2011, we only used data from Landsat 5; from 2011 to 2013, we chose Landsat 7 data; (iii) after 2013, we only used Landsat 8 data.

• We extracted the cloud-free pixels by using the cloud mask in the Landsat surface reflectance product.
• We used the cloud-free pixels’ reflectance to derive the NDVI. Before we calculated the NDVI, we inter-calibrated the red and near-infrared reflectance in Landsat 5 Thematic Mapper (TM), Landsat 7 ETM+, and Landsat 8 Operational Land Imager (OLI). The TM and ETM+ have similar band settings and research has shown that their NDVIs are consistent [15]. However, the reflective wavelength differences between ETM+ and OLI still affected the NDVI values [16]. Therefore, we used the coefficients in [16] to normalize the red and near-infrared band reflectance in OLI and ETM+. Next, we derived the NDVI from all the inter-calibrated cloud-free surface reflectance. After calculating the NDVI for each image, we generated the NDVI at each 16-day interval and set the highest NDVI at each pixel during this period as the 16-day interval NDVI.
• We performed NDVI gap-filling and wrote the quality control tag (QCtag) in the data quality layer. For pixels with NDVI observation, we wrote QCtag as 0 and kept the NDVI data. When the 16-day NDVI value was missing, we filled it by the linear interpolation method of two NDVI values of nearby dates. If the two nearby NDVI values were observed within 48 days, we defined this gap-filled NDVI as short-term gap-filled NDVI, and the QCtag was set as 1, which indicated that the reconstructed NDVI was based on short-term gap-filled data. If the two nearby NDVI values were more than 48 days apart, we defined this gap-filled NDVI as a long-term gap-filled NDVI, and the QCtag was set to 2 to indicate that the reconstructed NDVI was based on long-term gap-filled data. The data description of the gap-filled and reconstructed NDVI is shown in

  Table 1. Data quality tag of gap-filled and reconstructed NDVI.

  Quality Definition	Data Gap (x in days) between Two Nearest Observed NDVI in 16 Days Intervals	QCtag	NDVI-Filled Method
  Reconstructed NDVI based on observation	0	0	Not applied
  Reconstructed NDVI based on short-term gap-filling	0 < x ≤ 48	1	Linearly interpolated from the nearest two available NDVI
  Reconstructed NDVI based on long-term gap-filling	48 < x	2	Linearly interpolated from the nearest two available NDVI

  .
• We used the Savitzky–Golay filter with a window size = 3, an adaptation strength = 5 to smooth the NDVI data in the time series of each pixel. We were able to derive the time-series-reconstructed 30-meter-spatial-resolution NDVI of any given study area.

2.3. Generating GPP on GEE

We generated the 30-meter-spatial-resolution GPP at 16-day intervals on the GEE platform. The GEE is an open-access platform and provides access to Landsat 5, 7, and 8 data from 1984 to the present. Before we generated the GPP, we needed to prepare the time-series-reconstructed NDVI data, climatic data, and vegetation classification data. We summarized the climatic data at each 16-day interval. The PAR was transferred to the GEE platform from the short-wave radiation (SW) dataset from the National Centers for Environmental Prediction (NCEP) Climate Forecast System (CFS), calculated as:

$$\mathrm{PAR}=\frac{\sum {\mathrm{SW}}_{6\mathrm{h}}\times 0.45\times 6\times 3600}{1000000\times 16}$$

(6)

where the original SW_6h is in W m⁻². After transferring the PAR from equation (3), it was changed to MJ m⁻² day⁻¹. Air temperature and VPD were derived from the ERA5 dataset (https://www.ecmwf.int/en/forecasts/datasets/reanalysis-datasets/era5, accessed on 20 April 2022). The VPD was calculated from 2-m daily mean dewpoint temperature and 2-m air temperature at each 16-day interval [12]. We used a similar method to obtain a 16-day mean surface pressure. The CO₂ concentration at each 16-day interval was averaged from the nearest month’s CO₂ data in the dataset of Trends in Atmospheric Carbon Dioxide (https://gml.noaa.gov/ccgg/trends/data.html, accessed on 31 December 2021). We used the dynamic land-cover classification map (GLC-FCS30, Zhang et al., 2021 [17], https://data.casearth.cn/thematic/glc_fcs30?lang=en_US, accessed on 15 May 2022), which changed every five years from 1980 to 2020, with a fine classification system at 30-meter spatial resolution to classify the vegetation. The data sources are shown in

Table 2. Data sources for generating GPP on the GEE platform.

Data	Specific Name	Transferred Unit	Spatial Resolution	Data Source or Reference
NDVI	Normalized-difference vegetation index	Unitless	30 m	This study
Shortwave radiation	Downward short-wave radiation Flux surface 6-h average	MJ m−2 day−1	0.5°	NCEP
Surface pressure	Surface pressure	Pa	0.5°	NCEP
2-m air temperature	Mean 2-m air temperature	°C	0.25°	ERA5
2-m dewpoint temperature	Dewpoint 2-m temperature	°C	0.25°	ERA5
CO2	Global monthly mean CO2	ppm	/	https://gml.noaa.gov/ccgg/trends/gl_trend.html, accessed on 31 December 2021
Vegetation classification map	GLC-FCS30		30 m	[17]

.

Next, we optimized the model coefficients for the fine-spatial-resolution revised EC-LUE model. We used the reconstructed fine-resolution NDVI from Section 2.2, 16-day interval climatic data and the FLUXNET sites’ 16-day mean GPP from 2000 to 2014 [18] to optimize the model coefficients with the Markov Chain Monte Carlo approach in Huang et al. (2021) [19]. We quantified the GPP data quality by quality tag (NEE_VUT_MEAN_QC) in FLUXNET2015 dataset. When the 16-day mean NEE_VUT_MEAN_QC was less than 0.75, we defined the GPP at this 16-day point as gap-filled GPP; we defined the others as observed GPP. In the step of optimizing the model coefficients, we generated the GPP at 16-day intervals from 95 selected sites, in a similar way to Zheng et al. (2020) [3]. Next, in each site, we chose the odd years as the calibration dataset and the remaining half as the validation dataset. We optimized the coefficients with all the data in the calibration dataset at all the sites from nine vegetation types, including deciduous broadleaf forest (DBF), evergreen broadleaf forest (EBF), evergreen needleleaf forest (ENF), mixed forest (MF), grassland (GRA), savanna (SAV), shrubland (SHR), wetland (WET), and C3 cropland (CRO). The optimized coefficients of the fine spatial remote sensing product are shown in

Table 3. Model coefficients of the fine-spatial-resolution GPP from the GEE APP. Model coefficients are linked to Equations (1), (4) and (5).

Vegetation Type	εmax (g C m−2 MJ−1)	$\mathit{\vartheta }$	VPD0 (kPa)
DBF	3.26	71.09	1.16
EBF	2.50	49.42	0.94
ENF	2.89	45.69	0.81
MF	2.93	66.87	0.86
GRA	3.37	70.11	1.05
SAV	2.33	53.24	1.89
SHR	1.23	37.21	1.30
WET	2.94	76.79	1.28
CRO	4.50	64.00	1.50

.

2.4. Statistical Analysis

We compared the GPP in each product with the FLUXNET 16-day mean GPP at 193 sites of nine major vegetation types and derived the comparative results of GPP products at three spatial resolutions (i.e., 30-m, 500-m, 8-km). We used the determination coefficient (R²), root mean square error (RMSE), mean bias (bias), and relative predict error (RPE) to quantify the GPP estimation accuracy.

$${\mathrm{R}}^2=1-\frac{\sum {\left({\mathrm{GPP}}_f-{\mathrm{GPP}}_i\right)}^2}{\sum {\left({\mathrm{GPP}}_f-{\mathrm{GPP}}_{f,m}\right)}^2}$$

(7)

$$\mathrm{RMSE}=\sqrt{\frac{\sum {\left({\mathrm{GPP}}_f-{\mathrm{GPP}}_i\right)}^2}{n}}$$

(8)

$$\mathrm{bias}=\frac{\sum \left({\mathrm{GPP}}_i-{\mathrm{GPP}}_f\right)}{n}$$

(9)

$$\mathrm{RPE} \left(\%\right)=\frac{{\mathrm{GPP}}_{i,m}-{\mathrm{GPP}}_{f,m}}{{\mathrm{GPP}}_{f,m}}\times 100\%$$

(10)

The GPP_i, GPP_f, GPP_f,m, and GPP_i,m, n in Equations (7)–(10) represent the modeled GPP, GPP from flux towers, mean GPP data from flux towers, mean modeled GPP, and number of samples, respectively.

Subsequently, we chose some typical sites and compared their seasonal change in GPP with the GPP_FLUX. Next, we selected some typical regions to quantify the spatial heterogeneity of GPP_L, GPP_M, and GPP_A. We used all the pixels in the selected regions and calculated the standard deviation (STD, Equation (11)) to quantify the spatial heterogeneity. The coefficient of variation (CV, Equation (12)) of GPP during all the selected 16-day intervals, was used to characterize the temporal change in GPP at each pixel. The standard deviation and CV were calculated as:

$$\mathrm{STD}=\sqrt{\frac{\sum {\left({\mathrm{GPP}}_i-{\mathrm{GPP}}_{f,m}\right)}^2}{n}}$$

(11)

$$\mathrm{CV}=\mathrm{STD}/{\mathrm{GPP}}_{f,m}$$

(12)

Finally, to quantify the influence of missing Landsat data on fine-resolution GPP estimation, we classified the GPP into three groups of GPP_L driven by observed short-term gap-filled, and long-term gap-filled reconstructed NDVI: GPP_L0, GPP_L1, and GPP_L2, respectively. Next, we used the R², RMSE, and RPE of each group in nine vegetation types to evaluate the effect of missing remote sensing data on fine-resolution GPP modeling.

3. Results

3.1. Evaluation of Fine Resolution GPP

The fine-spatial-resolution GPP (GPP_L) correlated well with GPP measured from 193 eddy covariance towers (GPP_FLUX) (Figure 3A). Overall, the GPP_L had a coefficient of determination (R²) of around 0.64 and its RMSE was 2.34 g C m⁻² d⁻¹. For the single vegetation type, the GPP_L had a high level of agreement (R² > 0.5) with the GPP_FLUX for eight vegetation types, excluding the EBF (

Table 4. GPP validation for different vegetation types. Bold indicates the best statistic. The unit of RMSE and bias is g C m⁻² d⁻¹.

	Landsat R2	MODIS R2	AVHRR R2	Landsat RMSE	MODIS RMSE	AVHRR RMSE	Landsat bias	MODIS bias	AVHRR bias
CRO	0.53	0.37	0.40	4.84	5.59	5.47	−0.81	3.08	3.28
DBF	0.73	0.72	0.61	2.38	2.39	2.85	−0.31	−0.03	0.74
EBF	0.31	0.35	0.40	2.39	2.33	2.24	0.52	−0.07	0.32
ENF	0.66	0.71	0.41	1.80	1.67	2.38	−0.14	0.11	0.59
GRA	0.62	0.63	0.52	2.31	2.28	2.61	−0.08	0.41	0.80
MF	0.65	0.65	0.64	1.98	1.97	2.01	0.16	−0.53	−0.44
SAV	0.72	0.66	0.50	1.33	1.46	1.78	−0.02	0.69	0.78
SHR	0.70	0.68	0.51	0.97	1.01	1.24	0.11	−1.35	−1.62
WET	0.60	0.44	0.50	2.18	2.59	2.46	−0.01	−0.06	−0.17

). Moreover, the GPP_L had a low bias (<1 g C m⁻² d⁻¹) with the GPP_FLUX at all the selected sites among vegetation types. The GPP_L captured the seasonal profiles of the GPP_FLUX at the typical sites of each vegetation type (Figure 4). For the forest stands, the GPP_L tracked the seasonal change in GPP well in the selected sites. In particular, the GPP_L showed similar seasonal profiles to the EBF stand. At the non-forest stands, the GPP_L also captured the magnitude of the GPP_FLUX. The GPP_L tracked the GPP_FLUX well at the savanna site, except at the peak of the growing season (Figure 4H).

3.2. Comparison of GPP Estimates at Three Spatial Resolutions

In order to derive the intercomparison results from the three resolution GPP products, we compared 193 sites’ GPP_FLUX with three resolutions of GPP product in constant dates_. Compared to the other coarse-resolution GPP products (Figure 2), the GPP_L had higher R² and lower RMSE than the GPP_M (R² = 0.56, RMSE = 2.58 g C m⁻² d⁻¹) driven by the 500-meter-spatial-resolution MODIS data. The GPP_A driven by the 0.05°AVHRR data had the lowest accuracy, with an R² equal to 0.45 and an RMSE equal to 2.88 g C m⁻² d⁻¹. Moreover, the GPP_L showed the least bias (−0.09 g C m⁻² d⁻¹) compared to the GPP_M (0.21 g C m⁻² d⁻¹) and the GPP_A (0.60 g C m⁻² d⁻¹). Most of the observed and modeled GPP were close to the 1:1 line. The GPP_L showed higher R² than the GPP_M and GPP_A in five of the nine individual vegetation types (Table 4). The GPP_L significantly improved the GPP modeling, especially for cropland, savanna, and wetland, explaining more than 16%, 6%, and 16% of the variance in the coarse-spatial-resolution GPP, respectively. The GPP_L also corrected the high bias of the GPP (>1 g C m⁻² d⁻¹) for cropland and shrub in the GPP_M and GPP_A. The coarse-spatial-resolution GPP (GPP_M, GPP_A) also followed the seasonal change in the GPP but could not reproduce the low GPP of January (Figure 4C). For the dryland vegetation types (SHR, SAV), however, the GPP_M and GPP_A were not able to reproduce the seasonal cycle of the GPP_FLUX. The GPP_M and GPP_A showed overestimations of GPP in the Mediterranean dryland (Figure 4G)

The GPP_L captured the spatial and temporal change in GPP better than the coarse-resolution GPP in agroecosystems (Figure 5). The coarse-resolution GPP (GPP_M and GPP_A) showed less spatial variance than the fine-resolution GPP. The highest standard deviation of the GPP (GPP_STD) in the GPP_L occurred in August 2017 (4.561 g C m⁻² d⁻¹). At the same time, the GPP_M and GPP_A had a GPP_STD of 0.381 and 0.045 g C m⁻² d⁻¹, respectively (Figure 5B). Moreover, the GPP_L during the peak of the growing season (~17 g C m⁻² d⁻¹) was much higher than the GPP_M (~9 g C m⁻² d⁻¹) and GPP_A (~9 g C m⁻² d⁻¹). The GPP_L also showed that the annual mean GPP was higher in the irrigated regions (blue circles in Figure 5A,C,E) compared to the GPP_M. Overall, the GPP_L reproduced the GPP (0 g C m⁻² d⁻¹ at some northern pixels) of roads and non-planted areas (Figure 5C), but the GPP_M did not indicate these conditions (Figure 5E). The fine-spatial-resolution GPP also characterized the spatial heterogeneity of the seasonal change in GPP. The GPP_L showed that the overall CV of the GPP ranged from 0.75 to 1.23 in this region, but the GPP_M showed that the CV of the GPP ranged from 0.81 to 1.16. In the irrigated region, the highest CVs of the GPP_L and GPP_M were 1.2 and 1.15, respectively (Figure 5D,F). Moreover, the GPP_M showed almost no significant spatial difference in this region.

Dryland ecosystems include savanna, open shrubland, and deserts. In savanna ecosystems, the landscape includes random areas of grass and trees, which leads to high spatial heterogeneity in the GPP. Although the mean values of the GPP were similar to the GPP_L, GPP_M, and GPP_A at the research site (Figure 6), the spatial variance in the GPP was quite different. The GPP_L showed the highest GPP_STD, of around 3.378 g C m⁻² d⁻¹, while the GPP_STD in the GPP_M and GPP_A were 0.134 and 0.005 g C m⁻² d⁻¹, respectively (Figure 6D). The seasonal variations in the GPP showed different patterns (Figure 6C) of pixels representing trees (CV ~0.25) and pixels predominantly representing grass (CV~0.4). For the desert, the GPP_L captured the seasonal change in GPP in an oasis (Figure 7). However, since the GPP driven by MODIS data misclassified the vegetation type in this region, the GPP in this region was close to zero. The GPP_M had a maximum GPP of only around 0.002 g C m⁻² d⁻¹ in 2016, which suggested a significant underestimation compared to the GPP_L.

3.3. The Acquisition of Landsat Data Affects GPP Estimation

The data acquisition of Landsat NDVI strongly affected the fine-spatial-resolution GPP modeling accuracy. Figure 8A shows that the proportions of cloud-free NDVI observations in the total data acquisition ranged from 26% to 43% over all nine ecosystem types. The DBF, SHR, and WET had higher proportions (>40%) of cloud-free NDVI observations than the other vegetation types. A low proportion of cloud-free NDVI was found in CRO and EBF, and therefore, more than 56% of the periods of the NDVI values were reconstructed. In each vegetation type, the RPE in the GPP_L based on observed NDVI (GPP_L0) and the GPP_L based on the reconstructed NDVI from short-term gap-filling (GPP_L1) were similar (Figure 8B), ranging from -5% (ENF) to 22% (CRO). Conversely, the RPE in the GPP_L based on reconstructed NDVI from long-term gap-filling (GPP_L2) was higher than the RPE of the GPP_L0. Apart from the EBF, the lowest absolute RPE was 5% (DBF) and the highest was 49%, in the CRO. The differences between the GPP_L0 and the GPP_L2 were much greater (>30%) in the SHR and CRO than in the other vegetation types. The R² was much higher in the GPP_L0 than in the GPP_L2 in most vegetation types (Figure 8C). Meanwhile, the RMSE was similar among the GPP_L0, GPP_L1, and GPP_L2 (Figure 8D). This is because the mean value of the GPP was higher in the GPP_L0 when it was under clear sky, and the incident PAR was much higher than the GPP_L2 when it was under cloudy sky. Thus, the higher RMSE in the GPP_L0 was not linked to higher uncertainties in the GPP_L0 compared to the GPP_L2. Generally, our results showed that the GPP_L0 was less uncertain than the GPP_L2 for most of the vegetation types.

4. Discussion

4.1. Benefits of the One-Step Process of GPP

In this study, we provided a web APP for generating fine-spatial-resolution (30-m) GPP (Figure A1), which allows users to derive the GPP at any given study area around the world. This tool can be used to obtain and download the GPP on demand, providing one-step GPP processing. This APP has three advantages: high efficiency, workload reduction, and open access. First, we used this APP to generate 15 years of fine-spatial-resolution GPP from 193 sites in less than 3 h. Second, since this APP reduces the workload by eliminating data download and image preprocessing, it easily obtains spatiotemporally continuous GPP simulations. Third, since GEE is an open-access platform, this APP can be accessed by all users. For specific applications, users can use this platform to compare the fine-spatial-resolution GPP to other data. This helps cross-comparisons of GPP models, as well as the comparison of GPP to any given bio-physiological variable (e.g., soil water content, ozone concentration). Specifically, to our knowledge, this is the first APP to provide access to fine-spatial-resolution GPPs for any selected location around the world. Therefore, with this APP, users can easily evaluate GPP changes, especially in heterogeneous regions such as deserts, savanna, and cropland.

4.2. Fine-Spatial-Resolution Remote Sensing Data Improves GPP Estimations

The fine-spatial-resolution GPP product (GPP_L) showed good agreement with the GPP observed at the flux sites (Figure 3 and Figure 4; Table 4), which indicated that the GPP_L could characterize the spatial and temporal change in GPP better than the coarse-spatial-resolution GPP products. The GPP_L significantly improved the accuracy of the GPP simulations in the cropland and dryland ecosystems [20]. Our results suggested that the GPP_L provided more detailed information with which to delineate the GPP distribution (Figure 5, Figure 6 and Figure 7). Fine-resolution images can be used to characterize the spatial heterogeneity of roads, buildings, farmland in cropland ecosystems, and grass and trees in savanna ecosystems. The improved consideration of the seasonal change in GPP in savanna may improve our understanding of its dynamic carbon sequestration [20,21]. Furthermore, the GPP_L of evergreen broadleaf forest can represent the GPP decrease in the dry season better than coarse spatial GPP products (Figure 4C). Therefore, the GPP_L characterizes the complexity of GPP [10] and supports specific applications in spatially heterogeneous regions [4].

The GPP_L represented the magnitude of the GPP because we reparametrized the coefficients in the fine-spatial-resolution EC-LUE model. On the one hand, Nguy-Robertson et al. (2015) reported that at the Nebraska extension cropland experiment site (Figure 5) [22], the maximal daily GPP of soybean could reach 18 g C m⁻² d⁻¹, but both the GPP_M and the GPP_A significantly underestimated this condition. After optimizing the coefficients in the EC-LUE model, the peak of the GPP_L was 17.5 g C m⁻² d⁻¹, which was close to the observed GPP [8]. Thus, using the GPP_L to evaluate the amount of GPP in a cropland area can better characterize the growing condition of crops and precisely estimate their yield [23]. On the other hand, in desert regions, vegetation begins to grow after rainfall, but its distribution is sparse. The land-cover classification in coarse-spatial-resolution data may identify these areas as non-vegetated, so the GPP is set to zero. However, fine-spatial-resolution data can be used to detect very small changes in vegetation growth, as indicated by NDVI, in the form of pixels covered by vegetation [24]. Moreover, although the total amount of GPP in dryland may be lower, its expansion leads to substantial changes in GPP [25]. Therefore, the GPP_L product can precisely simulate the total amount of photosynthesis over a large region.

4.3. Challenges and Future Work

Even though the fine-spatial-resolution GPP can reflect the spatial heterogeneity of the GPP, its modeling accuracy is affected by the acquisition of remote sensing data. Unlike the coarse-spatial-resolution remote sensing data, which have more than two observations at any given location of the land surface each day, the Landsat series satellite produces sparse observations. On average, for each Landsat satellite had 23 observations at the same location each year, but due to cloud cover, 50% of the cloud-free data were missing [14,15,26]. Our results showed that the total amount of reconstructed NDVI based on the long-term gap-filling method is around 50% for most of the vegetation types (Figure 8A), which leads to high uncertainties in the GPP_L (Figure 8B). Similarly, a recent study showed that around 50% of Landsat data are cloud-covered, leading to uncertainties in crop yield estimations [27,28].

This indicates that the modeling accuracy of fine-spatial-resolution GPP is highly affected by the data gaps in the time series of fine-spatial-resolution remote sensing data. There are some solutions for this: first, different sources of fine-scale data can be harmonized, such as in the harmonization of Landsat and Sentinel-2 data [29], but this approach is limited by the shorter availability of Sentinel-2 data; second, researchers use data fusion methods with coarse-spatial-resolution remote sensing data [30,31], as in the case of MODIS-Landsat data generation, but this approach requires the recalibration of the reflectance at each channel [32,33]; third, researchers use a data assimilation method to generate Landsat time-series data with climate data [34], but this approach needs huge amounts of computing resources, so it is not available for large-scale evaluation on the cloud platform. Therefore, future work on mapping the GPP on a continental scale should apply fine-scale spatio-temporally reconstructed cloud-free remote sensing data.

5. Conclusions

We presented a web application for generating fine-spatial-resolution GPP driven by 30-m Landsat time-series data (i.e., GPP_L). The GPP_L precisely characterized the spatial heterogeneity and the seasonal variations in the GPP of complex landscapes, while the coarse-spatial-resolution remote sensing data could not reproduce these conditions. Significant improvements in GPP_L were found in agroecosystems and dryland ecosystems, such as shrubland and savanna, compared to the coarse-spatial-resolution GPP products. The EC-LUE model driven by fine spatial resolution with optimized coefficients also corrected the GPP estimation biases for the cropland and dryland, which could be used to reproduce the total amount of plant photosynthesis. Furthermore, the acquisition of the Landsat data strongly affected the estimation accuracy of the GPP_L. Therefore, future work should apply spatio-temporally continuous reconstructed fine-spatial-resolution remote sensing data to the GPP product. Furthermore, since the land-cover map on the GEE platform changes every five years, future work on GPP mapping should apply the dynamic land-cover map, which can capture the GPP changes caused by annual land cover changes. With our online APP on GEE, users can generate GPP at any given location and for any year since 1984 on demand. This will help with the development of further applications that investigate the relationship between GPP and other parameters.