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

: 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 ﬁne 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 Efﬁciency (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 ﬁne-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 − 2 d − 1 ) and bias ( − 0.09 g C m − 2 d − 1 ) than the coarse-spatial-resolution GPP products. In particular, the GPP L signiﬁcantly 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.


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 2 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-resolutionsatellite 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.

Revised EC-LUE Model
The EC-LUE model framework has the following form [12]: where FPAR is the fraction of absorbed photosynthetically active radiation calculated by fine-spatial-resolution normalized difference vegetation index (NDVI), calculated as: 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 −2 ); ε max (LUEmax, g C MJ −1 ) 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 2 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: 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: where VPD0 is biome-specific half-saturation coefficient of the VPD (in k Pa). The Cs is calculated as: with Ci as leaf internal CO 2 concentration [12]. The ϑ in Equation (5) is the CO 2 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.

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).

1.
We extracted the cloud-free pixels by using the cloud mask in the Landsat surface reflectance product.

2.
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.

3.
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.

4.
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.

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: where the original SW 6h is in W m −2 . After transferring the PAR from Equation (3), it was changed to MJ m −2 day −1 . 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 2 concentration at each 16-day interval was averaged from the nearest month's CO 2 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. 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. Table 3. Model coefficients of the fine-spatial-resolution GPP from the GEE APP. Model coefficients are linked to Equations (1), (4) and (5).

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 2 ), root mean square error (RMSE), mean bias (bias), and relative predict error (RPE) to quantify the GPP estimation accuracy.
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: 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 shortterm gap-filled, and long-term gap-filled reconstructed NDVI: GPP L0 , GPP L1 , and GPP L2 , respectively. Next, we used the R 2 , RMSE, and RPE of each group in nine vegetation types to evaluate the effect of missing remote sensing data on fine-resolution GPP modeling.

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 2 ) of around 0.64 and its RMSE was 2.34 g C m −2 d −1 . For the single vegetation type, the GPP L had a high level of agreement (R 2 > 0.5) with the GPP FLUX for eight vegetation types, excluding the EBF (Table 4). Moreover, the GPP L had a low bias (<1 g C m −2 d −1 ) 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).

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 2 and lower RMSE than the GPP M (R 2 = 0.56, RMSE = 2.58 g C m −2 d −1 ) 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 2 equal to 0.45 and an RMSE equal to 2.88 g C m −2 d −1 . Moreover, the GPP L showed the least bias (−0.09 g C m −2 d −1 ) compared to the GPP M (0.21 g C m −2 d −1 ) and the GPP A (0.60 g C m −2 d −1 ). Most of the observed and modeled GPP were close to the 1:1 line. The GPP L showed higher R 2 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 −2 d −1 ) 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 coarseresolution 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 −2 d −1 ). At the same time, the GPP M and GPP A had a GPP STD of 0.381 and 0.045 g C m −2 d −1 , respectively ( Figure 5B). Moreover, the GPP L during the peak of the growing season (~17 g C m −2 d −1 ) was much higher than the GPP M (~9 g C m −2 d −1 ) and GPP A (~9 g C m −2 d −1 ). 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 −2 d −1 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 −2 d −1 , while the GPP STD in the GPP M and GPP A were 0.134 and 0.005 g C m −2 d −1 , 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 −2 d −1 in 2016, which suggested a significant underestimation compared to the GPP L .

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 2 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.

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 onestep 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-spatialresolution 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.

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 (Figures 3 and 4; Table 4), which indicated that the GPP L could characterize the spatial and temporal change in GPP better than the coarse-spatialresolution 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 ( Figures 5-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 −2 d −1 , 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 −2 d −1 , 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.

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.

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.  Acknowledgments: This work used eddy covariance data acquired and shared by the FLUXNET community, including the following networks: AmeriFlux, AfriFlux, AsiaFlux, CarboAfrica, CarboEu-ropeIP, CarboItaly, CarboMont, ChinaFlux, Fluxnet-Canada, GreenGrass, ICOS, KoFlux, LBA, NECC, OzFlux-TERN, Swiss FluxNet, TCOS-Siberia, and USCCC. The ERA-Interim reanalysis data were provided by ECMWF and processed by LSCE. The FLUXNET eddy covariance data processing and harmonization was carried out by the European Fluxes Database Cluster, AmeriFlux Management Project, and the Fluxdata project of FLUXNET, with the support of CDIAC and ICOS Ecosystem Thematic Center, and the OzFlux, ChinaFlux, and AsiaFlux offices. This work was also supported by the Ministry of Education, Youth and Sports of the Czech Republic within the National Sustainability Programme I (NPU I), grant number LO1415 and by the project for national infrastructure support CzeCOS/ICOS, reg. no. LM2015061.

Conflicts of Interest:
The authors declare no conflict of interest.