#!/usr/bin/env python
u"""
sph_spline.py
Written by Tyler Sutterley (05/2022)
Interpolates a sparse grid over a sphere using spherical surface splines in
tension following Wessel and Becker (2008)
Adapted from P. Wessel, SOEST, U of Hawaii, April 2008
Uses Generalized Legendre Function algorithm from Spanier and Oldman
"An Atlas of Functions", 1987
CALLING SEQUENCE:
output = sph_spline(lon, lat, data, longitude, latitude, tension=0)
INPUTS:
lon: input longitude
lat: input latitude
data: input data
longitude: output longitude
latitude: output latitude
OUTPUTS:
output: interpolated data
OPTIONS:
tension: tension to use in interpolation (greater than 0)
PYTHON DEPENDENCIES:
numpy: Scientific Computing Tools For Python
https://numpy.org
scipy: Scientific Tools for Python
https://docs.scipy.org/doc/
cython: C-extensions for Python
http://cython.org/
REFERENCES:
Wessel, P. and J. M. Becker, 2008, Interpolation using a generalized
Green's function for a spherical surface spline in tension,
Geophysical Journal International, doi:10.1111/j.1365-246X.2008.03829.x
UPDATE HISTORY:
Updated 05/2022: updated docstrings to numpy documentation format
Updated 01/2022: added function docstrings
Updated 09/2017: using rcond=-1 in numpy least-squares algorithms
Updated 08/2016: using cythonized version of generalized Legendre function
treat case for no tension but x is equal to 1 within machine precision
Written 08/2016
"""
import numpy as np
import scipy.special
from spatial_interpolators.PvQv_C import PvQv_C
[docs]def sph_spline(lon, lat, data, longitude, latitude, tension=0.):
"""
Interpolates a sparse grid over a sphere using spherical
surface splines in tension
Parameters
----------
lon: float
input longitude
lat: float
input latitude
data: float
input data
longitude: float
output longitude
latitude: float
output latitude
tension: float, default 0.0
tension to use in interpolation
Returns
-------
output: float
interpolated data grid
References
----------
.. [Wessel2008] P. Wessel, and J. M. Becker,
"Interpolation using a generalized Green's function for a
spherical surface spline in tension,"
*Geophysical Journal International*, 174(1), 21--28, (2008).
`doi: 10.1111/j.1365-246X.2008.03829.x <https://doi.org/10.1111/j.1365-246X.2008.03829.x>`_
"""
#-- 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)
#-- 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')
if (tension < 0):
raise ValueError('TENSION must be greater than 0')
#-- 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)
sz = len(longitude.flatten())
#-- Find and remove mean from data
data_mean = data.mean()
data_range = data.max() - data.min()
#-- Normalize data
data_norm = (data - data_mean) / data_range
#-- compute linear system
N = len(data)
GG = np.zeros((N,N))
for i in range(N):
Rd = np.dot(np.transpose([xs,ys,zs]), np.array([xs[i],ys[i],zs[i]]))
#-- remove singleton dimensions and calculate spherical surface splines
GG[i,:] = SSST(Rd, P=tension)
#-- Compute model m for normalized data
m = np.linalg.lstsq(GG,data_norm,rcond=-1)[0]
#-- calculate output interpolated array (or matrix)
output = np.zeros((sz))
for j in range(sz):
Re = np.dot(np.transpose([xs,ys,zs]), np.array([XI[j],YI[j],ZI[j]]))
#-- remove singleton dimensions and calculate spherical surface splines
gg = SSST(Re, P=tension)
output[j] = data_mean + data_range*np.dot(gg, m)
#-- reshape output to original dimensions and return
if (np.ndim(longitude) > 1):
output = output.reshape(nlon,nlat)
return output
#-- SSST: Spherical Surface Spline in Tension
#-- Returns the Green's function for a spherical surface spline in tension,
#-- following Wessel and Becker [2008].
#-- If p == 0 or not given then use minimum curvature solution with dilogarithm
def SSST(x, P=0):
#-- floating point machine precision
eps = np.finfo(np.float).eps
if (P == 0):
#-- use dilogarithm (Spence's function) if using splines without tension
y = np.zeros_like(x)
if np.any(np.abs(x) < (1.0 - eps)):
k, = np.nonzero(np.abs(x) < (1.0 - eps))
y[k] = scipy.special.spence(0.5 - 0.5*x[k])
#-- Deal with special cases x == +/- 1
if np.any(((x + eps) >= 1.0) | ((x - eps) <= -1.0)):
k, = np.nonzero(((x + eps) >= 1.0) | ((x - eps) <= -1.0))
y[k] = scipy.special.spence(0.5 - 0.5*np.sign(x[k]))
else:
#-- if in tension
#-- calculate tension parameter
v = (-1.0 + np.lib.scimath.sqrt(1.0 - 4.0*P**2))/2.0
#-- Initialize output array
y = np.zeros_like(x, dtype=v.dtype)
A = np.pi/np.sin(v*np.pi)
#-- Where Pv solution works
if np.any(np.abs(x) < (1.0 - eps)):
k, = np.nonzero(np.abs(x) < (1.0 - eps))
y[k] = A*Pv(-x[k],v) - np.log(1.0 - x[k])
#-- Approximations where x is close to -1 or 1 using values from
#-- "An Atlas of Functions" by Spanier and Oldham, 1987 (590)
#-- Deal with special case x == -1
if np.any((x - eps) <= -1.0):
k, = np.nonzero((x - eps) <= -1.0)
y[k] = A - np.log(2.0)
#-- Deal with special case x == +1
if np.any((x + eps) >= 1.0):
k, = np.nonzero((x + eps) >= 1.0)
y[k] = np.pi*(1.0/np.tan(v*np.pi)) + 2.0*(np.euler_gamma +
scipy.special.psi(1.0+v)) - np.log(2.0)
#-- use only the real part (remove insignificant imaginary noise)
y = np.real(y)
return y
#-- Calculate Legendre function of the first kind for arbitrary degree v
def Pv(x,v):
P = np.zeros_like(x, dtype=v.dtype)
for i, val in enumerate(x):
if (val == -1):
p = np.inf
else:
#-- use compiled Cython version of PvQv (PvQv_C.so from PvQv_C.pyx)
p,q,k = PvQv_C(val, v)
# p,q,k = PvQv(val, v)
P[i] = p
return P
#-- Calculate generalized Legendre functions of arbitrary degree v
#-- Based on recipe in "An Atlas of Functions" by Spanier and Oldham, 1987 (589)
#-- Pv is the Legendre function of the first kind
#-- Qv is the Legendre function of the second kind
def PvQv(x, v):
iter = 0
if (x == -1):
P = -np.inf
Q = -np.inf
elif (x == +1):
P = 1.0
Q = np.inf
else:
#-- set a and R to 1
a = 1.0
R = 1.0
K = 4.0*np.sqrt(np.abs(v - v**2))
if ((np.abs(1 + v) + np.floor(1 + v.real)) == 0):
a = 1.0e99
v = -1.0 - v
#-- s and c = sin and cos of (pi*v/2.0)
s = np.sin(0.5*np.pi*v)
c = np.cos(0.5*np.pi*v)
w = (0.5 + v)**2
#-- if v is less than or equal to six (repeat until greater than six)
while (v.real <= 6.0):
v += 2
R = R*(v - 1.0)/v
#-- calculate X and g and update R
X = 1.0 / (4.0 + 4.0*v)
g = 1.0 + 5*X*(1.0 - 3.0*X*(0.35 + 6.1*X))
R = R*(1.0 - X*(1.0 - g*X/2))/np.sqrt(8.0*X)
#-- set g and u to 2.0*x
g = 2.0*x
u = 2.0*x
#-- set f and t to 1
f = 1.0
t = 1.0
#-- set k to 1/2
k = 0.5
#-- calculate new X
X = 1.0 + (1e8/(1.0 - x**2))
#-- update t
t = t*x**2 * (k**2.0 - w)/((k + 1.0)**2 - 0.25)
#-- add 1 to k
k += 1.0
#-- add t to f
f += t
#-- update u
u = u*x**2 * (k**2 - w)/((k + 1)**2 - 0.25)
#-- add 1 to k
k += 1.0
#-- add u to g
g += u
#-- if k is less than K and |Xt| is greater than |f|
#-- repeat previous set of operations until valid
while ((k < K) | (np.abs(X*t) > np.abs(f))):
iter += 1
t = t*x**2 * (k**2.0 - w) / ((k + 1.0)**2 - 0.25)
k += 1.0
f += t
u = u*x**2 * (k**2.0 - w) / ((k + 1.0)**2 - 0.25)
k += 1.0
g += u
#-- update f and g
f += (x**2*t/(1.0 - x**2))
g += (x**2*u/(1.0 - x**2))
#-- calculate generalized Legendre functions
P = ((s*g*R) + (c*f/R))/np.sqrt(np.pi)
Q = a*np.sqrt(np.pi)*((c*g*R) - (s*f/R))/2.0
#-- return P, Q and number of iterations
return (P, Q, iter)