A convolutional neural network for spatial downscaling of satellite-based solar-induced chlorophyll fluorescence (SIFnet)

Gross primary productivity (GPP) is the sum of leaf photosynthesis and represents a crucial component of the global carbon cycle. Space-borne estimates of GPP typically rely on observable quantities that co-vary with GPP such as vegetation indices using reflectance measurements (e.g., normalized difference vegetation index, NDVI, near-infrared reflectance of terrestrial vegetation, NIR_v, and kernel normalized difference vegetation index, kNDVI). Recent work has also utilized measurements of solar-induced chlorophyll fluorescence (SIF) as a proxy for GPP. However, these SIF measurements are typically coarse resolution, while many processes influencing GPP occur at fine spatial scales. Here, we develop a convolutional neural network (CNN), named SIFnet, that increases the resolution of SIF from the TROPOspheric Monitoring Instrument (TROPOMI) on board of the satellite Sentinel-5P by a factor of 10 to a spatial resolution of 500 m. SIFnet utilizes coarse SIF observations together with high-resolution auxiliary data. The auxiliary data used here may carry information related to GPP and SIF. We use training data from non-US regions between April 2018 until March 2021 and evaluate our CNN over the conterminous United States (CONUS). We show that SIFnet is able to increase the resolution of TROPOMI SIF by a factor of 10 with a r² and RMSE metrics of 0.92 and 0.17 mW m⁻² sr⁻¹ nm⁻¹, respectively. We further compare SIFnet against a recently developed downscaling approach and evaluate both methods against independent SIF measurements from Orbiting Carbon Observatory 2 and 3 (together OCO-2/3). SIFnet performs systematically better than the downscaling approach (r=0.78 for SIFnet, r=0.72 for downscaling), indicating that it is picking up on key features related to SIF and GPP. Examination of the feature importance in the neural network indicates a few key parameters and the spatial regions in which these parameters matter. Namely, the CNN finds low-resolution SIF data to be the most significant parameter with the NIR_v vegetation index as the second most important parameter. NIR_v consistently outperforms the recently proposed kNDVI vegetation index. Advantages and limitations of SIFnet are investigated and presented through a series of case studies across the United States. SIFnet represents a robust method to infer continuous, high-spatial-resolution SIF data.

1 Introduction

Photosynthesis represents the single largest CO₂ flux between the atmosphere and the biosphere. At the canopy level, the sum of all leaf photosynthesis is termed gross primary productivity (GPP), and accurate characterization of GPP represents a major uncertainty in the carbon cycle (Friedlingstein et al., 2019). Directly measuring GPP from remote sensing systems (e.g., satellites) is not presently possible. Instead, previous work has utilized stationary measurements of net ecosystem exchange (NEE) from flux towers that can be decomposed into GPP and respiration (e.g., Reichstein et al., 2005). Observable quantities from satellites (e.g., vegetation indices computed from reflectance data) are then related to GPP inferred from flux towers (e.g., Huete et al., 2006; Jung et al., 2019; Sims et al., 2006; Zeng et al., 2020) in light use efficiency (LUE) (Mahadevan et al., 2008) or machine learning models (Jung et al., 2019) to derive global estimates of GPP.

Vegetation indices such as the normalized difference vegetation index (NDVI) and near-infrared reflectance of terrestrial vegetation (NIR_v) combine two (or more) spectral bands with different absorption characteristics (Hanes, 2013) to infer quantities related to plant physiology and canopy structure. The MODIS instrument was launched on the Terra and Aqua satellites in 1999 and 2002, respectively. This instrument has proved particularly useful due, in part, to the long operational lifetime, and vegetation indices can be derived from the individual reflectance bands of MODIS. More recently launched satellites, like Sentinel-5P, carry instruments with the necessary signal-to-noise ratio and spectral resolution to retrieve solar-induced chlorophyll fluorescence (SIF). The electromagnetic signal SIF is emitted by chlorophylls during photosynthesis. SIF is emitted in the red–far-red wavelengths of 650–850 nm (Magney et al., 2020). It is a way, besides photochemistry and nonphotochemical quenching, for de-excitement of the chlorophylls (Turner et al., 2020; Köhler et al., 2018). Even though the link between chlorophyll fluorescence and photosynthesis is nonlinear at leaf and canopy scale, that does not hold for satellite scales, in which a linear relationship of SIF to GPP is frequently reported (e.g., Turner et al., 2021; Magney et al., 2020; Frankenberg et al., 2011; Joiner et al., 2011).

Vegetation indices (also termed greenness) can be regarded as a measure of photosynthetic capacity (Sellers, 1985), whereas SIF indicates photosynthetic activity. SIF has been shown to be a powerful proxy for estimating GPP (Magney et al., 2019; Turner et al., 2021), to capture the impact of drought on photosynthetic activities across different vegetation types (Shekhar et al., 2020a; Castro et al., 2020), and to assess the regional source of carbon emissions (Shekhar et al., 2020b).

Köhler et al. (2018) described the first retrievals of SIF from the TROPOspheric Monitoring Instrument (TROPOMI), the sole instrument on the Sentinel-5P satellite. The TROPOMI instrument has an equatorial crossing time of 13:30 local solar time and a 16 d orbit cycle. TROPOMI has a wide swath (2600 km across track) that allows for near-daily temporal resolution and a spatial resolution of 5.5 × 3.5 km. This was a substantial improvement to previous satellite instruments measuring SIF that were limited to 40 × 40 km spatial resolution (Joiner et al., 2013). Despite the higher spatial resolution of TROPOMI, there have been efforts to estimate SIF at finer spatial scales (e.g., Turner et al., 2020). This is motivated by the importance of fine-scale phenomena in the carbon cycle such as ecosystem fragmentation (e.g., Haddad et al., 2015).

