#!/usr/bin/env python
u"""
sph_radial_basis.py
Written by Tyler Sutterley (05/2022)
Interpolates data over a sphere using radial basis functions
with QR factorization option to eliminate ill-conditioning
CALLING SEQUENCE:
output = sph_radial_basis(lon, lat, data, longitude, latitude,
smooth=smooth, epsilon=epsilon, method='inverse')
INPUTS:
lon: input longitude
lat: input latitude
data: input data
longitude: output longitude
latitude: output latitude
OUTPUTS:
output: interpolated data
OPTIONS:
smooth: smoothing weights
epsilon: adjustable constant for distance functions
default is the mean Euclidean distance
method: radial basis function (** has option for QR factorization method)
multiquadric**
inverse_multiquadric** or inverse** (default)
inverse_quadratic**
gaussian**
linear
cubic
quintic
thin_plate: thin-plate spline
QR: use QR factorization algorithm of Fornberg (2007)
norm: distance function for radial basis functions (if not using QR)
euclidean: Euclidean Distance with distance_matrix (default)
GCD: Great-Circle Distance using n-vectors with angle_matrix
PYTHON DEPENDENCIES:
numpy: Scientific Computing Tools For Python
https://numpy.org
scipy: Scientific Tools for Python
https://docs.scipy.org/doc/
REFERENCES:
B Fornberg and C Piret, "A stable algorithm for flat radial basis functions
on a sphere." SIAM J. Sci. Comput. 30(1), 60-80 (2007)
B Fornberg, E Larsson, and N Flyer, "Stable Computations with Gaussian
Radial Basis Functions." SIAM J. Sci. Comput. 33(2), 869-892 (2011)
UPDATE HISTORY:
Updated 05/2022: updated docstrings to numpy documentation format
Updated 01/2022: added function docstrings
Updated 02/2019: compatibility updates for python3
Updated 09/2017: using rcond=-1 in numpy least-squares algorithms
Updated 08/2016: finished QR factorization method, added norm option
Forked 08/2016 from radial_basis.py for use over a sphere
Updated 08/2016: using format text within ValueError, edit constant vector
removed 3 dimensional option of radial basis (spherical)
changed hierarchical_radial_basis to compact_radial_basis using
compactly-supported radial basis functions and sparse matrices
added low-order polynomial option (previously used default constant)
Updated 01/2016: new hierarchical_radial_basis function
that first reduces to points within distance. added cutoff option
Updated 10/2014: added third dimension (spherical)
Written 08/2014
"""
from __future__ import print_function, division
import numpy as np
import scipy.special
from spatial_interpolators.legendre import legendre
[docs]def sph_radial_basis(lon, lat, data, longitude, latitude, smooth=0.,
epsilon=None, method='inverse', QR=False, norm='euclidean'):
"""
Interpolates a sparse grid over a sphere using radial basis
functions with QR factorization option
Parameters
----------
lon: float
input longitude
lat: float
input latitude
data: float
input data
longitude: float
output longitude
latitude: float
output latitude
smooth: float, default 0.0
smoothing weights
epsilon: float or NoneType, default None
adjustable constant for distance functions
method: str, default 'inverse'
compactly supported radial basis function
* ``'multiquadric'`` [#f1]_
* ``'inverse_multiquadric'`` [#f1]_ or ``'inverse'`` [#f1]_
* ``'inverse_quadratic'`` [#f1]_
* ``'gaussian'`` [#f1]_
* ``'linear'``
* ``'cubic'``
* ``'quintic'``
* ``'thin_plate'``
QR: bool, default False
use QR factorization algorithm of [Fornsberg2007]_
norm: str, default 'euclidean'
Distance function for radial basis functions
* ``'euclidean'``: Euclidean Distance with distance_matrix
* ``'GCD'``: Great-Circle Distance using n-vectors with angle_matrix
Returns
-------
output: float
interpolated data grid
References
----------
.. [Fornsberg2007] B. Fornberg and C. Piret,
"A stable algorithm for flat radial basis functions on a sphere,"
*SIAM Journal on Scientific Computing*, 30(1), 60--80, (2007).
`doi: 10.1137/060671991 <https://doi.org/10.1137/060671991>`_
.. [Fornsberg2011] B. Fornberg, E. Larsson, and N. Flyer,
"Stable Computations with Gaussian Radial Basis Functions,"
*SIAM Journal on Scientific Computing*, 33(2), 869--892, (2011).
`doi: 10.1137/09076756X <https://doi.org/10.1137/09076756X>`_
.. [#f1] has option for QR factorization method
"""
#-- remove singleton dimensions
lon = np.squeeze(lon)
lat = np.squeeze(lat)
data = np.squeeze(data)
longitude = np.squeeze(longitude)
latitude = np.squeeze(latitude)
#-- size of new matrix
if (np.ndim(longitude) > 1):
nlon,nlat = np.shape(longitude)
sz = np.int(nlon*nlat)
else:
sz = len(longitude)
#-- Check to make sure sizes of input arguments are correct and consistent
if (len(data) != len(lon)) | (len(data) != len(lat)):
raise Exception('Length of Longitude, Latitude, and Data must be equal')
if (np.shape(longitude) != np.shape(latitude)):
raise Exception('Size of output Longitude and Latitude must be equal')
#-- create python dictionary of radial basis function formulas
radial_basis_functions = {}
radial_basis_functions['multiquadric'] = multiquadric
radial_basis_functions['inverse_multiquadric'] = inverse_multiquadric
radial_basis_functions['inverse'] = inverse_multiquadric
radial_basis_functions['inverse_quadratic'] = inverse_quadratic
radial_basis_functions['gaussian'] = gaussian
radial_basis_functions['linear'] = linear
radial_basis_functions['cubic'] = cubic
radial_basis_functions['quintic'] = quintic
radial_basis_functions['thin_plate'] = thin_plate
#-- create python dictionary of radial basis function expansions
radial_expansions = {}
radial_expansions['multiquadric'] = multiquadratic_expansion
radial_expansions['inverse_multiquadric'] = inverse_multiquadric_expansion
radial_expansions['inverse'] = inverse_multiquadric_expansion
radial_expansions['inverse_quadratic'] = inverse_quadratic_expansion
radial_expansions['gaussian'] = gaussian_expansion
#-- check if formula name is listed
if method in radial_basis_functions.keys():
RBF = radial_basis_functions[method]
else:
raise ValueError("Method {0} not implemented".format(method))
#-- check if formula name is valid for QR factorization method
if QR and (method in radial_expansions.keys()):
expansion = radial_expansions[method]
elif QR and (method not in radial_expansions.keys()):
raise ValueError("{0} expansion not available with QR".format(method))
#-- create python dictionary of distance functions (if not using QR)
norm_functions = {}
norm_functions['euclidean'] = distance_matrix
norm_functions['GCD'] = angle_matrix
if norm in norm_functions:
norm_matrix = norm_functions[norm]
else:
raise ValueError("Distance Function {0} not implemented".format(norm))
#-- convert input lat and lon into cartesian X,Y,Z over unit sphere
phi = np.pi*lon/180.0
th = np.pi*(90.0 - lat)/180.0
xs = np.sin(th)*np.cos(phi)
ys = np.sin(th)*np.sin(phi)
zs = np.cos(th)
#-- convert output longitude and latitude into cartesian X,Y,Z
PHI = np.pi*longitude.flatten()/180.0
THETA = np.pi*(90.0 - latitude.flatten())/180.0
XI = np.sin(THETA)*np.cos(PHI)
YI = np.sin(THETA)*np.sin(PHI)
ZI = np.cos(THETA)
#-- Creation of data distance matrix (Euclidean or Great-Circle Distance)
#-- Data to Data
Rd = norm_matrix(np.array([xs, ys, zs]),np.array([xs, ys, zs]))
N,M = np.shape(Rd)
#-- if epsilon is not specified
if epsilon is None:
#-- calculate norm with mean distance
uix,uiy = np.nonzero(np.tri(N,M=M,k=-1))
epsilon = np.mean(Rd[uix,uiy])
#-- QR factorization algorithm of Fornberg (2007)
if QR:
#-- calculate radial basis functions using spherical harmonics
R,w = RBF_QR(th,phi,epsilon,data,expansion)
n_harm = np.sqrt(np.shape(R)[0]).astype(np.int)
#-- counter variable for filling spherical harmonic matrix
index = 0
#-- evaluation matrix E
E = np.zeros((sz,np.int(n_harm**2)))
for l in range(0,n_harm):
#-- Each loop adds a block of columns of degree l to E
E[:,index:2*l+index+1] = spherical_harmonic_matrix(l,THETA,PHI)
index += 2*l + 1
#-- calculate output interpolated array (or matrix)
output = np.dot(E,np.dot(R,w))
else:
#-- Calculation of the PHI Matrix with smoothing
PHI = np.zeros((N+1,M+1))
PHI[:N,:M] = RBF(epsilon, Rd) + np.eye(N,M=M)*smooth
#-- Augmentation of the PHI Matrix with a Constant Vector
PHI[:N,M] = np.ones((N))
PHI[N,:M] = np.ones((M))
#-- Computation of the Weights
DMAT = np.concatenate(([data,[0]]),axis=0)
w = np.linalg.lstsq(PHI,DMAT[:,np.newaxis],rcond=-1)[0]
#-- Computation of distance Matrix (Euclidean or Great-Circle Distance)
#-- Data to Mesh Points
Re = norm_matrix(np.array([XI,YI,ZI]),np.array([xs,ys,zs]))
#-- calculate radial basis function for data-to-mesh matrix
E = RBF(epsilon,Re)
#-- Augmentation of the Evaluation Matrix with a Constant Vector
P = np.ones((sz,1))
E = np.concatenate(([E, P]),axis=1)
#-- calculate output interpolated array (or matrix)
output = np.dot(E,w)
#-- reshape output to original dimensions and return
if (np.ndim(longitude) == 1):
return np.squeeze(output)
else:
return output.reshape(nlon,nlat)
#-- define radial basis function formulas
def multiquadric(epsilon, r):
#-- multiquadratic
f = np.sqrt((epsilon*r)**2 + 1.0)
return f
def multiquadratic_expansion(epsilon, mu):
c = -2.*np.pi*(2.*epsilon**2+1.+(mu+1./2.)*np.sqrt(1.+4.*epsilon**2)) / \
(mu + 1.0/2.0)/(mu + 3.0/2.0)/(mu - 1.0/2.0) * \
(2.0/(1.0 + np.sqrt(4.0*epsilon**2+1.0)))**(2.0*mu+1.0)
return c
def inverse_multiquadric(epsilon, r):
#-- inverse multiquadratic
f = 1.0/np.sqrt((epsilon*r)**2 + 1.0)
return f
def inverse_multiquadric_expansion(epsilon, mu):
c = 4.0*np.pi/(mu+1.0/2.0)*(2.0/(1.0+np.sqrt(4.0*epsilon**2+1.)))**(2*mu+1.)
return c
def inverse_quadratic(epsilon, r):
#-- inverse quadratic
f = 1.0/(1.0+(epsilon*r)**2)
return f
def inverse_quadratic_expansion(epsilon, mu):
c = 4.0*np.pi**(3.0/2.0)*scipy.special.factorial(mu) / \
scipy.special.gamma(mu + 3.0/2.0)/(1.0 + 4.0*epsilon**2)**(mu+1) * \
scipy.special.hyp2f1(mu+1,mu+1,2.*mu+2,4.*epsilon**2/(1.+4.*epsilon**2))
return c
def gaussian(epsilon, r):
#-- gaussian
f = np.exp(-(epsilon*r)**2)
return f
def gaussian_expansion(epsilon, mu):
c = 4.0*np.pi**(3.0/2.0)*np.exp(-2.0*epsilon**2) * \
scipy.special.iv(mu + 1.0/2.0, 2.0*epsilon**2)/epsilon**(2.0*mu + 1.0)
return c
def linear(epsilon, r):
#-- linear polynomial
return r
def cubic(epsilon, r):
#-- cubic polynomial
f = r**3
return f
def quintic(epsilon, r):
#-- quintic polynomial
f = r**5
return f
def thin_plate(epsilon, r):
#-- thin plate spline
f = r**2 * np.log(r)
#-- the spline is zero at zero
f[r == 0] = 0.0
return f
#-- calculate great-circle distance between between n-vectors
def angle_matrix(x,cntrs):
s,M = np.shape(x)
s,N = np.shape(cntrs)
A = np.zeros((M,N))
A[:,:] = np.arccos(np.dot(x.transpose(), cntrs))
A[np.isnan(A)] = 0.0
return A
#-- calculate Euclidean distances between points (default norm)
def distance_matrix(x,cntrs):
s,M = np.shape(x)
s,N = np.shape(cntrs)
#-- decompose Euclidean distance: (x-y)^2 = x^2 - 2xy + y^2
dx2=np.kron(np.ones((1,N)), np.sum(x * x, axis=0)[:,np.newaxis]).astype('f')
dxy=2.0*np.dot(x.transpose(), cntrs).astype('f')
dy2=np.kron(np.ones((M,1)), np.sum(cntrs * cntrs, axis=0)).astype('f')
D = np.sqrt(dx2 - dxy + dy2)
return D
#-- calculate spherical harmonics of degree l evaluated at (theta,phi)
def spherical_harmonic_matrix(l,theta,phi):
#-- calculate legendre polynomials
nth = len(theta)
Pl = legendre(l,np.cos(theta)).transpose()
#-- calculate degree dependent factors C and F
m = np.arange(0,l+1)#-- spherical harmonic orders up to degree l
C = np.sqrt((2.0*l + 1.0)/(4.0*np.pi))
F=np.sqrt(scipy.special.factorial(1+l-m-1)/scipy.special.factorial(1+l+m-1))
F=np.kron(np.ones((nth,1)), F[np.newaxis,:])
#-- calculate Euler's of spherical harmonic order multiplied by azimuth phi
mphi = np.exp(1j*np.dot(np.squeeze(phi)[:,np.newaxis],m[np.newaxis,:]))
#-- calculate spherical harmonics
Ylms = F*Pl[:,0:l+1]*mphi
#-- multiply by C and convert to reduced matrix (theta,Slm:Clm)
SPH = C*np.concatenate((np.imag(Ylms[:,:0:-1]),np.real(Ylms)), axis=1)
return SPH
#-- RBF interpolant with shape parameter epsilon through the node points
#-- (theta,phi) with function values f from Fornberg
#-- Outputs beta: the expansion coefficients of the interpolant with respect to
#-- the RBF_QR basis.
def RBF_QR(theta,phi,epsilon,data,RBF):
n = len(phi)
Y1 = np.zeros((n,n))
B1 = np.zeros((n,n))
#-- counter variable for filling spherical harmonic matrix
index = 0
#-- difference adding the next spherical harmonic degree
d = 0.0
#-- degree of the n_th spherical harmonic
l = 0
l_n = np.ceil(np.sqrt(n))-1
#-- floating point machine precision
eps = np.finfo(np.float).eps
while (d < -np.log10(eps)):
#-- create new variables for Y and B which will resize if (l > (l_n -1))
lmax = np.max([l_n,l])
Y = np.zeros((n,int((lmax+1)**2)))
Y[:,:index] = Y1[:,:index]
B = np.zeros((n,int((lmax+1)**2)))
B[:,:index] = B1[:,:index]
#-- Each loop adds a block of columns of SPH of degree l to Y and to B.
#-- Compute the spherical harmonics matrix
Y[:,index:2*l+index+1] = spherical_harmonic_matrix(l,theta,phi)
#-- Compute the expansion coefficients matrix
B[:,index:2*l+index+1] = Y[:,index:2*l+index+1]*RBF(epsilon,l)
B[:,index+l] = B[:,index+l]/2.0
#-- Truncation criterion
if (l > (l_n - 1)):
dN1 = np.linalg.norm(B[:,int(l_n**2):int((l_n+1)**2)], ord=np.inf)
dN2 = np.linalg.norm(B[:,int((l+1)**2)-1], ord=np.inf)
d = np.log10(dN1/dN2*epsilon**(2*(l_n-l)))
#-- copy B to B1 and Y to Y1
B1 = B.copy()
Y1 = Y.copy()
#-- Calculate column index of next block
index += 2*l+1
l += 1
#-- QR-factorization to find the RBF_QR basis
Q,R = np.linalg.qr(B)
#-- Introduce the powers of epsilon
X1 = np.kron(np.ones((n,1)), np.ceil(np.sqrt(np.arange(n,l**2))))
X2 = np.kron(np.ones((1,l**2-n)),
(np.ceil(np.sqrt(np.arange(1,n+1)))-1)[:,np.newaxis])
E = epsilon**(2.0*(X1 - X2))
#-- Solve the interpolation linear system
R_beta = np.transpose(E*np.linalg.lstsq(R[:n,:n], R[:n,n:], rcond=-1)[0])
R_new = np.concatenate((np.eye(n),R_beta),axis=0)
w = np.linalg.lstsq(np.dot(Y,R_new), data, rcond=-1)[0]
return (R_new,w)