Source code for irfpy.imacommon.imascipac_util

""" Utilities for science (raw count) matrix of IMA.

*Introduction*

Input matrix, :class:irfpy.imacommon.imascipac.CntMatrix object, will be processed
using pre-defined processing functions.  They include

- noise reduction
- interpolation

*How it works*

The input matrix (attribute matrix will be modified.
On the other hand, there is another matrix
without processing (raw data) remained as the name matrix_data.
The matrix has two attribute of the count rate data.

- matrix: Processed matrix
- matrix_data: Matrix in the original data file (non-processed data)

Functions dedicated to MEX/VEX are also implemented in :mod:irfpy.mima.scidata_util and
:mod:irfpy.vima.scidata_util.
"""
import numpy as _np

[docs]def noise_reduction_constant(cntmatrix, value=1):
""" Remove the counts below the constant value as noise.

:param cntmatrix: Count matrix. The attribute matrix is changed.
:type cntmatrix: :class:CntMatrix
:keyword value: The background value.
:returns: Nothing.

The counts in the matrix > value is unchanged.
The counts in the matrix <= value will be set to zero.

>>> import numpy as np
>>> matx = np.zeros([32, 16, 96, 16])
>>> matx[0, 0, 0, 0] = 1
>>> matx[0, 1, 2, 3] = 2
>>> matx[1, 2, 3, 4] = 4
>>> matx[:, :, 95, :] = -1
>>> import datetime
>>> from irfpy.imacommon.imascipac import CntMatrix
>>> from irfpy.imacommon import imamode
>>> cntmatrix = CntMatrix(datetime.datetime(1975, 1, 1, 0, 0, 0), matx, imamode.M24)
>>> print(cntmatrix.matrix.sum())
7.0
>>> print(cntmatrix.matrix[0, 0, 95, 0])
--
>>> noise_reduction_constant(cntmatrix)
>>> print(cntmatrix.matrix.sum())
6.0
>>> print(cntmatrix.matrix[0, 0, 95, 0])
--
>>> print(cntmatrix.matrix_data.sum())
7.0
"""
cntmatrix.matrix = _np.ma.where(cntmatrix.matrix <= value, 0, cntmatrix.matrix)

""" Add one count to all values >= 1 to counteract the one count reduction on board spacecraft.

:param cntmatrix: Count matrix. The attribute matrix is changed.
:type cntmatrix: :class:CntMatrix
:returns: Nothing.

The counts in the matrix == 0 is unchanged.
The counts in the matrix >= 1 will have one count added.

>>> import numpy as np
>>> matx = np.zeros([32, 16, 96, 16])
>>> matx[0, 0, 0, 0] = 0
>>> matx[0, 1, 2, 3] = 2
>>> matx[1, 2, 3, 4] = 4
>>> matx[:, :, 95, :] = -1    # To emulate the masked array.
>>> import datetime
>>> from irfpy.imacommon.imascipac import CntMatrix
>>> from irfpy.imacommon import imamode
>>> cntmatrix = CntMatrix(datetime.datetime(1975, 1, 1, 0, 0, 0), matx, imamode.M24)
>>> print(cntmatrix.matrix.max())
4.0
>>> print(cntmatrix.matrix[0, 0, 95, 0])
--
>>> print(cntmatrix.matrix[0, 0, 0, 0])
0.0
>>> print(cntmatrix.matrix.max())
5.0
>>> print(cntmatrix.matrix[0, 0, 95, 0])
--
>>> print(cntmatrix.matrix[0, 0, 0, 0])
0.0
>>> print(cntmatrix.matrix_data.max())     # matrix_data should not change.
4.0
"""
cntmatrix.matrix = _np.ma.where(cntmatrix.matrix >= 1, cntmatrix.matrix + 1, cntmatrix.matrix)

