""" 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
- one count addition
*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)
*See also*
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
>>> matx = np.ma.masked_less(matx, 0) # 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.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)
[docs]def one_count_addition(cntmatrix):
""" Add one count to all values >= 1 to counteract the one count reduction on board spacecraft.
codeauthor:: Moa Persson
: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.
>>> matx = np.ma.masked_less(matx, 0) # 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
>>> one_count_addition(cntmatrix)
>>> 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
>>> matx = np.ma.masked_less(matx, 0) # 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.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)
cntmatrix.matrix = _np.ma.masked_array(newmatrix, mask=cntmatrix.matrix_data.mask.copy())
[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
>>> matx = np.ma.masked_less(matx, 0) # 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)
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