"""
The following documents the pyfar coordinates class and functions for
coordinate conversion. More background information is given in
:py:mod:`coordinates concepts <pyfar._concepts.coordinates>`.
Available sampling schemes are listed at :py:mod:`~pyfar.samplings`.
"""
import numpy as np
from scipy.spatial import cKDTree
from scipy.spatial.transform import Rotation as sp_rot
import deepdiff
import re
from copy import deepcopy
import pyfar as pf
[docs]class Coordinates():
"""
Container class for storing, converting, rotating, querying, and plotting
3D coordinate systems.
"""
# Notes on the structure for developing -----------------------------------
#
# * All implemented cordinate systems are defined in a nested dictionary.
# The dictionary is returned by self._systems()
#
# * The current coordinate system is contained in self._system. It is
# generated by self._make_system()
#
# * A dictionary search is used to check if domain, coordinates, and units
# are valid
#
# * The coordinate points are stored in sefl._points and returned by the
# getter functions.
#
# To implement a new coordinate system -----------------------------------
#
# * Add the system to self._systems()
# * Add a converter to the module if it does not exist, e.g., cart2ellip()
# * Add class functions, e.g., self.set_ellip() and self.get_ellip() if
# they dont exist.
#
# NOTE: Coordinate systems must have unique coordinate names, unless the
# coordinates are identical, i.e., the spherical coordinate systems
# all contain the radius, becuase the definition is identical in all
# cases, but the cylindrical system has the coordinate radius_z
# because the definition differs from the sperical radius.
[docs] def __init__(self, points_1=None, points_2=None, points_3=None,
domain='cart', convention=None, unit=None,
weights=None, sh_order=None, comment=None):
"""
Create :py:func:`Coordinates` object with or without coordinate points.
The points that enter the Coordinates object are defined by the
`domain`, `convention`, and `unit` as illustrated in the
:py:mod:`coordinates concepts <pyfar._concepts.coordinates>`:
+--------------------+----------+------------+----------+----------+
| domain, convention | points_1 | points_2 | points_3 | unit |
+====================+==========+============+==========+==========+
| cart, right | x | y | z | met |
+--------------------+----------+------------+----------+----------+
| sph, top_colat | azimuth | colatitude | radius | rad, deg |
+--------------------+----------+------------+----------+----------+
| sph, top_elev | azimuth | elevation | radius | rad, deg |
+--------------------+----------+------------+----------+----------+
| sph, side | lateral | polar | radius | rad, deg |
+--------------------+----------+------------+----------+----------+
| sph, front | phi | theta | radius | rad, deg |
+--------------------+----------+------------+----------+----------+
| cyl, top | azimuth | z | radius_z | rad, deg |
+--------------------+----------+------------+----------+----------+
For more information run
>>> coords = Coordinates()
>>> coords.systems()
Parameters
----------
points_1 : array like, number
points for the first coordinate
points_2 : array like, number
points for the second coordinate
points_3 : array like, number
points for the third coordinate
domain : string
domain of the coordinate system
``'cart'``
Cartesian
``'sph'``
Spherical
``'cyl'``
Cylindrical
The default is ``'cart'``.
convention: string
coordinate convention (see above)
The default is ``'right'`` if domain is ``'cart'``,
``'top_colat'`` if domain is ``'sph'``, and ``'top'`` if domain is
``'cyl'``.
unit: string
unit of the coordinate system. By default the first available unit
is used, which is meters (``'met'``) for ``domain = 'cart'`` and
radians (``'rad'``) in all other cases (See above).
weights: array like, number, optional
sampling weights for the coordinate points. Must have same `size`
as the points points, i.e., if `points` have five entries, the
`weights` must also have five entries. The default is ``None``.
sh_order : int, optional
maximum spherical harmonic order of the sampling grid.
The default is ``None``.
comment : str, optional
comment about the stored coordinate points. The default is
``None``.
"""
if points_1 is None:
points_1 = []
if points_2 is None:
points_2 = []
if points_3 is None:
points_3 = []
# init emtpy object
super(Coordinates, self).__init__()
# set the coordinate system
self._system = self._make_system(domain, convention, unit)
# save coordinates to self
self._set_points(points_1, points_2, points_3, True)
# save meta data
self._set_weights(weights)
self._sh_order = sh_order
self._comment = comment
[docs] def set_cart(self, points_1, points_2, points_3,
convention='right', unit='met'):
"""
Enter coordinate points in cartesian coordinate systems.
The points that enter the Coordinates object are defined by the
`domain`, `convention`, and `unit`
+--------------------+----------+------------+----------+----------+
| domain, convention | points_1 | points_2 | points_3 | unit |
+====================+==========+============+==========+==========+
| cart, right | x | y | z | met |
+--------------------+----------+------------+----------+----------+
For more information run
>>> coords = Coordinates()
>>> coords.systems()
Parameters
----------
points_i: array like, number
points for the first, second, and third coordinate
convention : string, optional
convention in which the coordinate points are stored. The default
is ``'right'``.
unit : string, optional
unit in which the coordinate points are stored. The default is
``'met'`` for meters.
"""
# set the coordinate system
self._system = self._make_system('cart', convention, unit)
# save coordinates to self
self._set_points(points_1, points_2, points_3, True)
[docs] def get_cart(self, convention='right', unit='met', convert=False):
"""
Get coordinate points in cartesian coordinate systems.
The points that are returned are defined by the `domain`, `convention`,
and `unit`:
+--------------------+----------+------------+----------+----------+
| domain, convention | p[...,1] | p[...,1] | p[...,1] | units |
+====================+==========+============+==========+==========+
| cart, right | x | y | z | met |
+--------------------+----------+------------+----------+----------+
For more information run
>>> coords = Coordinates()
>>> coords.systems()
Parameters
----------
convention : string, optional
convention in which the coordinate points are stored. The default
is ``'right'``.
unit : string, optional
unit in which the coordinate points are stored. The default is
``'met'``.
convert : boolean, optional
if True, the internal representation of the samplings points will
be converted to the queried coordinate system. The default is
``False``, i.e., the internal presentation remains as it is.
Returns
-------
points : numpy array
coordinate points. ``points[...,0]`` holds the points for the first
coordinate, ``points[...,1]`` the points for the second, and
``points[...,2]`` the points for the third coordinate.
"""
# check if object is empty
if self.cshape == (0,):
raise ValueError('Object is empty.')
# make the new system
new_system = self._make_system('cart', convention, unit)
# return if system has not changed
if self._system == new_system:
return self._points
# convert to radians
pts = self._points.copy()
for nn, unit in enumerate(self._system['units']):
if unit == 'degrees':
pts[..., nn] = pts[..., nn] / 180 * np.pi
# convert to cartesian ...
# ... from spherical coordinate systems
if self._system['domain'] == 'sph':
if self._system['convention'] == 'top_colat':
x, y, z = sph2cart(pts[..., 0], pts[..., 1], pts[..., 2])
elif self._system['convention'] == 'top_elev':
x, y, z = sph2cart(pts[..., 0],
np.pi / 2 - pts[..., 1],
pts[..., 2])
elif self._system['convention'] == 'side':
x, z, y = sph2cart(pts[..., 1],
np.pi / 2 - pts[..., 0],
pts[..., 2])
elif self._system['convention'] == 'front':
y, z, x = sph2cart(pts[..., 0], pts[..., 1], pts[..., 2])
else:
# Can not be tested. Will only be raised if a coordinate system
# is not fully implemented.
raise ValueError(
(f"Conversion for {self._system['convention']} "
"is not implemented."))
# ... from cylindrical coordinate systems
elif self._system['domain'] == 'cyl':
if self._system['convention'] == 'top':
x, y, z = cyl2cart(pts[..., 0], pts[..., 1], pts[..., 2])
else:
# Can not be tested. Will only be raised if a coordinate system
# is not fully implemented.
raise ValueError(
(f"Conversion for {self._system['convention']} "
"is not implemented."))
else:
# Can not be tested. Will only be raised if a coordinate system
# is not fully implemented.
raise ValueError(
f"Conversion for {convention} is not implemented.")
# return points and convert internal state if desired
return self._return_system(x, y, z, new_system, convert)
[docs] def set_sph(self, points_1, points_2, points_3,
convention='top_colat', unit='rad'):
"""
Enter coordinate points in spherical coordinate systems.
The points that enter the Coordinates object are defined by the
`domain`, `convention`, and `unit`
+--------------------+----------+------------+----------+----------+
| domain, convention | points_1 | points_2 | points_3 | unit |
+====================+==========+============+==========+==========+
| sph, top_colat | azimuth | colatitude | radius | rad, deg |
+--------------------+----------+------------+----------+----------+
| sph, top_elev | azimuth | elevation | radius | rad, deg |
+--------------------+----------+------------+----------+----------+
| sph, side | lateral | polar | radius | rad, deg |
+--------------------+----------+------------+----------+----------+
| sph, front | phi | theta | radius | rad, deg |
+--------------------+----------+------------+----------+----------+
For more information run
>>> coords = Coordinates()
>>> coords.systems()
Parameters
----------
points_i: array like, number
points for the first, second, and third coordinate
convention : string, optional
convention in which the coordinate points are stored. The default
is ``'top_colat'``.
unit : string, optional
unit in which the coordinate points are stored. The default is
``'rad'``.
"""
# set the coordinate system
self._system = self._make_system('sph', convention, unit)
# save coordinates to self
self._set_points(points_1, points_2, points_3, True)
[docs] def get_sph(self, convention='top_colat', unit='rad', convert=False):
"""
Get coordinate points in spherical coordinate systems.
The points that are returned are defined by the `domain`,
`convention`, and `unit`:
+--------------------+----------+------------+----------+----------+
| domain, convention | p[...,1] | p[...,1] | p[...,1] | units |
+====================+==========+============+==========+==========+
| sph, top_colat | azimuth | colatitude | radius | rad, deg |
+--------------------+----------+------------+----------+----------+
| sph, top_elev | azimuth | elevation | radius | rad, deg |
+--------------------+----------+------------+----------+----------+
| sph, side | lateral | polar | radius | rad, deg |
+--------------------+----------+------------+----------+----------+
| sph, front | phi | theta | radius | rad, deg |
+--------------------+----------+------------+----------+----------+
For more information run
>>> coords = Coordinates()
>>> coords.systems()
Parameters
----------
convention : string, optional
convention in which the coordinate points are stored. The default
is ``'top_colat'``.
unit : string, optional
unit in which the coordinate points are stored. The default is
``'rad'``.
convert : boolean, optional
if True, the internal representation of the samplings points will
be converted to the queried coordinate system. The default is
``False``, i.e., the internal presentation remains as it is.
Returns
-------
points : numpy array
coordinate points. ``points[...,0]`` holds the points for the first
coordinate, ``points[...,1]`` the points for the second, and
``points[...,2]`` the points for the third coordinate.
"""
# check if object is empty
if self.cshape == (0,):
raise ValueError('Object is empty.')
# make the new system
new_system = self._make_system('sph', convention, unit)
# return if system has not changed
if new_system == self._system:
return self._points
# get cartesian system first
if not (self._system['domain'] == 'cart' and
self._system['convention'] == 'right'):
pts = self.get_cart('right', 'met')
# remove noise below eps
eps = np.finfo(np.float64).eps
pts[np.abs(pts) < eps] = 0
else:
pts = self._points.copy()
# convert to spherical...
# ... top polar systems
if convention[0:3] == 'top':
pts_1, pts_2, pts_3 = cart2sph(
pts[..., 0], pts[..., 1], pts[..., 2])
if convention == 'top_elev':
pts_2 = np.pi / 2 - pts_2
# ... side polar system
# (idea for simple converions from Robert Baumgartner and SOFA_API)
elif convention == 'side':
pts_2, pts_1, pts_3 = cart2sph(
pts[..., 0], pts[..., 2], -pts[..., 1])
# range angles
pts_1 = pts_1 - np.pi / 2
pts_2 = np.mod(pts_2 + np.pi / 2, 2 * np.pi) - np.pi / 2
# ... front polar system
elif convention == 'front':
pts_1, pts_2, pts_3 = cart2sph(
pts[..., 1], pts[..., 2], pts[..., 0])
else:
# Can not be tested. Will only be raised if a coordinate system
# is not fully implemented.
raise ValueError(
f"Conversion for {convention} is not implemented.")
# convert to degrees
if new_system['unit'] == 'deg':
pts_1 = pts_1 / np.pi * 180
pts_2 = pts_2 / np.pi * 180
# return points and convert internal state if desired
return self._return_system(pts_1, pts_2, pts_3, new_system, convert)
[docs] def set_cyl(self, points_1, points_2, points_3,
convention='top', unit='rad'):
"""
Enter coordinate points in cylindrical coordinate systems.
The points that enter the Coordinates object are defined by the
`domain`, `convention`, and `unit`
+--------------------+----------+------------+----------+----------+
| domain, convention | points_1 | points_2 | points_3 | unit |
+====================+==========+============+==========+==========+
| cyl, top | azimuth | z | radius_z | rad, deg |
+--------------------+----------+------------+----------+----------+
For more information run
>>> coords = Coordinates()
>>> coords.systems()
Parameters
----------
points_i: array like, number
points for the first, second, and third coordinate
convention : string, optional
convention in which the coordinate points are stored. The default
is ``'top'``.
unit : string, optional
unit in which the coordinate points are stored. The default is
``'rad'``.
"""
# set the coordinate system
self._system = self._make_system('cyl', convention, unit)
# save coordinates to self
self._set_points(points_1, points_2, points_3, True)
[docs] def get_cyl(self, convention='top', unit='rad', convert=False):
"""
Get coordinate points in cylindircal coordinate system.
The points that are returned are defined by the `domain`, `convention`,
and `unit`:
+--------------------+----------+------------+----------+----------+
| domain, convention | p[...,1] | p[...,1] | p[...,1] | units |
+====================+==========+============+==========+==========+
| cyl, top | azimuth | z | radius_z | rad, deg |
+--------------------+----------+------------+----------+----------+
For more information run
>>> coords = Coordinates()
>>> coords.systems()
Parameters
----------
convention : string, optional
convention in which the coordinate points are stored. The default
is ``'right'``.
unit : string, optional
unit in which the coordinate points are stored. The default is
``'met'``.
convert : boolean, optional
if True, the internal representation of the samplings points will
be converted to the queried coordinate system. The default is
False, i.e., the internal presentation remains as it is.
Returns
-------
points : numpy array
coordinate points. ``points[...,0]`` holds the points for the first
coordinate, ``points[...,1]`` the points for the second, and
``points[...,2]`` the points for the third coordinate.
"""
# check if object is empty
if self.cshape == (0,):
raise ValueError('Object is empty.')
# make the new system
new_system = self._make_system('cyl', convention, unit)
# return if system has not changed
if new_system == self._system:
return self._points
# convert to cartesian system first
if not (self._system['domain'] == 'cart' and
self._system['convention'] == 'right'):
pts = self.get_cart('right', 'met')
# remove noise below eps
eps = np.finfo(np.float64).eps
pts[np.abs(pts) < eps] = 0
else:
pts = self._points.copy()
# convert to cylindrical ...
# ... top systems
if convention == 'top':
pts_1, pts_2, pts_3 = cart2cyl(
pts[..., 0], pts[..., 1], pts[..., 2])
else:
# Can not be tested. Will only be raised if a coordinate system
# is not fully implemented.
raise ValueError(
f"Conversion for {convention} is not implemented.")
# convert to degrees
if new_system['unit'] == 'deg':
pts_1 = pts_1 / np.pi * 180
# return points and convert internal state if desired
return self._return_system(pts_1, pts_2, pts_3, new_system, convert)
@property
def weights(self):
"""Get sampling weights."""
return self._weights
@weights.setter
def weights(self, value):
"""Set sampling weights."""
self._set_weights(value)
@property
def sh_order(self):
"""Get the maximum spherical harmonic order."""
return self._sh_order
@sh_order.setter
def sh_order(self, value):
"""Set the maximum spherical harmonic order."""
self._sh_order = value
@property
def comment(self):
"""Get comment."""
return self._comment
@comment.setter
def comment(self, value):
"""Set comment."""
self._comment = value
@property
def cshape(self):
"""
Return channel shape.
The channel shape gives the shape of the coordinate points excluding
the last dimension, which is always 3.
"""
if self._points.size:
return self._points.shape[:-1]
else:
return (0,)
@property
def cdim(self):
"""
Return channel dimension.
The channel dimension gives the number of dimensions of the coordinate
points excluding the last dimension.
"""
if self._points.size:
return self._points.ndim - 1
else:
return 0
@property
def csize(self):
"""
Return channel size.
The channel size gives the number of points stored in the coordinates
object.
"""
return self._points.size // 3
[docs] def systems(self, show='all', brief=False):
"""
Print coordinate systems and their description on the console.
.. note::
All coordinate systems are described with respect to the right
handed cartesian system (``domain='cart'``, ``convention='right'``).
Distances are always specified in meters, while angles can be
radians or degrees (``unit='rad'`` or ``unit='deg'``).
Parameters
----------
show: string, optional
``'current'`` to list the current corrdinate system or ``'all'``
to list all coordinate systems. The default is ``'all'``.
brief : boolean, optional
Will only list the domains, conventions and units if True. The
default is ``False``.
Returns
-------
Prints to console.
"""
if show == 'current':
domain = self._system['domain']
convention = self._system['convention']
unit = self._system['unit']
elif show == 'all':
domain = convention = unit = 'all'
else:
raise ValueError("show must be 'current' or 'all'.")
# get coordinate systems
systems = self._systems()
# print information
domains = list(systems) if domain == 'all' else [domain]
if brief:
print('domain, convention, unit')
print('- - - - - - - - - - - - -')
for dd in domains:
conventions = list(systems[dd]) if convention == 'all' \
else [convention]
for cc in conventions:
# current coordinates
coords = systems[dd][cc]['coordinates']
# current units
if unit != 'all':
units = [units for units in systems[dd][cc]['units']
if unit == units[0][0:3]]
else:
units = systems[dd][cc]['units']
# key for unit
unit_key = [u[0][0:3] for u in units]
print(f"{dd}, {cc}, [{', '.join(unit_key)}]")
else:
for dd in domains:
conventions = \
list(systems[dd]) if convention == 'all' else [convention]
for cc in conventions:
# current coordinates
coords = systems[dd][cc]['coordinates']
# current units
if unit != 'all':
units = [units for units in systems[dd][cc]['units']
if unit == units[0][0:3]]
else:
units = systems[dd][cc]['units']
# key for unit
unit_key = [u[0][0:3] for u in units]
print("- - - - - - - - - - - - - - - - - "
"- - - - - - - - - - - - - - - - -")
print(f"domain: {dd}, convention: {cc}, unit: "
f"[{', '.join(unit_key)}]\n")
print(systems[dd][cc]['description_short'] + '\n')
print("Coordinates:")
for nn, coord in enumerate(coords):
cur_units = [u[nn] for u in units]
print(
f"points_{nn + 1}: {coord} ",
f"[{', '.join(cur_units)}]")
print('\n' + systems[dd][cc]['description'] + '\n\n')
[docs] def show(self, mask=None, **kwargs):
"""
Show a scatter plot of the coordinate points.
Parameters
----------
mask : boolean numpy array, None, optional
Plot points in red if ``mask==True``. The default is ``None``,
which the same color for all points.
kwargs : optional
keyword arguments are passed to ``matplotlib.pyplot.scatter()``.
If a mask is provided and the key `c` is contained in kwargs, it
will be overwritten.
Returns
-------
ax : matplotlib.axes._subplots.Axes3DSubplot
The axis used for the plot.
"""
if mask is None:
pf.plot.scatter(self, **kwargs)
else:
mask = np.asarray(mask)
assert mask.shape == self.cshape,\
"'mask.shape' must be self.cshape"
colors = np.full(mask.shape, pf.plot.color('b'))
colors[mask] = pf.plot.color('r')
pf.plot.scatter(self, c=colors.flatten(), **kwargs)
[docs] def find_nearest_k(self, points_1, points_2, points_3, k=1,
domain='cart', convention='right', unit='met',
show=False):
"""
Find the k nearest coordinates points.
Parameters
----------
points_i : array like, number
first, second and third coordinate of the points to which the
nearest neighbors are searched.
k : int, optional
Number of points to return. k must be > 0. The default is ``1``.
domain : string, optional
domain of the points. The default is ``'cart'``.
convention: string, optional
convention of points. The default is ``'right'``.
unit : string, optional
unit of the points. The default is ``'met'`` for meters.
show : bool, optional
show a plot of the coordinate points. The default is ``False``.
Returns
-------
index : numpy array of ints
The locations of the neighbors in the getter methods (e.g.,
``self.get_cart``). Dimension according to `distance` (see below).
Missing neighbors are indicated with ``csize``. Also see Notes
below.
mask : boolean numpy array
mask that contains ``True`` at the positions of the selected points
and ``False`` otherwise. Mask is of shape ``cshape``.
Notes
-----
``numpy.spatial.cKDTree`` is used for the search, which requires an
(N, 3) array. The coordinate points in self are thus reshaped to
(`csize`, 3) before they are passed to ``cKDTree``. The index that
is returned refers to the reshaped coordinate points. To access the
points for example use
>>> points_reshaped = self.get_cart().reshape((self.csize, 3))
>>> points_reshaped[index]
Examples
--------
Find frontal point from a spherical coordinate system
.. plot::
>>> import pyfar as pf
>>> coords = pf.samplings.sph_lebedev(sh_order=10)
>>> result = coords.find_nearest_k(1, 0, 0, show=True)
"""
# check the input
assert isinstance(k, int) and k > 0 and k <= self.csize,\
"k must be an integer > 0 and <= self.csize."
# get the points
_, index, mask = self._find_nearest(
points_1, points_2, points_3,
domain, convention, unit, show, k, 'k')
return index, mask
[docs] def find_nearest_cart(self, points_1, points_2, points_3, distance,
domain='cart', convention='right', unit='met',
show=False, atol=1e-15):
"""
Find coordinates within a certain distance in meters to query points.
Parameters
----------
points_i : array like, number
first, second and third coordinate of the points to which the
nearest neighbors are searched.
distance : number
Euclidean distance in meters in which the nearest points are
searched. Must be >= 0.
domain : string, optional
domain of the points. The default is ``'cart'``.
convention: string, optional
convention of points. The default is ``'right'``.
unit : string, optional
unit of the points. The default is ``'met'`` for meters.
show : bool, optional
show a plot of the coordinate points. The default is ``False``.
atol : float, optional
a tolerance that is added to `distance`. The default is`` 1e-15``.
Returns
-------
index : numpy array of ints
The locations of the neighbors in the getter methods (e.g.,
``get_cart``). Dimension as in :py:func:`~find_nearest_k`.
Missing neighbors are indicated with ``csize``. Also see Notes
below.
mask : boolean numpy array
mask that contains ``True`` at the positions of the selected points
and ``False`` otherwise. Mask is of shape ``cshape``.
Notes
-----
``numpy.spatial.cKDTree`` is used for the search, which requires an
(N, 3) array. The coordinate points in self are thus reshaped to
(`csize`, 3) before they are passed to ``cKDTree``. The index that
is returned refers to the reshaped coordinate points. To access the
points for example use
>>> points_reshaped = self.get_cart().reshape((self.csize, 3))
>>> points_reshaped[index]
Examples
--------
Find frontal points within a distance of 0.5 meters
.. plot::
>>> import pyfar as pf
>>> coords = pf.samplings.sph_lebedev(sh_order=10)
>>> result = coords.find_nearest_cart(1, 0, 0, 0.5, show=True)
"""
# check the input
assert distance >= 0, "distance must be >= 0"
# get the points
distance, index, mask = self._find_nearest(
points_1, points_2, points_3,
domain, convention, unit, show, distance, 'cart', atol)
return index, mask
[docs] def find_nearest_sph(self, points_1, points_2, points_3, distance,
domain='sph', convention='top_colat', unit='rad',
show=False, atol=1e-15):
"""
Find coordinates within certain angular distance to the query points.
Parameters
----------
points_i : array like, number
first, second and third coordinate of the points to which the
nearest neighbors are searched.
distance : number
Great circle distance in degrees in which the nearest points are
searched. Must be >= 0 and <= 180.
domain : string, optional
domain of the input points. The default is ``'sph'``.
convention: string, optional
convention of the input points. The default is ``'top_colat'``.
unit: string, optional
unit of the input points. The default is ``'rad'``.
show : bool, optional
show a plot of the coordinate points. The default is ``False``.
atol : float, optional
a tolerance that is added to `distance`. The default is ``1e-15``.
Returns
-------
index : numpy array of ints
The locations of the neighbors in the getter methods (e.g.,
``get_cart``). Dimension as in :py:func:`~find_nearest_k`.
Missing neighbors are indicated with ``csize``. Also see Notes
below.
mask : boolean numpy array
mask that contains ``True`` at the positions of the selected points
and ``False`` otherwise. Mask is of shape ``cshape``.
Notes
-----
``numpy.spatial.cKDTree`` is used for the search, which requires an
(N, 3) array. The coordinate points in self are thus reshaped to
(`csize`, 3) before they are passed to ``cKDTree``. The index that
is returned refers to the reshaped coordinate points. To access the
points for example use
``points_reshaped = points.get_sph().reshape((points.csize, 3))``
``points_reshaped[index]``
Examples
--------
Find top points within a distance of 45 degrees
.. plot::
>>> import pyfar as pf
>>> coords = pf.samplings.sph_lebedev(sh_order=10)
>>> result = coords.find_nearest_sph(0, 0, 1, 45, show=True)
"""
# check the input
assert distance >= 0 and distance <= 180, \
"distance must be >= 0 and <= 180."
# get radius and check for equality
radius = self.get_sph()[..., 2]
delta_radius = np.max(radius) - np.min(radius)
if delta_radius > 1e-15:
raise ValueError(
"find_nearest_sph only works if all points have the same \
radius. Differences are larger than 1e-15")
# get the points
distance, index, mask = self._find_nearest(
points_1, points_2, points_3,
domain, convention, unit, show, distance, 'sph', atol,
np.max(radius))
return index, mask
[docs] def find_slice(self, coordinate: str, unit: str, value, tol=0,
show=False, atol=1e-15):
"""
Find a slice of the coordinates points.
Parameters
----------
coordinate : str
coordinate for slicing.
unit : str
unit in which the value is passed
value : number
value of the coordinate around which the points are sliced.
tol : number, optional
tolerance for slicing. Points are sliced within the range
``[value-tol, value+tol]``. The default is ``0``.
show : bool, optional
show a plot of the coordinate points. The default is ``False``.
atol : number, optional
a tolerance that is added to `tol`. The default is ``1e-15``.
Returns
-------
index : tuple of numpy arrays
The indices of the selected points as a tuple of arrays. The length
of the tuple matches :py:func:`~cdim`. The length of each array
matches the number of selected points.
mask : boolean numpy array
mask that contains True at the positions of the selected points and
False otherwise. Mask is of shape self.cshape.
Notes
-----
`value` must be inside the range of the coordinate (see ``.systems``).
However, `value` +/- `tol` may exceed the range.
Examples
--------
Find horizontal slice of spherical coordinate system within a ring of
+/- 10 degrees
.. plot::
>>> import pyfar as pf
>>> coords = pf.samplings.sph_lebedev(sh_order=10)
>>> result = coords.find_slice('elevation', 'deg', 0, 5, show=True)
"""
# check if the coordinate and unit exist
domain, convention, index = self._exist_coordinate(coordinate, unit)
# get type and range of coordinate
c_info = self._systems()[domain][convention][coordinate]
# convert input to radians
value = value / 180 * np.pi if unit == 'deg' else value
tol = tol / 180 * np.pi if unit == 'deg' else tol
# check if value is within the range of coordinate
if c_info[0] in ["bound", "cyclic"]:
assert c_info[1][0] <= value <= c_info[1][1],\
f"'value' is {value} but must be in the range {c_info[1]}."
# get the search range
rng = [value - tol, value + tol]
# wrap range if coordinate is cyclic
if c_info[0] == 'cyclic':
low = c_info[1][0]
upp = c_info[1][1]
if rng[0] < c_info[1][0] - atol:
rng[0] = (rng[0] - low) % (upp - low) + low
if rng[1] > c_info[1][1] + atol:
rng[1] = (rng[1] - low) % (upp - low) + low
# get the coordinates
coords = eval(f"self.get_{domain}('{convention}')")
coords = coords[..., index]
# get the mask
if rng[0] <= rng[1]:
mask = (coords >= rng[0] - atol) & (coords <= rng[1] + atol)
else:
mask = (coords >= rng[0] - atol) | (coords <= rng[1] + atol)
# plot all and returned points
if show:
self.show(mask)
index = np.where(mask)
return index, mask
[docs] def rotate(self, rotation: str, value=None, degrees=True, inverse=False):
"""
Rotate points stored in the object around the origin of coordinates.
This is a wrapper for ``scipy.spatial.transform.Rotation`` (see this
class for more detailed information).
Parameters
----------
rotation : str
``'quat'``
rotation given by quaternions.
``'matrix'``
rotation given by matrixes.
``'rotvec'``
rotation using rotation vectors.
``'xyz'``
rotation using euler angles. Up to three letters. E.g., ``'x'``
will rotate about the x-axis only, while ``'xz'`` will rotate
about the x-axis and then about the z-axis. Use lower letters
for extrinsic rotations (rotations about the axes of the
original coordinate system xyz, which remains motionless) and
upper letters for intrinsic rotations (rotations about the axes
of the rotating coordinate system XYZ, solidary with the moving
body, which changes its orientation after each elemental
rotation).
value : number, array like
amount of rotation in the format specified by `rotation` (see
above).
degrees : bool, optional
pass angles in degrees if using ``'rotvec'`` or euler angles
(``'xyz'``). The default is ``True``. Use False to pass angles in
radians.
inverse : bool, optional
Apply inverse rotation. The default is ``False``.
Notes
-----
Points are converted to the cartesian right handed coordinate system
for the rotation.
Examples
--------
Get a coordinates object
>>> import pyfar as pf
>>> coordinates = pf.samplings.sph_gaussian(sh_order=3)
Rotate 45 degrees about the y-axis using
1. quaternions
>>> coordinates.rotate('quat', [0 , 0.38268343, 0 , 0.92387953])
2. a rotation matrix
>>> coordinates.rotate('matrix',
... [[ 0.70710678, 0 , 0.70710678],
... [ 0 , 1 , 0. ],
... [-0.70710678, 0 , 0.70710678]])
3. a rotation vector
>>> coordinates.rotate('rotvec', [0, 45, 0])
4. euler angles
>>> coordinates.rotate('XYZ', [0, 45, 0])
To see the result of the rotation use
>>> coordinates.show()
"""
# initialize rotation object
if rotation == 'quat':
rot = sp_rot.from_quat(value)
elif rotation == 'matrix':
rot = sp_rot.from_matrix(value)
elif rotation == 'rotvec':
if degrees:
value = np.asarray(value) / 180 * np.pi
rot = sp_rot.from_rotvec(value)
elif not bool(re.search('[^x-z]', rotation.lower())):
# only check if string contains xyz, everything else is checked in
# from_euler()
rot = sp_rot.from_euler(rotation, value, degrees)
else:
raise ValueError("rotation must be 'quat', 'matrix', 'rotvec', "
"or from ['x', 'y', 'z'] or ['X', 'Y', 'Z'] but "
f"is '{rotation}'")
# current shape
shape = self.cshape
# apply rotation
points = rot.apply(self.get_cart().reshape((self.csize, 3)), inverse)
# set points
self.set_cart(points[:, 0].reshape(shape),
points[:, 1].reshape(shape),
points[:, 2].reshape(shape))
[docs] def copy(self):
"""Return a deep copy of the Coordinates object."""
return deepcopy(self)
def _encode(self):
"""Return dictionary for the encoding."""
return self.copy().__dict__
@classmethod
def _decode(cls, obj_dict):
"""Decode object based on its respective ``_encode`` counterpart."""
obj = cls()
obj.__dict__.update(obj_dict)
return obj
@staticmethod
def _systems():
"""
Get class internal information about all coordinate systems.
Returns
-------
_systems : nested dictionary
List all available coordinate systems.
Key 0 - domain, e.g., 'cart'
Key 1 - convention, e.g., 'right'
Key 2a - 'short_description': string
Key 2b - 'coordinates': ['coordinate_1',
'coordinate_2',
'coordinate_3']
Key 2c - 'units': [['unit_1.1','unit_2.1','unit_3.1'],
...
['unit_1.N','unit_2.N','unit_3.N']]
Key 2d - 'description'
Key 2e - 'front': positive x (for debugging, meters and radians)
Key 2f - 'left' : positive y (for debugging, meters and radians)
Key 2g - 'back' : negative x (for debugging, meters and radians)
Key 2h - 'right': negative y (for debugging, meters and radians)
Key 2i - 'up' : positive z (for debugging, meters and radians)
Key 2j - 'down' : negative z (for debugging, meters and radians)
Key 2k,l,m - coordinate_1,2,3 : [type, [lower_lim, upper_lim]]
type can be 'unbound', 'bound', or 'cyclic'
"""
# define coordinate systems
_systems = {
"cart": {
"right": {
"description_short":
"Right handed cartesian coordinate system.",
"coordinates":
["x", "y", "z"],
"units":
[["meters", "meters", "meters"]],
"description":
"Right handed cartesian coordinate system with x,y, "
"and z in meters.",
"positive_x": [1, 0, 0],
"positive_y": [0, 1, 0],
"negative_x": [-1, 0, 0],
"negative_y": [0, -1, 0],
"positive_z": [0, 0, 1],
"negative_z": [0, 0, -1],
"x": ["unbound", [-np.inf, np.inf]],
"y": ["unbound", [-np.inf, np.inf]],
"z": ["unbound", [-np.inf, np.inf]]}
},
"sph": {
"top_colat": {
"description_short":
"Spherical coordinate system with North and South "
"Pole.",
"coordinates":
["azimuth", "colatitude", "radius"],
"units":
[["radians", "radians", "meters"],
["degrees", "degrees", "meters"]],
"description":
"The azimuth denotes the counter clockwise angle in "
"the x/y-plane with 0 pointing in positive x-"
"direction and pi/2 in positive y-direction. The "
"colatitude denotes the angle downwards from the z-"
"axis with 0 pointing in positve z-direction and pi "
"in negative z-direction. The azimuth and colatitude "
"can be in radians or degrees, the radius is always "
"in meters.",
"positive_x": [0, np.pi / 2, 1],
"positive_y": [np.pi / 2, np.pi / 2, 1],
"negative_x": [np.pi, np.pi / 2, 1],
"negative_y": [3 * np.pi / 2, np.pi / 2, 1],
"positive_z": [0, 0, 1],
"negative_z": [0, np.pi, 1],
"azimuth": ["cyclic", [0, 2 * np.pi]],
"colatitude": ["bound", [0, np.pi]],
"radius": ["bound", [0, np.inf]]},
"top_elev": {
"description_short":
"Spherical coordinate system with North and South "
"Pole. Conform with AES69-2015: AES standard for file "
"exchange - Spatial acoustic data file format (SOFA).",
"coordinates":
["azimuth", "elevation", "radius"],
"units":
[["radians", "radians", "meters"],
["degrees", "degrees", "meters"]],
"description":
"The azimuth denotes the counter clockwise angle in "
"the x/y-plane with 0 pointing in positive x-"
"direction and pi/2 in positive y-direction. The "
"elevation denotes the angle upwards and downwards "
"from the x/y-plane with pi/2 pointing at positive "
"z-direction and -pi/2 pointing in negative z-"
"direction. The azimuth and elevation can be in "
"radians or degrees, the radius is always in meters.",
"positive_x": [0, 0, 1],
"positive_y": [np.pi / 2, 0, 1],
"negative_x": [np.pi, 0, 1],
"negative_y": [3 * np.pi / 2, 0, 1],
"positive_z": [0, np.pi / 2, 1],
"negative_z": [0, -np.pi / 2, 1],
"azimuth": ["cyclic", [0, 2 * np.pi]],
"elevation": ["bound", [-np.pi / 2, np.pi / 2]],
"radius": ["bound", [0, np.inf]]},
"side": {
"description_short":
"Spherical coordinate system with poles on the "
"y-axis.",
"coordinates":
["lateral", "polar", "radius"],
"units":
[["radians", "radians", "meters"],
["degrees", "degrees", "meters"]],
"description":
"The lateral angle denotes the angle in the x/y-plane "
"with pi/2 pointing in positive y-direction and -pi/2 "
"in negative y-direction. The polar angle denotes the "
"angle in the x/z-plane with -pi/2 pointing in "
"negative z-direction, 0 in positive x-direction, "
"pi/2 in positive z-direction, pi in negative x-"
"direction. The polar and lateral angle can be in "
"radians and degree, the radius is always in meters.",
"positive_x": [0, 0, 1],
"positive_y": [np.pi / 2, 0, 1],
"negative_x": [0, np.pi, 1],
"negative_y": [-np.pi / 2, 0, 1],
"positive_z": [0, np.pi / 2, 1],
"negative_z": [0, -np.pi / 2, 1],
"lateral": ["bound", [-np.pi / 2, np.pi / 2]],
"polar": ["cyclic", [-np.pi / 2, np.pi * 3 / 2]],
"radius": ["bound", [0, np.inf]]},
"front": {
"description_short":
"Spherical coordinate system with poles on the x-axis."
" Conform with AES56-2008 (r2019): AES standard on "
"acoustics - Sound source modeling.",
"coordinates":
["phi", "theta", "radius"],
"units":
[["radians", "radians", "meters"],
["degrees", "degrees", "meters"]],
"description":
"Phi denotes the angle in the y/z-plane with 0 "
"pointing in positive y-direction, pi/2 in positive "
"z-direction, pi in negative y-direction, and 3*pi/2 "
"in negative z-direction. Theta denotes the angle "
"measured from the x-axis with 0 pointing in positve "
"x-direction and pi in negative x-direction. Phi and "
"theta can be in radians and degrees, the radius is "
"always in meters.",
"positive_x": [0, 0, 1],
"positive_y": [0, np.pi / 2, 1],
"negative_x": [0, np.pi, 1],
"negative_y": [np.pi, np.pi / 2, 1],
"positive_z": [np.pi / 2, np.pi / 2, 1],
"negative_z": [3 * np.pi / 2, np.pi / 2, 1],
"phi": ["cyclic", [0, 2 * np.pi]],
"theta": ["bound", [0, np.pi]],
"radius": ["bound", [0, np.inf]]}
},
"cyl": {
"top": {
"description_short":
"Cylindrical coordinate system along the z-axis.",
"coordinates":
["azimuth", "z", "radius_z"],
"units":
[["radians", "meters", "meters"],
["degrees", "meters", "meters"]],
"description":
"The azimuth denotes the counter clockwise angle in "
"the x/y-plane with 0 pointing in positive x-"
"direction and pi/2 in positive y-direction. The "
"heigt is given by z, and radius_z denotes the radius "
"measured orthogonal to the z-axis.",
"positive_x": [0, 0, 1],
"positive_y": [np.pi / 2, 0, 1],
"negative_x": [np.pi, 0, 1],
"negative_y": [3 * np.pi / 2, 0, 1],
"positive_z": [0, 1, 0],
"negative_z": [0, -1, 0],
"azimuth": ["cyclic", [0, 2 * np.pi]],
"z": ["unbound", [-np.inf, np.inf]],
"radius_z": ["bound", [0, np.inf]]}
}
}
return _systems
def _exist_system(self, domain=None, convention=None, unit=None):
"""
Check if a coordinate system exists and throw an error if it does not.
The coordinate systems are defined in self._systems.
Parameters
----------
domain : string
Sepcify the domain of the coordinate system, e.g., 'cart'.
convention : string
The convention of the coordinate system, e.g., 'top_colat'
units: string
The unit of the coordinate system (rad, deg, or met for radians,
degrees, or meters)
"""
if domain is None:
raise ValueError('The domain must be specified')
# get available coordinate systems
systems = self._systems()
# check if domain exists
assert domain in systems or domain is None, \
f"{domain} does not exist. Domain must be one of the follwing: "\
f"{', '.join(list(systems))}."
# check if convention exisits in domain
if convention is not None:
assert convention in systems[domain] or convention is None,\
f"{convention} does not exist in {domain}. Convention must "\
f"be one of the following: {', '.join(list(systems[domain]))}."
# check if units exist
if unit is not None:
# get first convention in domain
# (units are the same within a domain)
if convention is None:
convention = list(systems[domain])[0]
cur_units = [u[0][0:3] for u in
systems[domain][convention]['units']]
assert unit in cur_units, \
f"{unit} does not exist in {domain} convention "\
f"Unit must be one of the following: {', '.join(cur_units)}."
def _exist_coordinate(self, coordinate, unit):
"""
Check if coordinate and unit exist.
Returns domain and convention, and the index of coordinate if
coordinate and unit exists and raises a value error otherwise.
"""
# get all systems
systems = self._systems()
# find coordinate and unit in systems
for domain in systems:
for convention in systems[domain]:
if coordinate in systems[domain][convention]['coordinates']:
# get position of coordinate
index = systems[domain][convention]['coordinates'].\
index(coordinate)
# get possible units
units = [u[index][0:3] for u in
systems[domain][convention]['units']]
# return or raise ValueError
if unit in units:
return domain, convention, index
raise ValueError(
(f"'{coordinate}' in '{unit}' does not exist. See "
"self.systems() for a list of possible "
"coordinates and units"))
def _make_system(self, domain=None, convention=None, unit=None):
"""
Make and return class internal information about coordinate system.
"""
# check if coordinate system exists
self._exist_system(domain, convention, unit)
# get the new system
system = self._systems()
if convention is None:
convention = list(system[domain])[0]
system = system[domain][convention]
# get the units
if unit is not None:
units = [units for units in system['units']
if unit == units[0][0:3]]
units = units[0]
else:
units = system['units'][0]
unit = units[0][0:3]
# add class internal keys
system['domain'] = domain
system['convention'] = convention
system['unit'] = unit
system['units'] = units
return system
def _return_system(self, pts1, pts2, pts3, new_system, convert):
if convert:
# set the new system
self._system = new_system
# return points with conversion
self._set_points(pts1, pts2, pts3, True)
return self._points
else:
# return points without conversion
return self._set_points(pts1, pts2, pts3, system=new_system)
def _set_points(self, points_1, points_2, points_3,
convert=False, system=None):
"""
Check points and convert to matrix.
Parameters
----------
convert : boolean, optional
Set self._points if convert = True. Return points as
matrix otherwise. The fefault is False.
system: dict, optional
The coordinate system against which the range of the points are
checked as returned from self._make_system. If system = None
self._system is used.
Set self._points, which is an atleast_2d numpy array of shape
[L,M,...,N, 3].
"""
# get coordinate system against the points are checked
if system is None:
system = self._system
# cast to numpy array
pts = [np.atleast_2d(np.asarray(points_1, dtype=np.float64)),
np.atleast_2d(np.asarray(points_2, dtype=np.float64)),
np.atleast_2d(np.asarray(points_3, dtype=np.float64))]
# transpose
for nn, p in enumerate(pts):
if len(p.shape) == 2 and p.shape[0] == 1:
pts[nn] = np.transpose(p)
# shapes of non scalar entries
shapes = [p.shape for p in pts if p.shape != (1, 1)]
# check for equal shape
for nn in range(1, len(shapes)):
assert shapes[0] == shapes[nn], \
"points_1, points_2, and points_3 must be scalar or of the \
same shape."
# check the range of points
for nn, p in enumerate(pts):
# get type and range
c = system['coordinates'][nn]
c_type = system[c][0]
c_range = np.array(system[c][1])
# range to degrees
if system['units'][nn] == 'degrees':
c_range = np.round(c_range / np.pi * 180)
# check bounds (cyclic values could be wraped but this is safer)
if c_type in ['bound', 'cyclic']:
assert ((p >= c_range[0]) & (p <= c_range[1])).all(),\
f"Values of points_{nn+1} must be in the range \
{c_range}"
# repeat scalar entries if non-scalars exists
if len(shapes):
for nn, p in enumerate(pts):
if p.size == 1:
pts[nn] = np.tile(p, shapes[0])
# axis for concatenation
if shapes:
axis = len(shapes[0]) if shapes[0][-1] > 1 else len(shapes[0]) - 1
else:
axis = 1
# create axis for concatenation if it does not exist
for nn, p in enumerate(pts):
if p.ndim == axis:
pts[nn] = p[..., np.newaxis]
# concatenate
pts = np.concatenate((pts[0], pts[1], pts[2]), axis)
# remove noise below eps
eps = np.finfo(np.float64).eps
pts[np.abs(pts) < eps] = 0
if convert:
# save to class variable
self._points = pts
else:
return pts
def _set_weights(self, weights):
"""
Check and set sampling weights.
Set self._weights, which is an atleast_1d numpy array of shape
[L,M,...,N].
"""
# check input
if weights is None:
self._weights = weights
return
# cast to np.array
weights = np.asarray(weights, dtype=np.float64)
# reshape according to self._points
assert weights.size == self.csize,\
"weights must have same size as self.csize"
weights = weights.reshape(self.cshape)
# set class variable
self._weights = weights
def _find_nearest(self, points_1, points_2, points_3,
domain, convention, unit, show,
value, measure, atol=1e-15, radius=None):
# get KDTree
kdtree = self._make_kdtree()
# get target point in cartesian coordinates
coords = Coordinates(points_1, points_2, points_3,
domain, convention, unit)
points = coords.get_cart()
# querry nearest neighbors
points = points.flatten() if coords.csize == 1 else points
# get the points depending on measure and value
if measure == 'k':
# nearest points
distance, index = kdtree.query(points, k=value)
elif measure == 'cart':
# points within euclidean distance
index = kdtree.query_ball_point(points, value + atol)
distance = None
elif measure == 'sph':
# convert great circle to euclidean distance
x, y, z = sph2cart([0, value / 180 * np.pi],
[np.pi / 2, np.pi / 2],
[radius, radius])
value = np.sqrt((x[0] - x[1])**2
+ (y[0] - y[1])**2
+ (z[0] - z[1])**2)
# points within great circle distance
index = kdtree.query_ball_point(points, value + atol)
distance = None
# mask for scatter plot
mask = np.full((self.csize), False)
mask[index] = True
mask = mask.reshape(self.cshape)
# plot all and returned points
if show:
self.show(mask)
return distance, index, mask
def _make_kdtree(self):
"""Make a numpy KDTree for fast search of nearest points."""
xyz = self.get_cart()
kdtree = cKDTree(xyz.reshape((self.csize, 3)))
return kdtree
def __getitem__(self, index):
"""Return copied slice of Coordinates object at index."""
new = self.copy()
# slice points
new._points = np.atleast_2d(new._points[index])
# slice weights
if new._weights is not None:
new._weights = new._weights[index]
return new
def __array__(self):
"""Instances of Coordinates behave like `numpy.ndarray`, array_like."""
# make the new system
new_system = self._make_system('cart', convention='right', unit='met')
if self._system == new_system:
return self.get_cart()
# copy to avoid changing the coordinate system of the original object
return self.copy().get_cart()
def __repr__(self):
"""Get info about Coordinates object."""
# object type
if self.cshape != (0,):
obj = f"{self.cdim}D Coordinates object with {self.csize} points "\
f"of cshape {self.cshape}"
else:
obj = "Empty Coordinates object"
# coordinate convention
conv = "domain: {}, convention: {}, unit: {}".format(
self._system['domain'], self._system['convention'],
self._system['unit'])
# coordinates and units
coords = ["{} in {}".format(c, u) for c, u in
zip(self._system['coordinates'], self._system['units'])]
# join information
_repr = obj + "\n" + conv + "\n" + "coordinates: " + ", ".join(coords)
# check for sampling weights
if self._weights is None:
_repr += "\nDoes not contain sampling weights"
else:
_repr += "\nContains sampling weights"
# check for sh_order
if self._sh_order is not None:
_repr += f"\nSpherical harmonic order: {self._sh_order}"
# check for comment
if self._comment is not None:
_repr += f"\nComment: {self._comment}"
return _repr
def __eq__(self, other):
"""Check for equality of two objects."""
return not deepdiff.DeepDiff(self, other)
[docs]def cart2sph(x, y, z):
"""
Transforms from Cartesian to spherical coordinates.
Spherical coordinates follow the common convention in Physics/Mathematics.
The `colatitude` is measured downwards from the z-axis and is 0 at the
North Pole and pi at the South Pole. The `azimuth` is 0 at positive
x-direction and pi/2 at positive y-direction (counter clockwise rotation).
Cartesian coordinates follow the right hand rule.
.. math::
azimuth &= \\arctan(\\frac{y}{x}),
colatitude &= \\arccos(\\frac{z}{r}),
radius &= \\sqrt{x^2 + y^2 + z^2}
.. math::
0 < azimuth < 2 \\pi,
0 < colatitude < \\pi
Parameters
----------
x : numpy array, number
x values
y : numpy array, number
y values
z : numpy array, number
z values
Returns
-------
azimuth : numpy array, number
azimuth values
colatitude : numpy array, number
colatitude values
radius : numpy array, number
radii
Notes
-----
To ensure proper handling of the azimuth angle, the ``arctan2``
implementation from numpy is used.
"""
radius = np.sqrt(x**2 + y**2 + z**2)
z_div_r = np.where(radius != 0, z / radius, 0)
colatitude = np.arccos(z_div_r)
azimuth = np.mod(np.arctan2(y, x), 2 * np.pi)
return azimuth, colatitude, radius
[docs]def sph2cart(azimuth, colatitude, radius):
"""
Transforms from spherical to Cartesian coordinates.
Spherical coordinates follow the common convention in Physics/Mathematics.
The `colatitude` is measured downwards from the z-axis and is 0 at the
North Pole and pi at the South Pole. The `azimuth` is 0 at positive
x-direction and pi/2 at positive y-direction (counter clockwise rotation).
Cartesian coordinates follow the right hand rule.
.. math::
x &= radius \\cdot \\sin(colatitude) \\cdot \\cos(azimuth),
y &= radius \\cdot \\sin(colatitude) \\cdot \\sin(azimuth),
z &= radius \\cdot \\cos(colatitude)
.. math::
0 < azimuth < 2 \\pi
0 < colatitude < \\pi
Parameters
----------
azimuth : numpy array, number
azimuth values
colatitude : numpy array, number
colatitude values
radius : numpy array, number
radii
Returns
-------
x : numpy array, number
x values
y : numpy array, number
y values
z : numpy array, number
z vales
"""
r_sin_cola = radius * np.sin(colatitude)
x = r_sin_cola * np.cos(azimuth)
y = r_sin_cola * np.sin(azimuth)
z = radius * np.cos(colatitude)
return x, y, z
[docs]def cart2cyl(x, y, z):
"""
Transforms from Cartesian to cylindrical coordinates.
Cylindrical coordinates follow the convention that the `azimuth` is 0 at
positive x-direction and pi/2 at positive y-direction (counter clockwise
rotation). The `height` is identical to the z-coordinate and the `radius`
is measured orthogonal from the z-axis.
Cartesian coordinates follow the right hand rule.
.. math::
azimuth &= \\arctan(\\frac{y}{x}),
height &= z,
radius &= \\sqrt{x^2 + y^2},
.. math::
0 < azimuth < 2 \\pi
Parameters
----------
x : numpy array, number
x values
y : numpy array, number
y values
z : numpy array, number
z values
Returns
-------
azimuth : numpy array, number
azimuth values
height : numpy array, number
height values
radius : numpy array, number
radii
Notes
-----
To ensure proper handling of the azimuth angle, the ``arctan2``
implementation from numpy is used.
"""
azimuth = np.mod(np.arctan2(y, x), 2 * np.pi)
if isinstance(z, np.ndarray):
height = z.copy()
else:
height = z
radius = np.sqrt(x**2 + y**2)
return azimuth, height, radius
[docs]def cyl2cart(azimuth, height, radius):
"""
Transforms from cylindrical to Cartesian coordinates.
Cylindrical coordinates follow the convention that the `azimuth` is 0 at
positive x-direction and pi/2 at positive y-direction (counter clockwise
rotation). The `height` is identical to the z-coordinate and the `radius`
is measured orthogonal from the z-axis.
Cartesian coordinates follow the right hand rule.
.. math::
x &= radius \\cdot \\cos(azimuth),
y &= radius \\cdot \\sin(azimuth),
z &= height
.. math::
0 < azimuth < 2 \\pi
Parameters
----------
azimuth : numpy array, number
azimuth values
height : numpy array, number
height values
radius : numpy array, number
radii
Returns
-------
x : numpy array, number
x values
y : numpy array, number
y values
z : numpy array, number
z values
Notes
-----
To ensure proper handling of the azimuth angle, the ``arctan2``
implementation from numpy is used.
"""
x = radius * np.cos(azimuth)
y = radius * np.sin(azimuth)
if isinstance(height, np.ndarray):
z = height.copy()
else:
z = height
return x, y, z