"""In house numpy helpers.
.. codeauthor:: Yoshifumi Futaana
The aim of this module is to extend the numpy functionality.
.. autosummary::
lingrids
loggrids
min_where
max_where
shrink
squeeze_right
squeeze_left
isascending
isdescending
guess_upper
guess_lower
guess_interval
padlin
sort_periodic
unjump_periodic
span
broadcast
broadcast_r
rtp
lonlat2xyz
trim_epsilon
nearest_index
nearest_index_right
nearest_index_left
"""
import bisect as _bisect
from irfpy.util import exception as ex
import numpy as np
import numpy as _np
[docs]def min_where(arr, cond):
""" Conditional minimum.
Minimum is calculated but under the specific conditions.
For example, it can be used when you want to know the minimum of the given array
but above zero.
>>> arr = np.array([[0, 0, 7], [3, 2, 4]])
>>> print(min_where(arr, arr > 0))
2
:param arr: Numpy array.
:param cond: Condition. Should be shape shape as ``arr``.
:returns: Return the minimum of arr, but where cond is satisfied.
If all the cond is ``False``, ``nan`` is returned.
>>> print(min_where(arr, arr > 1000))
nan
"""
if arr.shape != cond.shape:
raise ValueError('Shape mismatch.')
arr_cond = np.ma.masked_array(arr, mask=(~cond))
minimum_value = arr_cond.min()
try:
minimum_value = float(minimum_value.filled(np.nan)) # In case the value minimum_value is masked.
except AttributeError as e:
pass
return minimum_value
[docs]def max_where(arr, cond):
""" Conditional maximum, similar to :func:`min_where`.
:param arr: Numpy array.
:param cond: Condition. Should be shape shape as ``arr``.
:returns: Return the minimum of arr, but where cond is satisfied.
If all the cond is ``False``, ``nan`` is returned.
>>> arr = np.array([[0, 0, 7], [3, 2, 4]])
>>> print(max_where(arr, arr < 4))
3
"""
if arr.shape != cond.shape:
raise ValueError('Shape mismatch.')
arr_cond = np.ma.masked_array(arr, mask=(~cond))
max_value = arr_cond.max()
try:
max_value = float(max_value.filled(np.nan)) # In case the value max_value is masked.
except AttributeError as e:
pass # This is nominal case. Usually max_value is float or int, not masked scalar.
return max_value
[docs]def shrink(array, newshape, operation='mean'):
""" The given array is shrinked in size with given operation.
:param array: Numpy array. Multidimensional.
:param newshape: New shape. The new shape must have the same number of array dimensions.
Number of each new axis should be the divisor of the corresponding array dimension.
:keyword operation: 'sum', 'mean', 'max', 'min', 'median'
.. note::
Keep in mind, "median" is not multidimensional value, but done dimension-by-dimension
from the last dimension.
Namely, the operation is like ``array2d.median(axis=1).median(axis=0)``.
>>> print(shrink(np.arange(24).reshape(4, 6), (2, 3), operation='mean')) # doctest: +NORMALIZE_WHITESPACE
[[ 3.5 5.5 7.5]
[15.5 17.5 19.5]]
>>> print(shrink(np.arange(24).reshape(4, 6), (2, 3), operation='max'))
[[ 7 9 11]
[19 21 23]]
"""
if operation not in ['sum', 'mean', 'max', 'min', 'median']:
raise ValueError("Operation ``{}`` is not supported.".format(operation))
if array.ndim != len(newshape):
raise ValueError("Dimensions should be kept (Original={}/New={})".format(
array.shape, newshape))
shape_num = [(d, c // d) for d, c in zip(newshape, array.shape)]
shape_num_index = [shape for num in shape_num for shape in num]
try:
array = array.reshape(shape_num_index)
except ValueError as e:
raise ValueError('Given shape is not compatible for shrinking. Original={}/New={}'.format(
str(array.shape), str(newshape)))
# Operation starts for each dimension.
for i in range(len(newshape)):
if operation == "sum":
array = array.sum(-1 * (i + 1))
elif operation == "mean":
array = array.mean(-1 * (i + 1))
elif operation == "max":
array = array.max(axis=-1 * (i + 1))
elif operation == "min":
array = array.min(axis=-1 * (i + 1))
elif operation == "median":
array = np.ma.median(array, axis=-1 * (i + 1))
else:
raise ValueError("Operation {} not supported.".format(operation))
return array
[docs]def squeeze_right(arr):
""" Similar to :func:`numpy:numpy.squeeze`, the given array is squeezed but only from right axes.
>>> arr = np.random.uniform(size=[1, 3, 5, 4, 2, 1])
>>> np.squeeze(arr).shape
(3, 5, 4, 2)
>>> squeeze_right(arr).shape
(1, 3, 5, 4, 2)
>>> print((arr[..., 0] == squeeze_right(arr)).all())
True
"""
shape = arr.shape
newshape = []
sq = True
for num in shape[::-1]:
if sq:
if num != 1:
sq = False
newshape.insert(0, num)
else:
pass # Ignore the dimension
else:
newshape.insert(0, num)
return arr.reshape(newshape)
[docs]def squeeze_left(arr):
""" Similar to :func:`numpy:numpy.squeeze`, the given array is squeezed but only from left axes.
>>> arr = np.random.uniform(size=[1, 3, 5, 4, 2, 1])
>>> np.squeeze(arr).shape
(3, 5, 4, 2)
>>> squeeze_left(arr).shape
(3, 5, 4, 2, 1)
>>> print((arr[0] == squeeze_left(arr)).all())
True
"""
shape = arr.shape
newshape = []
sq = True
for num in shape:
if sq:
if num != 1:
sq = False
newshape.append(num)
else:
pass # Ignore the dimension
else:
newshape.append(num)
return arr.reshape(newshape)
[docs]def padlin(vector, pad_width, iaxis, kwargs):
""" A helper function to use :func:`numpy:numpy.pad`.
Use :func:`numpy:numpy.pad` function with padding of constant inclination over the whole range.
>>> a = [1.0, 3, 7, 5, 2] # 5 element array. Average inclication is 0.25 = (2 -1 ) / 4
>>> b = np.pad(a, (2, 3), mode=padlin) # Specify the function, padlin, as a mode for np.pad.
>>> print(b) # doctest: +NORMALIZE_WHITESPACE
[0.5 0.75 1. 3. 7. 5. 2. 2.25 2.5 2.75]
In this example, the padding before is 2 element, which increments 0.25.
At the end, 3 elements are added, using 0.25 increments after the last value of 2.
"""
incl = (vector[-pad_width[1] - 1] - vector[pad_width[0]]) / (len(vector) - pad_width[0] - pad_width[1] - 1)
vector[:pad_width[0]] = -np.arange(pad_width[0], 0, -1) * incl + vector[pad_width[0]]
vector[-pad_width[1]:] = (np.arange(pad_width[1]) + 1) * incl + vector[-pad_width[1] - 1]
return vector
[docs]def isascending(array1d, allow_equal=True):
""" Check if the array is ascending order or not
:param array1d: Array to be checked.
:keyword allow_equal: Adjacent equal values are treated as ascending.
>>> arr = [1, 2, 3, 4, 5]
>>> print(isascending(arr))
True
>>> arr = [1, 2, -3, 4, 5]
>>> print(isascending(arr))
False
>>> arr = [1, 2, 2, 5]
>>> print(isascending(arr))
True
>>> arr = [1, 2, 2, 5]
>>> print(isascending(arr, allow_equal=False))
False
"""
array1d = np.copy(array1d)
if len(array1d) == 0:
return False
elif len(array1d) == 1:
return True
darr = array1d[1:] - array1d[:-1]
if allow_equal:
return darr.min() >= 0
else:
return darr.min() > 0
[docs]def isdescending(array1d, allow_equal=True):
""" Check if the array is ascending order or not
:param array1d: Array to be checked.
:keyword allow_equal: Adjacent equal values are treated as ascending.
>>> arr = [5, 4, 3, 2, 1]
>>> print(isdescending(arr))
True
>>> arr = [5, 4, -3, 2, 1]
>>> print(isdescending(arr))
False
>>> arr = [5, 2, 2, 1]
>>> print(isdescending(arr))
True
>>> arr = [5, 2, 2, 1]
>>> print(isdescending(arr, allow_equal=False))
False
"""
if len(array1d) == 0:
return False
elif len(array1d) == 1:
return True
array1d = np.copy(array1d)[::-1]
return isascending(array1d, allow_equal=allow_equal)
def _guess_upper_but_last(array1d):
""" Return the array of upper bound, except for the last one.
"""
arr_guess = np.zeros_like(array1d) - 1e20 # Set very strange value as default.
arr_guess[:-1] = (array1d[:-1] + array1d[1:]) / 2.0
return arr_guess
def _guess_lower_but_first(array1d):
arr_guess = np.zeros_like(array1d) - 1e20 # Set very strange value as default.
arr_guess[1:] = (array1d[:-1] + array1d[1:]) / 2.0
return arr_guess
[docs]def guess_interval(array1d, log=False, method="median"):
""" From the given 1-D array, a guessed interval array is returned.
The returned value is :func:`guess_upper` - :func:`guess_lower`.
"""
return guess_upper(array1d, log=log, method=method) - guess_lower(array1d, log=log, method=method)
[docs]def guess_upper(array1d, log=False, method="median"):
""" From the given 1D array, a guessed upper bound array is returned.
For the given ascending array, the upper bound is calculate from the mean with the next element.
For example, the array (2, 3, 5, 7, 11, 13, 17), the upper bound array is
(2.5, 4, 6, 9, 12, 15, *19*). All the elements but the last one is the mean with the next element.
The last one is guessed with the given *method*; by default, the median of the difference of the other element.
>>> arr = np.array([2., 3, 5, 7, 11, 13, 17])
>>> arr_upper = guess_upper(arr)
>>> print(arr_upper) # doctest: +NORMALIZE_WHITESPACE
[ 2.5 4. 6. 9. 12. 15. 18. ]
The array is assumed to be sorted ascendingly, but no check is done.
:param array1d: Ascending np.array.
:param log: Set *True* to calculate this in the log space.
:return: 1D array.
You can get this functionality in log space.
>>> arr = [1, 10, 100, 1000, 10000]
In this case, the upper value is the "Geometric mean".
>>> arr_upper = guess_upper(arr, log=True)
>>> print(arr_upper) # doctest: +NORMALIZE_WHITESPACE
[3.16227766e+00 3.16227766e+01 3.16227766e+02 3.16227766e+03 3.16227766e+04]
"""
arr = np.copy(array1d)
if log:
arr = np.log(arr)
### Get upper bound, but for the last element.
arr_guess = _guess_upper_but_last(arr)
meth = method.lower()
if meth in ('median',):
diffarr = arr_guess - arr
meandiff = np.median(diffarr[:-1])
arr_guess[-1] = arr[-1] + meandiff
### elif meth in ('average',): # will follow in the future..
else:
raise ex.IrfpyException('Only method="median" is supported yet')
if log:
arr_guess = np.exp(arr_guess)
return arr_guess
[docs]def guess_lower(array1d, log=False, method='median'):
""" From the given 1D array, a guessed lower bound array is returned.
See the similar function :func:`guess_upper` for details.
>>> arr = [2, 3, 5, 7, 11, 13, 17]
>>> lower = guess_lower(arr)
>>> print(lower) # doctest: +NORMALIZE_WHITESPACE
[ 1. 2.5 4. 6. 9. 12. 15. ]
>>> arr = [1, 10, 100, 1000, 10000]
>>> arr_lower = guess_lower(arr, log=True)
>>> print(arr_lower) # doctest: +NORMALIZE_WHITESPACE
[3.16227766e-01 3.16227766e+00 3.16227766e+01 3.16227766e+02 3.16227766e+03]
"""
arr = np.copy(array1d)
if log:
arr = np.log(arr)
### Get lower bound, but for the first element.
arr_guess = _guess_lower_but_first(arr)
meth = method.lower()
if meth in ('median',):
diffarr = arr_guess - arr
meandiff = np.median(diffarr[1:])
arr_guess[0] = arr[0] + meandiff
### elif meth in ('average',): # will follow in the future..
else:
raise ex.IrfpyException('Only method="median" is supported yet')
if log:
arr_guess = np.exp(arr_guess)
return arr_guess
[docs]def guess_bounds(array1d, log=False, method='median'):
""" From the given 1D array, a guessed bound array is returned.
From the given 1D array, a guessed bound array is returned.
Bound array is an array specifying the bounds of the elements of the original array.
For example, suppose we have an array ``[0, 1, 2, 3]``.
The corresponding bound array is ``[-0.5, 0.5, 1.5, 2.5, 3.5]``.
The 1st element (``0``) of the original array's bound is ``-0.5`` and ``0.5``.
The returned array must have 1 elements more than that of the original array.
:param array1d: Array-like.
:param log: If the estimation to be done in ``log`` space, use this option as *True*.
:param method: Method to guess. Currently ``median`` is only supported.
:return: A new array, with a size of (N+1), where N is the size of the original array (``array1d``).
>>> original_array = [0, 1, 2, 3]
>>> bound_array = guess_bounds(original_array)
>>> print(bound_array)
[-0.5, 0.5, 1.5, 2.5, 3.5]
"""
arr_u = guess_upper(array1d, log=log, method=method)
arr_l = guess_lower(array1d, log=log, method=method)
arr = arr_l.tolist() + [arr_u[-1]]
return arr
[docs]def sort_periodic(array1d, period):
""" From an array in an periodic system, ascending array is obtained.
:param array1d: An array.
:param period: The period.
:return: An array.
If you have an array [0, 100, 20, 120, 40] in an ascending system
with periodicity of 180. In this case, the returned array is
[0, 100, *200*, *300*, *400*].
The original element of "20" is smaller than "100", the one element before.
This is interpreted as "the value 20 is ascending from 100, but exceed
the upper bound so that the value is decrimented by period=180. Thus,
the ascending value should be 20+180=200".
>>> arr = [0, 100, 20, 120, 40]
>>> ascend_arr = sort_periodic(arr, 180)
>>> print(ascend_arr) # doctest: +NORMALIZE_WHITESPACE
[ 0 100 200 300 400]
"""
asc_arr = np.copy(array1d)
narr = asc_arr.shape[0]
for i in range(narr - 1):
if array1d[i + 1] < array1d[i]:
asc_arr[i + 1:] += period
return asc_arr
[docs]def unjump_periodic(array1d, period):
""" A given 1-d array's large gap (>period/2) is corrected.
*Use case*
If you have a series [179, -179, 179, -179] in a periodic system (-180 to 180, say longitude),
it could be more natural to consider it as [179, 181, 179, 181] or [-181, -179, -181, -179].
*Sample*
>>> arr = [179, -179, 179, -179]
>>> arr2 = unjump_periodic(arr, 360)
>>> print(arr2)
[179 181 179 181]
>>> arr = [179, -179, 179, np.nan, -179, 179, -179, np.nan, 179]
>>> print(unjump_periodic(arr, 360))
[ 179. 181. 179. nan -179. -181. -179. nan 179.]
>>> arr = [0, 180, 40]
>>> print(unjump_periodic(arr, 360))
[ 0 -180 -320]
"""
arr = np.array(array1d, copy=True)
diff = arr[1:] - arr[:-1]
while np.nanmax(diff) >= period / 2:
idx = np.where(diff >= period / 2)[0][-1] # The last one
arr[idx+1:] = arr[idx+1:] - period
diff = arr[1:] - arr[:-1]
while np.nanmin(diff) < -period / 2:
idx = np.where(diff < -period / 2)[0][-1]
arr[idx+1:] = arr[idx+1:] + period
diff = arr[1:] - arr[:-1]
if np.nanmax(diff) >= period / 2 or np.nanmin(diff) < -period / 2:
# I think this is not needed...
return unjump_periodic(arr, period)
return arr
[docs]def span(inarray):
""" Just return the tuple of (min, max) for the given 1D array.
:param inarray: np.ndarray
:returns: A tuple, (min, max).
Equivalent to (inarray.min(), inarray.max())
>>> arr = np.array([1, 2, 5, -3, 2, 8, 10, -2])
>>> print(span(arr))
(-3, 10)
"""
return inarray.min(), inarray.max()
[docs]def broadcast(array_in, shape_out):
""" Create the broadcasted array with the specific shape from the given array
.. note::
You can get the same results using ``numpy.broadcast_to`` function.
.. code-block:: python
import numpya as np
broadcasted_array = np.broadcast_to(array_in, shape_out)
>>> array_in = np.arange(10).reshape([5, 2])
>>> arr2 = np.broadcast_to(array_in, (4, 3, 5, 2)) # Equivalent to irfpy.util.nptools.broadcast
>>> print(arr2.shape)
(4, 3, 5, 2)
>>> print(span(arr2[:, :, 0, 0]))
(0, 0)
>>> print(span(arr2[:, :, 3, 1]))
(7, 7)
The array is broadcasted to the given shape_out.
:param array_in: Numpy array. Shape of (s0, s1, ..., sS).
:param shape_out: The shape to get. The shape_out must be (sX, sY, ..., s0, s1, ..., sS).
:returns: The broadcasted array.
>>> array_in = np.arange(10).reshape([5, 2]) # (5, 2) shaped array is prepared.
>>> arr2 = broadcast(array_in, (4, 3, 5, 2)) # Broadcast to (4, 3, 5, 2) array.
>>> print(arr2.shape)
(4, 3, 5, 2)
>>> print(span(arr2[:, :, 0, 0]))
(0, 0)
>>> print(span(arr2[:, :, 3, 1]))
(7, 7)
>>> broadcast(array_in, (5, 2, 3)) # doctest: +IGNORE_EXCEPTION_DETAIL
Traceback (most recent call last):
File "<stdin>", line 1, in ?
irfpy.util.exception.IrfpyException: 'Array cannot be broadcasted. Non-broadcastable. In:(5, 2) Out:(5, 2, 3)'
>>> broadcast(array_in, (3,)) # doctest: +IGNORE_EXCEPTION_DETAIL
Traceback (most recent call last):
File "<stdin>", line 1, in ?
irfpy.util.exception.IrfpyException: 'Array cannot be broadcasted. Lack of dimension. In:(5, 2) Out:(3,)'
"""
shape_in = array_in.shape
dimen_in = array_in.ndim
shape_out = tuple(shape_out)
dimen_out = len(shape_out)
### Out array dimension is smaller than the in, raise an exception
if dimen_in > dimen_out:
raise ex.IrfpyException('Array cannot be broadcasted. Lack of dimension. In:{} Out:{}'.
format(str(shape_in), str(shape_out)))
### Not broadcastable.
if shape_out[-dimen_in:] != shape_in:
raise ex.IrfpyException('Array cannot be broadcasted. Non-broadcastable. In:{} Out:{}'.
format(str(shape_in), str(shape_out)))
### Tile function
shape_tile = np.array(shape_out)
shape_tile[-dimen_in:] = 1 # All the "broadcasting indexes" should be 1.
array_out = np.tile(array_in, shape_tile)
return array_out
[docs]def broadcast_r(array_in, shape_out):
""" Broadcast, but in the revered direction
.. note::
You can get the same results by
.. code-block:: python
broadcasted_array = array_in[..., np.newaxis, np.newaxis, np.newaxis] + np.zeros(shape_out)
Here ``np.newaxis`` should repeat the proper times to make everythng consistent.
>>> array_in = np.arange(10).reshape([5, 2]) # arrayn_in is (5, 2) shape
>>> arr2 = array_in[..., np.newaxis, np.newaxis] + np.zeros([5, 2, 10, 3], dtype=int) # Equivalent to irfpy.util.nptools.broadcast_r
>>> print(arr2.shape)
(5, 2, 10, 3)
>>> print(span(arr2[0, 0, :, :]))
(0, 0)
>>> print(span(arr2[3, 1, :, :]))
(7, 7)
This looks 10 times faster than the broadcast_r method.
.. todo::
Change the algorithm using the above algorithm
For example, an array with the shape of (5, 2) will be broadcasted to
the array with the shape of (5, 2, 10, 4, 7). New_array[i, j, :, :, :] == Old_array[i, j].
:param array_in: np array with shape of (s0, s1, ...)
:param shape_out: The aimed shape. It should be (s0, s1, ..., sX, sY, ...)
:return: The broadcasted array from the "top" of the array.
:See also: :class:`broadcast`
>>> arrin = np.arange(10).reshape([5, 2]) # (5, 2) shaped array is prepared.
>>> arrout = broadcast_r(arrin, (5, 2, 10, 3)) # (5, 2, 10, 3) array will be returned.
>>> print(arrout.shape)
(5, 2, 10, 3)
>>> print(arrout[2, 1, :, :].max())
5
>>> print(arrout[2, 1, :, :].min())
5
>>> broadcast_r(arrin, (9, 5, 2)) # doctest: +IGNORE_EXCEPTION_DETAIL
Traceback (most recent call last):
File "<stdin>", line 1, in ?
irfpy.util.exception.IrfpyException: 'Array cannot be broadcasted. Non-broadcastable. In:(5, 2) Out:(9, 5, 2)'
"""
shape_in = array_in.shape
dimen_in = array_in.ndim
shape_out = tuple(shape_out)
dimen_out = len(shape_out)
### Out array dimension is smaller than the in, raise an exception
if dimen_in > dimen_out:
raise ex.IrfpyException('Array cannot be broadcasted. Lack of dimension. In:{} Out:{}'.
format(str(shape_in), str(shape_out)))
### Check if broadcast works or not.
if shape_out[:dimen_in] != shape_in:
raise ex.IrfpyException('Array cannot be broadcasted. Non-broadcastable. In:{} Out:{}'.
format(str(shape_in), str(shape_out)))
### Tile function. Shape to give np.tile is prepared.
# shape_out is (s0, s1, ..., s(di-1), sX, sY, ...)
shape_tile = np.roll(shape_out, -dimen_in) # (sX, sY, ..., s0, s1, ...)
shape_tile[-dimen_in:] = 1 # (sX, sY, ..., 1, 1, ...)
array_out = np.tile(array_in, shape_tile) # with shape of (sX, sY, ..., s0, s1, ...)
### Transpose.
transpose_index = np.arange(dimen_out) # (0, 1, 2, ..., do-1)
transpose_index = np.roll(transpose_index, dimen_in) # (..., do-1, 0, 1, 2, ...)
array_out = np.transpose(array_out, axes=transpose_index)
return array_out
[docs]def rtp(x, y, z, *, degrees: bool=False, primary_axis=None, secondary_axis=None, latitude: bool=False):
""" Cartesian x, y, and z is converted to r, |theta|, and |phi|.
:param x: X
:param y: Y
:param z: Z
:keyword primary_axis: The primary axis, default is (0, 0, 1).
:keyword secondary_axis: The secondary axis, default is (1, 0, 0).
:keyword degrees: If *True*, returned angles are in degrees. Default *False*, i.e., radians.
:keyword latitude: If *True*, returned polar angle is latitude 90 - |theta| (ranging from -90 to 90 deg)
:return: (r, |theta|, |phi|). |theta| and |phi| in radians (or degrees when ``degrees=True``.
r is the distance. |theta| is the polar angle, i.e.,
the angle from the primary axis (default, z).
|phi| is the azimuthal angle, the angle between the secondary axis and
the vector projected
to the plane perpendicular to the primary plane (default x-y plane).
:rtype: A tuple
I don't know why but the numpy or numpy.linalg package does not provide this.
(Let me know if there is ofiicial version).
It is very simple in math so can be easily implemented, but also good to have a function.
>>> print(rtp(0, 0, 1))
(1.0, 0.0, 0.0)
>>> r, t, p = rtp(1, 1, 1)
>>> print('{:.3f} {:.3f} {:.3f}'.format(r, t, p))
1.732 0.955 0.785
Array is allowed.
>>> r, t, p = rtp([1, 0, -1, 0], [1, 1, 1, 0], [0, 1, 1, 0]) # List can be allowed.
>>> print(r.shape)
(4,)
>>> print(t) # doctest: +NORMALIZE_WHITESPACE
[1.57079633 0.78539816 0.95531662 nan]
>>> print(p)
[0.78539816 1.57079633 2.35619449 0. ]
Multi-dimension.
>>> r, t, p = rtp([[1, 0],[ -1, 0]], [[1, 1], [1, 0]], [[0, 1], [1, 0]]) # List can be allowed.
>>> print(r.shape)
(2, 2)
>>> print(t) # doctest: +NORMALIZE_WHITESPACE
[[1.57079633 0.78539816]
[0.95531662 nan]]
>>> print(p)
[[0.78539816 1.57079633]
[2.35619449 0. ]]
You can specify the primary and secondary axis.
>>> r, t, p = rtp(1, 0, 0, primary_axis=[1, 0, 0], secondary_axis=[0, 1, 0], degrees=True)
>>> print(r, t, p)
1.0 0.0 0.0
>>> r, t, p = rtp(0, 2, 0, primary_axis=[1, 0, 0], secondary_axis=[0, 1, 0], degrees=True)
>>> print(r, t, p)
2.0 90.0 0.0
>>> r, t, p = rtp(0, 0, 3, primary_axis=[1, 0, 0], secondary_axis=[0, 1, 0], degrees=True)
>>> print(r, t, p)
3.0 90.0 90.0
>>> r, t, p = rtp(0, 2, 0, primary_axis=[1, 0, 0], secondary_axis=[0, 0, 1], degrees=True)
>>> print(r, t, p)
2.0 90.0 -90.0
>>> r, t, p = rtp([1, 0], [0, 0], [0, 3], primary_axis=[1, 0, 0], secondary_axis=[0, 1, 0], degrees=True)
>>> print(r, t, p) # doctest: +NORMALIZE_WHITESPACE
[1. 3.] [ 0. 90.] [ 0. 90.]
You can set *latitude* keyword to obtain the latitude.
>>> r, t, p = rtp(2, 0, 0, degrees=True)
>>> print(t)
90.0
>>> r, lat, lon = rtp(2, 0, 0, degrees=True, latitude=True)
>>> print(lat)
0.0
>>> r, lat, lon = rtp(2, 0, -2, degrees=True, latitude=True)
>>> print('{:.1f}'.format(lat))
-45.0
"""
if primary_axis is None:
primary_axis = _np.array([0, 0, 1])
if secondary_axis is None:
secondary_axis = _np.array([1, 0, 0])
prima = _np.array(primary_axis)
secon = _np.array(secondary_axis)
prima = prima / _np.linalg.norm(prima)
third = _np.cross(prima, secon)
third = third / _np.linalg.norm(third)
secon = _np.cross(third, prima)
secon = secon / _np.linalg.norm(secon)
vec = np.array([x, y, z]) # (3,) or (3, ...) shape
xx = _np.tensordot(secon, vec, axes=(0, 0))
yy = _np.tensordot(third, vec, axes=(0, 0))
zz = _np.tensordot(prima, vec, axes=(0, 0))
r = np.sqrt(np.array(xx ** 2 + yy ** 2 + zz ** 2))
if latitude:
theta = np.arcsin(np.clip(zz / r, -1, 1))
else:
theta = np.arccos(np.clip(zz / r, -1, 1))
phi = np.arctan2(yy, xx)
if degrees:
theta = _np.degrees(theta)
phi = _np.degrees(phi)
return (r, theta, phi)
[docs]def lonlat2xyz(longitude_deg, latitude_deg):
""" Simply convert longitude and latitude to array
:param longitude_deg: Longitude in degrees
:param latitude_deg: Latitude in degrees
:returns: Array with (..., 3) shape
>>> lon = 0
>>> lat = 0
>>> print(lonlat2xyz(lon, lat))
[1. 0. 0.]
>>> lon = [0, 0, 0, 90, 90, 90]
>>> lat = [-90, 0, 90, -90, 0, 90]
>>> xyz = lonlat2xyz(lon, lat)
>>> print(xyz.shape)
(6, 3)
>>> print(trim_epsilon(xyz))
[[ 0. 0. -1.]
[ 1. 0. 0.]
[ 0. 0. 1.]
[ 0. 0. -1.]
[ 0. 1. 0.]
[ 0. 0. 1.]]
"""
lonr = np.deg2rad(longitude_deg)
latr = np.deg2rad(latitude_deg)
x = np.cos(lonr) * np.cos(latr)
y = np.sin(lonr) * np.cos(latr)
z = np.sin(latr)
vec = np.array([x, y, z])
vec = np.moveaxis(vec, 0, -1)
return vec
[docs]def lingrids(x0, x1, nx):
""" Return the grids.
:param x0: x0
:param x1: x1
:param nx: Number of bins
:return: A tuple, ``(grid_c, grid_b, grid_d)``. grid_c is the center of the grid with (nx) shape,
and grid_b is the bound of the grid with (nx+1) shape. grid_d is the width of the bin with (nx) shape.
- ``grid_b[0] = x0 and grid_b[nx] = x1``.
- ``grid_c = (grid_b[1:] + grid_b[:-1]) / 2.``
- ``grid_b = np.linspace(x0, x1, nx+1)``.
- ``grid_d = (grid_b[1:] - grid_b[:-1])``.
>>> grid_c, grid_b, grid_d = lingrids(0, 10, 5)
>>> print(grid_b[0])
0.0
>>> print(grid_b[5])
10.0
>>> print(grid_c[0])
1.0
>>> print(grid_d[3]) # grid_d is always 2 in this case
2.0
"""
grid_b = np.linspace(x0, x1, nx + 1)
grid_c = (grid_b[1:] + grid_b[:-1]) / 2.
grid_d = grid_b[1:] - grid_b[:-1]
return grid_c, grid_b, grid_d
[docs]def loggrids(x0, x1, nx):
""" Return the grid.
:param x0: x0
:param x1: x1
:param nx: Number of bins
:return: A tuple, ``(grid_c, grid_b, grid_d)``. grid_c is the center of the grid with (nx) shape,
and grid_b is the bound of the grid with (nx+1) shape. grid_d is the width of the bin with (nx) shape.
- ``grid_b[0] = 10^x0 and grid_b[nx] = 10^x1``.
- ``grid_c = sqrt(grid_b[1:] * grid_b[:-1])``
- ``grid_b = np.linspace(x0, x1, nx+1)``.
- ``grid_d = (grid_b[1:] - grid_b[:-1])``.
>>> grid_c, grid_b, grid_d = loggrids(0, 10, 5)
>>> print(grid_b[0])
1.0
>>> print('{:.2e}'.format(grid_b[5]))
1.00e+10
>>> print('{:.2e}'.format(grid_c[0]))
1.00e+01
>>> print('{:.2e}'.format(grid_d[3])) # 10^8 - 10^6
9.90e+07
"""
grid_b = np.logspace(x0, x1, nx + 1)
grid_c = np.sqrt(grid_b[1:] * grid_b[:-1])
grid_d = grid_b[1:] - grid_b[:-1]
return grid_c, grid_b, grid_d
[docs]def trim_epsilon(arr, epsilon=1e-15):
r""" Trim the epsilon.
:param arr: A numpy array
:keyword epsilon: |epsilon|
:return: The input ``arr`` but with substituting zero (+0.0) for the components
where the absolute value is less than epsilon
.. math::
b_i = \left\{
\begin{array}{ll}
a_i & (|a_i| \ge \epsilon) \\
0 & (|a_i| < \epsilon)
\end{array}\right.
"""
return np.where(np.abs(arr) < epsilon, 0, arr)
[docs]def nearest_index(a, x):
""" Return the index where the x is nearest to in array a.
:param a: A list. It should be sorted in ascending order.
:param x: Data
:returns: The index of the data in the array ``a``
with the nearest data to ``x``.
In other words, the index ``i`` usually satisfies
``a[i] < x <= a[i+1]``.
>>> xlist = [100, 200, 400, 800, 1600]
>>> print(nearest_index(xlist, 50))
0
>>> print(nearest_index(xlist, 100))
0
>>> print(nearest_index(xlist, 149))
0
>>> print(nearest_index(xlist, 150))
1
>>> print(nearest_index(xlist, 151))
1
>>> print(nearest_index(xlist, 200))
1
>>> print(nearest_index(xlist, 300))
2
>>> print(nearest_index(xlist, 400))
2
>>> print(nearest_index(xlist, 1600))
4
>>> print(nearest_index(xlist, 1800))
4
"""
return nearest_index_right(a, x)
[docs]def nearest_index_right(xlist, x0):
""" Same as :func:`nearest_index`
"""
idxbi = _bisect.bisect(xlist, x0)
if idxbi == 0:
return idxbi
if idxbi == len(xlist):
return idxbi - 1
xm = x0 - xlist[idxbi - 1]
xp = xlist[idxbi] - x0
if xm < xp:
return idxbi - 1
else:
return idxbi
[docs]def nearest_index_left(xlist, x0):
""" Return the index where the x is nearest to in array a.
:param a: A list. It should be sorted in ascending order.
:param x: Data
:returns: The index of the data in the array ``a``
with the nearest data to ``x``.
In other words, the index ``i`` usually satisfies
``a[i] < x <= a[i+1]``.
>>> xlist = [100, 200, 400, 800, 1600]
>>> print(nearest_index_left(xlist, 50))
0
>>> print(nearest_index_left(xlist, 100))
0
>>> print(nearest_index_left(xlist, 149))
0
>>> print(nearest_index_left(xlist, 150))
0
>>> print(nearest_index_left(xlist, 151))
1
>>> print(nearest_index_left(xlist, 200))
1
>>> print(nearest_index_left(xlist, 300))
1
>>> print(nearest_index_left(xlist, 400))
2
>>> print(nearest_index_left(xlist, 1600))
4
>>> print(nearest_index_left(xlist, 1800))
4
"""
idxbi = _bisect.bisect(xlist, x0)
if idxbi == 0:
return idxbi
if idxbi == len(xlist):
return idxbi - 1
xm = x0 - xlist[idxbi - 1]
xp = xlist[idxbi] - x0
if xm <= xp:
return idxbi - 1
else:
return idxbi