#!/usr/bin/env python
u"""
interpolate.py
Written by Tyler Sutterley (12/2022)
Interpolators for spatial data
PYTHON DEPENDENCIES:
numpy: Scientific Computing Tools For Python
https://numpy.org
https://numpy.org/doc/stable/user/numpy-for-matlab-users.html
scipy: Scientific Tools for Python
https://docs.scipy.org/doc/
UPDATE HISTORY:
Written 12/2022
"""
from __future__ import annotations
import numpy as np
import scipy.spatial
import scipy.interpolate
import pyTMD.spatial
# PURPOSE: bilinear interpolation of input data to output data
[docs]def bilinear(
ilon: np.ndarray,
ilat: np.ndarray,
idata: np.ndarray,
lon: np.ndarray,
lat: np.ndarray,
fill_value: float = np.nan,
dtype: str | np.dtype = np.float64
):
"""
Bilinear interpolation of input data to output coordinates
Parameters
----------
ilon: np.ndarray
longitude of tidal model
ilat: np.ndarray
latitude of tidal model
idata: np.ndarray
tide model data
lat: np.ndarray
output latitude
lon: np.ndarray
output longitude
fill_value: float, default np.nan
invalid value
dtype: np.dtype, default np.float64
output data type
Returns
-------
data: np.ndarray
interpolated data
"""
# verify that input data is masked array
if not isinstance(idata, np.ma.MaskedArray):
idata = np.ma.array(idata)
idata.mask = np.zeros_like(idata, dtype=bool)
# find valid points (within bounds)
valid, = np.nonzero((lon >= ilon.min()) & (lon <= ilon.max()) &
(lat > ilat.min()) & (lat < ilat.max()))
# interpolate gridded data values to data
npts = len(lon)
# allocate to output interpolated data array
data = np.ma.zeros((npts), dtype=dtype, fill_value=fill_value)
data.mask = np.ones((npts), dtype=bool)
# initially set all data to fill value
data.data[:] = data.fill_value
# for each valid point
for i in valid:
# calculating the indices for the original grid
ix, = np.nonzero((ilon[0:-1] <= lon[i]) & (ilon[1:] > lon[i]))
iy, = np.nonzero((ilat[0:-1] <= lat[i]) & (ilat[1:] > lat[i]))
# corner data values for adjacent grid cells
IM = np.ma.zeros((4), fill_value=fill_value, dtype=dtype)
IM.mask = np.ones((4), dtype=bool)
# corner weight values for adjacent grid cells
WM = np.zeros((4))
# build data and weight arrays
for j,XI,YI in zip([0,1,2,3],[ix,ix+1,ix,ix+1],[iy,iy,iy+1,iy+1]):
IM.data[j], = idata.data[YI,XI].astype(dtype)
IM.mask[j], = idata.mask[YI,XI]
WM[3-j], = np.abs(lon[i]-ilon[XI])*np.abs(lat[i]-ilat[YI])
# if on corner value: use exact
if (np.isclose(lat[i],ilat[iy]) & np.isclose(lon[i],ilon[ix])):
data.data[i] = idata.data[iy,ix].astype(dtype)
data.mask[i] = idata.mask[iy,ix]
elif (np.isclose(lat[i],ilat[iy+1]) & np.isclose(lon[i],ilon[ix])):
data.data[i] = idata.data[iy+1,ix].astype(dtype)
data.mask[i] = idata.mask[iy+1,ix]
elif (np.isclose(lat[i],ilat[iy]) & np.isclose(lon[i],ilon[ix+1])):
data.data[i] = idata.data[iy,ix+1].astype(dtype)
data.mask[i] = idata.mask[iy,ix+1]
elif (np.isclose(lat[i],ilat[iy+1]) & np.isclose(lon[i],ilon[ix+1])):
data.data[i] = idata.data[iy+1,ix+1].astype(dtype)
data.mask[i] = idata.mask[iy+1,ix+1]
elif np.any(np.isfinite(IM) & (~IM.mask)):
# find valid indices for data summation and weight matrix
ii, = np.nonzero(np.isfinite(IM) & (~IM.mask))
# calculate interpolated value for i
data.data[i] = np.sum(WM[ii]*IM[ii])/np.sum(WM[ii])
data.mask[i] = np.all(IM.mask[ii])
# return interpolated values
return data
[docs]def spline(
ilon: np.ndarray,
ilat: np.ndarray,
idata: np.ndarray,
lon: np.ndarray,
lat: np.ndarray,
fill_value: float = None,
dtype: str | np.dtype = np.float64,
reducer=np.ceil,
**kwargs
):
"""
`Bivariate spline interpolation
<https://docs.scipy.org/doc/scipy/reference/generated/
scipy.interpolate.RectBivariateSpline.html>`_
of input data to output coordinates
Parameters
----------
ilon: np.ndarray
longitude of tidal model
ilat: np.ndarray
latitude of tidal model
idata: np.ndarray
tide model data
lat: np.ndarray
output latitude
lon: np.ndarray
output longitude
fill_value: float or NoneType, default None
invalid value
dtype: np.dtype, default np.float64
output data type
reducer: obj, default np.ceil
operation for converting mask to boolean
kx: int, default 1
degree of the bivariate spline in the x-dimension
ky: int, default 1
degree of the bivariate spline in the y-dimension
kwargs: dict
additional arguments for ``scipy.interpolate.RectBivariateSpline``
Returns
-------
data: np.ndarray
interpolated data
"""
# set default keyword arguments
kwargs.setdefault('kx', 1)
kwargs.setdefault('ky', 1)
# verify that input data is masked array
if not isinstance(idata, np.ma.MaskedArray):
idata = np.ma.array(idata)
idata.mask = np.zeros_like(idata, dtype=bool)
# interpolate gridded data values to data
npts = len(lon)
# allocate to output interpolated data array
data = np.ma.zeros((npts), dtype=dtype, fill_value=fill_value)
data.mask = np.ones((npts), dtype=bool)
# construct splines for input data and mask
if np.iscomplexobj(idata):
s1 = scipy.interpolate.RectBivariateSpline(ilon, ilat,
idata.data.real.T, **kwargs)
s2 = scipy.interpolate.RectBivariateSpline(ilon, ilat,
idata.data.imag.T, **kwargs)
s3 = scipy.interpolate.RectBivariateSpline(ilon, ilat,
idata.mask.T, **kwargs)
# evaluate the spline at input coordinates
data.data.real[:] = s1.ev(lon, lat)
data.data.imag[:] = s2.ev(lon, lat)
data.mask[:] = reducer(s3.ev(lon, lat)).astype(bool)
else:
s1 = scipy.interpolate.RectBivariateSpline(ilon, ilat,
idata.data.T, **kwargs)
s2 = scipy.interpolate.RectBivariateSpline(ilon, ilat,
idata.mask.T, **kwargs)
# evaluate the spline at input coordinates
data.data[:] = s1.ev(lon, lat).astype(dtype)
data.mask[:] = reducer(s2.ev(lon, lat)).astype(bool)
# return interpolated values
return data
[docs]def regulargrid(
ilon: np.ndarray,
ilat: np.ndarray,
idata: np.ndarray,
lon: np.ndarray,
lat: np.ndarray,
fill_value: float = None,
dtype: str | np.dtype = np.float64,
reducer=np.ceil,
**kwargs
):
"""
`Regular grid interpolation
<https://docs.scipy.org/doc/scipy/reference/generated/
scipy.interpolate.RegularGridInterpolator.html>`_
of input data to output coordinates
Parameters
----------
ilon: np.ndarray
longitude of tidal model
ilat: np.ndarray
latitude of tidal model
idata: np.ndarray
tide model data
lat: np.ndarray
output latitude
lon: np.ndarray
output longitude
fill_value: float or NoneType, default None
invalid value
dtype: np.dtype, default np.float64
output data type
reducer: obj, default np.ceil
operation for converting mask to boolean
bounds_error: bool, default False
raise Exception when values are requested outside domain
method: str, default 'linear'
Method of interpolation
- ``'linear'``
- ``'nearest'``
- ``'slinear'``
- ``'cubic'``
- ``'quintic'``
kwargs: dict
additional arguments for ``scipy.interpolate.RegularGridInterpolator``
Returns
-------
data: np.ndarray
interpolated data
"""
# set default keyword arguments
kwargs.setdefault('bounds_error', False)
kwargs.setdefault('method', 'linear')
# verify that input data is masked array
if not isinstance(idata, np.ma.MaskedArray):
idata = np.ma.array(idata)
idata.mask = np.zeros_like(idata, dtype=bool)
# interpolate gridded data values to data
npts = len(lon)
# allocate to output interpolated data array
data = np.ma.zeros((npts), dtype=dtype, fill_value=fill_value)
data.mask = np.ones((npts), dtype=bool)
# use scipy regular grid to interpolate values for a given method
r1 = scipy.interpolate.RegularGridInterpolator((ilat, ilon),
idata.data, fill_value=fill_value, **kwargs)
r2 = scipy.interpolate.RegularGridInterpolator((ilat, ilon),
idata.mask, fill_value=1, **kwargs)
# evaluate the interpolator at input coordinates
data.data[:] = r1.__call__(np.c_[lat, lon])
data.mask[:] = reducer(r2.__call__(np.c_[lat, lon])).astype(bool)
# return interpolated values
return data
# PURPOSE: Nearest-neighbor extrapolation of valid data to output data
# PURPOSE: calculate Euclidean distances between points
def _distance(c1: np.ndarray, c2: np.ndarray):
"""
Calculate Euclidean distances between points
Parameters
----------
c1: np.ndarray
first set of coordinates
c2: np.ndarray
second set of coordinates
Returns
-------
c: np.ndarray
Euclidean distance
"""
# decompose Euclidean distance: (x-y)^2 = x^2 - 2xy + y^2
dx2 = np.sum(c1**2)
dxy = np.dot(c1[np.newaxis,:], c2.T)
dy2 = np.sum(c2**2, axis=1)
# calculate Euclidean distance
D, = np.sqrt(dx2 - 2.0*dxy + dy2)
return D