Globally, 20 % and 70 % of the remaining forests are within a distance of 100 m and 1 km, respectively, from the forest edges, meaning that most of the forests are fragmented (Haddad et al., 2015). Reinmann and Hutyra (2017) show that the carbon uptake and storage of trees near the forest edge increase up to 13 ± 3 % and 10 ± 1 %, respectively. On the other hand most of our understanding about forest carbon fluxes comes from intact ecosystems, resulting in a mismatch between the ecosystems we are trying to quantify and the data we are using to do so (Smith et al., 2018). Higher-resolution estimates of photosynthetic activity might enable us to include fragmentation effects of ecosystems to global carbon cycle estimates or biosphere models like Jung et al. (2019), Wu et al. (2021), and Turner et al. (2021). Additionally, recent work has shown the importance of fine-scale variations in the urban biosphere on the overall carbon flux for a city (e.g., Miller et al., 2020).

There has been some recent work with the goal of increasing the resolution of existing global SIF estimates through downscaling methods (i.e., physics-based methods). For example, Turner et al. (2020, 2021) used NIR_v to partition SIF within a particular TROPOMI scene and oversampled it using a 16 d window afterwards, resulting in a daily 500 m SIF estimate over the conterminous United States (CONUS). Duveiller et al. (2020) downscaled GOME-2 satellite SIF from 0.5 to 0.05^∘ using a parameterization with a term for the fraction of absorbed photosynthetically active radiation (fAPAR), one for water stress, and one for heat stress based on MODIS data. Siegmann et al. (2021) used airborne data to downscale far-red SIF from canopy to leaf level.

Machine learning has also been used to create global, high-resolution SIF data sets. Li and Xiao (2019), Yu et al. (2019), and Zhang et al. (2018) used spectral bands from MODIS as input to neural networks that were trained with Orbiting Carbon Observatory 2 (OCO-2) SIF data to build global continuous SIF products at 0.05^∘ resolution. OCO-2 has a narrow swath, and, therefore, the networks are trained only in the regions where OCO-2 SIF is available by using MODIS data as input. After training, the global MODIS data are used as input to estimate SIF on a global scale. Gentine and Alemohammad (2018) use MODIS reflectance data as input and predict GOME-2 normalized by clear-sky irradiance. Multiplying that with a MODIS-derived photosynthetic active radiation product results in a MODIS only estimated SIF, termed RSIF. Zhang et al. (2021) trained a convolutional neural network (CNN) with MODIS data on the artificial GOSIF data set (Li and Xiao, 2019) at a resolution of 0.05^∘ and used the trained network and MODIS data at a resolution of 0.008^∘ to estimate SIF at 0.008^∘. The physics-based downscaling approach from Turner et al. (2020) can only consider one variable for weighting the SIF signal, while the machine-learning-based approaches in the literature can consider more than one variable – but many do not use SIF data as an input to their model, meaning that they estimate SIF based on reflectance data.

Here we build a convolutional neural network to obtain high-resolution SIF, named SIFnet. SIFnet increases the spatial resolution of TROPOMI SIF by considering coarse-resolution SIF with high-resolution auxiliary data as input. These auxiliary data consist of either proxies of SIF or photosynthetic drivers. SIFnet is trained using data with near-global coverage. Different model parameters (structure, input features, and scaling factors) are compared and evaluated. After training the model, the resolution of TROPOMI SIF is refined by a factor of 10 to a spatial resolution of 0.005^∘. This product is then compared against a recent downscaling method from the literature (Turner et al., 2020). Both high-resolution estimates are validated over CONUS against the independent SIF measurements of the OCO-2 and OCO-3 instruments (together OCO-2/3) (OCO-2 Science Team/Michael Gunson, 2020; OCO-3 Science Team/Michael Gunson, 2020).

2 Data sets

2.1 Data sources

The input data to the neural network are listed in Table 1. These diverse global data products are expected to capture a broad range of photosynthetic drivers. The table divides the data into time-varying or time-invariant and training or validation data. The native spatial resolution is shown in the last column of Table 1. In a first step, all data sets are aggregated to 0.05^∘ spatial resolution and 16 d time steps. In the case of a higher native spatial resolution, the data are regridded by computing the mean value that falls into the coarse-resolution grid cell. In the case of coarser resolutions than 0.05^∘, it is resampled to the common grid. Quality control flags and cloud filtering are applied when necessary.

Table 1 Data sets used in this work. ENF: evergreen needleleaf forest; EBF: evergreen broadleaf forest; DNF: deciduous needleleaf forest; DBF: deciduous broadleaf forest; MF: mixed forest; UF: unknown forest.

References: ¹ Köhler et al. (2018); ² Schaaf and Wang (2015); ³ Copernicus Climate Change Service (2017); ⁴ Entekhabi et al. (2010); ⁵ PySolar (2021); ⁶ Danielson and Gesch (2011); ⁷ Buchhorn et al. (2020); ⁸ Morreale et al. (2021); ⁹ OCO-2 Science Team/Michael Gunson (2020); ¹⁰ OCO-3 Science Team/Michael Gunson (2020).

MODIS measures the reflected radiance from the earth surface in seven different spectral bands covering the visible and infrared spectral region. Vegetation indices are computed by combining the near-infrared (where chlorophyll is non-absorbing) and the red band (where chlorophyll is highly absorbing) (Hanes, 2013). Specifically, the normalized difference vegetation index (NDVI) (Tucker, 1979), near-infrared vegetation index (NIR_v) (Badgley et al., 2017), the kernel NDVI (kNDVI) (Camps-Valls et al., 2021), and the enhanced vegetation index (EVI) (Huete et al., 2002) are computed as follows:

$$\begin{array}{}\text{(1)}& {\displaystyle }\mathrm{NDVI}={\displaystyle \frac{{\mathit{\rho }}_{\mathrm{NIR}}-{\mathit{\rho }}_{\mathrm{RED}}}{{\mathit{\rho }}_{\mathrm{NIR}}+{\mathit{\rho }}_{\mathrm{RED}}}},\text{(2)}& {\displaystyle }{\mathrm{NIR}}_v={\mathit{\rho }}_{\mathrm{NIR}}\cdot \mathrm{NDVI},\text{(3)}& {\displaystyle }\mathrm{kNDVI}=\tanh \left({\mathrm{NDVI}}^{\mathrm{2}}\right),\text{(4)}& {\displaystyle }\mathrm{EVI}=G\cdot {\displaystyle \frac{{\mathit{\rho }}_{\mathrm{NIR}}-{\mathit{\rho }}_{\mathrm{RED}}}{{\mathit{\rho }}_{\mathrm{NIR}}+C_{\mathrm{1}}\cdot {\mathit{\rho }}_{\mathrm{RED}}-C_{\mathrm{2}}\cdot {\mathit{\rho }}_{\mathrm{BLUE}}+L}}.\end{array}$$

EVI coefficients for MODIS are as follows: L=1, C₁=5, C₂=7.5, and G=2.5 (Huete et al., 2002). ρ_NIR is the near infrared band, ρ_RED is the red band, and ρ_BLUE is the blue band from the MODIS satellites.

Temperature and precipitation are taken from ERA5-Land data at the time step of interest and with a delay of one time step. Soil moisture has been shown to be a strong driver of global photosynthesis due, in part, to its impact on vapor-pressure deficit (Humphrey et al., 2021). Here we use the coarse-resolution NASA USDA Soil Moisture Active Passive (SMAP) soil moisture (Entekhabi et al., 2010) as a model input and explore its correlation with TROPOMI SIF. The cosine of the solar zenith angle (SZA) is a proxy for photosynthetically active radiation (PAR) under cloud-free conditions (Chen et al., 2020; Turner et al., 2021).

Time-invariant data sets consist of elevation data, fractional land cover classification, and forest fragmentation data. The land cover classification (Buchhorn et al., 2020) is resampled to 11 fractional classes. The forest fragmentation data consist of two bands and have a native resolution of 30 m. One band describes the share of forest within the grid cell (forest share) and the other how much of that forest is edge forest (defined as a maximum distance to an edge or other land cover type of 30 m). OCO-2 and OCO-3 have high spatial resolution (2.25 × 1.29 km) but small swaths (10 km) and a 16 d revisit time.

2.2 Covariation of input data sets with SIF

We are interested in understanding what these different data sets are telling us about SIF and also how they co-vary with each other. We compare all collected time variant data against TROPOMI SIF in the spatial and temporal domain. As a quantitative measure, we compute the Pearson correlation coefficient (r) (Benesty et al., 2009). Figure 1 shows a scatter comparison of SIF against the auxiliary data at the lowest resolution of the two corresponding sets. Negative SIF values (on the x axis in Fig. 1) are due to relatively high retrieval errors which scale with radiance levels (Köhler et al., 2018).

Figure 1 Scatter comparison of SIF to timely changing auxiliary data. The time span of measurements is from April 2018 to March 2021 at 16 d resolution. Longitude and latitude borders are from −180 to 180^∘ and −60 to 70^∘, respectively. The comparison resolution corresponds to the lowest resolution of the two corresponding products. For all MODIS data the resolution is 0.05^∘ and for precipitation, air temperature, surface soil moisture (ssm), and subsurface soil moisture (susm) 0.1^∘. To quantify the goodness of fit we compute the Pearson correlation coefficient (r) for each subplot (Benesty et al., 2009).

Figure 2 Pearson correlation coefficient of NIR_v and kNDVI to TROPOMI SIF. Data are compared at 0.05^∘ spatial resolution and in 16 d time steps starting in April 2018 until March 2021. The value per grid cell in (a) and (b) represents the Pearson correlation coefficient of the vegetation index to SIF in time. Panel (c) represents the difference in correlation of the vegetation indices to SIF.

Figure 2 shows spatial patterns of the Pearson correlation coefficients between both NIR_v and kNDVI to SIF. In both our spatial (Fig. 2) and temporal (Fig. 1) analyses, we find NIR_v is a better predictor for SIF than kNDVI, which contradicts the recent findings from Camps-Valls et al. (2021). However, Camps-Valls et al. (2021) used GOME-2 SIF instead of TROPOMI SIF.

The vegetation index NIR_v outperforms kNDVI in nearly all vegetated regions. Only central Asia, the Sahara, and very high latitudes show a better correlation of kNDVI with SIF. At the same time, these regions generally show a weaker correlation of vegetation indices with SIF.

3 Development and optimization of SIFnet

3.1 Training and optimization of the neural network

Convolutional neural networks (CNNs) are supervised machine learning methods that need matching feature and ground truth data pairs to compute the loss that is back propagated (Bishop, 2007). As such, we begin by coarsening SIF data to 0.5^∘ and use it with auxiliary data at 0.05^∘ as input to SIFnet, allowing us to estimate SIF at 0.05^∘. The model output is compared against the measured TROPOMI SIF at 0.05^∘. After optimizing the model it can resolve a scaling factor of 10 between coarse-resolution input SIF and model output SIF. Figure 3 visualizes this method. In the following step of estimating high-resolution SIF, the feature SIF data have a resolution of 0.05^∘ and auxiliary data of 0.005^∘, resulting in a model output of SIF at 0.005^∘.