[docs]def noise_reduction_isolated_counts(cntmatrix, value=1):
""" The counts isolated is considered as the noise.

:param cntmatrix: Count matrix. The attribute matrix is changed.
:type cntmatrix: :class:CntMatrix
:keyword value: The background value.
:returns: Nothing.

>>> import numpy as np
>>> matx = np.zeros([32, 16, 96, 16])
>>> matx[0, 0, 0, 0] = 1    # Isolated -> Removed
>>> matx[3, 5, 8, 10:12] = 1   # Not isolated.
>>> matx[0, 1, 2, 3] = 2   # Isolated, but more than value=1
>>> matx[1, 2, 3, 4] = 4   # Not isolated
>>> matx[2, 2, 3, 4] = 1   # Not isolated
>>> matx[:, :, 95, :] = -1  # Mask
>>> import datetime
>>> from irfpy.imacommon.imascipac import CntMatrix
>>> from irfpy.imacommon import imamode
>>> cntmatrix = CntMatrix(datetime.datetime(1975, 1, 1, 0, 0, 0), matx, imamode.M24)
>>> print(cntmatrix.matrix.sum())
10.0
>>> print(cntmatrix.matrix[0, 0, 95, 0])
--
>>> noise_reduction_isolated_counts(cntmatrix)
>>> print(cntmatrix.matrix.sum())
9.0
>>> print(cntmatrix.matrix[0, 0, 95, 0])
--
>>> print(cntmatrix.matrix[0, 0, 0, 0])   # Originally 1, but removed.
0.0
>>> print(cntmatrix.matrix[2, 2, 3, 4])   # Originally 1, but not isolated so that not removed.
1.0
>>> print(cntmatrix.matrix_data.sum())
10.0
"""
extended_matrix = _np.pad(cntmatrix.matrix.filled(0), 1, mode='constant', constant_values=0)    # Add each edges to 0.  Now the shape is (34, 18, 98, 18).

is_below = (cntmatrix.matrix <= value)
is_ax0_minu = (extended_matrix[:-2, 1:-1, 1:-1, 1:-1] == 0)
is_ax0_plus = (extended_matrix[2:, 1:-1, 1:-1, 1:-1] == 0)
is_ax1_minu = (extended_matrix[1:-1, :-2, 1:-1, 1:-1] == 0)
is_ax1_plus = (extended_matrix[1:-1, 2:, 1:-1, 1:-1] == 0)
is_ax2_minu = (extended_matrix[1:-1, 1:-1, :-2, 1:-1] == 0)
is_ax2_plus = (extended_matrix[1:-1, 1:-1, 2:, 1:-1] == 0)
is_ax3_minu = (extended_matrix[1:-1, 1:-1, 1:-1, :-2] == 0)
is_ax3_plus = (extended_matrix[1:-1, 1:-1, 1:-1, 2:] == 0)

### If all the above conditions are satisfied, the count is considered as noise.
is_noise = _np.array([is_below,
is_ax0_minu, is_ax0_plus,
is_ax1_minu, is_ax1_plus,
is_ax2_minu, is_ax2_plus,
is_ax3_minu, is_ax3_plus,
])     # (9, M32, A16, E32, P16) shape.
is_noise = is_noise.all(axis=0)

newmatrix = _np.where(is_noise, 0, cntmatrix.matrix)

[docs]def interpolate_zero_azimuth(cntmatrix):
r""" Interpolate the 0-azimuth channel from 1- and 15-aziuth channels.

Sometimes, azimuth channel 0's counts contain crosstalk over the other channels.
Thus, it might need correction.

Here the processing of the following is provided.

.. math::

C[iM, 0, iE, iP] = \frac{C[iM, 1, iE, iP] + C[iM, 15, iE, iP]}{2}

Only valid for 16 azimuth channel mode (e.g. mode.M24, mode.M25).

Here an emulating data isp repared.

>>> import numpy as np
>>> matx = np.random.randint(0, 50, [32, 16, 96, 16])   # Random counts.
>>> matx[30, 0, 10, 5] = 100
>>> matx[30, 1, 10, 5] = 30
>>> matx[30, 15, 10, 5] = 20
>>> matx[:, :, 95, :] = -1  # Mask
>>> import datetime
>>> from irfpy.imacommon.imascipac import CntMatrix
>>> from irfpy.imacommon import imamode
>>> cntmatrix = CntMatrix(datetime.datetime(1975, 1, 1, 0, 0, 0), matx, imamode.M24)

Now the processing.

>>> interpolate_zero_azimuth(cntmatrix)
>>> print(cntmatrix.matrix[30, 0, 10, 5])   # to be replaced to 25 (=(30 + 20) / 2)
25
>>> print((cntmatrix.matrix[:, 1:, :, :] == cntmatrix.matrix_data[:, 1:, :, :]).all())   # All matrix except for az=0 is the same.
True
"""

mode = cntmatrix.mode
if mode.azim != 16:
import warnings
warnings.warn('interpolate_zero_azimuth can only be applied to 16 azimuthal data.')
return

cntmatrix.matrix = cntmatrix.matrix.copy()   # Needed to avoid the backward propagation

cntmatrix.matrix[:, 0, :, :] = (cntmatrix.matrix[:, 15, :, :] + cntmatrix.matrix[:, 1, :, :]) / 2.0
`