''' IMA science data utility.
MEX/IMA data access and further utilities.
See also :ref:`tutorial_vima_scidata2`.
'''
import datetime
import copy as _copy
import logging as _logging
_logger = _logging.getLogger(__name__)
from irfpy.mima import SENSOR_NAME
from irfpy.mima import scidata_tree
from irfpy.mima import scidata
from irfpy.imacommon import imascipac
from irfpy.imacommon import scidata as cscidata
from irfpy.util import tdict
import irfpy.mima.energy
import irfpy.mima.fov
from irfpy.imacommon import imamode as mode
import numpy as np
import irfpy.util.energyspectrum
from irfpy.util import exception as ex
from irfpy.util import ringcache as _ringcache
from irfpy.mima import ghostbuster as _ghostbuster
_tree = scidata_tree.MimaTree()
_array2d_cache_emu = _ringcache.RingCache(5)
_array2d_cache_nonemu = _ringcache.RingCache(5)
def _arrange_matlab_file(fn, emulate_full=False):
""" Load the file, cache, and return the emulated/non-emulated data
"""
# Cache check
if emulate_full:
cache = _array2d_cache_emu
else:
cache = _array2d_cache_nonemu
_logger.debug('Cache: Size={}'.format(cache.size()))
if cache.hasKey(fn):
_logger.debug("Cache hit")
return cache.get(fn)
else:
_logger.debug('Cache not hit.')
# Load and process data
scimat = scidata.ImaCountFile(fn).mat
try:
arr2 = cscidata.matfile2timeseries(scimat, emulate_full=emulate_full)
except ex.IrfpyException as e:
_logger.error('irfpy meets an error during the processing of {}.'.format(fn))
_logger.error(' IrfpyException: {}'.format(str(e)))
raise e
# Cache the data
cache.add(fn, arr2)
return arr2
[docs]def getarray2d(t0, t1, emulate_full=False):
''' Return the :class:`irfpy.imacommon.TimeSeriesCntMatrix2D` object.
Return the dataset of MEX/IMA counts as
:class:`irfpy.imacommon.TimeSeriesCntMatrix2D` object.
corresponding to the time range ``t0`` and ``t1``.
>>> import datetime
>>> t0 = datetime.datetime(2014, 10, 19, 18, 0, 0)
>>> t1 = datetime.datetime(2014, 10, 19, 18, 15, 0)
>>> dataset = getarray2d(t0, t1)
>>> print(dataset)
<TimeSeriesCntMatrix2D:len=76 from 2014-10-19 17:59:49.653240 to 2014-10-19 18:14:49.747040>
:keyword emulate_full: Data is emulated to represent the full matrix,
namely (M32, A16, E96) size matrix.
For mode such as 28, where the original data is (M32, A8, E96) is
expanded by assuming the uniform count rate over the bins to the full size.
:returns: Time series of MEX/IMA counts for the given time period
in the :class:`irfpy.imacommon.TimeSeriesCntMatrix2D` object.
The returned object is sometimes pass to the :func:`irfpy.imacommon.imascipac.convert2to3`
function to get the time series 3-D count matrix, or, :func:`irfpy.imacommon.imascipac.timeseries2d_to_energytime_simple`
to get the E-t diagram type of matrix.
'''
ts = imascipac.TimeSeriesCntMatrix2D()
fns = _tree.getfiles(t0, t1)
arr = imascipac.TimeSeriesCntMatrix2D()
for fn in fns:
try:
arr2 = _arrange_matlab_file(fn, emulate_full=emulate_full)
arr = arr + arr2
except ex.IrfpyException as e:
_logger.warning('irfpy meets an error. The data block for the file {} is disregarded.'.format(fn))
_logger.warning(' IrfpyException: {}'.format(str(e)))
try:
arr = arr.clipped(t0, t1)
except tdict.DataNotInDbError as e:
arr = imascipac.TimeSeriesCntMatrix2D()
### Set the name of the sensor
for t, matrix in arr:
matrix.name = SENSOR_NAME
return _copy.deepcopy(arr)
[docs]def getarray3d(t0, t1, emulate_full=False):
""" TO get the time series of 3D array for IMA data.
:param t0: Start of the time.
:param t1: Stop of the time.
:keyword emulate_full: Emulate M24, the full mode.
:return: :class:`irfpy.imacommon.TimeSeriesCntMatrix` object.
You can get time series of 3D data as:
>>> import datetime
>>> t0 = datetime.datetime(2014, 10, 16, 18, 2, 0)
>>> t1 = datetime.datetime(2014, 10, 16, 18, 20, 0)
>>> dataset = getarray3d(t0, t1)
>>> print(len(dataset))
7
>>> t3, d3 = dataset[3] # 3rd data.
>>> print(t3.strftime('%FT%T.%f'))
2014-10-16T18:09:26.640680
>>> print(d3.matrix.max()) # Max of the counts
64.0
Here ``d3`` is the :class:`irfpy.imacommon.TimeSeriesCntMatrix` object.
Note that the first and last 3d data is NOT necessarily bounded.
This means that a part of elevation scans for those dataset could be masked.
Therefore, it is a user's responsibility to check the completeness of 3D data.
Below is the example. You may realize the first 10 elevation is masked.
The time obtained is the time of the returned data: in this case the time at
elevation index 10.
>>> t0, d0 = dataset[0]
>>> print(t0)
2014-10-16 18:01:50.609400
>>> print(d0.matrix[30, 0, 30, :]) # M=30, A=0, E=30, and all P is printed.
[-- -- -- -- -- -- -- -- -- -- 0.0 0.0 0.0 0.0 0.0 0.0]
"""
dat2d = getarray2d(t0, t1, emulate_full=emulate_full)
dat3d = imascipac.convert2to3(dat2d, emulate_full=emulate_full)
return _copy.deepcopy(dat3d)
[docs]def c2j_proton_simple(count_matrix3d):
"""The simplest way of calculation proton differential flux.
Look at the source code what was done.
:param count_matrix3d: An object of class:`irfpy.imacommon.CntMatrix`.
The unwanted counts (e.g. high mass signals, alpha signals, background) should
be masked or removed prior to call this function.
This object should be produced using the emulate mode enabled, if the operation
mode is not full (M24).
See :func:`getarray2d` and :func:`convert2to3` for emulation mode.
:returns: The matrix of the differential flux.
"""
### The data will be summed over in the mass
matrix = count_matrix3d.matrix # (M32, A16, E96, P16) array
matrix = matrix.sum(axis=0) # (A16, E96, P16) array now
### Obtain the energy step
etbl = irfpy.mima.energy.get_default_table_v5_late() # post 2009 only supported
### Obtain the g-factor for proton of PACC=4
from irfpy.mima.gfactor0 import GL_CR4
g_dphi_func = GL_CR4.Proton.PAC1() # G0 dPhi in cm2 rad eV/eV. This is a function.
g_dphi = g_dphi_func(etbl) # (96, ) shape. [cm2 rad eV/eV] For higher energy than 800 eV has valid numbers.
fov = irfpy.mima.fov.SimpleFOV()
thetatbl = np.deg2rad(90 - fov.elevation_angle(np.arange(16)))
sinthetatbl = np.sin(thetatbl)
dtheta = np.deg2rad(5.625)
dt = 192 / 96 / 16. # Assume 100% duty time... =125ms
c_ape = matrix.transpose([0, 2, 1]) # (A16, P16, E96) array
c2_ape = c_ape / g_dphi / etbl # (A16, P16, E96) array
c2_aep = c2_ape.transpose([0, 2, 1]) # (A16, E96, P16) array
c3_aep = c2_aep / sinthetatbl / dtheta # Now C / DA DOmega Effic, I think.
c4_aep = c3_aep / dt
return c4_aep # Indeed, differential flux I think...
[docs]def energy_spectrum_3d_nomass(count_matrix, mass=1.67262178e-27, negative=True):
"""
:param count_matrix: Count matrix (:class:`irfpy.imacommon.imascipac.CntMatrix`)
In principle, it works with any mode, but the only confirmed use is for Mode 24 (including emulated)
:keyword mass: The mass assumed. Default to be protons.
:keyword negative: If *True*, negative energy is considered as valid (Default) and confirms the energy
step always 96.
If *False*, negative energy bin is disregarded, and thus the energy step may change (according to the
energy table)
:return: Energy spectrum class
:rtype: :class:`irfpy.util.energyspectrum.DifferentialFluxGrid`
Limitation: Only for post-2009 data (due to the energy table being hard-coded).
>>> import datetime
>>> t0 = datetime.datetime(2014, 10, 18, 11, 30)
>>> t1 = datetime.datetime(2014, 10, 18, 11, 50)
>>> data2d = getarray2d(t0, t1)
>>> data3d = imascipac.convert2to3(data2d)
>>> print(len(data3d))
7
>>> cnt3d = data3d[4][1] # You can take a single 3D distribution.
>>> print(cnt3d.t)
2014-10-18 11:40:45.991700
>>> spec = energy_spectrum_3d_nomass(cnt3d)
>>> print('{:.1f}'.format(spec.density())) # One may get density!
10.3
Note the coordinate system is in IMA-fixed frame but different from the IMAS frame in SPICE.
In the system used in the returned energy spectrum is: z is along the symmetric axis.
.. todo::
Check the coordinate system
"""
logger = _logging.getLogger(__name__ + ".energy_spectrum_3d_nomass")
etbl = irfpy.mima.energy.get_default_table_v5_late(negative=False)
fov = irfpy.mima.fov.SimpleFOV()
phitbl = fov.azimuth_angle(np.arange(16), direction='velocity')
logger.info(phitbl)
thetatbl = 90 - fov.elevation_angle(np.arange(16), direction='velocity')
logger.info(thetatbl)
### Count3D matrix should be converted to the differential flux array.
dflx = c2j_proton_simple(count_matrix) ### This is the simplest way. Probably, it is not bad idea to get dflux directly from argument.
# Now in (A16, E96, P16) array
dflx = dflx.transpose([1, 2, 0]) # (E, P, A) order.
### Count3D matrix should be re-ordered to confirm all the index to sort ascendingly.
dflx = dflx[::-1, :, :] # Revert the index of energy
etbl = etbl[::-1] # Revert the index of energy
### Now instance the EnergySpectrum object.
energyspectrum = irfpy.util.energyspectrum.DifferentialFluxGrid(dflx, etbl, np.deg2rad(thetatbl), np.deg2rad(phitbl), mass=mass)
return energyspectrum
[docs]def get3d(t, emulate_full=False):
""" Get the 3-D data at the time of t0.
:param t: Time of observation
:param emulate_full: Emulate M24, the full mode, if *True*
:return: A tuple, (time_start, VEX/IMA_data), where time_start is the datetime object and
VEX/IMA_Data is an object of :class:`irf.imacommon.imascipac.CntMatrix`.
>>> t, data = get3d(datetime.datetime(2014, 10, 18, 11, 50))
>>> print(t)
2014-10-18 11:47:10.147960
>>> print(data)
<<class 'irfpy.imacommon.imascipac.CntMatrix'>(MEX/IMA)@2014-10-18T11:47:10.147960:MOD=24 >>24<<:CNTmax=25>
"""
dt = datetime.timedelta(minutes=10)
t0 = t - dt
t1 = t + dt
timeser3d = getarray3d(t0, t1, emulate_full=emulate_full)
tstart, data3d = timeser3d.get(t)
return tstart, _copy.deepcopy(data3d)
[docs]def get2d(t, emulate_full=False):
""" Get the 2-D data at a specific time.
:param t: Time of observation
:param emulate_full: Emulate M24 if *True*.
:return: A tuple, (time_start, VEX/IMA 2D data). The IMA 2D data is an object of :class:`irfpy.imacommon.imascipac.CntMatrix2D`.
>>> t, data = get2d(datetime.datetime(2014, 10, 18, 11, 50))
>>> print(t)
2014-10-18 11:49:58.147960
>>> print(data)
<<class 'irfpy.imacommon.imascipac.CntMatrix2D'>(MEX/IMA)@2014-10-18T11:49:58.147960:MOD=24 >>24<<:POL=14/14:CNTmax=2>
"""
dt = datetime.timedelta(minutes=10)
t0 = t - dt
t1 = t + dt
timeser2d = getarray2d(t0, t1, emulate_full=emulate_full)
tstart, data2d = timeser2d.get(t)
return tstart, _copy.deepcopy(data2d)
[docs]def interpret_as_fast24s(imacnt3d, emulate_full=False):
""" A given IMA count 3D matrix is interpreted as eight 3-D reduced matrix.
:param imacnt3d: Ima count.
:type imacnt3d: :class:`irfpy.imacommon.imascipac.CntMatrix`
:keyword emulate_full: If given, the matrix will be sized to (M32, A16, E96, P16).
:return: A :class:`irfpy.imacommon.imascipac.TimeseiresCntMatrix` containing eight 3-D reduced matrix.
Each component is :class:`irfpy.imacommon.imascipac.CntMatrix` with (M32, A16, E32, P6) matrix
if ``emulate_full`` is set to ``False``. The mode information is set to ``M40``.
If ``emulate_full`` is set to ``True``, then the matrix will be sized to (M32, A16, E96, P16),
but masked for the unused energy and polar bins. In this case the mode informatin is
set to ``M41``.
:raises ``irfpy.util.exception.IrfpyError``: If the given imacnt3d is not the full matrix (``M24``).
In this case, users should load the data (likely by :func:`get3d` or :func:`getarray3d`)
using ``emulate_full`` set to ``True``.
"""
if imacnt3d.mode != mode.M24:
raise ex.IrfpyError('To interpret the data as Fast mode (EEPROM13), the data should be loaded as mode 24. Use ``emulate`` option to decode.')
sp = imacnt3d.matrix # (M32,A16,E96,P16) array
sp = sp.reshape(32, 16, 3, 32, 16) # (M32,A16,(E0-32, E1-32, E2-32), P16)
sp = sp.transpose([0, 1, 3, 4, 2]) # (M32, A16, E32, P16, Es3)
sp = sp.reshape(32, 16, 32, 48) # (M32, A16, E32, P+Es48)
sp = sp.reshape(32, 16, 32, 8, 6) # (M32, A16, E32, T8, P6)
sp = sp.transpose([0, 1, 2, 4, 3]) # (M32, A16, E32, P6, T8)
tlist = [imacnt3d.t + i * datetime.timedelta(seconds=24) for i in range(8)]
timeseries = imascipac.TimeSeriesCntMatrix()
for i in range(8):
matrix = sp[..., i]
if matrix.mask.all():
continue
if emulate_full:
matrix2 = np.ma.zeros([32, 16, 96, 16]) - 1e20
matrix2[:, :, :32, :6] = matrix
matrix = np.ma.masked_less(matrix2, 0)
packet = imascipac.CntMatrix(tlist[i], matrix, mode.M41, mode_data=imacnt3d.mode_data, name=imacnt3d.name)
else:
packet = imascipac.CntMatrix(tlist[i], matrix, mode.M40, mode_data=imacnt3d.mode_data, name=imacnt3d.name)
timeseries.append(tlist[i], packet)
return timeseries
### Processing that is used for both MEX-VEX.
from irfpy.imacommon.imascipac_util import interpolate_zero_azimuth, noise_reduction_constant, noise_reduction_isolated_counts
[docs]def interpolate_missing_massring(cntmatrix):
"""Interpolate the disabled mass ring number.
Ring number 0, 4, 10, 22 are disabled for MEX/ASPERA-4 from the certain date (TBC).
After this method, counts in ring number 0 are always 0.
counts in ring numbers (4, 10, 22) are teh interpolation with neighboring ring numbers.
Only valid for Mass=32 mode.
"""
mode = cntmatrix.mode
if mode.mass != 32:
import warnings
warnings.warn('interpolate_zero_azimuth can only be applied to 32 mass-bin data.')
return
cntmatrix.matrix = cntmatrix.matrix.copy() # New object is needed not to propagate back.
cntmatrix.matrix[0] = 0
cntmatrix.matrix[4] = (cntmatrix.matrix[3] + cntmatrix.matrix[5]) / 2.0
cntmatrix.matrix[10] = (cntmatrix.matrix[9] + cntmatrix.matrix[11]) / 2.0
cntmatrix.matrix[22] = (cntmatrix.matrix[21] + cntmatrix.matrix[23]) / 2.0
[docs]def ghost_reduction_rule_based_3d_rev1(cntmatrix):
""" Remove ghost according to the rule defined by HN interpreted by YF.
.. todo::
This function takes too much time. Somewhere we should optimize that.
Ghost is removed from the MEX data.
:param cntmatrix: Count matrix. The ghost removal is done in place.
:type cntmatrix: :class:`irfpy.imascipac.CntMatrix`
:return: A matrix of ghost, which were removed.
:rtype: ``numpy.ndarray`` with a shape same as the ``cntmatrix.data``
>>> t, data = get3d(datetime.datetime(2016, 10, 1, 3, 52, 40))
>>> matrix_original = data.matrix.copy()
>>> ghost = ghost_reduction_rule_based_3d_rev1(data)
>>> print(ghost.shape)
(32, 16, 96, 16)
>>> print(ghost[..., 2].sum()) # Ghost count, removed
2.0
>>> (matrix_original - data.matrix)[..., 2].sum() # The difference between the original and the removed shall be the same as the total ghost count.
2.0
"""
ghost = np.zeros_like(cntmatrix.matrix) # MAEP
for po in range(cntmatrix.matrix.shape[3]): # Loop over polar axis
g = _ghost_reduction_rule_based_2d_rev1(cntmatrix.matrix[..., po], cntmatrix.hk.promsection, cntmatrix.hk.pacclevel)
ghost[..., po] = g
cntmatrix.matrix = cntmatrix.matrix - ghost
cntmatrix.log_ghost = {'method': 'ghost_reduction_rule_based_3d_rev1',
'matrix': ghost}
return ghost
[docs]def ghost_reduction_rule_based_2d_rev1(cntmatrix2d):
""" Remove ghost according to the rule defined by HN interpreted by YF.
Ghost is removed from the MEX data.
:param cntmatrix2d: Count matrix. The ghost removal is done in place.
:type cntmatrix2d: :class:`irfpy.imascipac.CntMatrix2D`
:return: A matrix of ghost, which were removed.
:rtype: ``numpy.ndarray`` with a shape same as the ``cntmatrix2d.data``, typically (32, 16, 96)
>>> t, data = get2d(datetime.datetime(2016, 10, 1, 3, 52, 40))
>>> print(data.polar)
[2, 2]
>>> matrix_original = data.matrix.copy()
>>> ghost = ghost_reduction_rule_based_2d_rev1(data)
>>> print(ghost.shape)
(32, 16, 96)
>>> print(ghost.sum()) # Ghost count, removed
2
>>> (matrix_original - data.matrix).sum() # The difference between the original and the removed shall be the same as the total ghost count.
2
"""
ghost = _ghost_reduction_rule_based_2d_rev1(cntmatrix2d.matrix, cntmatrix2d.hk.promsection, cntmatrix2d.hk.pacclevel)
_logger.warning(f'{id(cntmatrix2d.matrix)}, {cntmatrix2d.matrix.sum()}')
val = cntmatrix2d.matrix - ghost
_logger.warning(f'{id(cntmatrix2d.matrix)}, {cntmatrix2d.matrix.sum()}')
cntmatrix2d.matrix = val
_logger.warning(f'{id(cntmatrix2d.matrix)}, {cntmatrix2d.matrix.sum()}')
_logger.warning(f'{id(ghost)}, {ghost.sum()}')
cntmatrix2d.log_ghost = {'method': 'ghost_reduction_rule_based_2d_rev1',
'matrix': ghost}
return ghost
def _ghost_reduction_rule_based_2d_rev1(matrix, promsection, pacclevel):
ghost = np.zeros_like(matrix) # MAE
for az in range(matrix.shape[1]): # Loop over azimuth.
emmatrix = _ghostbuster.EMmatrix(matrix[:, az, :].T,
promsection, pacclevel) # Object containing (96,32) array
foreground, ghost1 = _ghostbuster.buster_rule1(emmatrix)
foreground, ghost2 = _ghostbuster.buster_rule2(emmatrix)
foreground, ghost3 = _ghostbuster.buster_rule3(emmatrix)
foreground, ghost4 = _ghostbuster.buster_rule4(emmatrix)
foreground, ghost5 = _ghostbuster.buster_rule5(emmatrix)
foreground, ghost6 = _ghostbuster.buster_rule6(emmatrix)
foreground, ghost7 = _ghostbuster.buster_rule7(emmatrix)
ghost_matrix1 = ghost1.get_matrix() # (96,32) array, same for below
ghost_matrix2 = ghost2.get_matrix()
ghost_matrix3 = ghost3.get_matrix()
ghost_matrix4 = ghost4.get_matrix()
ghost_matrix5 = ghost5.get_matrix()
ghost_matrix6 = ghost6.get_matrix()
ghost_matrix7 = ghost7.get_matrix()
_logger.debug('Buster1/2/3/4/5/6/7: {} {} {} {} {} {}'.format(ghost_matrix1.sum(),
ghost_matrix2.sum(),
ghost_matrix3.sum(),
ghost_matrix4.sum(),
ghost_matrix5.sum(),
ghost_matrix6.sum(),
ghost_matrix7.sum()))
ghost_matrix = np.nanmax([ghost_matrix1,
ghost_matrix2,
ghost_matrix3,
ghost_matrix4,
ghost_matrix5,
ghost_matrix6,
ghost_matrix7,
], axis=0) # (96,32) array
ghost[:, az, :] = ghost_matrix.T
return ghost