''' IMA science data, providing the count rate data access
It is a common part of IMA count rate science data access
in the Matlab format.
For nominal use, refer to :mod:`irfpy.mima.scidata`
and :mod:`irfpy.vima.scidata`.
'''
import logging as _logging
logger = _logging.getLogger(__name__)
import numpy as np
import irfpy.util.exception as ex
from irfpy.util import utc
from irfpy.imacommon import imamode
from irfpy.imacommon import imascipac
[docs]def reshape_mat3d(mat, npol=16):
''' Reshape the matlab contents into MAEPT order
The function reshapes the matrix to MAEPT order
(mass-azimuth-energy-polar-time).
Polar (elevaion) number should be specified, if
you know it is not 16 elevation mode.
:returns: (spectra, masslist, azimlist, polarlist, timelist).
spectra has (N, 32, 16, 96, 16) element, where N is the
number of 3-D data.
masslist, azimlist, polarlist has (N, 2, 16, 96, 16) shape.
timelist has (N, 16, 96, 16) shape.
'''
s, m, a, e, t = reshape_mat(mat) # Now, MAET order.
### s (32, 16, 96, T)
### m (2, 16, 96, T)
### a (2, 16, 96, T)
### e (2, 16, 96, T)
### T (16, 96, T)
splitidx = list(range(npol, s.shape[3], npol)) # 0-15, 16-31, ...
ss = np.array(np.split(s, splitidx, axis=3)) # N, 32, 16, 96, T/N
mm = np.array(np.split(m, splitidx, axis=3)) # N, 2, 16, 96, T/N
aa = np.array(np.split(a, splitidx, axis=3)) # N, 2, 16, 96, T/N
ee = np.array(np.split(e, splitidx, axis=3)) # N, 2, 16, 96, T/N
tt = np.array(np.split(t, splitidx, axis=2)) # N, 16, 96, T/N
ss = np.rollaxis(ss, 0, 5)
mm = np.rollaxis(mm, 0, 5)
aa = np.rollaxis(aa, 0, 5)
ee = np.rollaxis(ee, 0, 5)
tt = np.rollaxis(tt, 0, 4)
return ss, mm, aa, ee, tt
[docs]def reshape_mat(mat):
''' Reshape the matlab contents into MAET order.
The function reshape the matrix to MAET order.
:param mat: Matlab contents.
:returns: (spectra, masslist, azimlist, elevlist, timelist)
Each output has the shape of (32, 16, 96, T) for spectra;
(32, 16, T) for timelist;
and (2, 32, 16, T) for others. T is the number of time.
:param mat: Matlab format data file
:returns: Tuple of array. (spectra, mass, azimuth, polar, time).
'''
tl = mat['iontime'][0, :]
n = len(tl)
tl2 = tl.reshape(n // (1 * 32 * 16), 1, 32, 16) # Originally generated using parameter of 8-multiple. Leave it as 1-multiple, hohefully it will work "forever".
tl2 = np.rollaxis(tl2, 2)
tl2 = np.rollaxis(tl2, 3, 1)
tl2 = tl2.reshape((32, 16, n // (32 * 16)))
sp = mat['ionspectra']
sp2 = sp.reshape(96, n // (1 * 32 * 16), 1, 32, 16)
sp2 = np.rollaxis(sp2, 3)
sp2 = np.rollaxis(sp2, 4, 1)
sp2 = sp2.reshape((32, 16, 96, n // (32 * 16)))
el = mat['elev']
el2 = el.reshape(2, n // (1 * 32 * 16), 1, 32, 16)
el2 = np.rollaxis(el2, 3, 1)
el2 = np.rollaxis(el2, 4, 2)
el2 = el2.reshape((2, 32, 16, n // (32 * 16)))
ma = mat['masschannel']
ma2 = ma.reshape(2, n // (1 * 32 * 16), 1, 32, 16)
ma2 = np.rollaxis(ma2, 3, 1)
ma2 = np.rollaxis(ma2, 4, 2)
ma2 = ma2.reshape((2, 32, 16, n // (32 * 16)))
az = mat['azim']
az2 = az.reshape(2, n // (1 * 32 * 16), 1, 32, 16)
az2 = np.rollaxis(az2, 3, 1)
az2 = np.rollaxis(az2, 4, 2)
az2 = az2.reshape((2, 32, 16, n // (32 * 16)))
return (sp2, ma2, az2, el2, tl2)
[docs]def matfile2timeseries(matrix, emulate_full=False):
''' Getting the matrix data in Matlab, time series of matrix is returned.
:param matrix: Matlab data (dictionary) read via :func:`scipy.io.loadmat`,
or using higher-level way, e.g. :class:`irfpy.vima.scidata.ImaCountFile`
: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 2D matrix
:rtype: :class:`irfpy.imacommon.imascipac.TimeSeriesCntMatrix2D`
This is a function for developer.
>>> from irfpy.vima.scidata import get_ima_count_file_full
>>> import datetime
>>> matlab_data = get_ima_count_file_full(datetime.datetime(2008, 3, 3, 14, 55))
>>> timeseries2d = matfile2timeseries(matlab_data.mat)
'''
### The returned object
tseries = imascipac.TimeSeriesCntMatrix2D()
tl = matrix['iontime'][0, :] # (n,)
n = len(tl)
sp = matrix['ionspectra'] # (96, n)
el = matrix['elev'] # (2, n)
ma = matrix['masschannel'] # (2, n)
az = matrix['azim'] # (2, n)
logger.debug('SP shape = %s', str(sp.shape))
### HK. Embedded one, so the number of data should be (n,). (HK in HK file not considered.)
try:
hkmode = matrix['mode'][0] # (1, n) -> (n,)
logger.debug('HK shape = %s', str(hkmode.shape))
promsection = matrix['promsection'][0]
logger.debug('PROM shape = %s', str(promsection.shape))
pacclevel = matrix['pacclevel'][0]
logger.debug('PACC shape = %s', str(pacclevel.shape))
fifofill = matrix['fifo_fill'][0]
logger.debug('FIFO shape = %s', str(fifofill.shape))
except KeyError:
logger.error("The database looks old. See https://irfpy.irf.se/projects/aspera/cookbook/recipe_database_old_error.html")
raise
## Check if the data is already sorted.
tl2 = np.sort(tl)
if not (tl2 == tl).all():
raise ex.IrfpyException('The data given is not correctly ordered')
## tu is the unique time, ti is its start index.
tu, ti = np.unique(tl, return_index=True)
## End index.
ti2 = np.roll(ti, -1)
ti2[-1] = n
## Index list to be a list (iterable).
ti = ti.tolist()
ti2 = ti2.tolist()
for i0, i1 in zip(ti, ti2):
t0 = np.unique(tl[i0:i1])[0]
sp0 = np.array(sp[:, i0:i1], dtype=int)
el0 = np.array(el[:, i0:i1], dtype=int)
ma0 = np.array(ma[:, i0:i1], dtype=int)
az0 = np.array(az[:, i0:i1], dtype=int)
hk = imascipac.HK()
hk.mode = hkmode[i0]
hk.pacclevel = pacclevel[i0]
hk.promsection = promsection[i0]
hk.fifofill = fifofill[i0]
logger.debug('{} {} {} {} {}'.format(
t0, sp0.shape, el0.shape, ma0.shape, az0.shape))
### To guess the mode...
el0u0 = np.unique(np.sort(el0[0]))
el0u1 = np.unique(np.sort(el0[1]))
ma0u0 = np.unique(np.sort(ma0[0]))
ma0u1 = np.unique(np.sort(ma0[1]))
az0u0 = np.unique(np.sort(az0[0]))
az0u1 = np.unique(np.sort(az0[1]))
nmass = len(ma0u0)
if nmass != len(ma0u1):
logger.warning('\n'.join(('Data structure violation',
'Mass channel data corrupted?')))
logger.warning('Data for t={} ({}) is skipped'.format(
t0, utc.mat2dt(t0)))
continue
naz = len(az0u0)
if naz != len(az0u1):
logger.warning('\n'.join(('Data structure violation',
'Azimuth channel data corrupted?')))
logger.warning('Data for t={} ({}) is skipped'.format(
t0, utc.mat2dt(t0)))
continue
if not (len(el0u0) == len(el0u1) == 1):
logger.warning('\n'.join(('Data structure violation',
'Elevation channel data corrupted?')))
logger.warning('Data for t={} ({}) is skipped'.format(
t0, utc.mat2dt(t0)))
continue
numel = el0u1[0] - el0u0[0] + 1
if (numel, naz, nmass) == (1, 16, 32):
if sp0.shape != (96, 512,):
logger.warning('\n'.join(('Data structure looks M24',
'but inconsistent with the spectrum shape ({}).'.format(sp0.shape),
'Numel={numel}, naz={naz}, nmass={nmass}'.format(numel=numel, naz=naz, nmass=nmass))))
logger.warning('Data for t={} ({}) is skipped'.format(
t0, utc.mat2dt(t0)))
continue
mode = imamode.M24
elif (numel, naz, nmass) == (2, 16, 32):
if sp0.shape != (96, 512,):
logger.warning('\n'.join(('Data structure looks M25',
'but inconsistent with the spectrum shape ({}).'.format(sp0.shape),
'Numel={numel}, naz={naz}, nmass={nmass}'.format(numel=numel, naz=naz, nmass=nmass))))
logger.warning('Data for t={} ({}) is skipped'.format(
t0, utc.mat2dt(t0)))
continue
mode = imamode.M25
elif (numel, naz, nmass) == (8, 16, 32):
if sp0.shape != (96, 512,):
logger.warning('\n'.join(('Data structure looks M27',
'but inconsistent with the spectrum shape ({}).'.format(sp0.shape),
'Numel={numel}, naz={naz}, nmass={nmass}'.format(numel=numel, naz=naz, nmass=nmass))))
logger.warning('Data for t={} ({}) is skipped'.format(
t0, utc.mat2dt(t0)))
continue
mode = imamode.M27
elif (numel, naz, nmass) == (8, 8, 32):
if sp0.shape != (96, 256):
logger.warning('\n'.join(('Data structure looks M28',
'but inconsistent with the spectrum shape ({}).'.format(sp0.shape),
'Numel={numel}, naz={naz}, nmass={nmass}'.format(numel=numel, naz=naz, nmass=nmass))))
logger.warning('Data for t={} ({}) is skipped'.format(
t0, utc.mat2dt(t0)))
mode = imamode.M28
## elif...
else:
logger.warning('\n'.join(('Data structure unknown or not supported.',
'Numel={numel}, naz={naz}, nmass={nmass}'.format(numel=numel, naz=naz, nmass=nmass))))
logger.warning('Data for t={} ({}) is skipped'.format(
t0, utc.mat2dt(t0)))
continue
logger.debug('Guessed={}'.format(mode))
### Reshaping for each mode.
matx = None
if mode in (imamode.M24, imamode.M25, imamode.M27):
### (96, 512) to be (M32, A16, E96) order.`
if sp0.shape != (96, 512):
logger.warning('Wrong shape ({}) for mode {}.'.format(sp0.shape, mode))
logger.warning('Data for t={} ({}) is skipped'.format(
t0, utc.mat2dt(t0)))
continue
sp1 = sp0.reshape([96, 32, 16]) # EMA order.
sp2 = np.transpose(sp1, (1, 2, 0)) # MAE order.
ut0 = utc.mat2dt(t0)
scipac = imascipac.CntMatrix2D(ut0,
np.ma.masked_array(sp2),
[el0u0[0], el0u1[0]],
mode, hk=hk)
tseries.append(ut0, scipac)
### FOR M28, the processing depends on the emulation on or off.
elif mode in (imamode.M28,): # (96, 256) shape.
### (96, 256) shape to (M32, A8, E96) order.
if sp0.shape != (96, 256):
logger.warning('Wrong shape ({}) for mode {}.'.format(sp0.shape, mode))
logger.warning('Data for t={} ({}) is skipped'.format(
t0, utc.mat2dt(t0)))
continue
sp1 = sp0.reshape([96, 32, 8]) # EMA order.
sp2 = np.transpose(sp1, (1, 2, 0)) # MAE order.
ut0 = utc.mat2dt(t0)
if not emulate_full:
scipac = imascipac.CntMatrix2D(ut0,
np.ma.masked_array(sp2),
[el0u0[0], el0u1[0]],
mode, hk=hk)
tseries.append(ut0, scipac) # Just add the matrix into the array
else:
### Convert the sp2 into full spec 2D = Mode 27.
sp3 = np.ma.zeros([32, 16, 96]) # MAE order
sp3[:, ::2, :] = sp2 / 2.
sp3[:, 1::2, :] = sp2 / 2.
scipac = imascipac.CntMatrix2D(ut0,
np.ma.masked_array(sp3),
[el0u0[0], el0u1[0]],
imamode.M27, mode_data=imamode.M28, hk=hk)
tseries.append(ut0, scipac) # Just add the matrix into the array
return tseries
import unittest
import doctest
[docs]def doctests():
return unittest.TestSuite((
doctest.DocTestSuite(),
))
if __name__ == '__main__':
unittest.main(defaultTest='doctests')