Nonstationary Spatial Modeling of Massive Global Satellite Data

Huang Huang111 Statistics Program, King Abdullah University of Science and
Technology, Thuwal, 23955-6900, Saudi Arabia<EMAIL_ADDRESS>, Lewis
R. Blake222 Department of Applied Mathematics and Statistics, Colorado School
of Mines, Golden, CO, 80401, USA<EMAIL_ADDRESS><EMAIL_ADDRESS>,
Matthias Katzfuss333 Department of Statistics, Texas A&M University, College
Station, TX, 77843, USA<EMAIL_ADDRESS>, and Dorit M. Hammerling2

###### Abstract

Earth-observing satellite instruments obtain a massive number of observations
every day. For example, tens of millions of sea surface temperature (SST)
observations on a global scale are collected daily by the Moderate Resolution
Imaging Spectroradiometer (MODIS) instrument. Despite their size, such
datasets are incomplete and noisy, necessitating spatial statistical inference
to obtain complete, high-resolution fields with quantified uncertainties. Such
inference is challenging due to the high computational cost, the nonstationary
behavior of environmental processes on a global scale, and land barriers
affecting the dependence of SST. In this work, we develop a multi-resolution
approximation ($M$-RA ) of a Gaussian process (GP) whose nonstationary,
global covariance function is obtained using local fits. The $M$-RA requires
domain partitioning, which can be set up application-specifically. In the SST
case, we partition the domain purposefully to account for and weaken
dependence across land barriers. Our $M$-RA implementation is tailored to
distributed-memory computation in high-performance-computing environments. We
analyze a MODIS SST dataset consisting of more than 43 million observations,
to our knowledge the largest dataset ever analyzed using a probabilistic GP
model. We show that our nonstationary model based on local fits provides
substantially improved predictive performance relative to a stationary
approach.

Keywords: distributed computing, multi-resolution approximation, nonstationary
covariance, probabilistic prediction, sea surface temperature, Gaussian
process

##  1 Introduction

NASA’s large fleet of earth-observing satellites obtain a massive number of
observations each day. After pre-processing and retrieval, this results in
billions of geophysically calibrated and georeferenced swath measurements on a
global scale, referred to as Level-2 data by NASA. Measurements of a
particular variable are typically noisy and incomplete, leaving big gaps in
daily Level-2 maps, which prohibit the optimal use of these important data
sources and, in turn, prevent many scientific problems from being resolved.
Hence, the value and applicability of these data to answer important science
questions can be greatly improved if the true, complete spatial fields are
inferred optimally at high resolution, and uncertainties are rigorously
quantified.

Here we focus on sea-surface temperature (SST). As the oceans cover more than
70% of Earth’s surface and are critical for Earth’s life and ecosystems, SST
is essential for hurricane prediction, weather forecasting, and climate-change
research (O’Carroll et al.,, 2019). The Moderate Resolution Imaging
Spectroradiometer (MODIS) on the NASA Aqua satellite provides tens of millions
of Level-2 SST observations per day at a very high $1\textrm{km}^{2}$ spatial
resolution on a global scale. These observations exemplify the statistical
challenges associated with analyzing Level-2 satellite data: The datasets are
massive, global, and the underlying fields nonstationary. In addition, land
masses provide barriers, which can reduce spatial dependence at given
distances and thus provide further sources of nonstationarity.

Current practices for turning Level-2 satellite data into complete, gridded
data products, referred to as Level-3 data, are often relatively simple due to
convenience, computational limitations, or the need for very general products
that are useful to a broad user base. In many cases, including NASA’s official
Level-3 data products, data are gridded by simply computing averages over
coarse spatial and temporal grids. This approach is fast and intuitive but
also leads to a loss of fine-scale features and opaque uncertainties for
further use of the data, making it difficult or impossible to answer some
scientific questions.

Spatial statistical approaches provide a formal toolbox for analyzing spatial
data, enabling more accurate spatial predictions and rigorous assessment of
uncertainties. Based on a statistical model parameterizing the spatial
dependence structure, estimation of parameters and spatial prediction
(“kriging”) can be carried out (e.g., Cressie,, 1993). So-called spatial and
spatio-temporal hierarchical models (e.g., Cressie and Wikle,, 2011) are
essentially flexible extensions of kriging that can handle a variety of
different types of observations. However, most of these models are
computationally infeasible for very large satellite datasets. For oceanic data
analyses, optimal interpolation (Bretherton et al.,, 1976) is still the most
common algorithm for blending, say SST and data obtained by different
instruments (Chao et al.,, 2009). Optimal interpolation is essentially a
special case of kriging, but it has several drawbacks in its typical form. For
example, it is computationally infeasible for large datasets and lacks
suitable models for spatial nonstationarity, the latter of which is important
for many environmental variables, including SST, known to be smoother and less
variable offshore than nearshore (e.g., Chao et al.,, 2009).

Recent endeavors to scale spatial statistical approaches to large spatial data
were reviewed by Heaton et al., (2019), which included a numerical comparison
on a satellite dataset with about 0.15 million observations. Shirota et al.,
(2019) proposed a spatial process model as the sum of a low-rank Gaussian
predictive process (Banerjee et al.,, 2008) and a nearest-neighbor Gaussian
process (Datta et al.,, 2016) approximation of a stationary covariance using a
shared-memory parallelization scheme (OpenMP); their approach was applied to
remote-sensing data sets of around 17 million observations. Finley et al.,
(2020) implemented an R (R Core Team,, 2020) package for nearest-neighbor
Gaussian process approximations of stationary covariances with OpenMP and
showed an application example of around 38 million observations. Zammit-
Mangion and Rougier, (2020) considered nonstationarity modeling via a multi-
scale random process. They projected the process to a triangulation grid and
induced conditional independence so that the precision is sparse and inference
is feasible for large datasets. In the application, they sampled one million
observations from global SST satellite data and performed the inference using
shared-memory parallelism. Appel and Pebesma, (2020) proposed the
nonstationary $M$-RA framework in the spatio-temporal setting and applied
their R implementation to moderately large spatio-temporal SST datasets with
around 0.6 million observations.

Of the many SST products based on satellite data provided by the Group for
High Resolution Sea Surface Temperature (GHRSST), the multi-scale ultra-high-
resolution (MUR) SST (Chin et al.,, 2013) is possibly the most advanced
product currently available. It employs multi-resolution variational analysis
(MRVA), which uses wavelet-based multi-scale signal expansion to address the
irregularity in measurement locations and scale-dependent issues (Chin et
al.,, 1998). However, the MRVA uses basis functions of a fixed shape, which
for SST can result in over-smoothing in coastal areas and under-smoothing off-
shore. Also, additional approximations are needed to accommodate the
irregularly spaced data and the uncertainty estimates, which are only point-
wise.

Here, we provide an efficient method for analyzing global Level-2 satellite
data based on the multi-resolution approximation ($M$-RA) for spatial
processes first introduced in Katzfuss, (2017) and Katzfuss and Gong, (2020),
which uses a large number of basis functions to capture spatial variation at
all scales while being highly computationally efficient. In contrast to
related multi-resolution wavelet models (e.g., the MRVA), the $M$-RA is
directly applicable to irregularly spaced observations, allows proper
probabilistic inference, and the basis functions can adjust flexibly to
spatially-varying dependence structure.

