Interactive comment on “ Mapping of satellite Earth observations using moving window block kriging ”

Abstract. Global gridded maps (a.k.a. Level 3 products) of Earth system properties observed by satellites are central to understanding the spatiotemporal variability of these properties. They also typically serve either as inputs into biogeochemical models or as independent data for evaluating such models. Spatial binning is a common method for generating contiguous maps, but this approach results in a loss of information, especially when the measurement noise is low relative to the degree of spatiotemporal variability. Such "binned" fields typically also lack a quantitative measure of uncertainty. Geostatistical mapping has previously been shown to make higher spatiotemporal resolution maps possible, and also provides a measure uncertainty associated with the gridded products. This study proposes a flexible moving window block kriging method that can be used as a tool for creating high spatiotemporal resolution maps from satellite data. It relies only on the assumption that the observed physical quantity exhibits spatial correlation that can be inferred from the observations. The method has several innovations relative to previously applied methods: (1) it provides flexibility in the spatial resolution of the contiguous maps, (2) it is applicable for physical quantities with varying spatiotemporal coverage (i.e., density of measurements) by utilizing a more general and versatile data sampling approach, and (3) it provides rigorous assessments of the uncertainty associated with the gridded products. The method is demonstrated by creating Level 3 products from observations of column-integrated carbon dioxide (XCO2) from the GOSAT (Greenhouse Gases Observing Satellite) satellite, and solar induced fluorescence (SIF) from the GOME-2 (Global Ozone Monitoring Experiment-2) instrument.


Introduction
Satellite measurements of Earth's surface and atmospheric quantities have enormous benefits for Earth system science due to their global coverage and near-real-time availability.They provide key constraints for developing models representing our understanding of the functioning of the Earth system.However, due to orbit geometries and geophysical limitations, a uniform or contiguous global coverage of these observations in space and/or time is not possible.This necessitates creation of contiguous maps for obtaining measurements at unsampled times and locations for understanding overall patterns, driving biogeochemical or physical models, and/or validating model predictions.Due to their widespread utility, global gridded maps are often part of the standard suite of satellite data products and are often termed "Level 3" data (e.g., NASA, 2014).
In the case of column-integrated carbon dioxide (XCO 2 ) and solar induced fluorescence (SIF) observations, the two illustrative applications that will be used in this work, gridded products have been used, for example, to evaluate the representation of water stress in models of photosynthesis (Lee et al., 2013), to assess the performance of a terrestrial biosphere model in representing global CO 2 distributions (Hammerling et al., 2012b), and to constrain a model to assess the relative roles of variations in atmospheric transport and carbon exchange in explaining atmospheric CO 2 variability over the Amazon (Parazoo et al., 2013).The generation of Level 3 products is also often part of the standard processing sequence of observations (e.g., GOSAT Project, 2014;CO 2 DAAD, 2014).
Presently, "binning" is the most widespread method for creating such contiguous maps of satellite data.Such binning typically involves computing the mean of the observa-tions that fall within a grid cell (a.k.a."bin") of an appropriate geographic size and time window (for applications of binning in the context of satellite retrievals of atmospheric concentration of carbon dioxide see Kulawik et al., 2010, andCrevoisier et al., 2009).However, this simplicity comes with some limitations such as (1) the mean being computed from a different number of measurements across grid-cells, (2) the inability to take into account any redundancy among nearby observations in computing the mean, and (3) the lack of gap filling properties for grid cells that may contain no observations for a given time window.
The methodological deficiencies of binning can be overcome by using kriging, a geostatistical interpolation approach that takes into account the spatial and/or temporal correlation in the observations.Kriging is a best linear unbiased estimator, with the various implementations of ordinary kriging relying on the assumption of intrinsic stationarity.More typically, a covariance function is used to represent spatial correlation, and second-order stationarity is assumed, i.e., that the mean is constant and the covariance is only a function of the distance between observations (for kriging see Chiles and Delfiner, 2012).Because the mean and covariance of Earth system observations vary substantially, the kriging tools need to be modified to reflect this nonstationarity.One such method is moving window kriging, in which kriging is performed locally and the covariance parameters are determined locally within pre-specified spatial and/or temporal subdomains (e.g., Haas, 1990).The ability of the moving window kriging to reflect local uncertainty has been emphasized to be the most important advantage over kriging methods relying on the global covariance models (e.g., Harris et al., 2010;Walter et al., 2001;Van Tooren and Haas, 1993).Due to this advantage, the moving window kriging has been previously used for creating contiguous maps of satellite remote sensing observations of column-averaged CO 2 (XCO 2 ) (e.g., Hammerling et al., 2012a and b).
This work proposes a further development of the moving window kriging method for application with satellite observations of Earth system properties.Whereas Hammerling et al. (2012a, b) used ordinary kriging as the basis for obtaining estimates at the spatial support (i.e., resolution or spatial footprint size) of observations, we propose a moving window block kriging method that can yield estimates at any resolution equal to or greater than that of the observations (for discussions on change of support in the context of remote sensing see Atkinson and Curran, 1995, Collins and Woodcock, 1999, and Braverman, 2011).The main advantages of the proposed tools are that they make it possible to (1) select the spatial support/resolution of the mapped quantities, (2) handle large volumes of data by developing subsampling techniques that can make moving window block kriging computationally feasible for a large number of satellite measurements, and (3) provide rigorous assessments of the uncertainty associated with the contiguous maps.

