## 2 Materials and Methods
### Modeling Multivariate Spatial Dependence between SIF and XCO\({}_{2}\)
Finally, note that since \(\mathbf{s}\) and \(\mathbf{u}\) are locations on Earth's surface, we avoid Euclidean distance and let \(||\mathbf{s}-\mathbf{u}||\) be the chordal distance between \(\mathbf{s}\) and \(\mathbf{u}\) (e.g., [46]). To estimate the parameters of the full bivariate Matern covariance model, we use the standardized residuals, \(\mathbf{R}_{i}^{*}=(\mathbf{R}_{i1}^{*},\ldots,\mathbf{R}_{in_{i}}^{*})^{T}\), \(i=1,2\), given by Equation (6). We first compute empirical semivariograms and cross-semivariograms, and then we fit the parameters of the covariance and cross-covariance functions simultaneously using a multivariate extension of the weighted-least-squares approach developed in [47]. This process is described in Appendix A.1; a check for the validity of the fitted parameters is also established there. Generating the Spatially Contiguous coSIF Data Product with Quantified Uncertainties Based on Cokriging  
From the noisy and relatively sparse spatial data \(\mathbf{Z}_{1}\) (SIF) and \(\mathbf{Z}_{2}\) (XCO\({}_{2}\)), we seek the optimal spatial prediction, \(\widehat{Y}_{1}(\cdot)\), and its root-mean-squared prediction error, \(RMSPE\left(\widehat{Y}_{1}(\cdot)\right)\), for the latent SIF process, \(Y_{1}(\cdot)\), at each of the 0.05-degree land-based grid cell locations \(\{\mathbf{s}_{1},\ldots,\mathbf{s}_{m}\}\), \(m=591,\!119\), in the domain \(D\) spanning North America, shown in Figure 3. In generating our SIF data product, we produce spatial predictions at the 0.05-degree resolution for two main reasons: (1) 0.05 degrees is emerging as the standard for Level 3 SIF data products (e.g., [20; 23; 25; 26; 27]); (2) predictions are at the same resolution as the pre-processed data (Section 2.2). Since the spatial covariance model is formulated at the 0.05-degree resolution, upscaling to coarser resolutions is straightforward (e.g., [48]). The practice of leveraging interdependence between multiple related processes to optimally predict one of those processes is known in geostatistics as cokriging (e.g., [32]). Recall from Equation (2) that we model the spatial process \(Y_{1}(\cdot)\) as the combination of large-scale variation and any remaining spatial variation through the decomposition \(Y_{1}(\mathbf{s})=\mu_{1}(\mathbf{s})+\psi_{1}(\mathbf{s})\), for \(\mathbf{s}\in D\). Having standardized the SIF and XCO\({}_{2}\) data as in Equation (6), our task in cokriging is to optimally predict \(\psi_{1}^{*}(\cdot)=(Y_{1}(\cdot)-\widehat{\mu}_{1}(\cdot)-\hat{\mu}_{R,1})/ \hat{S}_{R,1}\), which is the standardized version of \(\psi_{1}(\cdot)\) defined in Equation (3). From predictions of \(\psi_{1}^{*}(\cdot)\), we can obtain predictions of \(Y_{1}(\cdot)\) on the original scale by reversing the standardization with a simple linear transformation. As in Section 2.3, information about the behavior of \(\psi_{1}^{*}(\cdot)\) (for SIF) and \(\psi_{2}^{*}(\cdot)\) (for XCO\({}_{2}\)) is taken from the standardized spatial residuals \(\mathbf{R}_{1}^{*}\) and \(\mathbf{R}_{2}^{*}\) described in Equation (6). For a given set of months, \(\mathbf{R}_{1}^{*}\) and \(\mathbf{R}_{2}^{*}\) have a combined size on the order of 50,000 elements (see Figure 2), so that traditional cokriging is computationally intractable due to a required matrix inversion that is shown in Equation (11) below. However, spatial prediction is a relatively local operation ([40], Section 3.4), so for a given prediction location \(\mathbf{s}_{0}\), we use a local neighborhood defined by the 150 values whose corresponding locations are closest in chordal distance to \(\mathbf{s}_{0}\), respectively, in \(\mathbf{R}_{1}^{*}\) and \(\mathbf{R}_{2}^{*}\). Let \(\mathbf{R}_{01}^{*}\) and \(\mathbf{R}_{02}^{*}\) be these 150-dimensional vectors of local residuals, and define \(\mathbf{R}_{0}^{*}=(\mathbf{R}_{01}^{*\tau},\,\mathbf{R}_{02}^{*\tau})^{T}\).