In this work, we extend this basic $M$-RA approach in several ways. We combine
the $M$-RA with a nonstationary, global covariance function; specifically, we
employ a nonstationary Matérn covariance in $\mathbb{R}^{3}$ with spatially
varying parameters (Paciorek and Schervish,, 2006) on the 2-sphere embedded in
$\mathbb{R}^{3}$, with parameters varying according to local fits of isotropic
Matérn covariances. We also develop an ocean-specific scheme, where we
partition boundaries and place the basis-function for the $M$-RA in such a way
that we can specifically account for, and weaken, dependence across land
barriers emulating the dependence structure of SST across different ocean
basins. On the computational side, we provide and describe an implementation
of the $M$-RA tailored to distributed-memory computation in high-performance-
computing environments that can split the computing task efficiently between
many computational nodes; further details of this implementation can be found
in several technical notes (Huang et al.,, 2019; Blake et al.,, 2019, 2021).

The resulting method is well suited to the analysis of Level-2 satellite data,
enabling the nonstationary analysis of massive, global data, while fully
quantifying (even joint) uncertainties and capturing dependence at all spatial
scales. Our method is applicable to a broad array of Level-2 satellite data
across all NASA Earth science focus areas, where it will allow scientific
questions to be answered more accurately and precisely.

In this work, we focus on one important application: providing a new high-
resolution SST product based on Level-2 MODIS data with more than 43 million
observations. On the Cheyenne supercomputer at the National Center for
Atmospheric Research, we obtain the full posterior predictive distribution of
SST on a fine grid in a computationally highly efficient manner, as well as
with high statistical efficiency. Our flexible nonstationary model captures
the spatially varying dependence structure, which leads to substantially
improved prediction relative to the stationary model (Huang et al.,, 2019;
Blake et al.,, 2019, 2021) in terms of both accuracy and uncertainty
quantification. Thus, our approach is highly useful for oceanographic research
that requires fine-scale joint distributions. For example, the magnitude of
gradients in ultra-high-resolution SST products is strongly related to the
intensity of the upwelling of cold, nutrient-rich water from depth (Vazquez-
Cuervo et al.,, 2017), and thus accurate gradient calculations allow proper
investigation of changes in upwelling and the impact on coastal fisheries due
to climatic changes. To our knowledge, our MODIS SST dataset is the largest
ever analyzed using a probabilistic Gaussian process model, especially with a
nonstationary covariance structure.

The remainder of the manuscript is organized as follows. Section 2 introduces
the nonstationary $M$-RA model and our distributed-computing implementation.
Section 3 discusses the SST data, important aspects related to performing
predictions, the nonstationary prediction results with the improvement over a
stationary model shown by several assessment metrics, and the complete and
high-resolution SST product. Section 4 concludes and points out future
research directions.

##  2 Modeling

###  2.1 A Brief Review of the $M$-RA Model

In this section, we provide a brief review of the $M$-RA method. For more
details on the $M$-RA model, the reader is directed to Katzfuss, (2017).
Assume that $y(\cdot)\sim\mathcal{GP}\big{(}0,C(\cdot,\cdot)\big{)}$ is a
Gaussian process (GP) with covariance function $C(\cdot,\cdot)$. Direct
inference based on $n$ observations of this GP is computationally infeasible
for large $n$. Similar to wavelet approaches, the $M$-RA process is specified
as a linear combination of basis functions at multiple levels of spatial
resolution, which can capture spatial structure from very fine to very large
scales. The $M$-RA requires a hierarchical partitioning of the domain
$\mathcal{D}$, in which we recursively split each region into $J$ (e.g.,
$J=2$) subregions, up to some level $M$ (e.g., $M=10$). Within each region, we
specify a grid consisting of a small number $r$ (e.g., $r=30$) of locations
(called knots). These knots will form the centers of basis functions, whose
shape is then determined iteratively by the predictive process (Banerjee et
al.,, 2008). The basis functions at each resolution are approximately optimal,
in that the predictive process based on $r$ knots is the Nyström approximation
of the first $r$ terms in the Karhunen–Loève expansion of the process to be
approximated (Sang and Huang,, 2012). More precisely, the $M$-RA process is a
combination of the basis functions at all $M$ resolutions:

$y_{M\textrm{-RA}}(\mathbf{s})=\sum_{m=1}^{M}\mathbf{b}_{m}(\mathbf{s})^{\prime}\bm{\eta}_{m},\quad\mathbf{s}\in\mathcal{D},$
(1)

where for each resolution $m=1,\ldots,M$, $\mathbf{b}_{m}(\cdot)$ is a vector
of $p_{m}r$ spatial basis functions, $\bm{\eta}_{m}\stackrel{{\scriptstyle
ind.}}{{\sim}}\mathcal{N}_{p_{m}r}(\mathbf{0},\bm{\Lambda}_{m})$ is a vector
of corresponding basis-function weights, and $p_{m}$ is the number of regions
at resolution $m$. Figure 1 illustrates the $M$-RA in a simple toy example.

Figure 1: A toy example of basis functions $\mathbf{b}_{m}(\cdot)$ in the
$M$-RA on a one-dimensional spatial domain $\mathcal{D}=[0,1]$ with $M=4$
levels, $J=2$ subregions per region (vertical lines), and $r=1$ knot per
region.

One important property of the $M$-RA is that both the basis functions
$\mathbf{b}_{m}(\cdot)$ and the covariance matrix of the weights,
$\bm{\Lambda}_{m}$, are automatically chosen to optimally approximate the true
process $y(\cdot)$ and its covariance function $C(\cdot,\cdot)$. Because
$C(\cdot,\cdot)$ can be any covariance function suitable for the application
or data of interest, the basis-function model is easily interpretable and
highly flexible (e.g., allowing the dependence structure of the process of
interest to vary over the spatial domain). Another property, crucial for
computational feasibility, is that for a larger $m$, the number of basis
functions in $\mathbf{b}_{m}(\cdot)$ increases while the support of each
function decreases. This is done in such a way that any location
$\mathbf{s}\in\mathcal{D}$ is only covered by $r$ basis functions per
resolution (i.e., $\mathbf{b}_{m}(\mathbf{s})$ has only $r$ nonzero elements),
and thus $Mr$ basis functions in total. Similarly, $\bm{\Lambda}_{m}$ is a
block-diagonal matrix of increasing dimension for increasing $m$, but the
blocks along the diagonal are only of size $r\times r$. This block-sparse form
ensures that all computations involving the $M$-RA are highly scalable for
massive datasets.

The $M$-RA can also be viewed as a special case of an ordered conditional
approximation (Vecchia,, 1988), resulting in a block-sparse inverse Cholesky
factor (Katzfuss and Guinness,, 2021; Jurek and Katzfuss,, 2021).

###  2.2 Inference via Distributed Computing

In general, $M$-RA inference proceeds in two stages. First, the prior
covariance matrices $\bm{\Lambda}_{m}$ and the basis functions
$\mathbf{b}_{m}(\cdot)$ (at knot/observation locations) are computed
sequentially for $m=1,2,\ldots,M$. Second, the posterior distributions of the
basis-function weights $\bm{\eta}_{m}$ given the data are computed
sequentially in the opposite direction, for $m=M,M-1,\ldots,1$.