Figure 3 CNN model structure and training and estimation method. Yellow and red blocks are convolutional and ReLU layers, respectively. Notation of convolutional layers: k(X1,X2): kernel sizes are X1,X2; chY: number of channels is Y. For training, the data are upscaled. We input auxiliary data at the target resolution and SIF data at a factor of 10 coarser.

Figure 3 shows our chosen CNN model structure for SIFnet. The model consists of convolutional and rectified linear unit (ReLU) layers that are arranged in a sequence. After the first convolutional block there is a residual connection that skips one ReLU and two convolutional layers. Convolutional kernel sizes are either (3,3) or (1,1). This structure is adapted from the literature findings from, for example, Lim et al. (2017). Further, several model structures (Sect. S4.3 in the Supplement) with a different amount of layers, channels, or residual blocks are compared. The chosen model structure represents the best trade-off between complexity and performance. The input feature collinearity and principal component analysis (PCA) presented in Sect. S3 show that some input features have high correlations with each other. A total of 9 out of the 19 PCs in the time variant and 13 out of the 15 PCs in the time-invariant data carry above 99 % of the variance. This suggests that fewer channels should be used in the CNN layers than the feature dimension (because some variables are similar). Therefore the number of channels in the first layer of SIFnet reduces the complexity from 34 to 16 channels (Fig. 3). More complex model structures did not result in a notably improved loss metrics (Sect. S4.3).

For training SIFnet we use 3 years of data (April 2018–March 2021) in 16 d time steps. The study regions are shown in Sect. S2. There are five folds used as training data: two folds over Asia, one over Europe, one over the southern part of Africa, and one over South America. Our validation region is North America (Fig. S4). The hyperparameter tuning is done by training the model on the five folds and computing the loss of the validation data. The parameters are optimized to minimize the loss of the validation data set. Due to computational reasons and the size of the data set, we do not apply a cross validation in the optimization process. The final product consists of high-resolution SIF at 0.005^∘ and is validated against independent SIF measurements of the instruments OCO-2 and OCO-3.

We center and scale each feature individually by subtracting the mean and normalizing by the standard deviation. For data augmentation of the training data, we use random crops and random flips. Each day of one fold has a matrix size of 1200 × 900 pixels. We analyze 69 d in 16 d steps over 3 years. For each input during the training process we randomly crop a matrix with a size of 100 × 100 pixels. As some areas have a large fraction of missing values (e.g., due to water or clouds), we only use cropped matrices that consist of > 80 % valid pixels in the SIF product. Further, we randomly flip vertically and horizontally, both with a probability of 0.5. These data augmentation methods provide us with a huge database that should avoid overfitting the network parameters. During training, all missing values in the data are set to zero. That mainly affects water regions as the share of missing values in the SIF data used is 91.2 % caused by water. In case there is a missing value in the SIF training sample, all feature values of this pixel are also set to zero to ensure the network does not learn false relationships between the predictors and the target variable (that also applies to vegetated regions). For the MODIS bands we applied the quality index value 0 (best quality only). This filtering also removes pixels that include clouds. To ensure a high coverage we interpolated in time for MODIS. Further, training and test folds are selected based on coverage; i.e., the regions near the Equator (between ±22.5^∘) are not included in the data cubes as MODIS reflectance is sensitive to clouds which appear frequently in these regions (compare Figs.  S2 and S4). All static variables have full coverage on land. ERA5 data has full spatial and temporal coverage. We did not apply any further quality filtering on SMAP soil moisture data. The data are provided as a level 3 product on Google Earth Engine.

Our individual loss function is comprised of two loss terms. We use the mean squared error (MSE) loss in combination with the structural dissimilarity index (DSSIM). The DSSIM is the countermeasure of the structural similarity index (SSIM): DSSIM = 1 − SSIM​​​​​​​ (Brunet et al., 2011). Therefore, we are not only optimizing the overall deviation of the estimated SIF to the measured SIF but also the structural patterns. Section S4.4 shows the benefit of including both MSE and SSIM terms in the loss function. Equation (5) shows our loss function:

$$\begin{array}{}\text{(5)}& \begin{array}{rl}\mathcal{L}=& a\cdot \mathrm{MSE}+b\cdot \mathrm{DSSIM}\\ =& a\cdot {\displaystyle \frac{\mathrm{1}}{n}}\cdot \sum _{i=\mathrm{1}}^n{\left(y_i-\stackrel{\mathrm{̃}}{y_i}\right)}^{\mathrm{2}}\\ & +b\cdot \left(\mathrm{1}-{\displaystyle \frac{\mathrm{2}\cdot {\mathit{\mu }}_Y\cdot {\mathit{\mu }}_{\stackrel{\mathrm{̃}}{Y}}\cdot \left(\mathrm{2}\cdot {\mathit{\sigma }}_{Y\stackrel{\mathrm{̃}}{Y}}+c_{\mathrm{2}}\right)}{\left({\mathit{\mu }}_Y^{\mathrm{2}}+{\mathit{\mu }}_{\stackrel{\mathrm{̃}}{Y}}^{\mathrm{2}}+c_{\mathrm{1}}\right)\cdot \left({\mathit{\sigma }}_Y^{\mathrm{2}}\cdot {\mathit{\sigma }}_{\stackrel{\mathrm{̃}}{Y}}^{\mathrm{2}}+c_{\mathrm{2}}\right)}}\right),\end{array}\end{array}$$

