Matching high resolution satellite data and flux tower footprints improves their agreement in photosynthesis estimates

Mapping canopy photosynthesis in both high spatial and temporal resolution is essential for carbon cycle monitoring in heterogeneous areas. However, well established satellites in sun-synchronous orbits such as Sentinel-2, Landsat and MODIS can only provide either high spatial or high temporal resolution but not both. Recently established CubeSat satellite constellations have created an opportunity to overcome this resolution trade-off. In particular, Planet Fusion allows full utilization of the CubeSat data resolution and coverage while maintaining high radiometric quality. In this study, we used the Planet Fusion surface reflectance product to calculate daily, 3-m resolution, gap-free maps of the nearinfrared radiation reflected from vegetation (NIRvP). We then evaluated the performance of these NIRvP maps for estimating canopy photosynthesis by comparing with data from a flux tower network in Sacramento-San Joaquin Delta, California, USA. Overall, NIRvP maps captured temporal variations in canopy photosynthesis of individual sites, despite changes in water extent in the wetlands and frequent mowing in the crop fields. When combining data from all sites, however, we found that robust agreement between NIRvP maps and canopy photosynthesis could only be achieved when matching NIRvP maps to the flux tower footprints. In this case of matched footprints, NIRvP maps showed considerably better performance than in situ NIRvP in estimating canopy photosynthesis both for daily sum and data around the time of satellite overpass (R = 0.78 vs. 0.60, for maps vs. in situ for the satellite overpass time case). This difference in performance was mostly due to the higher degree of consistency in slopes of NIRvP-canopy photosynthesis relationships across the study sites for flux tower footprintmatched maps. Our results show the importance of matching satellite observations to the flux tower footprint and demonstrate the potential of CubeSat constellation imagery to monitor canopy photosynthesis remotely at high spatio-temporal resolution.