When the data size is large, even though the $M$-RA framework hierarchically
imposes low-rank representation and block-independence structure within each
subregion at each resolution (which eases the computational burden and memory
storage requirement), performing predictions is still intractable if only one
computer is used. Therefore, a distributed-computing implementation is a must
when building a high-resolution product from massive satellite observations.
For ease of presentation, we explain the distributed computing scheme
conceptually on a one-dimensional domain with the same setting as in Figure 1,
where $M=4$ and $J=2$. Additional technical details are provided in a tech
note (Huang et al.,, 2019). We assume a cluster of three nodes is available,
where each node has its own memory and computing cores, and communication
among nodes is possible via networks. In this example, we have eight regions
at the finest level. We assign similar numbers of regions at the finest level
to each node and subsequently all ancestral regions at coarser levels. Figure
2 shows the working regions for each node.

Figure 2: Working regions for each node, indicated in blue.

Several prior quantities related to the covariance among knots need to be
computed and stored. These computations go from $m=1$ to 4, and at each level
$m$, all computations in different regions are independent and can be executed
in parallel. However, the computations for posterior distributions given the
data from $m=3$ to 1 are more complicated. Figure 3 summarizes different
computation scenarios when obtaining the posterior.

Figure 3: Summary of posterior computations. At each computing level, green
indicates working regions at the previous level; red with a cross indicates
regions where this node handled all children, so no synchronization is needed;
red regions with “S” or “W” require synchronization, where “S” represents the
“supervisor” that receives data and “W” is the “worker” that sends data.

When $m=3$, relevant quantities in the first region for Node 1 (red with a
cross) can be obtained without the need to communicate with other nodes,
because both its children at level 4 were handled in Node 1. The second region
for Node 1 needs synchronization with Node 2, because some of its children
were handled by Node 1 and some by Node 2. Node 2 sends previously computed
quantities for its second child at level 4 to Node 1, so Node 1 has all the
needed information to compute the posterior for the second region (red with a
“S”). After the synchronization, the second region (red with a “W”) for Node 2
is eliminated, because Node 1 has all its information, and it is no longer
green when $m=2$ for Node 2. Similar computations are executed until we finish
all the computations at level 1.

For predictions, we then need to proceed from Level 1 to $M$, again
calculating relevant quantities, and the prediction distributions for new
locations are obtained at level $M$. Computations for different regions at the
same level are executed in parallel and only depend on ancestral regions.

We use Message Passing Interface (MPI) to operate the synchronization in C++.
We also use OpenMP to execute parallel computations of different regions at
the same level for further computing acceleration.

###  2.3 Nonstationarity Modeling

Previous work has shown the excellent prediction performance of the $M$-RA
model with a stationary model for very large datasets (Huang et al.,, 2019;
Blake et al.,, 2019, 2021). However, for massive datasets on the globe, the
stationarity assumption typically does not hold and flexible nonstationary
models are required. Pioneering works have investigated the extension of the
$M$-RA to account for nonstationarity. For example, Benedetti et al., (2021)
proposed a mixture $M$-RA model where a shrinking prior is induced as one
component for the $M$-RA basis function weights in a Bayesian framework.

In this work, we directly extend the underlying covariance model to be
nonstationary and provide high-performance computing implementation for large
datasets. We build the nonstationary covariance model based on the approach
proposed by Paciorek and Schervish, (2006). The closed form for the
nonstationary correlation function $R^{NS}(\mathbf{s},\mathbf{s}^{\prime})$ on
$\mathbb{R}^{p}$, $p\in\mathbbm{N}_{+}$ is,

$R^{NS}(\mathbf{s},\mathbf{s}^{\prime})=\lvert\bm{\Sigma}_{\mathbf{s}}\rvert^{\frac{1}{4}}\lvert\bm{\Sigma}_{\mathbf{s}^{\prime}}\rvert^{\frac{1}{4}}\left|\frac{\bm{\Sigma}_{\mathbf{s}}+\bm{\Sigma}_{\mathbf{s}^{\prime}}}{2}\right|^{-\frac{1}{2}}R^{S}\left(\mathbf{s}-\mathbf{s}^{\prime};\dfrac{\bm{\Sigma}_{\mathbf{s}}+\bm{\Sigma}_{\mathbf{s}^{\prime}}}{2}\right),$
(2)

where $\bm{\Sigma}_{\mathbf{s}}$ is the covariance matrix of a Gaussian kernel
centered at $\mathbf{s}$, and $R^{S}(\mathbf{h};\bm{\Sigma})$ is a stationary
correlation function in Euclidean space. We consider a stationary Matérn
correlation function
$R^{S}(\mathbf{h};\bm{\Sigma})=\big{(}\sqrt{2\nu\mathbf{h}^{\top}\bm{\Sigma}^{-1}\mathbf{h}}\big{)}^{\nu}K_{\nu}\big{(}\sqrt{2\nu\mathbf{h}^{\top}\bm{\Sigma}^{-1}\mathbf{h}}\big{)}/\\{\Gamma(\nu)2^{\nu-1}\\}$,
where $K_{\nu}(\cdot)$ is the modified Bessel function of the second kind with
order $\nu$. It is also possible to allow the smoothness parameter $\nu$ to
vary spatially. However, identifiability issues are noticed when
$\bm{\Sigma}_{\mathbf{s}}$ and $\nu$ are both spatially varying (Anderes and
Stein,, 2011).

With this approach, we can model our data as points on the sphere using
chordal distance to yield a valid nonstationary covariance on the sphere,
which means $\mathbf{s}\in\mathcal{D}=\mathbb{R}^{3}$. Applying spatially
varying partial sill parameters $\sigma^{2}(\mathbf{s})$, nugget parameters
$\tau^{2}(\mathbf{s})$, and range parameters $\beta(\mathbf{s})$ determining
the isotropic kernel
$\bm{\Sigma}_{\mathbf{s}}=\beta^{2}(\mathbf{s})\mathbf{I}_{3}$, we have all
the parameters
$\bm{\theta}(\mathbf{s})=\\{\sigma^{2}(\mathbf{s}),\beta(\mathbf{s}),\tau^{2}(\mathbf{s}),\nu\\}$
that need to be inferred, and the nonstationary covariance model is

$\begin{array}[]{rcl}C(\mathbf{s},\mathbf{s}^{\prime};\bm{\theta}(\mathbf{s}))&=&\sigma(\mathbf{s})\sigma(\mathbf{s}^{\prime})\left(\dfrac{2\beta(\mathbf{s})\beta(\mathbf{s}^{\prime})}{\beta^{2}(\mathbf{s})+\beta^{2}(\mathbf{s}^{\prime})}\right)^{3/2}\left(\dfrac{2\sqrt{\nu}\|\mathbf{s}-\mathbf{s}^{\prime}\|}{\sqrt{\beta^{2}(\mathbf{s})+\beta^{2}(\mathbf{s}^{\prime})}}\right)^{\nu}\times\\\
&&K_{\nu}\left(\dfrac{2\sqrt{\nu}\|\mathbf{s}-\mathbf{s}^{\prime}\|}{\sqrt{\beta^{2}(\mathbf{s})+\beta^{2}(\mathbf{s}^{\prime})}}\right)/\\{\Gamma(\nu)2^{\nu-1}\\}+\tau^{2}(\mathbf{s})\mathbbm{1}_{\\{\mathbf{s}=\mathbf{s}^{\prime}\\}},\end{array}$
(3)