where n is the number of data points, y_i is the data point i in measured (target variable) SIF, $\stackrel{\mathrm{̃}}{y_i}$ is the data point i in estimated SIF, Y values are all data points of measured (target variable) SIF, $\stackrel{\mathrm{̃}}{Y}$ values are all data points of estimated SIF, μ_Y is the mean of Y, ${\mathit{\mu }}_{\stackrel{\mathrm{̃}}{Y}}$ is the mean of $\stackrel{\mathrm{̃}}{Y}$, σ_Y is the variance of Y, ${\mathit{\sigma }}_{\stackrel{\mathrm{̃}}{Y}}$ is the variance of $\stackrel{\mathrm{̃}}{Y}$, ${\mathit{\sigma }}_{Y\stackrel{\mathrm{̃}}{Y}}$ is the covariance of Y and $\stackrel{\mathrm{̃}}{Y}$, and c₁=(k₁L)² and c₂=(k₂L)² are variables for stabilization with L= 2 bit px${}^{-\mathrm{1}}-\mathrm{1}$, k₁=0.01, and k₂=0.03. The parameters a and b define the weights of the overall loss of the two individual losses. The overall model performance did not show a notable sensitivity to different a and b values. To approximately keep the individual losses in the same order of magnitude, we set a=1 and b=0.3. DSSIM is in the range of 0 to 1, with 0 meaning structurally similar and 1 structurally dissimilar.

We use the optuna library for the hyperparameter tuning of the learning rate, weight decay, and epoch of the CNN (Akiba et al., 2019). Here, a tree-structured Parzen estimator sampler suggests the parameters of the next trial which is based on a Gaussian mixture model. Section S4.1 provides more details on this hyperparameter tuning.

3.2 Results of model optimization

Figure 4 summarizes the results of the optimized model. We observe an overall r² of 0.92, SSIM of 0.87, and RMSE of 0.17 mW m⁻² sr⁻¹ nm⁻¹ between the estimated SIF from SIFnet and retrieved SIF from TROPOMI at 0.05^∘ (Fig. 4e). SSIM is calculated by comparing the average SIF signal of the 3 years under investigation. Figure 4e shows the three metrics for each month of the year. We observe the lowest r² values in January, February, and March. These are associated with low SIF values and, consequently, lower signal-to-noise ratios which drive the decreased performance. SSIM also indicates reduced performance during this time period. RMSE values are correlated with overall productivity with the lowest RMSE in winter; this is expected as this metric depends on the magnitude of the signal.

Figure 4 Test set results of CNN training at 0.05^∘. Panel (a) shows low-resolution SIF that is used as model input, (b) shows the estimated SIF at 0.05^∘ by SIFnet, (c) shows the measured TROPOMI SIF at 0.05^∘ from Köhler et al. (2018), (d) shows the scatter comparison between TROPOMI SIF and the SIFnet estimate at 0.05^∘, and (e) shows for each investigated month the metrics r², SSIM, and RMSE. Metrics are calculated at 16 d resolution and averaged to monthly values afterwards.

3.3 Which features drive SIFnet?

We are particularly interested in understanding which features drive our neural net. Here we evaluate the feature importance using the permutation feature importance method (Breiman, 2001; Fisher et al., 2019; Gregorutti et al., 2015, 2017) with our North American validation data at a target resolution of 0.05^∘. The method first computes the RMSE including all input features (RMSE_orig.). We then apply the following three steps.

1. Shuffle all pixels of one input feature randomly in time and space.

2. Compute the new RMSE of the estimation (RMSE${}_{\mathrm{F},\mathrm{shuf}.}$).

3. Compare the shuffled RMSE to the original: $d_{\mathrm{F},\mathrm{shuf}.}={\mathrm{RMSE}}_{\mathrm{F},\mathrm{shuf}.}/{\mathrm{RMSE}}_{\mathrm{orig}.}$

Figure 5 Feature importance. Panel (a) shows the total RMSE of the permuted feature divided by the RMSE without feature permutation, and (b) shows the RMSE of each pixel with permuted features divided by the RMSE without feature permutation. Some input variables are clustered, and all variables of that class are permuted at the same time. ρ_MODIS: all seven MODIS bands; LULC: all 11 land cover classes; other VIs: kNDVI, NDVI, and EVI; Mereor.: temperature, precipitation, temperature with 16 d delay, and precipitation with 16 d delay; SM: surface soil moisture and subsurface soil moisture.

Figure 5 shows the feature importance of clustered input classes and individual features to the overall estimation. Multiple applications of the feature permutation yielded negligible differences in feature importance. Figure 5a shows the RMSE share of shuffled data to the RMSE of unshuffled data. SIFnet finds low-resolution SIF (SIF_LR) to be the most important input variable, followed by the vegetation index NIR_v and the cosine of the solar zenith angle (cos (SZA)). All other variables do not contribute notably to the model output. This result strengthens our findings from Figs. 1 and 2 that NIR_v is better correlated with SIF from TROPOMI than kNDVI. Further, our feature importance is in line with Dechant et al. (2022) in which they find a high correlation of SIF with NIR_v multiplied with photosynthetic active radiation (PAR), of which the cos (SZA) can be used as a proxy. The CNN is a data-driven method and is not restricted by LUE terms. Although SM and meteorology (air temperature and precipitation) play a key role for photosynthesis, we find that they are not important to our model output. This does not necessarily imply that SIF is not linked to these parameters. This can be explained by the following. (1) The variables SM and those from ERA-5 are at coarser resolution than the actual model output of the training phase which is at 0.05^∘ (10 000 and 11 132 m for SM and ERA-5, respectively). Therefore each pixel at the resolution of 0.05^∘ does not have its unique value for SM or ERA-5, but multiple cells can be within one SM or ERA-5 pixel. (2) Not only do the auxiliary data of the model estimate higher-resolution SIF, but they are computed together with coarse-resolution SIF. Therefore, events like heat stress that impact a bigger area than the actual model output might be represented in the coarse-resolution SIF. (3) We have aggregated the data used to 16 d time steps. LUE parameters influencing SIF might have a bigger impact on the estimation at higher temporal resolutions.

