An anisotropic model for global climate data

TL;DR

This paper introduces a simple method to create axially symmetric Gaussian processes on the sphere, effectively capturing the directional anisotropy observed in global climate data for improved geostatistical modeling.

Contribution

It proposes a novel, elementary approach to model anisotropic effects in global climate data using Gaussian processes on the sphere.

Findings

Effective modeling of directional anisotropy in climate data

Simplified approach for Gaussian process construction on the sphere

Potential for improved geostatistical analysis

Abstract

We present a new, elementary way to obtain axially symmetric Gaussian processes on the sphere, in order to accommodate for the directional anisotropy of global climate data in geostatistical analysis.

Equations

$$\hat{Y}(x^{\ast })=\mathbb{E}(Y(x^{\ast })|Y(x_1)=y_1,\cdots ,Y(x_n)=y_n).$$

$$Var(Y(x^{\ast })|Y(x_1)=y_1,\cdots ,Y(x_n)=y_n),$$

$$\forall {\lambda }_1,\cdots ,{\lambda }_n\in \mathbb{R},\forall x_1\cdots x_n\in X,\sum _{i,j=1}^n{\lambda }_i{\lambda }_{j}F(x_i,x_j)\ge 0.$$

$$d(x,y)={\cos }^{-1}(\sin {\phi }_1\cdot \sin {\phi }_2+\cos {\phi }_1\cdot \cos {\phi }_2\cdot \cos ({\theta }_2-{\theta }_1)).$$

$$K_X(x,y)=F({\theta }_x-{\theta }_y\mathrm{mod}2\pi ,{\phi }_x,{\phi }_y).$$

$$F({\theta }_x-{\theta }_y\mathrm{mod}2\pi ,{\phi }_x,{\phi }_y)=F({\theta }_x-{\theta }_y\mathrm{mod}2\pi ,{\phi }_y,{\phi }_x).$$

$$K(x,y)=\sigma e^{-\frac{|x_1-y_1|}{r_1}}\cdot e^{-\frac{|x_2-y_2|}{r_2}}$$

$$K(x,y)=\sigma e^{-\frac{|{\phi }_x-{\phi }_y|}{r_{\phi }}}\cdot e^{-\frac{|{\theta }_x-{\theta }_y|}{r_{\theta }}},$$

$$K(x,y)=\sigma e^{-\frac{d(x,y)}{r_{iso}}}\cdot e^{-\frac{|{\phi }_x-{\phi }_y|}{r_{\phi }}}.$$

$$K(x,y)=K_{iso}(x,y)\cdot K_{\phi }({\phi }_x,{\phi }_y)$$

Full text

The following article has been accepted to the Statistical models for climate data invited specialized session at the 2019 conference of the Italian Statistics Society ”SIS2019”.

1 Introduction

Following the increasing amount of measurements and computing power available, geostatistics has evolved in the last decades to consider massive and global scale datasets. This adaptation presents new challenges, such as taking into account the spherical nature of Earth, and the considerable non-stationarity of global datasets, while maintaining computing feasibility.

We focus in this article on Gaussian process based geostatistics for climate data. Isotropic models on the sphere, which assume that the data has the same statistical behavior in every direction, are well understood and the literature gives rich classes of covariances (see for example Gneiting_2013 and the reviews Jeong_al_2017 and Porcu_al_2018 ).

However, climate data exhibits various nonstationary effects at the global scale. In particular, it is clear that the variability of climate variables is higher in the north-south direction, than in the east-west. Jones Jones_1963 notices this issue and defines axially symmetric Gaussian processes on the sphere, which are stationary in longitude only. He also characterizes their decomposition on the spherical harmonic basis, which Stein Stein_2007 truncates to a finite order to carry on a statistical analysis of a Total Ozone dataset at a global scale. While very flexible, Stein’s model involves the estimation of a large number of coefficients, and still doesn’t fit the local behaviour of the Ozone dataset. This calls for explicit axially symmetric covariances that relies on a small number of parameters. Such covariance functions are obtained through partial differentiation of isotropic processes by Jun and Stein Jun_Stein_2007 ; Jun_Stein_2008 , and further studied by Hitczenko and Stein Hitczenko_Stein_2012 . Recently Porcu et al. Porcu_al_2019 proposed to modify variograms of isotropic covariances to obtain axially symmetric analogues, while Huang et al. Huang_al_2012 consider products of separated covariances on latitudes and longitudes.