where $\mathbbm{1}_{\\{\\}}$ denotes the indicator function (see Appendix for
derivation).

To be able to calculate the nonstationary covariance between any pair of
locations shown in Equation (2), estimates of $\bm{\theta}(\mathbf{s})$ are
required for any location $\mathbf{s}\in\mathcal{D}$. We propose to compute
local estimates on a regular grid; these estimates can then be smoothed by
regressing them on a set of basis functions to obtain a smooth map
$\bm{\theta}(\mathbf{s})$. More specifically, after specifying a grid over the
study area, a local estimate of $\bm{\theta}(\mathbf{s}^{\ast})$ is obtained
at each grid point $\mathbf{s}^{\ast}$ by the following procedure (which can
be carried out in parallel):

  1. 1.

We build two boxes centered at $\mathbf{s}^{\ast}$ with different sizes,
$B_{1}(\mathbf{s}^{\ast})$ and $B_{2}(\mathbf{s}^{\ast})$, where
$B_{2}(\mathbf{s}^{\ast})$ is larger than $B_{1}(\mathbf{s}^{\ast})$.

  2. 2.

We randomly choose up to $N_{s}$ of the observations in
$B_{1}(\mathbf{s}^{\ast})$ to obtain information on the local short-range
dependence, and $N_{b}$ observations in $B_{2}(\mathbf{s}^{\ast})\backslash
B_{1}(\mathbf{s}^{\ast})$ that contain information on the long-range
dependence.

  3. 3.

We find the maximum likelihood estimates of the partial sill $\sigma^{2}$,
range $\beta$, smoothness $\nu$, and nugget $\tau^{2}$ in the stationary
Matérn covariance function
$\mathcal{M}(\mathbf{s}-\mathbf{s}^{\prime})=\sigma^{2}\big{(}\sqrt{2\nu}\|\mathbf{s}-\mathbf{s}^{\prime}\|/\beta\big{)}^{\nu}K_{\nu}\big{(}\sqrt{2\nu}\|\mathbf{s}-\mathbf{s}^{\prime}\|/\beta\big{)}/\\{\Gamma(\nu)2^{\nu-1}\\}+\tau^{2}\mathbbm{1}_{\\{\mathbf{s}=\mathbf{s}^{\prime}\\}}$
based on the $N_{s}+N_{b}$ chosen observations. We then assign the obtained
parameter estimates to the lattice point $\mathbf{s}^{\ast}$.

Once estimates are obtained for each grid point, we smooth the local parameter
estimates and extend them to locations not on the grid by regressing the
estimates onto a set of Wendland basis functions. We place $N_{w}$ compactly
supported Wendland basis functions centered at $\tilde{\mathbf{s}}$ with the
length of support $\ell$ as follows,

$\textrm{WL}(\mathbf{s};\tilde{\mathbf{s}})=\frac{1}{3}\bigg{(}1-\frac{\|\mathbf{s}-\tilde{\mathbf{s}}\|}{\ell}\bigg{)}^{6}\bigg{(}35\Big{(}\frac{\|\mathbf{s}-\tilde{\mathbf{s}}\|}{\ell}\Big{)}^{2}+18\frac{\|\mathbf{s}-\tilde{\mathbf{s}}\|}{\ell}+3\bigg{)}\times\mathbbm{1}_{\\{\|\mathbf{s}-\tilde{\mathbf{s}}\|<\ell\\}},$

We place $\tilde{\mathbf{s}}$ in the study area according to an icosahedral
strategy, which ensures equivalent spacing among $\tilde{\mathbf{s}}$ across
the domain. For each parameter estimation, we can add a Lasso penalty to avoid
overfitting, and the value of $N_{w}$, the penalty parameter for the Lasso
regularization, and the support length $\ell$ can be optimized via cross-
validation with respect to the mean squared prediction error (MSPE).

##  3 MODIS Sea Surface Temperature Analysis

###  3.1 Data Description

We use GHRSST Level-2P Global Sea Surface Skin Temperature data from the MODIS
instrument on the NASA Aqua satellite (GDS2), which are available from NASA’s
Physical Oceanography Distributed Active Archive Center (PODAAC) at
https://podaac.jpl.nasa.gov/dataset/MODIS_A-JPL-L2P-v2019.0. Both daytime and
nighttime measurements collected at a $1\textrm{km}^{2}$ resolution are
provided. For exploratory purposes, we analyze the daytime and nighttime data
of days 1 through 7 of January, April, October, and December to inspect
seasonal SST fluctuations and determine when the highest number of
observations are available. We notice that the data in these dates tend to
have more observations in the nighttime than daytime. Examples of the
nighttime SST are shown in Figure S1 in the Supplementary Material. The raw
data contain a quality level flag from 0 to 5, indicating the accuracy of the
recorded observations. The User Guide in the data product suggests that
observations with quality levels 2 or greater are usable, so we eliminate all
the observations with quality levels less than 2.

We select the date with the largest amount of available data, which is
nighttime, April 7, 2019, as the primary data set of consideration for the
remainder of our analysis, illustrated in Figure 4. Nighttime SST from April 1
to 6, 2019, will be used as additional reference data sources for the
covariance modeling.

Figure 4: Global nighttime SST on April 7, 2019, with $n=45{,}057{,}134$
observations.

###  3.2 Removing Mean Structure

We denote the SST at location $\mathbf{s}\in\mathcal{D}$ by $z(\mathbf{s})$.
We model $z(\mathbf{s})=\mu(\mathbf{s})+y(\mathbf{s})$, where
$\mu(\mathbf{s})$ is the mean structure to be estimated so that the resulting
residual random field $y(\mathbf{s})$ can be assumed a zero-mean Gaussian
process. We tested a large suite of trend models, including (local) polynomial
regressions, thin-plate splines, and cubic smoothing splines. Of these models,
the cubic smoothing spline class of models as a function of latitude yields
the best balance between accurate estimation of the mean and over-fitting.
Ten-fold cross-validation for models with the number of cubic smoothing spline
basis functions ranging from five to fifteen is used to select the optimal
model. In each fold in the cross-validation, we bin the training data into 180
latitudinal bins and compute the mean SST of each bin. We fit with a number of
cubic smoothing splines to these binned data and predict the SST on the
validation locations. Figure S2 in the Supplementary Material shows the
results of the cross-validation error versus the selected number of basis
functions. We select 11 cubic smoothing spline functions for their relatively
low degrees of freedom and the ability to accurately approximate the mean
function $\mu(\mathbf{s})$. Figure 5 depicts the fit using the 11 selected
basis functions. Because the observations are very sparse for latitudes
greater than $60^{\circ}$, we choose our study area as
$[180^{\circ}\textrm{W},180^{\circ}\textrm{E}]\times[60^{\circ}\textrm{S},60^{\circ}\textrm{N}]$,
leading to 43,802,698 observations in total. Figure 6 shows the residual
spatial field $y(\mathbf{s})$ in the study area after the mean trend removal.
More details about the residuals are shown in Figure S3 in the Supplementary
Material.

Figure 5: Averaged binned SST in 180 latitudinal bins (black dots) with the
cubic smoothing spline fit superimposed (orange line). Figure 6: SST residuals
on April 7, 2019. The removed mean is estimated by cubic smoothing splines
using 11 basis functions.