Figure 5b shows the spatial feature importance over the validation set in North America for the four most important features. We observe that SIF_LR has the biggest impact in the eastern US, which corresponds strongly to the high vs. low productivity regions in the US. NIR_v is a strong predictor in the southeastern US and in shrub regions in the western US. The contribution of cos (SZA) is highest at high latitudes and weakens at lower latitudes. NIR_v is found to be less predictive of SIF at high latitudes. The land mask is the fourth most important input feature and contributes most in shrub regions. These four features consistently stand out as the strongest predictors. Other inputs such as fragmentation and soil moisture were not found to be strong predictors here. In Sect. S4.6 we test higher scaling factors between low- and high-resolution SIF. Even with scaling factors of 20 and 50 low-resolution SIF stays the most and second most important input feature, respectively.

We also examined different combinations of inputs such as directly including the MODIS bands as opposed to vegetation indices derived from MODIS bands (see Fig. S11). Low-resolution SIF remains the most important feature, followed by the NIR band ρ_NIR. The land cover products increase in relevance. Interestingly, when low-resolution SIF is omitted as input for the model, we observe contrasting results to Fig. 5 in which NIR_v is no longer a leading predictor. We find that cos (SZA), ρ_NIR, kNDVI, and NDVI are the four most important features in this case (see Fig. S12). This may result from the collinearity between input features or suggests that another combination of ρ_NIR and ρ_RED is better correlated to SIF than, for example, NIR_v or kNDVI. This finding was robust to multiple optimizations and permutations.

3.4 Comparison of SIFnet to downscaled SIF

Figure 6 shows the 0.005^∘ SIF estimated by SIFnet and downscaled SIF from Turner et al. (2020). The difference between the two SIF estimates can be seen in Fig. 6c. SIFnet predicts lower SIF in the western US drylands and higher SIF over forested regions in the eastern US. This prediction of lower SIF in drylands is interesting because Turner et al. (2020) resorted to an ad hoc bias correction in these regions due to a low signal-to-noise ratio. Recent work from Wang et al. (2022) concluded that SIF and NIR_v capture complementary events in western US drylands as a proxy for GPP and that the linear correlation of SIF to GPP was substantially lower in these regions compared to other vegetation types. We also observe systematic differences in the predicted SIF in urban areas. These regions are further evaluated in Fig. S17. Notably, the SIFnet estimate is systematically lower than the downscaling estimate in most urban regions examined here, Seattle being a notable exception. Both SIFnet and the downscaling approach allocate SIF to large urban parks and green spaces, but SIFnet predicts little-to-no SIF over the rest of the urban area. In particular, SIFnet estimates nearly zero SIF in the urban core of Los Angeles and San Francisco. SIFnet and the downscaling method predict comparable SIF as we move away from the urban core. Fine-scale features in the urban region are visible in both SIF estimates such as the Schiller Woods in Chicago (42.0^∘ N, 87.8^∘ W).

Figure 6 SIFnet estimated SIF at 0.005^∘ for CONUS and its comparison to downscaled SIF. Panel (a) shows the SIFnet estimated SIF at 0.005^∘, (b) shows the downscaled SIF from Turner et al. (2020), and (c) shows the difference between SIFnet and downscaled SIF. Negative values imply a higher SIFnet SIF and positive values a higher downscaled SIF value.

4 Validation against OCO-2/3 SIF

The differences in SIF predicted from SIFnet and the downscaled SIF beg the question: which is correct? Here we evaluate both SIF products against independent SIF observations from OCO-2 and OCO-3 (OCO-2 Science Team/Michael Gunson, 2020; OCO-3 Science Team/Michael Gunson, 2020). These instruments have higher spatial resolution than TROPOMI and, as such, can be used to evaluate the high-resolution patterns predicted by both SIFnet and the downscaling approach. Specifically, OCO-2 and OCO-3 have nadir footprint sizes of 2.25 × 1.29 km. However, OCO-2 and OCO-3 do not provide full spatial coverage. They observe narrow swaths that are ∼ 10 km across-track. OCO-3 also provides a scanning mode to observe urban areas. Here, we use quality-checked OCO-2 data from April 2018 until March 2021 and OCO-3 data from July 2019 until March 2021. To compare the ungridded OCO-2 and OCO-3 data against the SIF estimated from TROPOMI, we compute the weighted average of all 0.005^∘ grid cells that fall within the bounds of an OCO footprint. Here, the TROPOMI estimates are subsampled to 0.0005^∘ (approx. 50 m at Equator), and the mean value is computed for all values which fall into the OCO footprint. For a quantitative comparison between OCO-2/3 and the SIFnet and downscaled estimate, the metrics r, r², and RMSE are computed.

The high-resolution SIF estimates from SIFnet, the downscaling, and OCO are instantaneous SIF measurements taken at a specific time of the day, while the time of TROPOMI observations can differ substantially. Here, we compute the daily average SIF by scaling with the cosine of the SZA (Frankenberg et al., 2011):

