Consider a bivariate spatial process Y(s)≡(Y1(s),Y2(s))′, s∈G, where we refer to G⊂R2 as the geographic domain. The cross-covariance functions are Cij,G(s,u)=cov(Yi(s),Yj(u));i,j=1,2;s,u∈G, where {Cij,G(⋅,⋅)} are nonstationary and asymmetric (cross-)covariance functions. Under warping, they can expressed as simpler, usually stationary and symmetric (cross-)covariance functions, {Cij,D(⋅,⋅)}, on a deformed space D.
First, consider two warping functions f1(⋅),f2(⋅) such that f1,f2:G→D. In this case, Cij,G(s,u)=Cij,D(fi(s),fj(u))=Coij,D(fi(s)−fj(u));i,j=1,2;s,u∈G, where {Coij,D(⋅)} are stationary and symmetric (cross-)covariance functions on the warped domain D.
In practice, any asymmetry present is likely to be simple and dominated by global shifts and rotations. To model asymmetry, we thus propose expressing each fi(⋅),i=1,2, as a composition of a shared warping function f(⋅), and an aligning function gi(⋅). That is, for i,j=1,2, we let Cij,G(s,u)=Cij,D(f∘gi(s),f∘gj(u))=Coij,D(f∘gi(s)−f∘gj(u));s,u∈G. We can choose the aligning functions gi(⋅),i=1,2, to be simple transformations that are commonly used to align spatial fields and which can include translations and rotations [@wiens2019surface], and the shared warping functions f(⋅) to be composed of deep injective warping functions, which are able to introduce complex nonstationary behaviour (Zammit-Mangion et al. 2019).