Our approach consists in considering the product of an isotropic covariance with an additional covariance on latitudes. Compared to Huang et al.,Huang_al_2012 , our method has the advantage of giving Gaussian processes that are continuous over the whole sphere, including the poles.

In Section 2, we give a description of two climatic datasets, which motivate our work and allow performance evaluation of models. Section 3 recalls some generalities on Gaussian processes on the sphere and introduce our model. In Section 4, we describe our ongoing work to compare model performances and properties to the existing literature.

2 Climatic datasets

In this section, we describe two important sources of climatic data, which motivates our study.

The radiosonde observations database from NOAA333US National Oceanic and Atmospheric Administration ESRL444Earth System Research Laboratory gathers worldwide radiosonde observations, that are considered as an important atmospheric observation standard. We propose to use it to assess model performances on real data.

We also consider the reanalysis ERA-Interim archive from ECMWF555European Centre for Medium-Range Weather Forecasts, that does not present the drawbacks of real data (missing data, inhomogeneous spatial covering) but rather constitutes a resource whose coherence and completeness allows to extract tailor-fitted datasets to assess a virtually unlimited variety of hypothesis.

2.1 Radiosonde observations

As mentioned above, the NOAA ESRL database gathers worldwide radiosonde observations. There are around 900 ground stations performing radiosonde observations in the world. Ascending balloons equipped with measurement devices are launched, radio-transmitting the observations of meteorological variables about every few seconds to the ground.

These measurements are stored at a restricted number of altitudes levels, corresponding to pressure levels: there are 22 mandatory levels (from the ground up to 1 hPa) common to every sounding, to which are added a varying number of significant levels (28 on the average), that depend on the launch and correspond to important variations in measured temperature or dew point depression. Missing data at the highest mandatory levels are frequent as the balloon explodes between 30 and 35km of altitude.

There are typically two launches per day in a station, occurring at times close to 00 and 12 UTC, but additional and missing launches occur. Data is available publicly and can be downloaded from the NOAA ESRL servers.

2.2 The ERA-Interim reanalysis

Produced by the ECMWF, the ERA-Interim archive is a reanalysis of global atmospheric data from 1979, updated in real-time to nowadays.

This model output is built from the input of numerous meteorological datasets of various nature (mostly satellites measurements, but also radiosondes, land, boat, planes,… measurements), which are assimilated and extrapolated following a 4-dimensional variational analysis.

Era-Interim provides values for a considerable number of meteorological variables. Atmospheric variables are given at a vertical resolution of sixty pressure levels (from the ground to 0.1 hPa) on the horizontal reduced Gaussian grid N128, which is regular in latitude but has a variable longitudinal resolution to maintain a ground resolution of approximately 79km (see details in ERA_grid ). Time resolution is 6h. Surface variables are given with the same horizontal display and a time resolution of 3 hours. In both cases it is also possible to download data which have been extrapolated on regular latitude/longitude grids at various resolutions.

Advantages of this archive consist in its completeness, its coherence in the methodology and in the physical consistency of the data. From this, it is possible to tailor-fit datasets to a virtually unlimited number of applications.

Specifications of the archive are given in ERAreport , while the reanalysis process is described in ERAarticle . Data is available publicly and can be downloaded from the ECMWF servers.

3 Spatial Modeling

We briefly recall generalities about Gaussian process based geostatistics (for which we refer to Cressie Cressie_book_1992 for more details) and Gaussian processes on the sphere before to introduce our new model.