###  3.3 $M$-RA Domain Decomposition

To perform $M$-RA inference for our global SST data, we hierarchically
decompose the spatial domain using GIS software (Environmental Systems
Research Institute,, 2011) in a manner such that we consider realistic spatial
dependence informed by physical boundaries and oceanographic knowledge.
Broadly speaking, we place knots and partitions in a way that weakens
dependence across land barriers, by placing partition boundaries along (i.e.,
on top of) land, and avoiding knot placement at lower resolutions too close to
land barriers. We choose $M=15$ levels in total. For level $m=2,\ldots,M-1$,
each region in level $m-1$ is split into $J=2$ subregions as constructed
regions in level $m$. We place $r=49$ knots in each region for level $m<M$. We
use global landmasses accurate to 10 km resolution based on the World Geodetic
System 1984 (WGS84) coordinate system to create the ocean map for the study
area between latitudes $60^{\circ}$S and $60^{\circ}$N. We only place the
knots on the open ocean, which, as customary for SST products, excludes
smaller bodies of water (i.e., Hudson Bay, the Mediterranean Sea, the Baltic
Sea, the Black Sea, the Caspian Sea, the Red Sea, and the Persian Gulf).
Placement of knots for levels $m=1,\ldots,10$ and region boundaries for
$m=2,\ldots,10$ are manually performed in the GIS software. In the ocean,
knots are placed strategically near partition boundaries for both the current
and the next level, so that the introduced approximation errors near
boundaries are small. Figure 7 shows the region boundaries and knots for the
first four levels.

A $m=1$

B $m=2$

C $m=3$

D $m=4$

Figure 7:  First four levels of the hierarchical domain decomposition. Knots
and partitions at a given level are shown in the same color. In each plot, we
show the knots at the current level and the boundaries of the current and
coarser partitions. Partitions are not placed on land as we only consider SST
on open water and therefore land masses are not part of the domain.

For levels $m\geq 11$, the region areas are small (the maximum region area at
this level is less than $7.34\times 10^{5}$ km2), and most of the large-scale
spatial correlation structure has been adequately captured in coarser levels.
Therefore, for simplicity, we split regions perpendicular to the longer
spatial dimension at the coordinate mean and place knots randomly in each
subregion. To capture the very fine-scale variation, the observation locations
are used as knots in the finest level $m=M=15$. To avoid slow computation for
regions that happen to contain a large number of observations, we do not fix
$J=2$ but keep partitioning the regions in $m=14$ until all the obtained
subregions, which are regions in $m=15$, contain fewer than 2,000
observations; hence, the number of subregions $J$ varies between regions at
$m=14$, with $J=2^{k}$ for some $k\in\mathbb{N}$ for each region.

###  3.4 Assessing Nonstationarity

To specify the nonstationary covariance model, we create a $2^{\circ}\times
2^{\circ}$ grid in the study area
$[180^{\circ}\textrm{W},180^{\circ}\textrm{E}]\times[60^{\circ}\textrm{S},60^{\circ}\textrm{N}]$
and obtain the local estimates of $\bm{\theta}(\mathbf{s}^{\ast})$ via the
procedure described in Section 2.3. Grid points on lands are removed. We
denote by $\texttt{lon}(\mathbf{s}^{\ast})$ and
$\texttt{lat}(\mathbf{s}^{\ast})$ the associated longitude and latitude for
$\mathbf{s}^{\ast}$. At each grid point $\mathbf{s}^{\ast}$,
$B_{1}(\mathbf{s}^{\ast})$ and $B_{2}(\mathbf{s}^{\ast})$ are chosen as
$B_{1}(\mathbf{s}^{\ast})=[\texttt{lon}(\mathbf{s}^{\ast})-2^{\circ},\texttt{lon}(\mathbf{s}^{\ast})+2^{\circ}]\times[\texttt{lat}(\mathbf{s}^{\ast})-2^{\circ},\texttt{lat}(\mathbf{s}^{\ast})+2^{\circ}]$
and
$B_{2}(\mathbf{s}^{\ast})=[\texttt{lon}(\mathbf{s}^{\ast})-20^{\circ},\texttt{lon}(\mathbf{s}^{\ast})+20^{\circ}]\times[\texttt{lat}(\mathbf{s}^{\ast})-20^{\circ},\texttt{lat}(\mathbf{s}^{\ast})+20^{\circ}]$.
We choose $N_{s}=800$ and $N_{b}=100$ for the observations in
$B_{1}(\mathbf{s}^{\ast})$ and $B_{2}(\mathbf{s}^{\ast})\backslash
B_{1}(\mathbf{s}^{\ast})$, respectively. We observe noisy estimates if the
number of chosen observations in a small box $B_{1}(\mathbf{s}^{\ast})$ is
insufficient due to a gap in the satellite data. Therefore, we eliminate all
the estimates at lattice points $\mathbf{s}^{\ast}$ where the number of
available observations in $B_{1}(\mathbf{s}^{\ast})$ is smaller than 800.
Figure S4 in the Supplementary Material shows the local estimates of the
smoothness parameter from April 1 to 7, 2019. We observe that the variability
of the smoothness estimates is small, which meets the requirement that a fixed
smoothness parameter is used in the nonstationary model (2). As the majority
of the estimates are around 0.5, we fix $\nu=0.5$, and thus $R^{S}(\cdot)$ in
Equation (2) is the exponential function, which brings computational benefits
by avoiding the expensive evaluation of the Bessel function. The nonstationary
covariance function now becomes (see the derivation in the Appendix):

$C(\mathbf{s},\mathbf{s}^{\prime})=\sigma(\mathbf{s})\sigma(\mathbf{s}^{\prime})\left(\dfrac{2\beta(\mathbf{s})\beta(\mathbf{s}^{\prime})}{\beta^{2}(\mathbf{s})+\beta^{2}(\mathbf{s}^{\prime})}\right)^{3/2}\exp\left(-\dfrac{\|\mathbf{s}-\mathbf{s}^{\prime}\|}{\sqrt{\frac{1}{2}\big{(}\beta^{2}(\mathbf{s})+\beta^{2}(\mathbf{s}^{\prime})\big{)}}}\right)+\tau^{2}(\mathbf{s})\mathbbm{1}_{\\{\mathbf{s}=\mathbf{s}^{\prime}\\}}.$
(4)

We then locally estimate the parameters $\sigma(\cdot)$, $\beta(\cdot),$ and
$\tau(\cdot)$ at lattice points again with the exponential covariance
function. Figure S5 in the Supplementary Material depicts all the local
estimates $\hat{\bm{\theta}}(\mathbf{s}^{\ast})$ from April 1 to 7, 2019 for
the nonstationary covariance function (4). We observe that the patterns during
the week are consistent, suggesting that the spatial covariance structure
remains similar during the short time.