$$\begin{array}{}\text{(6)}& \mathrm{Daily}\mathrm{SIF}(x,y)=\mathrm{SIF}({\mathit{\tau }}_{\mathrm{s}},x,y)\cdot {\displaystyle \frac{{\int }_{{\mathit{\tau }}_{\mathrm{o}}}^{{\mathit{\tau }}_{\mathrm{f}}}\cos \left[\mathrm{SZA}\right(\mathit{\tau },x,y\left)\right]\mathrm{d}\mathit{\tau }}{\cos \left[\mathrm{SZA}\right({\mathit{\tau }}_{\mathrm{s}},x,y\left)\right]}},\end{array}$$

where Daily SIF(x,y) is the daily integrated SIF estimate, $\mathrm{SIF}({\mathit{\tau }}_{\mathrm{s}},x,y)$ is the instantaneous SIF at the individual measurement time, SZA is the solar zenith angle, τ_s is the time of the satellite measurement, τ_o is the time of sunrise, and τ_f is the time of sunset. This implicitly assumes that both PAR and SIF scale with cos (SZA) under cloud-free conditions, and we neglect Rayleigh scattering, as well as gas absorption. Although this approach neglects several water or light conditions, it provides our best estimate of daily SIF and enables comparison between multiple SIF products with different measurement times (Turner et al., 2021; Köhler et al., 2018; Frankenberg et al., 2011). The method is equivalent to the daily correction scheme for OCO-2, OCO-3, and TROPOMI (OCO-2 Science Team/Michael Gunson, 2020; OCO-3 Science Team/Michael Gunson, 2020; Frankenberg et al., 2011). Additionally, we performed a sensitivity study in which we trained SIFnet using daily corrected SIF and found the results to be generally insensitive to the use of instantaneous vs. daily corrected SIF (see Fig. S18). Following this, we chose to apply the daily correction after deriving the high-resolution SIF.

Figure 7 Validation of SIFnet and downscaled SIF to OCO-2 and OCO-3 SIF over CONUS. Comparison from April 2018 until March 2021 in 16 d time steps. Daily OCO-2 and OCO-3 data are assigned to the closest 16 d time step. Panel (a) shows the gridded correlation of the two products against the combined data of OCO-2 and OCO-3. We first compute the SIF data from SIFnet and downscaling estimate that falls into the OCO footprint. Then we assign every OCO footprint to the closest grid point on the 1^∘ grid dependent on the center location of that footprint and compute the Pearson correlation coefficient. Panel (b) shows the scatter comparison of the weighted average of all grid cells on the 0.005^∘ estimated SIFnet SIF (ours) and downscaled SIF (Turner et al., 2020) that fall into the OCO-2 or OCO-3 footprint.

Figure 7 shows a comparison of both SIFnet and the downscaled SIF to OCO-2 and OCO-3. Specifically, Fig. 7a shows the correlation of SIFnet and the downscaled estimate with OCO-2/3 for every 1^∘ pixel over CONUS. Both SIFnet and the downscaled SIF generally show good agreement with r in excess of 0.7 for most of the high-productivity regions. We observe weaker correlations in the western drylands due, in part, to a lower signal-to-noise ratio. Overall, we find SIFnet to perform systematically better than the downscaled SIF, as shown in the difference plot. Figure 7b summarizes these spatial patterns in a scatterplot comparison. SIFnet again shows better performance than the downscaled SIF against OCO-2, OCO-3, and OCO-2/3. The Pearson correlation coefficient r is 0.78 and 0.72 for the SIFnet and downscaled estimate, respectively, when comparing to all OCO data (right column in Fig. 7b). The generally high RMSE indicates different scales and variability in the data sets.

Deviations between TROPOMI and OCO-2/3 also appear at a grid of 0.05^∘ (Fig. S19). The r² coefficient is 0.61 and 0.62 between TROPOMI and OCO-2 and OCO-3 SIF, respectively. Indeed, one might expect better correlations here as both present SIF at 740 nm. However, as pointed out in Köhler et al. (2018), the uncertainty of both TROPOMI and OCO-2 SIF is expected to lead to a certain spread between the data sets. In addition, we do not account for differences in acquisition times and viewing–illumination geometry, which can lead to additional uncertainties in this comparison. For reference, when comparing single footprints of TROPOMI SIF to aggregated OCO-2 SIF for June 2018 globally, Köhler et al. (2018) found a r² of 0.67; only additional aggregation leads to a r² of 0.88. The mean deviation of TROPOMI SIF to OCO-2 SIF is close to the average standard deviation of TROPOMI SIF (0.4 mW m⁻² sr⁻¹ nm⁻¹). In our analysis, from the 16 d product from TROPOMI SIF for April 2018 until March 2021 at 0.05^∘, we observe an average error in the TROPOMI SIF of 0.43 mW m⁻² sr⁻¹ nm⁻¹ for the CONUS. That error is close to the RMSE between instantaneous TROPOMI SIF and instantaneous OCO-2 SIF (0.37 mW m⁻² sr⁻¹ nm⁻¹). To compare TROPOMI and OCO-2/3 SIF we aggregate the OCO-2/3 footprints to the same grid as our TROPOMI data (0.05^∘). As we aggregate multiple OCO-2 or OCO-3 footprints to match one TROPOMI grid cell at 0.05^∘, the certainty of the OCO measurements increases, and therefore the RMSE between TROPOMI and OCO SIF decreases.

Figure 8 SIFnet and downscaled SIF, the difference between these, and the difference in correlation to OCO-2 and OCO-3 for four urban regions. The first column shows the SIFnet estimate, the second the downscaled SIF from Turner et al. (2020), the third the difference between SIFnet and downscaled SIF, and the last the difference in correlation of the high-resolution SIF estimates to combined OCO-2 and OCO-3 data on a 0.02^∘ grid multiplied with the ℒ¹ norm between the SIFnet and downscaled estimate.

