Source code for spharpy.interpolate

"""Interpolation module for spherical data."""

from scipy import interpolate as spinterpolate
import numpy as np




class SmoothSphereBivariateSpline(spinterpolate.SmoothSphereBivariateSpline):
    """Smooth bivariate spline approximation in spherical coordinates.
    The implementation uses the method proposed by Dierckx [#]_.

    Parameters
    ----------
    sampling : :py:class:`pyfar.Coordinates`
        Coordinates object containing the positions for which the data
        is sampled
    data : array, float
        Array containing the data at the sampling positions. Has to be
        real-valued.
    w : array, float
        Weighting coefficients. The default is ``None``, which sets the
        weighting coefficients to ``1``.
    s : float, 1e-4
        Smoothing factor > 0
        Positive smoothing factor defined for estimation condition:
        ``sum((w(i)*(r(i) - s(theta(i), phi(i))))**2, axis=0) <= s``
        Default ``s=len(w)`` which should be a good value if ``1/w[i]`` is
        an estimate of the standard deviation of ``r[i]``. The default
        value is ``1e-4``
    eps : float, 1e-16
        Sets the threshold used to determine the effective rank of the
        underlying linear system of equations. Values smaller than this
        threshold are treated as zero during the fitting process.
        Depends on the used data type and numerical precision. The default
        is ``1e-16``.

    Note
    ----
    This is a wrapper for scipy's SmoothSphereBivariateSpline only using
    Coordinates objects as container for the azimuth and elevation angles.
    For detailed information see scipy's documentation.

    References
    ----------
    .. [#] P. Dierckx, “Algorithms for smoothing data on the sphere with
           tensor product splines,” p. 24, 1984.

    Examples
    --------
    >>> n_max = 10
    >>> sampling = spharpy.samplings.equalarea(
    ...     n_max, n_points=500, condition_num=np.inf)
    >>> Y = spharpy.spherical.spherical_harmonic_basis_real(n_max, sampling)
    >>> y_vec = spharpy.spherical.spherical_harmonic_basis_real(
    ...     n_max, Coordinates(1, 0, 0))
    >>> data = Y @ y_vec.T
    >>> data = np.sin(sampling.elevation)*np.sin(2*sampling.azimuth)
    >>> interp_grid = spharpy.samplings.hyperinterpolation(35)
    >>> data_grid = np.sin(interp_grid.elevation)*np.sin(2*interp_grid.azimuth)
    >>> interpolator = interpolate.SmoothSphereBivariateSpline(
    ...     sampling, data, s=1e-4)
    ...
    >>> interp_data = interpolator(interp_grid)

    """ # noqa: 501

    def __init__(self, sampling, data, w=None, s=1e-4, eps=1e-16):
        theta = sampling.colatitude
        phi = sampling.azimuth
        if np.any(np.iscomplex(data)):
            raise ValueError("Complex data is not supported.")
        super().__init__(theta, phi, data, w, s, eps)



    def __call__(self, interp_grid, dtheta=0, dphi=0):
        """Evaluate the spline on a new sampling grid.

        Parameters
        ----------
        interp_grid : :py:class:`pyfar.Coordinates`
            Coordinates object containing a new set of points for which data
            is to be interpolated.
        dtheta : int, optional
            Order of theta derivative. Default is ``0``.
        dphi : int, optional
            Order of phi derivative. Default is ``0``.

        """ # noqa: 501
        theta = interp_grid.colatitude
        phi = interp_grid.azimuth
        return super().__call__(
            theta, phi, dtheta=dtheta, dphi=dphi, grid=False)




    def get_coeffs(self):
        """Get the coefficients of the spline.

        Returns
        -------
        coeffs : ndarray, float
            The coefficients of the spline. The shape is
            (n_coeffs_theta, n_coeffs_phi, n_coefficients).
            The coefficients are ordered by increasing degree and order.
        """
        return super().get_coeffs()




    def get_knots(self):
        """Get the knots of the spline.

        Returns
        -------
        knots : tuple of ndarray, float
            The knots of the spline. The first element is the knots in theta
            direction and the second element is the knots in phi direction.
        """
        return super().get_knots()




    def get_residual(self):
        """Get the residual of the spline.

        Returns
        -------
        residual : float
            The residual of the spline, which is the sum of squared errors
            between the data and the spline evaluation.
        """
        return super().get_residual()