  Rapid developments in satellite remote-sensing technology have enabled the
collection of geospatial data on a global scale, hence increasing the need for
covariance functions that can capture spatial dependence on spherical domains.
We propose a general method of constructing nonstationary, locally anisotropic
covariance functions on the sphere based on covariance functions in R^3. We
also provide theorems that specify the conditions under which the resulting
correlation function is isotropic or axially symmetric. For large datasets on
the sphere commonly seen in modern applications, the Vecchia approximation is
used to achieve higher scalability on statistical inference. The importance of
flexible covariance structures is demonstrated numerically using simulated data
and a precipitation dataset.