Introduction
Being able to estimate terrestrial gross primary productivity (GPP) accurately with remote sensing techniques is important for monitoring climate change effects, vegetation responses to environmental extremes and agricultural applications, among other uses. Apart from the chosen approach to estimate GPP from remote sensing variables, the spatial and temporal resolution of the data are important. High spatial resolution satellite data are needed for agricultural applications and monitoring of heterogeneous natural ecosystems. For example, small-scale agricultural fields, which can be fine scale heterogeneous ecosystems, require spatial GPP details to distinguish boundaries between fields and to understand subfield variation in crop growth (Duveiller and Defourny, 2010;Houborg and McCabe, 2016;Kimm et al., 2020). High temporal resolution is also important within short periods during transition season or extreme environmental events such as droughts.
High spatio-temporal resolution of satellite imagery is not only important to monitor vegetation, it is also crucial for the validation of remote sensing-based GPP products with flux tower measurements (Chen et al., 2011;Ran et al., 2016). While the footprint of eddy covariance flux towers covers a relatively large area (more than 1 km 2 for 80 % cumulative contribution of fluxes) (Chen et al., 2009;Chen et al., 2012), its location, size and shape change continuously, driven by surface roughness and meteorological factors such as wind direction and speed (Kljun et al., 2015;Prabha et al., 2008;Schmid, 1997). Although footprint variations can affect the feasibility to match flux tower GPP with satellite products in any type of ecosystem due to spatial heterogeneity within the flux tower footprint (Giannico et al., 2018), it is a key factor in heterogeneous landscapes as contributions to the tower GPP can come from different land cover types (Ran et al., 2016). Due to the limitations in the spatial resolution of available satellite imagery, this aspect of flux tower footprint shape, size and dynamics has mostly been circumvented by selecting only flux tower sites within large patches of relatively homogeneous ecosystems (Liang et al., 2019;Zhang et al., 2020). Even in such cases of consistent land cover within the flux tower footprint, however, the assumption of spatial homogeneity of GPP within the flux tower footprint might not be justified. An important aspect related to this is that even for tower-based, nearsurface sensing techniques, covering the entire flux tower footprint with tower-based optical sensors is virtually impossible even when a bi-hemispheric viewing geometry is used (Gamon, 2015;Liu et al., 2017;Marcolla and Cescatti, 2018).
Remotely monitoring GPP at both high spatial and temporal resolution has been challenging due to the trade-off between spatial and temporal resolutions in observations from sun-synchronous satellites. Most publicly available satellite-derived surface reflectance products, therefore, have either coarse spatial resolution (e.g., 250 m -1 km) with high revisit frequency (e.g., one day for Moderate Resolution Imaging Spectrometer; MODIS) (Justice et al., 1998) or fine spatial resolution (e.g., 10 m -30 m) with revisit frequency (e.g., up to 16 days for Sentinel 2 and Landsat 8) (Claverie et al., 2018;Drusch et al., 2012;Roy et al., 2014), which can effectively result in much lower temporal resolution in many areas due to frequent cloud cover. While many image fusion techniques have been developed (Zhu et al., 2018) and rather successfully applied to partly overcome the trade-off between spatial and temporal resolutions, these approaches still have important limitations. For example, forcing data with a coarse spatial resolution can result in bias against ground measurements, especially in heterogeneous landscapes (Kong et al., 2021).
A recent solution to spatiotemporal trade-offs in satellite observations have been satellite constellations, which consist of a high number of, typically small satellites (i.e., CubeSat; (https://www.cubesat.org/)). CubseSat consists of multiple cubic units (1 unit = 10 cm × 10 cm × 10 cm) but weights, less than 1.33kg unit -1 . On the other hand, it introduces another challenge of crosscalibration between a large number of sensors. The Cubesat-enabled Spatio-Temporal Enhancement Method (CESTEM) overcomes the challenges of low radiometric quality and cross-sensor inconsistency among CubeSat images by radiometric normalization using satellite images from rigorously calibrated sensors (Helder et al., 2020;Houborg and McCabe, 2018a). Indeed, CESTEM outperformed other image fusion products in capturing spatial and temporal variation in an in situ NDVI dataset for a heterogeneous rice paddy landscape (Kong et al., 2021) and has shown great potential for agricultural applications (Aragon et al., 2018;Aragon et al., 2021;Houborg and McCabe, 2018b).
Moreover, Planet Fusion, which is based on CESTEM algorithm, conducts both inter-sensor radiometric harmonization and gap-filling process to deliver daily cloud-free 4-band (blue, green, red, and near infrared; NIR) surface reflectance data with 3 m resolution (Planet Fusion Team, 2021).
Recently, NIR radiation reflected from vegetation has shown promising results in terms of estimating GPP at the field scale of a diverse range of vegetation types Dechant et al., 2022;Dechant et al., 2020). NIR reflectance from vegetation is approximated using the NIRv index as the product of the normalized difference vegetation index (NDVI) and NIR reflectance (Badgley et al., 2017). It accurately estimates monthly and annual GPP variations over globally distributed flux tower sites (Badgley et al., 2019). However, reflectance does not include information regarding the amount of incoming photosynthetically active radiation (PAR), which is a dominant factor driving daily GPP variation. Therefore, NIR radiation reflected from vegetation (NIRvP), which is the product of NIRv and PAR, has been proposed as a structural proxy for GPP estimation (Dechant et al., 2022;Dechant et al., 2020;Wu et al., 2020) and has been widely applied at hourly to daily timescales for site-level and larger scale of GPP estimation since then (Baldocchi et al., , 2022Dechant et al., 2020;Jiang et al., 2021;Liu et al., 2020). Baldocchi et al. (2020) found that NIRvP from in situ spectral sensors tracked diurnal to seasonal variations of tower-based GPP well for individual sites, but that the slopes in the NIRvP-GPP regressions varied considerably across sites and years. As the study of Baldocchi et al. (2020) included spatially heterogeneous sites such as wetlands with considerable variation in water extent and cropland with frequent mowing, the variation in regression slopes could at least partly be due to the footprint mismatch between spectral sensors and eddy covariance systems on the flux towers. This aspect is explored in detail in our study by generating daily, 3 m GPP maps from CubeSat NIRvP and considering different scenarios to link the satellite-based GPP to the tower observations. In particular, we attempted to answer the following scientific questions: Does matching daily flux tower footprints with high resolution (3 m) NIRvP improve correlation to GPP compared to a larger fixed area around the tower and in situ spectral measurements?
(ii) How do the spatial and temporal resolution of satellite data impact the relationship between NIRvP and GPP when matching the flux tower footprint?

Study sites
To evaluate Planet Fusion-based GPP estimates, , we selected study wetland and crop sites that had in situ measure measurements of NIRvP in the Sacramento-San Joaquin River delta in California, USA ( Figure 1), which are registered in AmeriFlux and FLUXNET (Baldocchi et al., 2001;Pastorello et al., 2020). The five study sites included three restored wetlands (US-Tw4, US-Myb, and US-Snf) and two crop fields (US-Bi1 and US-Bi2 which are alfalfa and corn, respectively), located on Sherman, Bouldin, and Twitchell Islands (Table 1). Ideal conditions for eddy covariance measurements include a homogenous and flat landscape and few disturbances, which is physically impossible in the real world (Chu et al., 2021). For example, the wetland sites had a complex mosaic of water and vegetation. The alfalfa site underwent repeated cuttings and regrowth over an annual course. The corn site had a short growing season, and was plowed and flooded for a period during the year Moreover, the summer growing season is mostly cloud free, hence a suitable venue for evaluating satellite products.US-Bi1  (Eichelmann et al., 2018;Valach et al., 2021). US-Snf (Kusak et al., 2020) (Pasture) is a wetland, more specifically peatland (Kasak et al., 2021) and is also used as a pasture to graze cattle throughout the seasons. The peak growing season over pasture site was around December to June. All study sites are registered with the AmeriFlux Network with publicly available data (Hemes et al., 2019).  we applied a neural network procedure (Moffat et al., 2007) to fill the gaps in 30-min NEE time series and then partitioned NEE into ecosystem assimilation and ecosystem respiration flux densities using the nighttime approach (Reichstein et al., 2005). A detailed description about flux data processing, gapfilling and partitioning can be found in Eichelmann et al. (2018) and Knox et al. (2015). .
In situ NIRv data were derived from NDVI sensors (Decagon SRS-Ni, Pullman, WA, USA) that measured incident and reflected radiation in the red (630 nm; full-width, half-maximum, 50 nm) and NIR (800 nm; full-width, half-maximum, 40 nm) spectral bands. NDVI sensors were mounted on booms extended from the flux tower (Table 1). We measured bi-hemispheric reflectance by using hemispherical view for both upward and downward looking sensors to increase the area of the footprint viewed by the NDVI sensors . Although the field of view of the NDVI sensors with bi-hemispheric view include the tower structure, this is expected to be a small, and temporally constant bias. In particular, the NDVI of the tower structure elements is expected to be very low, effectively decreasing the contribution to NIRv. Planet Fusion (Section 2.3) provides a bidirectional reflectance factor that differs from bi-hemispheric reflectance; therefore, we tested the impacts of the different reflectance quantities on NIRv over two study sites using satellite-derived bidirectional reflectance distribution function (BRDF) parameters and the result demonstrated a negligible difference ( Figure A1). Moreover, radiometric calibration had been done for each NDVI sensor before and after the growing season but slow drift year after year in NDVI sensor over the site (US-Tw4)was not a scope of uncertainty in this study as we only used single year data. At each tower, in situ incoming PAR data (μmol m −2 s −1 ) were measured by quantum sensors (PAR-Lite or PQS1, Kipp & Zonen, Delft, Netherlands).
To align the time between continuous in situ measurements and Planet Fusion data that rely on morning to around noon overpass time of satellites such as Sentinel 2 (around 1015hh -1030hh) and Landsat 8 (around 1015hh -1030hh), and MODIS nadir BRDF-adjusted reflectance (local solar noon), we averaged in situ measurements (i.e., GPP, PAR, and NIRv, hereafter NIRv in situ ) between 1030hh. and 1230hh to align it to Planet Fusion data. For daily GPP and in situ daily NIRvP (daily NIRvP in situ ), we summed half-hourly GPP and half-hourly in situ NIRvP over the whole day (0000hh -2400hh) to obtain daily sum values, respectively.

Planet Fusion NIRvP
NIRvP was calculated on the basis of Planet Fusion (PF) surface reflectance data (Planet Fusion Team, 2021 in theory, we assumed the time window between 1030hh and 1230hh when aligning it to in situ data. MODIS affects the daily variation of Planet Fusion more than Sentinel 2 and Landsat. Specifically, MCD43A4 products (Schaaf et al., 2002), which are forcing data for Planet Fusion, use the data from both MODIS terra (overpass at 1030hh) and MODIS aqua (overpass at 1330hh).
where ρNIR and ρRed are reflectance in the NIR and red regions, respectively  Kobayashi and Iwabuchi (2008)) with an artificial neural network. The daily PAR BESS data showed a strong linear relationship to daily PAR data with little bias (R 2 = 0.97, relative bias = -0.8%) ( Figure A2).

Flux footprint model
Flux footprints for each of the sites were calculated with the model from Kormann and Meixner (2001).
This model outperformed other commonly-used flux footprint models according to a CO2 gas release experiment (Kumari et al., 2020;Rey-Sanchez et al., in prep). Flux footprints were calculated for those half hours between 1030hh -1230hh, and later aggregated at the daily level to produce 50, 60, 70, and 80% footprint contours. An 80% footprint contour indicates that 80% of the sources/sinks contributing to the signal detected by the eddy-covariance tower are located within the indicated area. Multiple contours were calculated to obtain a visual inspection of the exponential decay of the footprint past the 50 % contour line. The 80 % threshold selection is supported by Chu et al. (2021), who found that footprint area beyond the 80% contour had a minimal influence on the final calculations of footprint representativeness in a study of monthly footprint-weight maps.. To obtain highly accurate footprint contours we calculated aerodynamic canopy height using the algorithm of Pennypacker and Baldocchi (2016), which has been shown to track detailed trends in the growth of wetland canopies at our study site (Kasak et al., 2020) and in other sites (Chu et al., 2018). Roughness length and displacement height were calculated as 0.1, and 0.66 of canopy height, respectively.

Evaluation
We evaluated the performance of the Planet Fusion -derived NIRvP product for GPP estimation by comparing it with in situ measurements from both NDVI sensor and eddy covariance system data (Chu et al., 2021). For that, we compared satellite-derived NIRv (and NIRvP) data with in situ measurements using several metrics. The coefficient of determination (R 2 ) of the linear regression between NIRv(P) and GPP, which equals the square of the Pearson correlation coefficient, was used as the main metric to characterize the strength of the linear relationship and thus the GPP estimation performance of NIRv(P).
The linear regression slope indicates the variation in the relationship of NIRv PF (or NIRvP PF ) to GPP for different ecosystem types. Root mean square error (RMSE), bias, and relative bias were calculated as follows: Eq. (6) where A is the satellite product (i.e., MODIS, HLS, or Planet Fusion), B is in situ NIRv, and E is the mean operator.
First, we directly evaluated Planet Fusion NIRv against NIRv calculated from in situ sensors, which can be considered a direct ground validation. We estimated the footprints of in situ NDVI sensors based on their heights and locations (Table 1). Because NDVI sensors measure hemispherical irradiance, we approximate the 80% footprint of in situ NIRv from the area of 80% of upwelling irradiance, which can be estimated by considering the sensor height above the ground (Table 1)  We compared the Planet Fusion NIRv product within the estimated footprint of the NDVI sensors (NIRv NDVIsensor PF ) and in situ NIRv measurements. We extracted HLS pixels covering NDVI sensor footprints (Figure 3 and A3) and extracted Planet Fusion NIRv data within the inHLS footprint (NIRv inHSL HLS ) to quantify spatial variation in NIRv.
Second, we evaluated Planet Fusion NIRv (and NIRvP) for GPP estimation by comparing it to GPP from the eddy covariance systems. The flux footprint indicates the source area of trace gases detected by an EC system and therefore can be used for pixel selection during CubeSat GPP evaluation. In this study, we used cumulative eddy covariance footprints up to 80% at satellite overpass time for each day (1030hh -1230hh) (ECfootprint). When extracting Planet Fusion pixels within the ECfootprint, we weighted the pixel values based on footprint contribution. We used an 80% accumulated footprint area for the weighting. We split the weighting factors from 0-50%, 50-60%, 60-70% and 70-80% (Eq. (7)). Eq. (7) Where k is each pixel number in the total number(N) of pixels within the footprint area, CLD % is D % footprint contour line.
We also chose a fixed 100 m by 100 m area centered on the flux towers ( Figure 3; Table 2), which is greater than the footprint of in situ NDVI sensor and partly covers the flux tower footprint, . Such fixed footprints have been applied in previous studies because of the use of coarse spatial resolution satellite products (Heinsch et al., 2006;Kim et al., 2006;Verma et al., 2015). We extracted Planet Fusion pixels from the fixed-area footprint (100 m x 100 m)centered on the towers (NIRv 100m PF ) and the estimated footprint area at satellite overpass time for each day (NIRv ECfootprint PF Additionally, we used annual accumulated footprint over Alfalfa (Bi1), Corn (Bi2), and Freshwater wetland (Tw4) ( Figure A3) to evaluate the impact of temporal resolution of Planet Fusion NIRvP on GPP estimation. The comparison of Planet Fusion pixels within annual accumulated footprints to GPP is in ( Figure A9). We also resampled Planet Fusion data to a 30 m resolution using the nearest-neighbor interpolation method without antialiasing, which produced the results closest to in situ measurements when resampling Planet Fusion-like data from 3 to 30 m resolution (Kong et al., 2021). Next, we extracted pixels within ECfootprint to test the effect of spatial resolution on GPP estimation performance by comparison with the original 3 m resolution (Li et al., 2008   To evaluate the time series of NIRvP data against tower GPP data, we used nonparametric singular spectrum analysis, which decomposes and reconstructs time-series data to detect trends, remove longterm trends, and emphasize short-term variation (Ghil et al., 2002;Mahecha et al., 2007). A time series can be regarded as a collection of additive components (e.g., trends, regular oscillations, and noise (Ghil et al., 2002)); therefore, we evaluated both trends (with low-frequency, smooth components) and oscillations (with high-frequency coherence). Singular spectrum analysis in this study consisted of two stages: decomposition and reconstruction. During the decomposition stage, we formed a time window with length L, which is half of the study period (L = N/2, where N is the length of a flux measurement with a unit of days), for each site and slid it along the time series to construct the L × K (K = N-L + 1) Hankel matrix X. Then, we applied singular value decomposition to X to extract the eigenvalues and eigenvectors of XX T . The leading singular value decomposition component (indicated by the largest eigenvalue) typically corresponds to the time-series trend (Alexandrov, 2009); therefore, we reconstructed the trend by inverting the projection from the leading component to a time series with length N. Subsequently, high-frequency (detrended) components were computed as the difference between the original time series and the trend. We conducted singular spectrum analysis using the R package Rssa v1.0.2 (Golyandina et al., 2018).

Results
This section reports the performance of Planet Fusion derived NIRv and NIRvP across different footprint types against in situ NIRv, NIRvP, and GPP.  Figure 5). Thus, the overall performance of NIRvP MODIS in comparison with NIRvP (R 2 = 0.51) was worse than the overall performance of NIRvP NDVIsensor PF (R 2 = 0.79) and NIRvP HLS (R 2 = 0.73) against NIRvP in situ , which had more consistent slopes ( Figure 5,   (Table A5).  (Table A5).
The performance of NIRvP PF in estimating GPP at each site individually was only clearly improved for one (US-Myb, Palustrine wetland) out of the five sites when the satellite footprint was matched to the ECfootprint (Figure 7,  (Figure 7; Table A6).
The slopes of linear regressions between NIRv PF and NIRvP PF and GPP showed large differences for the NDVIsensor and 100 m footprints but converged when Planet Fusion data was matched with the footprint type of ECfootprint ( Figure 6). More specifically, the variability in linear regression slopes of both NIRv and NIRvP with GPP showed the following ranking (in order of decreasing variability): Planet Fusion NIRv (NIRvP) with different footprint types (ECfootprint > 100m > NDVIsensor) > NIRv (NIRvP ) (Figure 7). In the cropland sites, the linear regression slopes of NIRv and NIRvP with GPP were larger at US-Bi2 (Corn, a C4 plant) than at US-Bi1 (Alfalfa, a C3 plant). These trends were also observed in both NIRv PF and NIRvP PF for all footprint types ( Figure 6;  (Figure 7). However, for some sites, and especially for the 100 m fixed footprint type, NIRvP performed slightly worse than NIRv.  (Table A6).  (Table A6).    (Table A6).

Discussion
Overall, we found that matching Planet Fusion NIRvP to the flux tower footprint considerably improved the agreement with flux tower-based GPP across sites compared to using a fixed footprint (100 m, a hectare area around tower; NDVIsensor, the footprint of in situ NDVI sensor) even when compared to the in situ NIRvP observations (Figure 6 and 7). The effects of matching the flux tower footprint on the regression slope between NIRvP and GPP was largest for the two wetland sites and a corn site (Figure 7). For daily GPP estimation from Planet Fusion data, the different approaches showed comparable performances (Figure 10). NIRvP considerably outperformed NIRv for GPP estimation around satellite overpass time, especially when pooling the data from all sites and using Planet Fusion NIRvP. A more detailed discussion on the different aspects is presented in the subsections below.

Flux tower footprint matching and effects of spatial and temporal resolution on GPP estimation
The high spatial resolution of Planet Fusion NIRvP (NIRvP PF ) (3 m) allowed us to match the footprints of in situ NDVI sensors and eddy covariance systems with Planet Fusion pixels, which led to considerable improvement in the agreement between NIRvP PF and flux tower GPP compared to either a fixed footprint (100 m by 100m) or the relatively small area covered by in situ NDVI sensors (NDVIsensor footprint) (Figure 3, 6, and 7). In fact, the performance of Planet Fusion NIRvP in daily flux tower footprints (NIRvP ECfootprint PF ) was even better than that of in situ sensors with bi-hemispheric view, which did not fully cover the EC footprint (Figure 3). This is a remarkable finding as it implies that the uncertainties associated with satellite observations of NIRvP are smaller than the impacts of the footprint mismatch between in situ optical sensors and the flux tower footprint.
An important aspect of our study was to characterize how the footprint mismatch affected the NIRvP-GPP relationships and to identify the factors explaining such patterns. We found that the effects of footprint mismatch had relatively small impacts on the NIRvP-GPP correlation for individual sites but large impacts on the data combined from all sites (Figure 6 and 7). This was due to inconsistent NIRvP-GPP regression slopes when the flux tower footprint was not matched (Figure 7). These results indicate that while the temporal GPP dynamics are relatively consistent in different parts of the ecosystems, the canopy structure factors such as leaf area index, leaf angles and clumping to which NIRvP is sensitive (Dechant et al., 2022;Dechant et al., 2020) can vary considerably within the flux tower footprints. Freshwater wetland (US-Tw4) and Palustrine wetland (US-Myb) showed the highest sensitivity of regression slopes to the choice of footprint schemes (Figure 6 and 10). Although one might first suspect the impact of the water fraction or the spatial patterns of the water (Matthes et al., 2014) rather than vegetation canopy structure in the footprint of optical sensors to explain the different slopes in wetland ecosystems, NIRvP is largely insensitive to water background as both the NDVI factor and NIR reflectance have low values for water (Chen et al., 2018;Weiss and Crabtree, 2011). Therefore, the most likely explanation for the large sensitivity of wetland sites to the footprint type is the spatial heterogeneity in canopy structure. Apart from insufficient spatial coverage (Gamon, 2015), observations from in situ spectral sensors can also be particularly biased as they typically include the tower structure and its surroundings which are not representative of the flux tower footprint as vegetation near the tower is cleared for the safety of the tower structure. The high consistency between NIRvP and GPP patterns at each site during the study period (Figure 11 and S3) is likely due to the strong impact of environmental drivers such as PAR, temperature and humidity.
We found that the NIRvP-GPP relationships indeed improved by footprint matching rather than other factors that could affect the relationship. As NIRvP is only a proxy for GPP but not a direct observation, the comparison of satellite NIRvP with ground GPP is not a direct ground validation and might be affected by other limitations of NIRvP to estimate GPP, in addition to other uncertainties related to satellite vs. ground observations. Therefore, we used two different ways to exclude such uncertainties as alternative factors explaining our findings of improved NIRvP-GPP relationships when matching Planet Fusion with flux tower footprints. First, we compared entirely in situ-based NIRvP-GPP relationships with relationships where the Planet Fusion area matched the footprint of the in situ NDVI sensor. As these results showed rather similar patterns in slope differences between sites ( Figure   6 and 7), the main factor for different slopes indeed seems to be the footprint size and site location.
Second, we also conducted a direct ground-validation analysis for NIRv and NIRvP and found strong agreement between Planet Fusion and in situ NIRv and NIRvP (Figure 4 and 5). This further confirms that the Planet Fusion -based satellite products are reliable and that the patterns in the results of NIRvP-GPP relationships for different footprints are indeed due to the area covered rather than other sources of uncertainty. Although we found the footprint mismatch to be the dominant factor explaining NIRvP-GPP slope differences between sites, the photosynthetic pathway of the vegetation also needs to be taken into account. Among the cropland sites, the slope of Planet Fusion NIRvP against GPP in C4 plants (corn in US-Bi2) was generally higher than that in C3 plants (alfalfa in US-Bi1) ( Figure 6). This higher slopes of C4 compared to C3 plants is consistent with previous studies although we found a smaller difference (Badgley et al., 2019;Baldocchi et al., 2020;Dechant et al., 2022). We found that the gap-filling of Planet Fusion data had the effect of increasing the spread in NIRvP-GPP regression slopes ( Figure 5). This effect, however, was small and mostly due to a single site (US-Myb; Palustrine wetland) and also cannot explain the differences between different footprint types as the gap-filled data was used for all of them.
The high spatial resolution of Planet Fusion improved the performance of NIRvP-based GPP estimation more compared to the high temporal resolution. To quantify the effects of spatial and temporal resolution of satellite imagery on the NIRvP-based GPP estimation results, we conducted further analyses. First, we evaluated if the footprint matching could still be effective at the coarser spatial resolution of 30 m instead of 3 m. Although the performance of 30 m data was still relatively robust, better performance (R 2 increase ~ 0.06) was found for 3 m at some of the sites, especially at the wetland and corn sites (Figure 9; Table A6), where the water surface extent can change at scales below 30 m (Halabisky et al., 2016). In the direct ground validation for NIRv and NIRvP, we found that the coarse 500 m MODIS pixels had poor performance with large slope differences between the sites (Figure 4 and 5) due to low proportion of footprint area within a pixel (McCombs et al., 2019), while the finer 30 m HLS and 3 m Planet Fusion pixels showed robust performance. Nevertheless, we found substantial land surface heterogeneity even within the 30 m satellite pixels ( Figure A8) indicating the advantage of the Planet Fusion high-resolution imagery which cannot be replaced by fusion products based on MODIS-HLS without compromising the performance (Kong et al., 2021). Second, we evaluated the impact of temporal resolution by comparing the results for daily matching of the EC footprint with that of matching the annual footprint ( Figure A9). We found similar performance in the two cases. Although there is daily variation in the footprints at each site, the annual footprints are good representatives of the area measured by the flux tower because of a strongly dominant wind direction related to the typical meteorological conditions in this area on an annual time scale ( Figure A3; Figure   S8). These findings are, therefore, specific to the sites we used and bigger differences between annual and daily footprint matching are expected for other sites with more variable wind directions (Kim et al., 2006).
The near-daily revisit frequency of CubeSats (Roy et al., 2021) allowed us to track canopy dynamics in detail. For the study period in each site, overall, CubeSats provided 49 -56% of daily data when gapfilled data is removed, which is a higher daily retrieval rate compared to HLS (22 -28%) (Figure 4; Table A5). The near-daily Planet Fusion acquisition frequency enabled the detection of surface changes like mowing events ( Figure S3 and S4) (Kong et al., 2021). Nonetheless, cloud-induced data gaps can still present a significant obstacle for uninterrupted land surface monitoring. Over the study sites, the Planet Fusion -based gap-filling process enabled daily tracking of photosynthesis ( Figure 5; Table A5).
Even though the data-gaps were not negligible throughout the study period ( Figure A8), Planet Fusion data were in general successfully gap-filled ( Figure 4) in part due to relatively short temporal gaps (Kong et al., 2021;Luo et al., 2018). While the gap-filling had considerable impacts on NIRv, the performance of Planet Fusion NIRvP compared to in situ NIRvP was robust even for gap-filled data ( Figure 4 and 5). Daily gap-filled NIRvP ECfootprint PF responded to GPP drops caused by drastic surface changes (e.g., mowing) (Figure 8; Figure S3) which agrees with previous findings that Planet Fusionbased evapotranspiration maps are able to capture evaporation fluxes on a day-to-day basis (Aragon et al., 2021). Further improvements in the capacity to track diurnal canopy dynamics may involve spatiotemporal image fusion between CubeSat and geostationary satellites that scan the same area in every 5-15 minutes (Khan et al., 2021) or integration of ECOSTRESS data that measure thermal infrared with 70 m resolution every 1-5 days (Fisher et al., 2020;Li et al., 2021).

Roles of radiation component in GPP mapping
At the time of satellite overpass, we found NIRvP PF estimated GPP better than NIRv PF when pooling data from all sites and also for most sites individually (Figure 6, 7 and 8). When NIRv was first proposed, it successfully explained the variation of GPP at the monthly time scale (Badgley et al., 2017).
As the variations of daily NIRv and, above all, PAR values are mostly averaged out at monthly time scales, NIRv alone can estimate monthly GPP with reasonable performance. However, NIRv alone raises an issue for GPP estimation at the time of satellite overpass and daily scales as PAR plays an important role (Dechant et al., 2022;Dechant et al., 2020). Although PAR variations are most important for estimating diurnal variations of GPP and NIRvP, PAR still plays an important role at the time of satellite overpass due to considerable seasonal PAR variations ( Figure A2 and S3). In fact, when examining our results closely, a stronger tendency towards non-linear relationships to GPP can be observed for NIRv compared to NIRvP when considering sites individually (Figure 7). Somewhat similar non-linear effects have also been reported before for statistical GPP estimation from reflectance measured by the hyperspectral sensor without including PAR information (Dechant et al., 2019).
However, another potential explanation for the non-linear effects at rather low values could be the nonzero NIRv values for senescent or dead vegetation outside of the main growing season when GPP is essentially zero ( Figure S3). As PAR is low during those times, NIRvP has a somewhat dampened nonlinearity compared to NIRv. Assuming that background effects are indeed partly causing the nonlinearities, using soil-adjusted NIRv may further improve the results (Jiang et al., 2021). Overall, the better performance of NIRvP ECfootprint PF for estimating GPP (Figure 6; Table A6) is consistent with both a mechanistic understanding of GPP drivers and the previous literature on the matter.
While some aspects of the differences in performance of NIRv compared to NIRvP for GPP estimation (Figure 6 and 7) might appear counter-intuitive at first sight, a closer examination of the results can explain those patterns. In particular, the difference between NIRv and NIRvP for GPP estimation was considerably smaller for in situ data than for Planet Fusion data and for the latter it was considerably larger for pooled data from all sites than for individual sites. Furthermore, the slope consistency and R 2 values between sites showed similar relative sensitivities to footprints for NIRv PF and for NIRvP PF although the R 2 values differed considerably in absolute values. All these results can be understood by taking into account the effects of gap-filling on NIRv PF , as the gap-filling led to decrease in agreement with in situ NIRv ( NIRv ) ( Figure A8) but this larger disagreement essentially disappeared for NIRvP (Figure 4 and 5). This compensatory effect of PAR may also suggest a critical control of PAR on regulating seasonal GPP variations at the time of satellite overpass. It should be noted that the effect of gap-filling on the performance of NIRv PF was mostly reflected in the correlation to GPP at individual sites rather than in the slope consistency between sites (Figure 4 and 5).
Daily PAR maps enable NIRv PF maps to conduct spatial upscaling of daily GPP estimation from point eddy covariance sites to landscape scale. Currently, estimates of GPP at eddy covariance sites cover only a small area of the terrestrial surface, with limited numbers of sites that are strongly biased geographically (Ciais et al., 2014;Schimel et al., 2015). We found that daily NIRvP ECfootprint PF , which was temporally upscaled from satellite overpassing time, was strongly correlated to daily GPP ( Figure   10) as GPP at a specific time of the day with a cosine correction generally shows a strong linear relationship with daily GPP (Ryu et al., 2012;Zhang et al., 2018). In other words, the relationship between NIRvP ECfootprint PF and GPP around the overpass time resulted in good performance of daily NIRvP ECfootprint PF in estimating daily GPP ( Figure 10). Furthermore, NIRvP includes information regarding canopy structure and the amount of incoming radiation absorbed and scattered by vegetation Dechant et al., 2022;Dechant et al., 2020). Assuming that canopy structure changes are small within a day, the multiplication of daily PAR by Planet Fusion-derived NIRv will robustly estimate daily GPP. Consequently, NIRv ECfootprint PF × daily PAR BESS performed similarly to daily NIRvP ECfootprint PF in estimating daily GPP (Figure 10 and 11). Accordingly, Planet Fusion -derived NIRv combined with daily PAR maps enables us to generate daily GPP maps with a high degree of fidelity ( Figure 10), which is consistent with the finding of previous studies at different scales (Jiang et al., 2021).

Limitations and perspectives
Our findings indicate that further tests of the slopes of the relationships between NIRvP and GPP are required. The study period for each site varied from 6 months to 1 year (Table 1), which may be sufficient for testing seasonal variation; however, it is insufficient for testing whether these slopes are stable over long periods. Several studies have reported strong correlations between NIRvP and GPP Dechant et al., 2022;Dechant et al., 2020), but the consistencies of their slopes over multiple years, which could include severe stress such as drought, have not been evaluated. Our study sites (cropland and wetlands) are well-known for low water stress. While Baldocchi et al. (2020) reported a strong and rather consistent relationship between NIRvP and GPP at an annual grassland site when pooling data from twelve years, the stability of slopes was not investigated. Also, the longer-term robustness of the relationships between NIRvP and GPP could differ between ecosystems. The accumulation of Planet Fusion data over longer periods will allow us to capture interannual variation in canopy GPP and test if the promising findings of Baldocchi et al. (2020) and Badgley et al. (2019) also hold at finer spatio-temporal scales across ecosystems.
Limitations are expected in the performance of NIRvP to capture shorter-term responses of GPP to vegetation stress as NIRvP is entirely based on canopy structure and radiation information and does not include physiological vegetation signals. In particular, NIRv closely approximates the product of the fraction of absorbed PAR (fPAR) times the fraction of photons that escape from the canopy (fesc) (Zeng et al., 2019). Therefore, NIRvP well captures APAR, the product of fPAR and PAR, which is a strong driver of GPP. Although there is convincing evidence that fesc can capture seasonal variations in light use efficiency (LUE) in crops in the absence of strong environmental stress (Dechant et al., 2020;Liu et al., 2020), fesc might not capture the faster physiological responses at the onset of droughts as fesc is driven by canopy structure variables such as leaf area index, clumping and leaf inclination (Badgley et al., 2017;Qiao et al., 2019;Zeng et al., 2019), which typically respond more slowly than leaf physiology. Therefore, remote sensing variables that carry physiological signals related to LUE such as the photochemical reflectance index (PRI) and physiological emission yield of sun-induced chlorophyll fluorescence (ΦF) should be considered to refine NIRvP-based estimates of GPP in situations of strong environmental stress (Dechant et al., 2022;Kimm et al., 2021;Wang et al., 2020). While SIF contains both the structural information of NIRvP and the physiological information of ΦF, the variations of ΦF tend to be so small that separating them from the structural component with the help of NIRvP might be necessary to use this information for improved GPP estimation (Dechant et al., 2022). contribute to biases in our results ( Figure S1). Regarding the temporal patterns, Planet Fusion data will be associated with larger uncertainties during more extended gap periods ( Figure A8; Figure 4). In those periods, NIRvP tracked GPP better than NIRv did, which could be explained by the role of day-to-day variations of PAR in predicting GPP. Moreover, the study sites are mostly cloud-free during the growing season, which is an ideal case, other sites might have broken clouds with shadows. In that case, the spatial mismatch between satellite-derived PAR and NIRv might degrade the performance of NIRvPbased GPP estimation as PAR with coarse spatial resolution cannot provide the radiation information separately between cloudy and cloud-free areas. In future works on cloudy areas, PAR derived from geostationary satellites with a high spatial resolution (Zhang et al., 2021) could be an option although the spatial resolution is still much coarser than the Planet Fusion imagery. Regarding in situ measurements, continuous eddy covariance measurements can be relatively free from outliers compared to snapshot satellite imagery but sampling errors can also lead to outliers ( Figure 6, Figure S4) even after integrating samples at 20 Hz to 30-min data (Moffat et al., 2007;Richardson et al., 2006). At a daily time scale, uncertainty (one standard deviation) in tower-derived GPP is typically assumed to be 15-20% of the measurements (Falge et al., 2002;Hagen et al., 2006;Richardson et al., 2006). Moreover, the assumptions made during the partitioning of GPP from NEE could induce uncertainty (Reichstein et al., 2005). Theoretically, both PAR and GPP should be zero in nighttime but gap-filled GPP have both negative and positive values in the nighttime due to an unideal temperature function of nighttime NEE.
Therefore, the daily (0000hh-2400hh) GPP values are larger than daytime (0800hh-1800hh) GPP ( Figure S9). Since we targeted to estimate daily GPP, which is needed for annual accumulated GPP, NIRvP-daily GPP relationship will be biased when we force nighttime GPP to 0 or exclude nighttime GPP in daily GPP ( Figure S9). In the case of US-Snf (Pasture), the peak growing season over US-Snf (Pasture) was around December to June but the data were used from June to December in this study (Table 1). Hence, when vegetations are less active, the presence of cattle in US-Snf (Pasture), which can affect carbon dioxide fluxes (Baldocchi et al., 2012;Detto et al., 2010), might influence the relationship between NIRvP and GPP. Regarding footprint matching between Planet Fusion NIRvP and GPP, we weighted Planet Fusion pixels within the flux footprints based on the contour lines rather than the spatially fully explicit footprint contributions per pixel. Therefore, we also encourage future studies to develop pixel-based footprint model for calculating the footprint contribution per pixel.
Matching the footprint between Planet Fusion-derived NIRvP and GPP measurements enables us to use the full potential of flux tower data for spatial upscaling. Previous studies have already pointed out the importance of the spatial heterogeneity within flux tower footprints as the primary source of uncertainty in spatial upscaling (Giannico et al., 2018;Ran et al., 2016)  (2021)) would be to use NIRv as proxy for GPP, assuming sufficiently constant PAR in the selected area. To characterize the spatial heterogeneity, one could use semi-variogram (Kim et al., 2006), which is a suitable geostatistical method for analyzing spatial patterns. Such an analysis could well lead to unexpected results as sites that were previously characterized as homogeneous landscapes (e.g., cropland) may show considerable spatial heterogeneities ( Figure A8). Apart from the footprint characterization that might also be of general interest to the eddy covariance community, Planet Fusion NIRvP will allow us to also include sites with heterogeneous footprints and landcover for spatial upscaling, which have not previously been feasible as the traditional satellite products with coarser spatio-temporal resolution had to rely on a selection of homogeneous sites to avoid biased results (Chen et al., 2011;Ran et al., 2016). Although the footprint matching also improved temporal correlation of NIRvP against GPP for each site (Figure 8), it was not as much as the convergence of slopes across the sites (Figure 7). NIRvP-based GPP estimation varied within the sites depending on the vegetation types and canopy structures but the slope of NIRvP-based daily GPP against in situ GPP over the study period were similar ( Figure 10) because climate factors were the dominant drivers. Nevertheless, variations of NIRvP-GPP regression slopes could appear in very dense forest, ecosystems with high species diversity, and vertically heterogeneous canopies (Ishii et al., 2004), which were not tested in this study. In particular, evergreen needleleaf trees stay green during winter although their photosynthesis is nearly nil, which can lead to seasonally changing slopes between NIRvP and GPP (Kim et al., 2021).
Our results highlight the potential of CubeSat-based imagery for daily canopy GPP estimation.
Previous studies have applied CubeSat datasets to monitor phenology, leaf area index, NDVI, and evapotranspiration (Aragon et al., 2018;Aragon et al., 2021;Houborg and McCabe, 2018b;Hwang et al., 2020;Kimm et al., 2020). The results of our study suggest that CubeSat imagery can be used to obtain daily canopy photosynthesis data at fine spatial and temporal resolution, thus overcoming Incorporation of inexpensive spectral sensors into flux towers that can monitor NIRv (Garrity et al., 2010;Kim et al., 2019;Ryu et al., 2010) and using existing sensor networks to increase ground observation of NIRv across the globe (Carrer et al., 2021) will be a useful strategy for evaluating and improving CubeSat-based fine-resolution GPP mapping.

Conclusions
In this study, we evaluated the performance of satellite-derived NIRvP as a proxy for GPP by comparing it with in situ GPP from a flux tower network over a heterogeneous landscape consisting of mixed wetlands and croplands. We found that the best agreement between flux tower GPP and NIRvP calculated from the Planet Fusion-based product was achieved when matching the Planet Fusion imagery to the tower flux footprints. In fact, the agreement of Planet Fusion imagery matched to the flux tower footprint was even considerably better than the agreement between in situ NIRvP and tower GPP. The improvement for GPP estimation due to the flux tower footprint matching was evidenced by the higher across-site consistency in linear regression slopes between Planet Fusion NIRvP and GPP compared to either fixed footprint Planet Fusion NIRvP or in situ NIRvP. This indicates a large impact of the mismatch between flux tower footprint and the area covered by optical sensors due to spatial heterogeneity within the temporally varying flux tower footprint. We found the largest effects of such mismatches in wetland sites that have a temporally varying fraction of vegetation in a given area. The performance of Planet Fusion-based NIRvP was robust across sites and times in the growing season both for GPP estimation at the time of satellite overpass and for daily summed values. Accordingly, Planet Fusion NIRvP serves as a viable metric for creating daily, high spatial resolution GPP maps at the landscape scale which could prove very useful for precision agriculture or the monitoring of heterogeneous and dynamic natural ecosystems such as wetlands. We encourage future efforts to test the proposed approach in a wide variety of ecosystems during water-limited or light saturation periods.
Overall, our results demonstrated significant advantages of the high spatio-temporal resolution Planet Fusion products for daily canopy GPP estimation and for upscaling flux data for model validation and parameterization. More generally, our findings indicate the large potential of sensor data fusion for high fidelity vegetation monitoring at unprecedented spatial and temporal resolutions.     ) with GPP measured in study sites (colored circles) around satellite overpassing time. Red lines indicate linear regression model slopes for all sites. R 2 is the coefficient of determination, p-value indicates the significance of the linear regression, and n is the number of samples used in the linear regression model. We assumed that the shape of the cumulative footprint over US-Bi1 (Alfalfa) is similar to that over the nearby cropland site (i.e., US-Bi2, Corn), and we considered the inter-annual variations of cumulative footprints were small. The annual cumulative footprint over US-Tw4 (Freshwater wetland) (2016) and US-Bi2 (2017) were used ( Figure A3). NIRv and NIRvP are given in unitless and unit of μmol m -2 s -1 , respectively.