3.1 Gaussian process geostatistics

The classical approach in geostatistics, given observed data $y_{1},\cdots,y_{n}\in\mathbb{R}$ at locations $x_{1},\cdots,x_{n}$ in some space $S$ consists in considering a Gaussian process $(Y_{x})_{x\in S}$ indexed by $S$, and estimating the value of the variable of interest at a new location $x^{*}$ by the Kriging estimator

[TABLE]

The quality of this prediction critically relies on the choice of the Gaussian process $Y$, whose statistical behavior should fit as much as possible the data of interest. The classical approach consists in selecting a Gaussian process from a parametric class $Y_{\theta}$, where $\theta$ is a vector of parameters, which is typically done by maximizing the likelihood. The key asset of this approach is that prediction comes with a statistical model of the data, which allows in particular the computation of the conditional variance

[TABLE]

that quantifies the prediction uncertainty at the location $x^{*}$.

Due to the $O(n^{3})$ numerical complexity of both Kriging prediction and likelihood evaluation, this approach proves itself computationally heavy – if not unfeasible – on massive datasets. There exists an important literature on adaptations to large datasets. See for example Heaton_al_2018 for a comparison of methods a single case study.

However in the majority of methods, the primary ingredient remain Gaussian processes, and finding classes of Gaussian processes that suit practical applications constitute an entire field of research.

Let us recall that the statistical properties of a Gaussian process are entirely characterized by its expectation and covariance functions. As the choice of a a mean function is free from any constrains, in this article all Gaussian fields will be assumed centered. In this setting, considering a Gaussian process is equivalent to considering a covariance function.

In contrast with the mean function, not every function is admissible as a covariance. Let us recall that a function $F:X\times X\rightarrow\mathbb{R}$ is the covariance of a Gaussian process $(X_{x})_{x\in X}$ if and only if it is symmetric and positive definite, that is

[TABLE]

In this case we will often say that $F$ is a valid covariance.

The finding of classes of covariance functions that give Gaussian properties which fit data from application is the first, necessary step of Gaussian process modeling.

3.2 Gaussian processes on the sphere

In order to model variables on the Earth, we consider Gaussian processes indexed by the sphere $\mathbb{S}^{2}$, which we will take of radius $1$. To a given point $x$ on the sphere are associated its longitude $\theta_{x}\in[-\pi,\pi]$ and latitude $\varphi_{x}\in[-\pi/2,\pi/2]$, in radians. Observe that $\theta_{x}$ and $\varphi_{x}$ are well defined for every $x$ distinct from the north or south poles, for which $\varphi_{x}$ is $\pi/2$ (resp. $-\pi/2$) but $\theta_{x}$ can take any value. This singularity of the spherical coordinates has some practical consequences when it comes to define a covariance on the whole sphere, as we will see in Section 3.3.

Given a Gaussian random field $X_{x\in\mathbb{S}^{2}}$ and $x,y\in\mathbb{S}^{2}$, we will write $K_{X}(x,y)$ or $K_{X}(\theta_{x},\theta_{y},\varphi_{x},\varphi_{y})$ for its covariance function.

As in the Euclidean case, a simplifying assumption is to consider stationary models.A A Gaussian process on the sphere is called isotropic if $K_{X}(x,y)=F(d(x,y))$, where $d(x,y)$ is the great circle distance on the sphere, given by the expression:

[TABLE]

We refer to GneitingGneiting_2013 , Porcu et al. Porcu_al_2018 , Jeong et al. Jeong_al_2017 and references within for a state of the art on isotropic Gaussian processes.

However at the global scale, climatic data exhibits important nonstationarity effects: in particular, it is clear that climatic variables are typically correlated at a shorter range in the latitudinal direction than in the longitudinal. Jones Jones_1963 proposes to address this issue and defines axially symmetric Gaussian processes, which are stationary only in longitude variable, that is to say their covariance verifies

