irfpy.util.nptools

In house numpy helpers.

Code author: Yoshifumi Futaana

The aim of this module is to extend the numpy functionality.

lingrids(x0, x1, nx)	Return the grids.
loggrids(x0, x1, nx)	Return the grid.
min_where(arr, cond)	Conditional minimum.
max_where(arr, cond)	Conditional maximum, similar to min_where().
shrink(array, newshape[, operation])	The given array is shrinked in size with given operation.
squeeze_right(arr)	Similar to numpy:numpy.squeeze(), the given array is squeezed but only from right axes.
squeeze_left(arr)	Similar to numpy:numpy.squeeze(), the given array is squeezed but only from left axes.
isascending(array1d[, allow_equal])	Check if the array is ascending order or not
isdescending(array1d[, allow_equal])	Check if the array is ascending order or not
guess_upper(array1d[, log, method])	From the given 1D array, a guessed upper bound array is returned.
guess_lower(array1d[, log, method])	From the given 1D array, a guessed lower bound array is returned.
guess_interval(array1d[, log, method])	From the given 1-D array, a guessed interval array is returned.
padlin(vector, pad_width, iaxis, kwargs)	A helper function to use numpy:numpy.pad().
sort_periodic(array1d, period)	From an array in an periodic system, ascending array is obtained.
unjump_periodic(array1d, period)	A given 1-d array's large gap (>period/2) is corrected.
span(inarray)	Just return the tuple of (min, max) for the given 1D array.
broadcast(array_in, shape_out)	Create the broadcasted array with the specific shape from the given array
broadcast_r(array_in, shape_out)	Broadcast, but in the revered direction
rtp(x, y, z, *[, degrees, primary_axis, ...])	Cartesian x, y, and z is converted to r, \(\theta\), and \(\phi\).
lonlat2xyz(longitude_deg, latitude_deg)	Simply convert longitude and latitude to array
trim_epsilon(arr[, epsilon])	Trim the epsilon.
nearest_index(a, x)	Return the index where the x is nearest to in array a.
nearest_index_right(xlist, x0)	Same as nearest_index()
nearest_index_left(xlist, x0)	Return the index where the x is nearest to in array a.

irfpy.util.nptools.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

Parameters:

• arr – Numpy array.

• 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

irfpy.util.nptools.max_where(arr, cond)

Conditional maximum, similar to min_where().

Parameters:

• arr – Numpy array.

• 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

irfpy.util.nptools.shrink(array, newshape, operation='mean')

The given array is shrinked in size with given operation.

Parameters:

• array – Numpy array. Multidimensional.

• 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.

• 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'))   
[[ 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]]

irfpy.util.nptools.squeeze_right(arr)

Similar to 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

irfpy.util.nptools.squeeze_left(arr)

Similar to 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

irfpy.util.nptools.padlin(vector, pad_width, iaxis, kwargs)

A helper function to use numpy:numpy.pad().

Use 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)   
[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.

irfpy.util.nptools.isascending(array1d, allow_equal=True)

Check if the array is ascending order or not

Parameters:

• array1d – Array to be checked.

• 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

irfpy.util.nptools.isdescending(array1d, allow_equal=True)

Check if the array is ascending order or not

Parameters:

• array1d – Array to be checked.

• 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

irfpy.util.nptools.guess_interval(array1d, log=False, method='median')

From the given 1-D array, a guessed interval array is returned.

The returned value is guess_upper() - guess_lower().

irfpy.util.nptools.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)   
[  2.5   4.    6.    9.   12.   15.   18. ]

The array is assumed to be sorted ascendingly, but no check is done.

Parameters:

• array1d – Ascending np.array.

• log – Set True to calculate this in the log space.

Returns:

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)    
[3.16227766e+00   3.16227766e+01   3.16227766e+02   3.16227766e+03  3.16227766e+04]

irfpy.util.nptools.guess_lower(array1d, log=False, method='median')

From the given 1D array, a guessed lower bound array is returned.

See the similar function guess_upper() for details.

>>> arr = [2, 3, 5, 7, 11, 13, 17]
>>> lower = guess_lower(arr)
>>> print(lower)   
[ 1.    2.5   4.    6.    9.   12.   15. ]

>>> arr = [1, 10, 100, 1000, 10000]
>>> arr_lower = guess_lower(arr, log=True)
>>> print(arr_lower)    
[3.16227766e-01   3.16227766e+00   3.16227766e+01   3.16227766e+02  3.16227766e+03]

irfpy.util.nptools.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.

Parameters:

• array1d – Array-like.

• log – If the estimation to be done in log space, use this option as True.

• method – Method to guess. Currently median is only supported.

Returns:

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]

irfpy.util.nptools.sort_periodic(array1d, period)

From an array in an periodic system, ascending array is obtained.

Parameters:

• array1d – An array.

• period – The period.

Returns:

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)   
[  0 100 200 300 400]

irfpy.util.nptools.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]

irfpy.util.nptools.span(inarray)

Just return the tuple of (min, max) for the given 1D array.

Parameters:

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)

irfpy.util.nptools.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.

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.

Parameters:

• array_in – Numpy array. Shape of (s0, s1, …, sS).

• 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))   
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,))   
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,)'

irfpy.util.nptools.broadcast_r(array_in, shape_out)

Broadcast, but in the revered direction

Note

You can get the same results by

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].

Parameters:

• array_in – np array with shape of (s0, s1, …)

• shape_out – The aimed shape. It should be (s0, s1, …, sX, sY, …)

Returns:

The broadcasted array from the “top” of the array.

See also:

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))   
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)'

irfpy.util.nptools.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\).

Parameters:

• x – X

• y – Y

• z – Z

• primary_axis – The primary axis, default is (0, 0, 1).

• secondary_axis – The secondary axis, default is (1, 0, 0).

• degrees – If True, returned angles are in degrees. Default False, i.e., radians.

• latitude – If True, returned polar angle is latitude 90 - \(\theta\) (ranging from -90 to 90 deg)

Returns:

(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).

Return type:

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)   
[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)   
[[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)   
[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

irfpy.util.nptools.lonlat2xyz(longitude_deg, latitude_deg)

Simply convert longitude and latitude to array

Parameters:

• longitude_deg – Longitude in degrees

• 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.]]

irfpy.util.nptools.lingrids(x0, x1, nx)

Return the grids.

Parameters:

• x0 – x0

• x1 – x1

• nx – Number of bins

Returns:

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

irfpy.util.nptools.loggrids(x0, x1, nx)

Return the grid.

Parameters:

• x0 – x0

• x1 – x1

• nx – Number of bins

Returns:

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

irfpy.util.nptools.trim_epsilon(arr, epsilon=1e-15)

Trim the epsilon.

Parameters:

• arr – A numpy array

• epsilon – \(\epsilon\)

Returns:

The input arr but with substituting zero (+0.0) for the components where the absolute value is less than epsilon

\[\begin{split}b_i = \left\{ \begin{array}{ll} a_i & (|a_i| \ge \epsilon) \\ 0 & (|a_i| < \epsilon) \end{array}\right.\end{split}\]

irfpy.util.nptools.nearest_index(a, x)

Return the index where the x is nearest to in array a.

Parameters:

• a – A list. It should be sorted in ascending order.

• 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

irfpy.util.nptools.nearest_index_right(xlist, x0)

Same as nearest_index()

irfpy.util.nptools.nearest_index_left(xlist, x0)

Return the index where the x is nearest to in array a.

Parameters:

• a – A list. It should be sorted in ascending order.

• 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