Methods
The proposed approach builds on the work of Hammerling et al. (2012a, b), with the goal of increasing the applicability and the flexibility of the nonstationary local kriging approach presented therein.The main innovations are twofold.The first is to allow for flexibility in the spatial support of the estimates (i.e., the spatial resolution at which the mapping is conducted).The second is to provide a general approach for subsampling available observations in a manner that (i) captures the local correlation structure in the vicinity of each estimation grid cell and (ii) makes the statistical mapping approach computationally feasible in the case of applications with a very large number of observations.The mapping proceeds in three steps for each grid cell and each estimation time on a regular grid, in order to create a contiguous map of the satellite observations.These steps are outlined in the subsections below and include subsampling of the observations, characterization of the local spatial covariance structure, and interpolation at the desired spatial resolution.In Sect.3, the new mapping approach is applied to two prototypical examples of satellite observations, namely observations of column-integrated concentration of atmospheric CO 2 concentrations (XCO 2 ) and observations of surface solar induced fluorescence (SIF), measured by the GOSAT (Greenhouse Gases Observing Satellite) satellite, and by the GOME (Global Ozone Monitoring Experiment) instrument, respectively.

Subsampling of observations
The goal of the subsampling strategy is to preferentially sample observations in the vicinity of a given estimation grid cell, such that both the characterization of the local spatial covariance structure and the ultimate mapped estimate and its associated uncertainty are representative of local variability.This is accomplished by selecting the total number of observations to be used, N, where N is selected to be large enough to yield a representative sample but small enough to make mapping computationally feasible on a given computational platform.For the applications presented in Sect.3, N = 500 and N = 1000 for the XCO 2 and SIF mapping, respectively.
N observations are selected for each estimation grid cell by assigning a relative selection probability to each observation based on that observation's separation distance from the centroid of the grid cell.This selection probability could be application specific, but for the applications presented here we selected where P s is the relative probability of a given observation being selected, and h is the great circle distance between the location x i of an observation and the centroid x j of the estimation grid cell: Geosci.Model Dev., 8, 3311-3319, 2015 www.geosci-model-dev.net/8/3311/2015/where r is the radius of the Earth and ϕ i and λ i are the latitude and longitude of location x i .
The form of P s in Eq. ( 1) ensures that a comparable number of observations is selected within any equal-area concentric band around an estimation grid cell, thereby also ensuring that observations that are at close distances to one another are preferentially close to the estimation location.This is a desirable feature because observations that are close to one another define the shape of the variogram at short separation distances (Sect.2.2), and the variogram should reflect variability in the vicinity of the estimation grid cell.Different forms of P s could also be used, for example if more/fewer observations along a given direction were desirable in order to better represent expected correlations along a given direction.
In previous work (Alkhaled et al., 2008;Hammerling et al., 2012a, b), a fixed application-specific window size was instead defined within which all available observations were used, together with a user-defined fraction of observations outside of the window.The window size was based in part on expected scales of variability in the satellite observations.The updated approach presented here reduces the number of user-selected parameters and explicitly provides a mechanism for ensuring the computational feasibility of mapping in the case of very large data sets, such as the SIF example examined here.