[TABLE]

Furthermore an axially symmetric Gaussian process is said to be latitudinally reversible if

[TABLE]

3.3 A new class of axially symmetric covariances

In the Euclidean case, a natural idea to obtain anisotropic covariances is to consider the product of covariances with respect to different coordinates. For example considering for $x,y\in\mathbb{R}^{2},$

[TABLE]

with distinct scale parameters $r_{1}$ and $r_{2}$ yields a model that is correlated at a longer range in a chosen direction. On the sphere one might consider in a similar way products of covariances in latitude and longitude, such as for $x,y\in\mathbb{S}^{2},$

[TABLE]

but in this case, since the longitude coordinate can take any value at the poles of the sphere (see Section 3.2), the covariance is not well defined when $x$ or $y$ is one of the two poles. Moreover, there is no way to chose a value for $K$ at these points to obtain a continuous kernel, since when $x$ (or $y$) closes to one of the pole, the value taken by $K(x,y)$ depend on the way $x$ approaches the pole.

Discontinuity of the covariance kernel yields Gaussian fields that are discontinuous in $L^{2}$ and whose realizations are not almost surely continuous. This is unwanted in most applications where variables exhibits a continuous behavior over the sphere, such as climatic applications. In the case of (8) we obtain an expected singular behavior at the poles.

To address this issue, we consider instead of (8) a covariance

[TABLE]

Notice that since $(x,y)\mapsto d(x,y)$ and $x\mapsto\varphi_{x}$ are continuous function on the sphere, this time $K$ is continuous everywhere. The term $e^{-\frac{|\varphi_{x}-\varphi_{y}|}{r_{\varphi}}}$ can be seen as an additional decorrelation to the isotropic covariance $\sigma e^{-\frac{d(x,y)}{r_{iso}}}$ in the latitudinal direction.

Extending this idea to kernels other than the exponential, we obtain the following class of covariances.

Theorem 3.1

Let $K_{iso}$ be an isotropic covariance on the sphere and $K_{\varphi}$ be a covariance on $[-\pi,\pi]$.

The kernel defined by

[TABLE]

is a latitudinally reversible, axially symmetric covariance on the sphere.

Furthermore, if $K_{iso}$ and $K_{\varphi}$ are continuous, $K$ is continuous on the whole sphere, and as such, a Gaussian field with covariance $K$ is continuous in $L^{2}$ sense, and has almost surely continuous trajectories.

Proof

As product of two valid covariances, $K$ is a valid covariance (see Schur product theorem in Zhang Zhang_book_2005 ). Axial symmetry of $K$ is a direct consequence of the isotropy of $K_{iso}$ and the fact that $K_{\varphi}$ does not depend on $\theta_{x}$ and $\theta_{y}$.

In contrast with Jun and Stein Jun_Stein_2007 ; Jun_Stein_2008 , whose approach seems to give richer classes of covariances, allowing for complex interactions between variables, our method permits to easily modify the classical isotropic covariances to obtain axially symmetric analogues. Porcu et al. Porcu_al_2019 attain a similar goal, but focusing on the variograms. In comparison with Huang et al. Huang_al_2012 , our covariances present the advantage to be continuous on the whole sphere.

4 Perspectives: performances and properties of axially symmetric models

In an ongoing work, we assess the performances of our new class of axially symmetric covariances, compared to the propositions of Jun and Stein Jun_Stein_2007 ; Jun_Stein_2008 , Huang Huang_al_2012 , Porcu et al. Porcu_al_2019 . We take advantage of the two data sources that were presented in Section 2 to assess their performances in prediction, training on real data and assessing prediction with a global covering.

We expect that the simple expression (10) of our covariances allows us to study their differentiability at the origin in every direction, allowing for Gaussian processes with homogeneous trajectory regularity, or at the contrary for models which exhibit a regularity that depends on the direction of the trajectory.