Figure 8 presents a detailed comparison of SIFnet and the downscaled SIF in four US cities. The first column shows the SIFnet estimate, the second the downscaled SIF from Turner et al. (2020), the third the difference between SIFnet and downscaled SIF, and the last column the difference in correlation of the high-resolution SIF estimates to combined OCO-2 and OCO-3 data on a 0.02^∘ grid multiplied by the ℒ¹ norm between the SIFnet and the downscaled SIF. The column on the right highlights both regions where the differences in predicted SIF are large and which product is performing better. As such, the right column will show white in areas where the difference in predicted SIF is small or the correlation with OCO is similar. While we observe large differences in predicted SIF for the urban areas (column 3), we do not find one product to perform systematically better in urban areas. This likely indicates the complexity in the SIF signal arising from urban areas. Additionally, urban areas make up a small fraction of the overall land mass and, as such, do not represent a large share of the training data in SIFnet. These factors likely contribute to the heterogeneous performance observed in the right column of Fig. 8.

Figure 9 SIFnet SIF, downscaled SIF, MODIS NIR_v, and a Google Earth cut-out for a part of Chicago. Left panel shows the SIFnet estimate, second panel shows the downscaled estimate from Turner et al. (2020), third panel shows MODIS NIR_v for Chicago, and last panel shows the Google Earth cut-out (Google LLC, 2021). For panels 1–3 the average data for April 2018 until March 2021 is shown.

However, there are some notable successes of SIFnet in urban areas that can be mapped directly to features in the urban area. Figure 9 shows both SIFnet and the downscaled SIF along with NIR_v from MODIS and a true color image of Chicago. A feature clearly stands out in both the downscaled SIF and NIR_v image. This is a region with missing NIR_v and effectively no downscaled SIF. However, SIFnet does not show a strong gradient here. This region corresponds to the Chicago airport. In the MODIS NIR_v image it is visible that there are no valid data available for that region for the 3 investigated years. The downscaling method from Turner et al. (2020) relies only on NIR_v in the weighting function. If there are no data available for the region, they are interpolated in space and time. Here, it shows that the method seems to fail in urban regions where no MODIS NIR_v signal is available. SIFnet handles this region better and seems to rely on other auxiliary data if there is no MODIS NIR_v available. In Fig. 8 it is also visible that the SIFnet estimate correlates better with the OCO-X data than the downscaled SIF for the region of the Chicago airport.

5 Conclusions

Here, we develop a convolutional neural network (CNN) model named SIFnet to increase the resolution of TROPOMI SIF by a factor of 10. The novelty of our method consists of using coarse-resolution SIF measurements together with high-resolution auxiliary data as model input to estimate high-resolution SIF. After optimization and hyperparameter tuning of SIFnet, the estimated SIF at 500 m resolution yields an r² and RMSE of 0.92 and 0.17, respectively, when compared against validation data (Fig. 4). We further compare the output of SIFnet against a recently developed downscaling method to estimate high-resolution SIF (Turner et al., 2020) and evaluate both methods against independent observations from the Orbiting Carbon Observatory 2 and 3 (OCO-2/3). SIFnet is found to perform systematically better than the downscaling approach when compared against independent measurements. Through interpretable machine learning methods, we identify the key features that SIFnet utilizes to accurately predict high-resolution spatial patterns of SIF. We find that SIFnet relies heavily on the low-resolution SIF feature (SIF_LR) and the vegetation index NIR_v (Fig. 5).

SIFnet is a multi-layer CNN that increases the spatial resolution of the TROPOMI SIF by a factor of 10. Our model uses auxiliary data sets related to gross primary productivity and SIF as inputs and yields a high-resolution SIF estimate. The model is trained using three years of data from Asia, Europe, Africa, and South America. North America is used as the validation data set. Our loss function is comprised of two terms: the mean squared error and the structural dissimilarity index. The combination of these two terms improved the performance of our model.

SIFnet was further compared to the recent downscaled SIF product developed by Turner et al. (2020). The two high-resolution estimates showed pronounced differences across the western US drylands. This difference is particularly interesting because these drylands tend to be low-productivity regions and traditionally have been difficult for SIF to accurately capture due to the low signal-to-noise ratio. Both high-resolution SIF estimates were compared to independent observations from OCO-2/3. SIFnet performed systematically better than the downscaled SIF (r=0.78 for SIFnet, r=0.72 for downscaling). SIFnet and the downscaling method also yielded differences in urban regions. However, there was substantial heterogeneity in the performance of SIFnet and downscaling in urban areas. One product did not perform systematically better than the other within urban areas. The mixed results in urban areas likely relates to both the complexity of the photosynthetic activity in urban areas, as well as the lack of training data, as urban areas represent a small fraction of the total landmass.

We adapted techniques from the area of interpretable machine learning to assess the key features driving SIFnet. Specifically, we conducted random permutations to input data sets and assessed the impact on the resulting RMSE. From this, we found that SIFnet relies most heavily on the low-resolution SIF feature (SIF_LR). The second most important factor is the MODIS vegetation index NIR_v. NIR_v is also found to outperform the recently proposed kNDVI vegetation index, in contrast to Camps-Valls et al. (2021). The interpretable machine learning approach also allowed us to identify spatial regions of importance for the different parameters. Interestingly, SIFnet relies more heavily on NIR_v in the western drylands where the SIF signal-to-noise ratio is low. This implies that SIFnet is picking up on key physics that lead to the improved performance relative to the downscaling method. Overall, SIFnet represents a robust method to infer continuous high-spatial-resolution information about processes related to gross primary productivity.