1 Introduction
Geostatistical methods for spatial and spatiotemporal data are in great demand from fields such as earth and climate sciences, epidemiology and agriculture. A comprehensive overview can be found in cressie2011statistics. Spatially stationary processes are commonly used as models in geostatistical applications, but often the assumption of stationarity and isotropic covariance functions are difficult to hold in real applications; see for example, guttorp1994space, le2006statistical, damian2001bayesian. A recent review of methods that allow nonhomogenous covariance models is schmidt2020flexible.
We propose a semiparametric method of nonstationary spatial covariance function that expands upon the work of sampson1992nonparametric. Data observed in a set
2 Spatial Deformation Model (SDM)
Let be spatial locations with for all in a geostatistical domain and with for all Moreover, let such that for all and are bijective, differentiable functions, for
Let be a Gaussian random field with spatial covariance function For any pair of spatial sites we have with given by
where and is a stationary, isotropic correlation function with parameters We remark that this model may not be identifiable, since is invariant for shifts and rigid rotations in Moreover, scaling by a constant will produce the same if is also scaled by
Nevertheless the model can be used for Kriging and spatial interpolation without issues.
For example, consider and given by the swirl transformation shown in Figure 1. We sampled from a spatiotemporal process on a regular grid of points, with and separable covariance function
where and A realization of a random field on a fine mesh grid is also shown in Figure 1, obtained with the RandomFields package (schlather2015analysis).
3 Spatial Deformation Estimation
Let
then is the semivariogram of the random field . In particular, for a stationary, isotropic
where is a positive nondecreasing function of . Note valid functions are conditionally negative definite (cressie1993statistics).
If there are ways to obtain a sample covariance matrix, such as temporal replicates, then we can construct a sample variogram as
If the variogram is isotropic on a set of artificial coordinates , then it must be have form
where Therefore, the variogram entries can be seen as a dispersion metric, with
3.1 Nonmetric multidimensional scaling
sampson1992nonparametric first step is to consider a nonmetric multidimensional scaling approach (hereafter, nMDS; see kruskal1964multidimensional; mardia1979multivariate; cox2000multidimensional). From a set of dispersions we seek a monotone transformation such that
and
and an artificial set of coordinates with interpoint Euclidean distances that minimizes the stress criterion
Here we need to impose constraints on such that the function is conditionally positive definite. kruskal1964multidimensional only finds an isotonic transformation, so sampson1992nonparametric adapted the optimization algorithm as follows:

Set k = 1. Find an initial configuration using, for example, Kruskal’s isotonic nMDS.

Fit a variogram model to the data using the coordinates obtaining

Obtain via MDS performed on

Set and return to step 2 until convergence.
3.2 Spline smoothing
The second step in sampson1992nonparametric approach is to find an approximation for such that This can be done if we approximate the functions with, for example, thinplate splines (wahba1980some). Thinplate splines are the minimizers of the variational problem
The solution for a fixed is
where and are estimating by plugin.
3.3 Tensor product of Bsplines
We propose a different approach for the spline smoothing by making use of low rank approximation (finite approximation) for the deformation functions. Specifically, let and assume we have a collection of points be a sequence of points in . Define a map such that
(1) 
where the functions , form a partition of unity.
Particularly, suppose that are the Bsplines basis functions and there is a integer number such that a map can be represented by
(2) 
where and Let for some index set such that then the functions are well approximated by
where are fixed positive integers and are Bspline basis functions (see, e.g., ramsay2005functional).
To guarantee that the functions do not fold onto themselves, we must guarantee that is locally invertible and differentiable. In fact, a diffeomorphism is desirable (see, e.g., perrin1999modelling). A necessary condition is that everywhere in the domain,
(3) 
but in practice such constraint is difficult to implement, requiring the evaluation of the Jacobian for every pair We will show in Section 4 that such condition can be translated into a condtion onto the coefficients for Note however that it is enough to ensure since a change in signs would imply a discontinuous Jacobian.
3.4 Simultaneous estimation of covariance function and deformation
The covariance parameters
as well as the mean vector
are estimated by maximizing the profile loglikelihood. Without loss of generality set Thus our proposed Bspline approach can be estimated as
Given solve the optimization problem
where corresponds to the optimization target seen as a function of only, with parameters fixed.

Given obtain such that
s.t. where are the row vector of basis functions and is the Kronecker product, and for some choice of such that if

