## Data and Methods
### Data Use  
The main remote sensing data used in this study were level-1 geometrically corrected (L1G) Landsat 5, Landsat 7, and Landsat 8 data acquired on 2 February 1995, 21 March 2003, and 27 March 2014, respectively, which were obtained from the United States Geological Survey (USGS) via the EarthExplorer database ([[https://earthexplorer.usgs.gov](https://earthexplorer.usgs.gov)]([https://earthexplorer.usgs.gov](https://earthexplorer.usgs.gov))). However, clouds and cloud shadows obstructed the view of the land surface during dataset acquisition in the study area. Clouds and their shadows are one of the main problems when using Landsat data in the tropics [47; 48]. The brightening effect of the clouds and the darkening effect of cloud shadows create significant sources of noise in the Landsat data and cause problems in cloud detection, which is an initial step in most analyses, including error estimation for vegetation indexes [49]. Clouds and their shadows covered 1.77%, 2.87% and 0.13% of the entire area of Denpasar for Landsat 5 on 2 February 1995, Landsat 7 on 21 March 2003, and Landsat 8 on 27 March 2014. Most clouds and cloud shadows cover the agriculture and grass, settlement, and forest land use types. The area of clouds and cloud shadows for the three main Landsat datasets was determined using the automated cloud-cover assessment (ACCA) algorithm [50]. The application of the ACCA algorithm requires green, red, near-infrared (NIR), shortwave infrared (SWIR), and thermal infrared (TIR) spectral bands in surface reflectance form. To replace the data lost due to clouds and cloud shadows, Landsat data acquired adjacent to the main data are employed to eliminate the impact of pixels lost in GPP estimation. As a consideration before using these data to replace the lost data, first, the Landsat data must have the date of acquisition closest to the primary data to avoid land surface changes due to agricultural processes. Second, the land surface in the main data and the secondary data used as a patch for lost pixels must have similar spectral characteristics, especially for unchanged land surfaces such as built-up areas. Third, the Landsat data used as the missing pixel cover must be cloud- and cloud shadow-free in the area to be patched. Based on these three considerations, a Landsat 5 acquisition on 26 May 1995, a Landsat 7 acquisition on 28 March 2003, and a Landsat 8 acquisition on 11 March 2014 were used as patches for pixels lost in 1995, 2003 and 2014, respectively. However, only 1.11% of the 2.87% pixels lost in 2003 are covered by secondary data, and the Landsat 7 acquisition on 22 April 2003 was utilized to cover the nearly 1.76% missing pixels that remained. The spectral similarity between the main data and secondary data is derived by correlating the remote sensing indices that are used to estimate GPP from VPM. The two remote sensing indices are the enhanced vegetation index (EVI) and the land surface water index (LSWI). EVI was proposed based on a feedback-based approach that incorporates both background adjustment and atmospheric resistance concepts into the normalized difference vegetation index (NDVI) [51]. The EVI employs an additional blue band for atmospheric correction and variable soil and canopy background reflectance [52] and is calculated from blue, red, and NIR bands. The LSWI, which is proposed by [PERSON] et al. [53], is the normalized difference between NIR and SWIR spectral bands. Because the SWIR spectral band is sensitive to vegetation water content and soil moisture, a combination of the SWIR and NIR spectral bands has been utilized to derive water-sensitive vegetation indices. For instance, when the leaf water content or soil moisture increases, SWIR absorption increases and SWIR reflectance decreases, resulting in an increase in the LSWI value [40].