Characterization of spatial covariance
The characterization of the local covariance structure of the observations around each estimation grid cell, based on the subsampled observations, proceeds as described in Hammerling et al. (2012a, Sect. 2.1), except that (1) all possible pairs of observations are included in the formulation of the raw variogram and the nugget-effect variance, representative of the retrieval/measurement errors, is not spatially uniform.The reader is referred to that earlier publication for additional details.
Briefly, for each estimation grid cell, a raw variogram is calculated based on the subsampled observations: where γ is the raw variogram value for a given pair of observations y(x i ) and y(x j ), and h is the great circle distance between the locations (x i and x j ) of these observations, as defined in Eq. (2).A parametric function, the theoretical variogram, is fitted to the raw variogram using non-linear least squares.For the prototypical applications presented here, an exponential variogram function with a nugget effect was used, because it yields a valid covariance function on a sphere (Huang et al., 2011), provided a good match to the known physical characteristics of the observations, and fit the observed variability well: where σ 2 and l are the variance and correlation length of the quantity being mapped, and σ 2 nug is the nugget variance, typically representative of measurement and retrieval errors in the case of satellite observations.The nugget component can be either prescribed (as in the XCO 2 example in Sect.3) or estimated (as in the SIF example in Sect.3), depending on the availability of information about measurement and retrieval errors.
The variogram parameters can be used to define a corresponding local spatial covariance structure for the mapped quantity (XCO 2 or SIF, in the prototypical examples presented here).For the variogram function in Eq. ( 4) this becomes The nugget effect is correspondingly used to define the covariance structure of the measurement and retrieval errors:

Mapping using moving window block kriging
Ordinary kriging, a minimum variance linear unbiased mapping method for spatial data, was used in Hammerling et al. (2012a, b) to create contiguous maps of XCO 2 .In this approach, the spatial support (i.e., footprint) of the estimates corresponds to that of the observations.Although the mapping can be performed at any spatial interval (e.g., once per 1 • × 1 • grid cell), the estimates remain representative of the variability at the scale of the observations.
Here, we instead use block kriging (e.g., Webster, 2000), an approach that yields estimates that represent an average within a specified area.This makes it possible to disassociate the native footprint of the observations from the resolution of the mapped product, thereby making it possible to create contiguous maps at any desired spatial resolution equivalent to or greater than the size of the observation footprints.As with moving window ordinary kriging, block kriging provides an optimal estimate of the quantity being mapped (XCO 2 and SIF, in the prototypical examples presented here) for each estimation location, based on the subsampled observations (Sect.2.1) and the local covariance structure (Sect.2.2), together with a rigorous assessment of the uncertainty associated with the estimate.
The linear system of equations that is solved to obtain the N weights λ assigned to the subsampled observations for a given estimation grid cell is where Q is a N × N covariance matrix among the N observations with individual entries as defined in Eq. ( 5), R is a N × N diagonal measurement and retrieval error covariance matrix among the N observations as defined in Eq. ( 6), 1 is a N× 1 unity vector, T denotes the vector transpose operation, and q A is an N× 1 vector of the spatial covariances between the estimation grid cell and the N observation locations, defined as where q A,i is the covariance between the grid cell and observation i, and q h i,j is defined as in Eq. ( 5) based on the distance h i,j between observation i and n regularly spaced locations within the grid cell.In general, the larger the n the better the representation of the area (i.e., grid cell) to observation covariance.For practical purposes, in the applications presented here, n is defined based on the relative footprint of the observations compared to that of the estimation grid cells.
The system in Eq. ( 7) is solved for λ and the Lagrange multiplier ν.These parameters are then used to define the estimate (ẑ) and estimation uncertainty variance (σ 2 ẑ ) for the grid cell as where y is the N× 1 vector of subsampled observations, and σ AA is the variance of the mapped quantity (XCO 2 or SIF, in the prototypical examples presented here) at the resolution of the estimation grid cell, defined as where q h j,k is defined as in Eq. ( 5) based on the distance h j,k between any combination of the n regularly spaced locations within the grid cell defined previously.

