# Copyright (c) 2022 CNES
#
# All rights reserved. Use of this source code is governed by a
# BSD-style license that can be found in the LICENSE file.
"""
RTree spatial index
-------------------
"""
from typing import Optional, Tuple
import numpy as np
from . import core, geodetic
[docs]class RTree:
"""R*Tree spatial index for geodetic scalar values.
Args:
system: WGS of the coordinate system used to transform equatorial
spherical positions (longitudes, latitudes, altitude) into ECEF
coordinates. If not set the geodetic system used is WGS-84.
Default to ``None``.
dtype: Data type of the instance to create.
ndims: The number of dimensions of the tree. This dimension must be
at least equal to 3 to store the ECEF coordinates of the points.
Default to ``3``.
"""
def __init__(self,
system: Optional[geodetic.Spheroid] = None,
dtype: Optional[np.dtype] = None,
ndims: int = 3):
"""Initialize a new R*Tree."""
dtype = dtype or np.dtype("float64")
if ndims < 3:
raise ValueError("ndims must be >= 3")
if dtype == np.dtype("float64"):
self._instance = getattr(core, f"RTree{ndims}DFloat64")(system)
elif dtype == np.dtype("float32"):
self._instance = getattr(core, f"RTree{ndims}DFloat32")(system)
else:
raise ValueError(f"dtype {dtype} not handled by the object")
self.dtype = dtype
[docs] def bounds(
self
) -> Tuple[Tuple[float, float, float], Tuple[float, float, float]]:
"""Returns the box able to contain all values stored in the container.
Returns:
A tuple that contains the coordinates of the minimum and maximum
corners of the box able to contain all values stored in the
container or an empty tuple if there are no values in the container.
"""
return self._instance.bounds()
[docs] def clear(self) -> None:
"""Removes all values stored in the container."""
return self._instance.clear()
[docs] def __len__(self):
"""Returns the number of values stored in the tree."""
return self._instance.__len__()
[docs] def __bool__(self):
"""Returns true if the tree is not empty."""
return self._instance.__bool__()
[docs] def packing(self, coordinates: np.ndarray, values: np.ndarray) -> None:
"""The tree is created using packing algorithm (The old data is erased
before construction.)
Args:
coordinates: a matrix ``(n, ndims)`` where ``n`` is the number of
observations and ``ndims`` is the number of coordinates in
order: longitude and latitude in degrees, altitude in meters and
then the other coordinates defined in Euclidean space if
``dims`` > 3. If the shape of the matrix is ``(n, ndims)`` then
the method considers the altitude constant and equal to zero.
values: An array of size ``(n)`` containing the values associated
with the coordinates provided.
"""
self._instance.packing(coordinates, values)
[docs] def insert(self, coordinates: np.ndarray, values: np.ndarray) -> None:
"""Insert new data into the search tree.
Args:
coordinates: a matrix ``(n, ndims)`` where ``n`` is the number of
observations and ``ndims`` is the number of coordinates in
order: longitude and latitude in degrees, altitude in meters and
then the other coordinates defined in Euclidean space if
``dims`` > 3. If the shape of the matrix is ``(n, ndims)`` then
the method considers the altitude constant and equal to zero.
values: An array of size ``(n)`` containing the values associated
with the coordinates provided.
"""
self._instance.insert(coordinates, values)
[docs] def query(self,
coordinates: np.ndarray,
k: Optional[int] = 4,
within: Optional[bool] = True,
num_threads: Optional[int] = 0) -> Tuple[np.ndarray, np.ndarray]:
"""Search for the nearest K nearest neighbors of a given point.
Args:
coordinates: a matrix ``(n, ndims)`` where ``n`` is the number of
observations and ``ndims`` is the number of coordinates in
order: longitude and latitude in degrees, altitude in meters and
then the other coordinates defined in Euclidean space if
``dims`` > 3. If the shape of the matrix is ``(n, ndims)`` then
the method considers the altitude constant and equal to zero.
k: The number of nearest neighbors to be searched. Defaults to
``4``.
within: If true, the method ensures that the neighbors found are
located within the point of interest. Defaults to ``false``.
num_threads: The number of threads to use for the computation. If 0
all CPUs are used. If 1 is given, no parallel computing code is
used at all, which is useful for debugging. Defaults to ``0``.
Returns:
A tuple containing a matrix describing for each provided position,
the distance, in meters, between the provided position and the found
neighbors and a matrix containing the value of the different
neighbors found for all provided positions.
If no neighbors are found, the distance and the value are set to
``-1``.
"""
return self._instance.query(coordinates, k, within, num_threads)
[docs] def inverse_distance_weighting(
self,
coordinates: np.ndarray,
radius: Optional[float] = None,
k: Optional[int] = 9,
p: Optional[int] = 2,
within: Optional[bool] = True,
num_threads: Optional[int] = 0) -> Tuple[np.ndarray, np.ndarray]:
"""Interpolation of the value at the requested position by inverse
distance weighting method.
Args:
coordinates: a matrix ``(n, ndims)`` where ``n`` is the number of
observations and ``ndims`` is the number of coordinates in
order: longitude and latitude in degrees, altitude in meters and
then the other coordinates defined in Euclidean space if
``dims`` > 3. If the shape of the matrix is ``(n, ndims)`` then
the method considers the altitude constant and equal to zero.
radius: The maximum radius of the search (m). Defaults The maximum
distance between two points.
k: The number of nearest neighbors to be used for calculating the
interpolated value. Defaults to ``9``.
p: The power parameters. Defaults to ``2``.
within: If true, the method ensures that the neighbors found are
located around the point of interest. In other words, this
parameter ensures that the calculated values will not be
extrapolated. Defaults to ``true``.
num_threads: The number of threads to use for the computation. If 0
all CPUs are used. If 1 is given, no parallel computing code is
used at all, which is useful for debugging. Defaults to ``0``.
Returns:
The interpolated value and the number of neighbors used in
the calculation.
"""
return self._instance.inverse_distance_weighting(
coordinates, radius, k, p, within, num_threads)
[docs] def radial_basis_function(
self,
coordinates: np.ndarray,
radius: Optional[float] = None,
k: Optional[int] = 9,
rbf: Optional[str] = None,
epsilon: Optional[float] = None,
smooth: Optional[float] = 0,
within: Optional[bool] = True,
num_threads: Optional[int] = 0) -> Tuple[np.ndarray, np.ndarray]:
"""Interpolation of the value at the requested position by radial basis
function interpolation.
Args:
coordinates: a matrix ``(n, ndims)`` where ``n`` is the number of
observations and ``ndims`` is the number of coordinates in
order: longitude and latitude in degrees, altitude in meters and
then the other coordinates defined in Euclidean space if
``dims`` > 3. If the shape of the matrix is ``(n, ndims)`` then
the method considers the altitude constant and equal to zero.
radius: The maximum radius of the search (m). Defaults The maximum
distance between two points.
k: The number of nearest neighbors to be used for calculating the
interpolated value. Defaults to ``9``.
rbf: The radial basis function, based on the radius, :math:`r`
given by the distance between points. This parameter can take
one of the following values:
* ``cubic``: :math:`\\varphi(r) = r^3`
* ``gaussian``: :math:`\\varphi(r) = e^{-(\\dfrac{r}
{\\varepsilon})^2}`
* ``inverse_multiquadric``: :math:`\\varphi(r) = \\dfrac{1}
{\\sqrt{1+(\\dfrac{r}{\\varepsilon})^2}}`
* ``linear``: :math:`\\varphi(r) = r`
* ``multiquadric``: :math:`\\varphi(r) = \\sqrt{1+(
\\dfrac{r}{\\varepsilon})^2}`
* ``thin_plate``: :math:`\\varphi(r) = r^2 \\ln(r)`
Default to ``multiquadric``
epsilon: adjustable constant for gaussian or multiquadrics
functions. Default to the average distance between nodes.
smooth: values greater than zero increase the smoothness of the
approximation. Default to 0 (interpolation).
within: If true, the method ensures that the neighbors found are
located around the point of interest. In other words, this
parameter ensures that the calculated values will not be
extrapolated. Defaults to ``true``.
num_threads: The number of threads to use for the computation. If 0
all CPUs are used. If 1 is given, no parallel computing code is
used at all, which is useful for debugging. Defaults to ``0``.
Returns:
The interpolated value and the number of neighbors used in the
calculation.
"""
adjustable = ['gaussian', 'inverse_multiquadric', 'multiquadric']
non_adjustable = ['cubic', 'linear', 'thin_plate']
rbf = rbf or adjustable[-1]
if epsilon is not None and rbf in non_adjustable:
raise ValueError(
f"epsilon must be None for {', '.join(non_adjustable)} RBF")
if rbf not in adjustable + non_adjustable:
raise ValueError(f"Radial basis function {rbf!r} is not defined")
rbf = "".join(item.capitalize() for item in rbf.split("_"))
return self._instance.radial_basis_function(
coordinates, radius, k, getattr(core.RadialBasisFunction, rbf),
epsilon, smooth, within, num_threads)
[docs] def window_function(
self,
coordinates: np.ndarray,
radius: float,
k: Optional[int] = 9,
wf: Optional[str] = None,
arg: Optional[float] = None,
within: Optional[bool] = True,
num_threads: Optional[int] = 0) -> Tuple[np.ndarray, np.ndarray]:
"""Interpolation of the value at the requested position by window
function.
The interpolated value will be equal to the expression:
.. math::
\\frac{\\sum_{i=1}^{k} \\omega(d_i,r)x_i}
{\\sum_{i=1}^{k} \\omega(d_i,r)}
where :math:`d_i` is the distance between the point of interest and
the :math:`i`-th neighbor, :math:`r` is the radius of the search,
:math:`x_i` is the value of the :math:`i`-th neighbor, and
:math:`\\omega(d_i,r)` is weight calculated by the window function
describe above.
Args:
coordinates: a matrix ``(n, ndims)`` where ``n`` is the number of
observations and ``ndims`` is the number of coordinates in
order: longitude and latitude in degrees, altitude in meters and
then the other coordinates defined in Euclidean space if
``dims`` > 3. If the shape of the matrix is ``(n, ndims)`` then
the method considers the altitude constant and equal to zero.
radius: The maximum radius of the search (m).
k: The number of nearest neighbors to be used for calculating the
interpolated value. Defaults to ``9``.
wf: The window function, based on the distance the distance between
points (:math:`d`) and the radius (:math:`r`). This parameter
can take one of the following values:
* ``blackman``: :math:`w(d) = 0.42659 - 0.49656 \\cos(
\\frac{\\pi (d + r)}{r}) + 0.076849 \\cos(
\\frac{2 \\pi (d + r)}{r})`
* ``blackman_harris``: :math:`w(d) = 0.35875 - 0.48829
\\cos(\\frac{\\pi (d + r)}{r}) + 0.14128
\\cos(\\frac{2 \\pi (d + r)}{r}) - 0.01168
\\cos(\\frac{3 \\pi (d + r)}{r})`
* ``boxcar``: :math:`w(d) = 1`
* ``flat_top``: :math:`w(d) = 0.21557895 -
0.41663158 \\cos(\\frac{\\pi (d + r)}{r}) +
0.277263158 \\cos(\\frac{2 \\pi (d + r)}{r}) -
0.083578947 \\cos(\\frac{3 \\pi (d + r)}{r}) +
0.006947368 \\cos(\\frac{4 \\pi (d + r)}{r})`
* ``lanczos``: :math:`w(d) = \\left\\{\\begin{array}{ll}
sinc(\\frac{d}{r}) \\times sinc(\\frac{d}{arg \\times r}),
& d \\le arg \\times r \\\\ 0,
& d \\gt arg \\times r \\end{array} \\right\\}`
* ``gaussian``: :math:`w(d) = e^{ -\\frac{1}{2}\\left(
\\frac{d}{\\sigma}\\right)^2 }`
* ``hamming``: :math:`w(d) = 0.53836 - 0.46164
\\cos(\\frac{\\pi (d + r)}{r})`
* ``nuttall``: :math:`w(d) = 0.3635819 - 0.4891775
\\cos(\\frac{\\pi (d + r)}{r}) + 0.1365995
\\cos(\\frac{2 \\pi (d + r)}{r})`
* ``parzen``: :math:`w(d) = \\left\\{ \\begin{array}{ll} 1 - 6
\\left(\\frac{2*d}{2*r}\\right)^2
\\left(1 - \\frac{2*d}{2*r}\\right),
& d \\le \\frac{2r + arg}{4} \\\\
2\\left(1 - \\frac{2*d}{2*r}\\right)^3
& \\frac{2r + arg}{2} \\le d \\lt \\frac{2r +arg}{4}
\\end{array} \\right\\}`
* ``parzen_swot``: :math:`w(d) = \\left\\{\\begin{array}{ll}
1 - 6\\left(\\frac{2 * d}{2 * r}\\right)^2
+ 6\\left(1 - \\frac{2 * d}{2 * r}\\right), &
d \\le \\frac{2r}{4} \\\\
2\\left(1 - \\frac{2 * d}{2 * r}\\right)^3 &
\\frac{2r}{2} \\ge d \\gt \\frac{2r}{4} \\end{array}
\\right\\}`
arg: The optional argument of the window function. Defaults to
``1`` for ``lanczos``, to ``0`` for ``parzen`` and for all
other functions is ``None``.
within: If true, the method ensures that the neighbors found are
located around the point of interest. In other words, this
parameter ensures that the calculated values will not be
extrapolated. Defaults to ``true``.
num_threads: The number of threads to use for the computation. If 0
all CPUs are used. If 1 is given, no parallel computing code is
used at all, which is useful for debugging. Defaults to ``0``.
Returns:
The interpolated value and the number of neighbors used in the
calculation.
"""
wf = wf or "blackman"
if wf not in [
"blackman",
"blackman_harris",
"boxcar",
"flattop",
"gaussian",
"hamming",
"lanczos",
"nuttall",
"parzen",
"parzen_swot",
]:
raise ValueError(f"Window function {wf!r} is not defined")
if wf in ["gaussian", "lanczos", "parzen"]:
if arg is None:
defaults = dict(gaussian=None, lanczos=1, parzen=0)
arg = defaults[wf]
if wf == "lanczos" and arg < 1: # type: ignore
raise ValueError(
f"The argument of the function {wf!r} must be "
"greater than 1")
if wf == "parzen" and arg < 0: # type: ignore
raise ValueError(
f"The argument of the function {wf!r} must be "
"greater than 0")
if wf == "gaussian" and arg is None:
raise ValueError(
f"The argument of the function {wf!r} must be "
"specified")
else:
if arg is not None:
raise ValueError(f"The function {wf!r} does not support the "
"optional argument")
wf = "".join(item.capitalize() for item in wf.split("_"))
return self._instance.window_function(coordinates, radius, k,
getattr(core.WindowFunction, wf),
arg, within, num_threads)
[docs] def __getstate__(self) -> Tuple:
"""Return the state of the object for pickling purposes.
Returns:
The state of the object for pickling purposes.
"""
return (self.dtype, self._instance.__getstate__())
[docs] def __setstate__(self, state: Tuple):
"""Set the state of the object from pickling.
Args:
state: The state of the object for pickling purposes.
"""
if len(state) != 2:
raise ValueError("invalid state")
_class = RTree(None, state[0])
self.dtype = _class.dtype
_class._instance.__setstate__(state[1])
self._instance = _class._instance