Repeat steps 1 and 2 until convergence.
We remark that the matrices are sparse when using Bsplines, therefore solutions to step 2 are scalable. To evaluate the inequality consider the following: write
or simply where is a vector of Bspline basis functions evaluated at , similarly is a vector of Bspline basis evaluated at , and is a matrix of spline coefficients for the tensor product approximation of the th function . Similarly,
is the partial derivative of with respect to and if we write then This allows us to see that
so the determinant of the Jacobian as a function of is an inner product of
, weighted by the skewSymmetric matrix
Ensuring the inequality for all values of remains a difficult task, but we can chose a set of basis for which does not depend on4 Constrained spatial deformation estimation
We will use an approach similar to musse2001topology. Consider Bsplines of degree 1 on . Assume that there are equally spaced inner knots, where . Since the knots are equally spaced, they can be written as where In this case, the Bspline bases are given by
with derivatives
Consider where and is between and (where and ). Then there are only 4 bases that evaluate to nonzero values, indexed by , , and , so
and therefore, looking only at the nonzero pairs , we have
where
note is therefore proportional to
The above equations describe a plane in with coefficients depending on . Now, since and , where (and similarly for ), we have four restrictions to consider:

When and ,

When and ,

When and ,

When and ,
This collection of constraints, for and together imply in a nonfolding deformation map, and can be enforced with constrained optimization routines. We have employed svanberg2002class constrained optimization algorithm, avaliable in the nloptr package (johnson2020nloptr). The code is available as an R package in https://github.com/guiludwig/bsplinedef.
5 Simulation study
To evaluate the performance of the algorithm, we conducted a simulation study based on the swirl function shown in Figure 1. Each sample has spatial points on a regular grid in the geographical domain We simulated from a Gaussian random field with mean function and covariance function with parameters and First, consider a single realization of the spatial random field. The estimated deformation maps are shown in Figure 2. We have obtained the constrained Bspline deformations (hereafter, bdef) with and basis functions. The estimated deformation function does not fold, even though the true deformation function is difficult to be recovered. The case when performs better than or
indicating that there are features in the deformation map that need a large number of degrees of freedom to be estimated. On the other hand, the functions using
sampson1992nonparametric (SG) method have smoothing parameters and set to match the bdef approach. They start showing folding at and cannot recover deformation maps that require a number of degrees of freedom larger than Note that the SG maps were stretched or shrunk to fit the plot area, but the bdef maps did not require this step.The comparison of estimated covariance matrices allow us to overlook the identifiability issues with rotations, shifts and scaling of the estimated maps. In Figure 3 we show a scatterplot of the upperdiagonal entries of the estimated covariance matrices for the data, versus the true covariance matrices. We remark that the bdef method shows no apparent bias and becomes more accurate as the number of degrees of freedom increase. On the other hand, SG has good performance at a small number of degrees of freedom (more smoothing), but tends to overestimate matrix entries as the number of degrees of freedom for the thinplate spline increases.
6 Case study: Rainfall data in southeastern Brazil
The dataset we use to illustrate our method comes from meteorological surveys conducted by INMET – Instituto Nacional de Meteorologia, Brazil. The measurements are made at every 15 minutes, and daily accumulated values are made available. Following rozante2010combining, we grouped the rainfall data in periods of 10 days.
Since we seek temporally stationary data, we decided to restrain the data collection to 20180101 to 20180330 (rainfall season), for a total of 9 time periods (period 1: 20180101 to 20180110, …, period 9: 20180322 to 20180331).
We selected the 50 stations of the southeastern region of Brazil that had complete observations available during the period. The stations are shown in the left panel of Figure 4. Data can be obtained at http://www.inmet.gov.br/portal/index.php?r=bdmep/bdmep.
In the right panel of Figure 4, we show the estimated deformation map, using Bspline basis functions. The estimated deformation map reveals topographical features shown in the left panel of the same Figure. Weather stations located in northwestern flat lands of Minas Gerais are treated as closer to each other than stations near the rough coast of Rio de Janeiro and São Paulo states, where elevation changes are abrupt.
In Figure 5 we perform conditional simulation of the Gaussian random field (Kriging) for the 10day periods starting in 2018–01–01 and 2018–01–11, showing the resulting prediction maps.
Acknowledgments
Ronaldo Dias was supported by the FAPESP foundation, grant 2018/046549. Guilherme Ludwig was supported by FAPESP grants 2018/046549 and 2019/035170, as well as FAEPEX grant 2253/19.
References
Appendix A Derivation of the nonfolding constraint
We can derive the constraint presented in Section 4 as follows: Observe that
or
Further simplifying leads to
Comments
There are no comments yet.