As explained in Section 2.3, we smooth the local parameter estimates by
regressing the estimates onto a set of Wendland basis functions. To compensate
for data vacancy in some areas due to gaps on April 7, we stack all the
parameter estimates from April 1 to 7. With exploratory analysis, we find that
logarithmic transformations to the parameter estimates (which must be
positive) make them closer to normality, so we perform the regression on the
logarithmic scale. We place Wendland basis functions centered at
$\tilde{\mathbf{s}}$ in the study area
$[-180^{\circ}\textrm{W},-180^{\circ}\textrm{E}]\times[60^{\circ}\textrm{S},60^{\circ}\textrm{N}]$
and remove $\tilde{\mathbf{s}}$ falling on land, so that basis functions only
obtain their maximum on the ocean. By cross-validation of $N_{w}$ and $\ell$,
we obtain 1399, 1239, and 1337 non-zero basis-function coefficients for the
partial sill, range, and nugget, respectively, that lead to optimal smoothing
results. The center locations of the selected basis functions are shown in
Figure 8(A). Figure 8 also illustrates the spatially varying parameter
estimate map built from the Wendland basis functions.

Figure 8: (A): The Wendland basis function center locations with non-zero
coefficients. (B)–(D): Spatially varying sill $\sigma^{2}(\mathbf{s})$, range
$\beta(\mathbf{s})$ in km, and nugget $\tau^{2}(\mathbf{s})$ estimates in the
exponential covariance function built from 1399, 1239, 1337 Wendland basis
functions, respectively.

These spatially varying parameters are used in the covariance function (4) to
be approximated by the $M$-RA .

###  3.5 Prediction Performance Assessment

A common feature in satellite datasets is large gaps that exist due to the
inability to obtain measurements under certain atmospheric conditions such as
cloud cover. To provide a complete high-resolution data product, the ability
to fill large gaps with accurate (probabilistic) predictions is of great
importance. Therefore, we randomly place large gaps and hold out the
observations within the gaps to assess prediction performance. More precisely,
we randomly sample a location (lon, lat) in the ocean and build a gap centered
at the location with size $10^{\circ}\times 10^{\circ}$ as
$B=[\textrm{lon}-5^{\circ},\textrm{lon}+5^{\circ}]\times[\textrm{lat}-5^{\circ},\textrm{lat}+5^{\circ}]$.
We hold out all the observations in $B$ as testing data $D^{\textrm{test}}$
and use the rest as training data $D^{\textrm{train}}$, i.e.,
$D^{\textrm{test}}=\\{y(\mathbf{s}_{i}):\big{(}\texttt{lon}(\mathbf{s}_{i}),\texttt{lat}(\mathbf{s}_{i})\big{)}\in
B\\}$ and
$D^{\textrm{train}}=\\{y(\mathbf{s}_{i}):\big{(}\texttt{lon}(\mathbf{s}_{i}),\texttt{lat}(\mathbf{s}_{i})\big{)}\in
B^{c}\\}$. We repeat this procedure 100 times to obtain datasets for 100
independent experiments:
$D^{\textrm{train}}_{1},D^{\textrm{test}}_{1},\ldots,D^{\textrm{train}}_{100},D^{\textrm{test}}_{100}$.
For a more uniform assessment of prediction performance across the 100
experiments, in each experiment $k=1,\ldots,100$, we make sure the number of
observations in $B_{k}$ is greater than 50,000 (by resampling the center point
until the criterion is met) and randomly choose 50,000 samples in building the
testing dataset. Figure 9 shows one example, the training and testing dataset
for Gap 1, $D^{\textrm{train}}_{1},D^{\textrm{test}}_{1}$.

Figure 9: Training and testing data (red box) for Gap 1,
$D^{\textrm{train}}_{1},D^{\textrm{test}}_{1}$.

In each experiment $k$, we perform the $M$-RA with domain configurations shown
in Section 3.3 to obtain predictive distributions of all locations in
$D^{\textrm{test}}_{k}$ conditional on $D^{\textrm{train}}_{k}$ with the mean
structure described in Section 3.2 and covariance model (4). We use the MSPE
to assess the point prediction performance and the log-score and Continuous
Ranked Probability Score (CRPS) to assess the probabilistic prediction
performance. The log-score and CRPS are strictly proper scoring rules for the
assessment of probabilistic prediction quality (the scores attain the minimum
when the realized observations are exactly from the predictive distribution;
Gneiting and Raftery,, 2007). For Gaussian predictive distributions with
prediction mean $\hat{y}(\mathbf{s})$ and standard deviation
$\hat{\varsigma}(\mathbf{s})$ at location $\mathbf{s}$, the log-score and CRPS
have the form as follows (e.g., Gneiting and Katzfuss,, 2014),

$\begin{array}[]{rcl}\textrm{Log-
score}(\mathbf{s})&=&\big{[}\log(2\pi)+{\\{y(\mathbf{s})-\hat{y}(\mathbf{s})\\}^{2}}/{\hat{\varsigma}^{2}(\mathbf{s})}\big{]}/2+\log\big{(}\hat{\varsigma}(\mathbf{s})\big{)},\\\
\textrm{CRPS}(\mathbf{s})&=&\\{y(\mathbf{s})-\hat{y}(\mathbf{s})\\}\Big{[}2\Phi\Big{(}\dfrac{y(\mathbf{s})-\hat{y}(\mathbf{s})}{\hat{\varsigma}(\mathbf{s})}\Big{)}-1\Big{]}+2\hat{\varsigma}\phi\Big{(}\dfrac{y(\mathbf{s})-\hat{y}(\mathbf{s})}{\hat{\varsigma}(\mathbf{s})}\Big{)}-\dfrac{\hat{\varsigma}(\mathbf{s})}{\sqrt{\pi}},\end{array}$

where $\Phi(\cdot)$ and $\phi(\cdot)$ are the standard normal cumulative
distribution function and probability density function, respectively. In each
experiment, the log-score and CRPS reported are the mean of the corresponding
metrics over all the testing locations.

For comparison, we also show the results from the stationary $M$-RA model. For
the stationary model, we use the $M$-RA first to estimate the four stationary
covariance model parameters via maximum likelihood estimation and then obtain
predictions based on the estimated stationary $M$-RA model. The MSPE, log-
score, and CPRS from the stationary and nonstationary models in the 100
experiments are summarized in Figure 10 and Table 1. We observe that
incorporation of nonstationarity yields more accurate point prediction and
substantially improves the predictive distribution. We also provide an
intuitive way to illustrate the improvement of the probabilistic forecast in
Figure 11, where we calculate the proportion of testing observations falling
into the interval obtained by the Gaussian predictive distribution with the
nominal level. Though slightly inflating the variance, the nonstationary model
has confidence interval coverage substantially closer to the theoretical
values than the stationary model for all levels.

Figure 10: Histogram of MSPE, log-score, and CPRS from the stationary (filled with black dots) and nonstationary (solid blue) models in the 100 experiments. Figure 11: Confidence interval coverage from the stationary (black dashed line) and nonstationary (dotted blue line) models in the 100 experiments. The solid magenta line shows the theoretical coverage from the nominal level. Table 1: Average of MSPE, log-score, and CPRS of the corresponding metrics over the 100 experiments for the stationary and nonstationary models, respectively. “#gaps nonstationary better” provides the number of experiments where the nonstationary model improves upon the stationary model with respect to the corresponding metric.  Model | MSPE | Log score | CRPS
---|---|---|---
Stationary | 1.08 | 1.57 | 0.60
Nonstationary | 1.01 | 1.21 | 0.49
#gaps nonstationary better | 62 | 96 | 92

###  3.6 High-Resolution Product