Example applications
The mapping approach described in Sect. 2 is demonstrated using two prototypical examples of satellite observations: (1) observations of column-integrated concentrations of atmospheric CO 2 (XCO 2 ) from the GOSAT satellite, and (2) observations of surface solar-induced fluorescence (SIF) from the GOME-2 instrument.These applications differ in the spatial footprint (i.e., support) of the observations (nadir footprint of about 10.5 km diameter at sea level (Kuze et al., 2009) and 40 km × 80 km (Joiner et al., 2013), respectively), the volume of available data (approximately 2 × 10 3 and 2 × 10 5 observations per week, respectively), the timescales of variability, and the degree of spatial variability and nonstationary in the observed quantity.

Global land XCO 2 fields observed by GOSAT
The Japanese GOSAT (e.g., Kuze et al., 2009) was launched in 2009 and is the first satellite dedicated to global greenhouse gas monitoring, including CO 2 and CH 4 .GOSAT flies in a polar, sun-synchronous orbit with a 3-day repeat cycle and an approximately 13:00 LT overpass time.GOSAT XCO 2 data are being used to examine a number of questions in carbon cycle science, including comparing observed and modeled XCO 2 fields (Hammerling et al., 2012b), quantifying sources and sinks of CO 2 (e.g., Deng et al., 2014;Basu et al., 2013Basu et al., , 2014;;Chevallier et al., 2014;Takagi et al., 2014), detecting perturbations in the carbon cycle (Guerlet et al., 2013), and interpreting seasonal changes in the carbon balance (Parazoo et al., 2013).Measurements of XCO 2 (a.k.a."Level 2" data) are derived using a number of retrieval algorithms, among them NASA's Atmospheric CO 2 Observations from Space (ACOS) algorithm (e.g., O'Dell et al., 2012;Crisp et al., 2012).Filtered and bias-corrected data from the most up to date version of this algorithm (ACOS v3.4 release 3) are used here to demonstrate the mapping approach presented in Sect. 2. Approximately 900 successful retrievals are available per 3-day repeat cycle, with the majority of observations being over land.These data have substantial retrieval uncertainties (e.g., O'Dell et al., 2012) and include large gaps (e.g., Fig. 1).These features prevent the application of simple spatial and temporal binning techniques for generating XCO 2 maps at spatiotemporal scales that are directly useful for addressing existing uncertainties in carbon cycle science.
The approach described in Sect. 2 is used to create continuous maps, a.k.a.Level 3 data, based on XCO 2 observations obtained over two repeat cycles, namely 2-7 August 2009 (Fig. 1).A 6-day period is used to balance the competing goals of including as many observations as possible, while avoiding time periods over which the XCO 2 field itself would change substantially (see discussion in Hammerling et al., 2012a).Maps of XCO 2 and associated uncertainties are created at native (Fig. 2a, b) and 1 • × 1 • (Fig. 2c, d) resolutions, in order to examine and demonstrate the impact of  resolution on mapping uncertainty.Targeting different resolutions is made possible by the use of the moving window block kriging approach presented here.N = 500 subsampled observations are used per estimation location.These maps can also be compared to those presented for an equivalent period in Hammerling (2012b, Auxiliary Figs. 2 and 3), with methodological differences as described in Sect. 2 and representative of the estimated XCO 2 at the native resolution of sounding footprints (nadir footprint ∼ 10.5 km diameter) with estimates at 1 • × 1.25 • intervals.Results show that, because of the information content of the sparse observations, the estimated fields (Fig. 2a, c) are similar at native and 1 • × 1 • resolutions; however, estimating directly at the coarser 1 • × 1 • resolution yields lower uncertainties as observations become more informative for spatially averaged quantities (Fig. 3).

Global land solar-induced fluorescence fields observed by GOME-2
A series of recent studies has demonstrated the potential use of satellite observations of SIF for understanding and quantifying photosynthetic CO 2 uptake at large scales, using data from the GOSAT satellite (e.g., Joiner et al., 2011;Frankenberg et al., 2011;Guanter et al., 2012;Joiner et al., 2012;Lee et al., 2013;Frankenberg et al., 2012), the SCIAMACHY (SCanning Imaging Absorption spectroMeter for Atmospheric CHartographY) instrument on board Envisat (e.g., Joiner et al., 2012), the GOME-2 instrument on board MetOp-A (Meteorological operational satellite-A; e.g., Joiner et al., 2013), and the Orbiting Carbon Observatory (OCO-2) (e.g., Frankenberg et al., 2014).Satellite measurements of fluorescence can be used with land surface models to improve the representation of GPP (gross primary production) and to understand the GPP response to environmental stress (e.g., Lee et al., 2013).Among available data sets, GOME-2 provides the highest spatial and temporal density of data.Until now, studies of SIF have relied on spatially and temporally binned average observations at monthly or coarser timescales and 1 • or coarser spatial scales (e.g., Fig. 4).The coarse spatial and temporal scales were used to overcome, through the use of simple averaging, spatial gaps in observations and the relatively high uncertainties associated with individual retrievals.One of the limitations of such an approach is that it inherently discards information about SIF variability at fine spatial and temporal scales, which is important for understanding the impact of transient effects such as changes in phenology and water availability (Lee et al., 2013), and developing biospheric models that can represent these effects correctly.A second limitation is the lack of a direct and robust quantification of the uncertainty associated with the mapped products, complicating uncertainty analysis in subsequent applications using the data.
As a second demonstration of the mapping approach proposed here, we use SIF GOME-2 V.14 data (Joiner et al., 2013) with the approach described in Sect. 2 to create contiguous maps of SIF at a single spatial resolution (1 • × 1 • ) but at multiple temporal resolutions.The examination of shorter time periods was selected in order to more directly respond to scientific opportunities in the use of SIF data and to complement the spatial-resolution-focused demonstration of Sect.3.1.Maps of SIF and associated uncertainties are created at 1-, 6-, and 31-day temporal resolutions in August, 2009 (Fig. 5), where August 2009 was chosen for convenience to correspond with the XCO 2 application presented in Sect.3.1.N = 1000 subsampled observations are used per estimation location.The monthly map can also be compared to the monthly binned map presented in Fig. 4.
Results show that the proposed approach can leverage nearby observations to create realistic contiguous maps even at 1-day resolution (Fig. 5a, b), although, as expected, uncertainties are reduced (Fig. 5d) at coarse temporal resolutions, just as was seen for coarser spatial resolutions in the XCO 2 application.In fact, the regions with little or no SIF data for the 1-day application are clearly visible as highuncertainty bands in Fig. 5b, and a user could explicitly decide whether such uncertainties are acceptable or too high for a given scientific application.When maps are intended to be used to drive and/or validate biogeochemical models, having the ability to choose a desirable balance between temporal resolution and mapping, uncertainty presents a considerable advantage.
Ideally, the temporal resolution at which maps are obtained is as fine as possible so as to capture the dynamics of the observed physical quantity, in this case SIF.The choice of optimal temporal resolution thus, in general, defines a trade-off between having sufficient observations for adequate spatial coverage while minimizing the impact of temporal variability in the quantity being examined (Hammerling et al., 2012a).In Fig. 6 it is apparent that the presented approach makes it possible to obtain maps at temporal resolutions much higher than the monthly (or coarser) res- olution of current binned products.As expected, the more abundant observations available at 6-day temporal resolution (Fig. 6d) lead to decreased estimation uncertainty compared to those at 1-day resolution (Fig. 6b).However, at monthly temporal resolutions (Fig. 6e, f) the temporal variability in SIF over a 31-day period increases the discrepancy among (spatially) nearby observations, leading to increased uncertainties at coarse timescales.This effect is apparent in comparing Fig. 6d and f, as uncertainty increases over, for example, eastern South America.A similar trade-off was also noted in selecting mapping timescales for XCO 2 (Hammerling et al., 2012a) and further speaks to the advantage of being able to select a mapping timescale based on scientific need and uncertainty tolerance, as is possible with the approach presented here.ping method is tested in terms of the accuracy of the best estimates (the difference between estimates and withheld observations) and the accuracy of the uncertainty bounds (the degree to which the reported uncertainties capture the difference between estimates and withheld observations) and bias (the mean difference between estimates and withheld observations).
The accuracy of the maps at daily temporal resolution and native spatial resolution is evaluated using the mean absolute difference (MAD) and the root mean squared difference (RMSD) between the mapped estimates and observations held back in leave-one-out cross-validation (Table 1).Although an absolute target value for these accuracy metrics is not available, it is interesting to note that the MAD and RMSD are comparable to the reported measurement uncertainty in both satellite data sets (0.77 ppm for GOSAT XCO 2 , 0.55 mW m −2 sr −1 nm −1 for GOME-2 SIF).We also compare the GOSAT XCO 2 values to those obtained from by applying the method developed in Hammerling et al. (2012a), which yielded a MAD of 0.86 ppm and a RMSD of 1.20 ppm, demonstrating comparable performance, but with the additional benefits provided by the new method as described in Sect. 2.
Estimation uncertainties reflect the locations and number of observations surrounding the estimation location, the degree of spatial variability in the mapped field in the vicinity of the estimation location, and the spatiotemporal support of the estimates.The accuracy of the uncertainties obtained from the mapping method is evaluated by quantifying the reliability with which the uncertainty bounds associated with the estimates capture the values of the withheld observations.Specifically, we calculate the percentage of estimation locations where the withheld observations fall outside of the 1, 2, and 3 estimation standard deviation (σz ) uncertainty bounds.For independent, normally distributed data, these percentages should be approximately 32, 5 and 0.3 %, respectively.Although these assumptions do not hold here, these values still provide a general indication of expected performance.
For both applications, the percentage of observations falling outside of the uncertainty bounds is lower than would be expected for normally distributed data (Table 1), showing good mapping accuracy.These percentages are very similar when the analysis is repeated using the method developed by Hammerling et al. (2012a).The lower percentages are due to the fact that observations are not normally distributed.
Finally, the bias of the developed method is quantified using the mean difference between estimates and the withheld observations in the leave-one-out cross-validation.Theoretically, mean difference should approach zero as the number of cross-validation points increases if the method provides perfectly unbiased estimates.The mean difference for both applications (Table 1) was several orders of magnitude lower than the observed spatial gradients in the mapped quantities (e.g., Figs. 1 and 4) and was not statistically significant (p > > 0.05: p = 0.86 for GOSAT XCO 2 ; p = 0.63 for GOME-2 SIF).The approach therefore yields unbiased estimates.

Conclusions
In this study we propose a flexible moving window block kriging method that can be used as a tool for creating high spatiotemporal resolution maps from satellite data.The method can be applied in a stand-alone mode, or as a part of broader satellite data processing package.The resulting maps can also be incorporated into biogeochemical and physical models of the Earth system.The approach relies only on the assumption that the observed physical quantity exhibits spatial correlation that can be inferred from the observations.The method has several advantages over previously applied methods: (1) it allows for the creation of contiguous maps at varying spatiotemporal resolution, (2) it can be applied for creating contiguous maps for physical quantities with varying spatiotemporal coverage (a.k.a.density of measurements), and (3) it provides assessments of the uncertainty of interpolated values.The approach emphasizes the use of local covariance structures in predictions by an arbitrary selection of the sampling function, limiting the radius around estimation locations and adjusting the number of sampled points to a fraction of available measurements.The approach also limits the number of partially subjective ancillary parameters required, making it applicable across a variety of applications.
The method was demonstrated by creating Level 3 products from two data sets with considerably different spatiotemporal properties.While the GOSAT XCO 2 observations were relatively sparse, the GOME-2 SIF data had a much higher spatiotemporal density.In the case of GOSAT XCO 2 , the effects of making predictions at different spatial supports (i.e., resolutions) were analyzed, showing that a decrease in the resolution slightly affects estimates ("smoothing" effect) and more significantly estimation uncertainties (reduced uncertainties at coarser resolution).In the case of GOME-2 SIF, the focus was kept on the effect of different aggregation time periods by creating maps at higher temporal resolutions.This example demonstrated the importance of being able to select a mapping timescale based on scientific need and uncertainty tolerance as optimal temporal resolution results from a trade-off between having sufficient observations for adequate spatial coverage, while minimizing the impact of temporal variability in the quantity being examined.In this it was shown that even daily Level 3 maps could be successfully created by the proposed method.For both data sets, the method was shown to yield precise, accurate, and unbiased estimates.The results clearly indicate that contiguous maps can be created at different spatial resolutions for time periods shorter than achievable by binning/averaging.
The resulting maps can be used to support the development of improved models of the Earth system, both by serving as driver data and validation data for such models.

Table 1 .
Cross-validation results of GOSAT XCO 2 and GOME-2 SIF data sets, including mean absolute difference, root mean squared difference, percent of observations lying outside of 1, 2, and 3 standard deviations (σ ẑ) of the mapping uncertainty and mean difference.