The high-resolution SST product for the entire study area using the
nonstationary $M$-RA model is shown in Figure 12. The complete map of
prediction means resolves SST variation even on very fine scales and does not
show obvious artifacts on boundaries between partitions, implying a good
domain decomposition scheme in our analysis. The prediction standard
deviations reflect the missing data pattern and show nonstationary structures
across the globe (e.g., large values in the ocean near the south of the
Labrador Sea). The $M$-RA prediction is performed on 25 nodes on the Cheyenne
supercomputer (Computational and Information Systems Laboratory,, 2019), each
of which has 36 2.3-GHz Intel Xeon cores, and takes 35 seconds (each node
launches one MPI process and 36 OpenMP threads).

Figure 12: High-resolution SST product using the nonstationary $M$-RA model.

##  4 Summary and Future work

Addressing the urgent need to provide statistically rigorous full-coverage
global data products based on exceedingly massive satellite observations, we
have developed a nonstationary version and distributed-memory implementation
of the $M$-RA modeling framework. We employ this method, following a careful
analysis of the mean and nonstationary covariance structure, to provide a
daily global SST product based on 43 million observations. The prediction
performance of the nonstationary model is substantially better compared to the
stationary model and highlights the need to incorporate nonstationarity when
dealing with global geophysical data. Additional future satellite instruments
with increasingly high resolution are expected to strengthen the need for the
kind of methodology provided here.

One particularly attractive feature of a statistically rigorous approach is
the possibility of extending the SST product into the form of joint
probability distributions, which in turn enables us to resolve SST gradients
with proper uncertainties at any location in the ocean. This is important for
the study of atmosphere-ocean coupling, and it allows for answering scientific
questions regarding, for example, the effect of climate change on upwelling
and the implications on coastal fisheries. We plan on conducting such studies
in collaboration with oceanographers in future research.

## Acknowledgments

MK was partially supported by National Science Foundation (NSF) Grants
DMS–1521676, DMS–1654083, and DMS–1953005, and by the National Aeronautics and
Space Administration (80NM0018F0527). We would like to thank Samuel H. Blake
for assistance with GIS. We would like to acknowledge high-performance
computing support from Cheyenne (doi:10.5065/D6RX99HX) provided by NCAR’s
Computational and Information Systems Laboratory, sponsored by the National
Science Foundation.

## Appendix

We derive the form for the nonstationary exponential covariance function given
in Equation (3). We let
$\bm{\Sigma}_{\mathbf{s}}=\beta^{2}(\mathbf{s})\mathbf{I}_{3},\bm{\Sigma}_{\mathbf{s}^{\prime}}=\beta^{2}(\mathbf{s}^{\prime})\mathbf{I}_{3}$
and plug into Equation (2) and we have:

$R^{NS}(\mathbf{s},\mathbf{s}^{\prime})=\left(\beta(\mathbf{s})^{6}\right)^{\frac{1}{4}}\left(\beta(\mathbf{s}^{\prime})^{6}\right)^{\frac{1}{4}}\left(\frac{2}{\beta^{2}(\mathbf{s})+\beta^{2}(\mathbf{s}^{\prime})}\right)^{{3}/{2}}R^{S}\left(\mathbf{s}-\mathbf{s}^{\prime};\dfrac{\beta^{2}(\mathbf{s})+\beta^{2}(\mathbf{s}^{\prime})}{2}\mathbf{I}_{3}\right).$

Considering
$R^{S}(\mathbf{h};\bm{\Sigma})=\big{(}\sqrt{2\nu\mathbf{h}^{\top}\bm{\Sigma}^{-1}\mathbf{h}}\big{)}^{\nu}K_{\nu}\big{(}\sqrt{2\nu\mathbf{h}^{\top}\bm{\Sigma}^{-1}\mathbf{h}}\big{)}/\\{\Gamma(\nu)2^{\nu-1}\\}$,
we have

$(\mathbf{s}-\mathbf{s}^{\prime})^{\top}\left(\dfrac{\beta^{2}(\mathbf{s})+\beta^{2}(\mathbf{s}^{\prime})}{2}\mathbf{I}_{3}\right)^{-1}(\mathbf{s}-\mathbf{s}^{\prime})=\dfrac{2\|\mathbf{s}-\mathbf{s}^{\prime}\|^{2}}{\beta^{2}(\mathbf{s})+\beta^{2}(\mathbf{s}^{\prime})}$

and subsequently

$R^{NS}(\mathbf{s},\mathbf{s}^{\prime})=\left(\dfrac{2\beta(\mathbf{s})\beta(\mathbf{s}^{\prime})}{\beta^{2}(\mathbf{s})+\beta^{2}(\mathbf{s}^{\prime})}\right)^{3/2}\left(\dfrac{2\sqrt{\nu}\|\mathbf{s}-\mathbf{s}^{\prime}\|}{\sqrt{\beta^{2}(\mathbf{s})+\beta^{2}(\mathbf{s}^{\prime})}}\right)^{\nu}K_{\nu}\left(\dfrac{2\sqrt{\nu}\|\mathbf{s}-\mathbf{s}^{\prime}\|}{\sqrt{\beta^{2}(\mathbf{s})+\beta^{2}(\mathbf{s}^{\prime})}}\right).$

The mapping
$(\mathbf{s},\mathbf{s}^{\prime})\mapsto\sigma(\mathbf{s})\sigma(\mathbf{s}^{\prime})$
is positive definite (e.g., Berg et al.,, 1984). Products of positive definite
kernels are positive definite. After adding an independent white noise with
variance $\tau^{2}(\mathbf{s})$, we get the valid nonstationary exponential
covariance function given in Equation (3).

## References

  * Anderes and Stein, (2011) Anderes, E. B. and Stein, M. L. (2011).  Local likelihood estimation for nonstationary random fields.  Journal of Multivariate Analysis, 102(3):506–520.
  * Appel and Pebesma, (2020) Appel, M. and Pebesma, E. (2020).  Spatiotemporal multi-resolution approximations for analyzing global environmental data.  Spatial Statistics, 38:100465.
  * Banerjee et al., (2008) Banerjee, S., Gelfand, A. E., Finley, A. O., and Sang, H. (2008).  Gaussian predictive process models for large spatial data sets.  Journal of the Royal Statistical Society, Series B, 70(4):825–848.
  * Benedetti et al., (2021) Benedetti, M. H., Berrocal, V. J., and Narisetty, N. N. (2021).  Identifying regions of inhomogeneities in spatial processes via an m-ra and mixture priors.  Biometrics.
  * Berg et al., (1984) Berg, C., Christensen, J. P. R., and Ressel, P. (1984).  Harmonic analysis on semigroups: theory of positive definite and related functions, volume 100.  Springer.
  * Blake et al., (2019) Blake, L. R., Huang, H., Vanderwende, B., and Hammerling, D. (2019).  A shallow-tree multi-resolution approximation for distributed and high-performance computing systems.  NCAR Technical Notes NCAR/TN-559+STR.
  * Blake et al., (2021) Blake, L. R., Huang, H., Vanderwende, B., and Hammerling, D. (2021).  The deep-tree approach: An improved parallel matlab implementation of the multi-resolution approximation for massive spatial data on high-performance computing systems.  NCAR Technical Notes NCAR/TN-565+STR.
  * Bretherton et al., (1976) Bretherton, F. P., Davis, R. E., and Fandry, C. (1976).  A technique for objective analysis and design of oceanographic experiments applied to MODE-73.  Deep Sea Research and Oceanographic Abstracts, 23(7):559–582.
  * Chao et al., (2009) Chao, Y., Li, Z., Farrara, J. D., and Hung, P. (2009).  Blending sea surface temperatures from multiple satellites and in situ observations for coastal oceans.  Journal of Atmospheric and Oceanic Technology, 26(7):1415–1426.
  * Chin et al., (1998) Chin, T., Milliff, R. F., and Large, W. (1998).  Basin-scale, high-wavenumber sea-surface wind fields from a multiresolution analysis of scatterometer data.  Journal of Atmospheric and Oceanic Technology, 15:741–763.
  * Chin et al., (2013) Chin, T. M., Vazquez, J., and Armstrong, E. M. (2013).  A multi-scale, high-resolution analysis of global sea surface temperature.  Algorithm Theoretical Basis Document, version 1.3.
  * Computational and Information Systems Laboratory, (2019) Computational and Information Systems Laboratory (2019).  Cheyenne: HPE/SGI ICE XA System (NCAR Community Computing).  Boulder, CO: National Center for Atmospheric Research.  doi:10.5065/D6RX99HX.
  * Cressie, (1993) Cressie, N. (1993).  Statistics for Spatial Data, revised edition.  John Wiley & Sons, New York, NY.
  * Cressie and Wikle, (2011) Cressie, N. and Wikle, C. K. (2011).  Statistics for Spatio-Temporal Data.  Wiley, Hoboken, NJ.
  * Datta et al., (2016) Datta, A., Banerjee, S., Finley, A. O., and Gelfand, A. E. (2016).  Hierarchical nearest-neighbor gaussian process models for large geostatistical datasets.  Journal of the American Statistical Association, 111(514):800–812.
  * Environmental Systems Research Institute, (2011) Environmental Systems Research Institute (2011).  ArcGIS Desktop: Release 10.
  * Finley et al., (2020) Finley, A. O., Datta, A., and Banerjee, S. (2020).  spNNGP R package for nearest neighbor gaussian process models.  arXiv preprint arXiv:2001.09111.
  * Gneiting and Katzfuss, (2014) Gneiting, T. and Katzfuss, M. (2014).  Probabilistic forecasting.  Annual Review of Statistics and Its Application, 1:125–151.
  * Gneiting and Raftery, (2007) Gneiting, T. and Raftery, A. E. (2007).  Strictly proper scoring rules, prediction, and estimation.  Journal of the American statistical Association, 102(477):359–378.
  * Heaton et al., (2019) Heaton, M. J., Datta, A., Finley, A. O., Furrer, R., Guinness, J., Guhaniyogi, R., Gerber, F., Gramacy, R. B., Hammerling, D., Katzfuss, M., Lindgren, F., Nychka, D. W., Sun, F., and Zammit-Mangion, A. (2019).  A case study competition among methods for analyzing large spatial data.  Journal of Agricultural, Biological, and Environmental Statistics, 24(3):398–425.
  * Huang et al., (2019) Huang, H., Blake, L. R., and Hammerling, D. M. (2019).  Pushing the limit: A hybrid parallel implementation of the multi-resolution approximation for massive data.  NCAR Technical Notes NCAR/TN-558+STR.
  * Jurek and Katzfuss, (2021) Jurek, M. and Katzfuss, M. (2021).  Multi-resolution filters for massive spatio-temporal data.  Journal of Computational and Graphical Statistics.
  * Katzfuss, (2017) Katzfuss, M. (2017).  A multi-resolution approximation for massive spatial datasets.  Journal of the American Statistical Association, 112(517):201–214.
  * Katzfuss and Gong, (2020) Katzfuss, M. and Gong, W. (2020).  A class of multi-resolution approximations for large spatial datasets.  Statistica Sinica, 30(4):2203–2226.
  * Katzfuss and Guinness, (2021) Katzfuss, M. and Guinness, J. (2021).  A general framework for Vecchia approximations of Gaussian processes.  Statistical Science, 36(1):124–141.
  * O’Carroll et al., (2019) O’Carroll, A. G., Armstrong, E. M., Beggs, H. M., Bouali, M., Casey, K. S., Corlett, G. K., Dash, P., Donlon, C. J., Gentemann, C. L., Høyer, J. L., et al. (2019).  Observational needs of sea surface temperature.  Frontiers in Marine Science, 6:420.
  * Paciorek and Schervish, (2006) Paciorek, C. and Schervish, M. (2006).  Spatial modelling using a new class of nonstationary covariance functions.  Environmetrics, 17(5):483–506.
  * R Core Team, (2020) R Core Team (2020).  R: A Language and Environment for Statistical Computing.  R Foundation for Statistical Computing, Vienna, Austria.
  * Sang and Huang, (2012) Sang, H. and Huang, J. Z. (2012).  A full scale approximation of covariance functions.  Journal of the Royal Statistical Society, Series B, 74(1):111–132.
  * Shirota et al., (2019) Shirota, S., Finley, A. O., Cook, B. D., and Banerjee, S. (2019).  Conjugate nearest neighbor gaussian process models for efficient statistical interpolation of large spatial data.  arXiv preprint arXiv:1907.10109.
  * Vazquez-Cuervo et al., (2017) Vazquez-Cuervo, J., Torres, B., Dewitte, B., Menemenlis, D., Chin, T. M., and Armstrong, E. M. (2017).  Relationship between SST gradients and upwelling off Peru and Chile: Model/data analysis.  International Journal of Remote Sensing, 38(23):6599–6622.
  * Vecchia, (1988) Vecchia, A. (1988).  Estimation and model identification for continuous spatial processes.  Journal of the Royal Statistical Society, Series B, 50(2):297–312.
  * Zammit-Mangion and Rougier, (2020) Zammit-Mangion, A. and Rougier, J. (2020).  Multi-scale process modelling and distributed computation for spatial data.  Statistics and Computing, 30(6):1609–1627.

Supplementary Material

Figure S1: Global nighttime SST for the first, fourth, and seventh days of
January, April, July and October 2019. Dates are presented in YYYY-MM-DD
format. Rows correspond to months and columns to days. The number of all
observations available for each day ($n$) is presented in the titles. Scale
ranges from 275 to 305 Kelvin. Figure S2: Cross validation error versus number
of basis functions (see Section 3.2). We select 11 cubic smoothing spline
functions for their relatively low degrees of freedom and small error. Figure
S3: Boxplots of the nighttime SST residuals using 11 cubic smoothing splines
aggregated into 180 latitudinal bins (see Section 3.2). Bin whiskers are
extended to the range of observations. The blue vertical dashed line is placed
at the bin containing the equator and the horizontal magenta dashed line is
placed at zero. The majority of residuals center at zero, and we observe
different variability of residuals at different latitudes, indicating a
nonstationary covariance model for the residuals is needed. Figure S4: Local
estimates of the smoothness parameter from April 1 to 7, 2019. Most of them
are relatively stable around 0.5, supporting the assumption of exponential
covariance being used (see Section 3.4). Figure S5: Locally estimated
parameters of the exponential covariance function for the nighttime SST
residuals from April 1 to 7, 2019. Consistent patterns are observed across the
week, motivating us to use all of them to obtain smoothed estimates by the
Wendland basis functions (see Section 3.4).