''' MEX/IMA science data, providing the count rate data access
.. note::
This is quick and dirty implementation; It was copied from :mod:`irfpy.vima.scidata`.
.. note::
The latest appraoch of implementing scidata accessor is in :mod:`scidata_util`.
Consider using that module.
:class:`ImaCountFile` is the main class.
This corresponds to a IMA matlab formatted data file.
:func:`get_ima_count_file` will return the data file contents
corresponding to the time specified.
A simple way of getting data is as follows.
You may use the :func:`get_ima_count_file` function.
>>> t = datetime.datetime(2006, 12, 5, 7, 50)
>>> imafile = get_ima_count_file(t)
>>> if imafile is not None:
... mat = imafile.mat # Matlab data dictionary.
... print(os.path.basename(imafile.file)) # doctest: +SKIP
20061205070000.mat
If no data is found, you will get None as returned.
>>> t = datetime.datetime(2006, 12, 5, 7)
>>> imafile = get_ima_count_file(t)
>>> print(imafile)
None
How to handle the instance ``mat``?
See an independent document (prepared for VEX but can be applied to MEX)
:ref:`tutorial_vima_scidata`.
The data path should be specified by
[mima] imamatbase in ``.irfpyrc``.
'''
import os
import datetime
import logging
import warnings
import numpy as np
import scipy.io
from irfpy.util.irfpyrc import Rc
from irfpy.util import exception as _ex
from irfpy.util.utc import num2date
from irfpy.asperacommon.excepts import DbPathNotFoundError
from irfpy.asperacommon.excepts import DatafileNotFoundError
import irfpy.imacommon.scidata as imascidata
_logger = logging.getLogger(__name__)
[docs]class ImaCountFile:
''' MEX/IMA data file.
>>> imacnt = ImaCountFile(ImaCountFile.get_sample_filename())
'''
rc = Rc()
[docs] @classmethod
def get_sample_filename(cls):
p = cls.get_base_path()
fn = os.path.join(p, '2010', '11', '09', 'DDS_IMA_20101109T140000.mat')
logging.debug(fn)
return fn
[docs] @classmethod
def get_uribase_path(cls):
import warnings
warnings.warn('get_uribase_path is no longer supported. Use get_base_path', DeprecationWarning)
return cls.get_base_path()
[docs] @classmethod
def get_uri(cls, yr, mo, dy, hr):
warnings.warn('get_uri is no longer supported. Use get_filename', DeprecationWarning)
return cls.get_filename(yr, mo, dy, hr)
[docs] @classmethod
def get_base_path(cls):
# cls.rc = Rc()
path = cls.rc.get('mima', 'imamatbase')
if path is None:
raise DbPathNotFoundError(
'No entry of [mima]/imamatbase in irfpyrc.')
return path
[docs] @classmethod
def get_filename(cls, yr, mo, dy, hr):
ttry = datetime.datetime(yr, mo, dy, hr, 0, 0)
subpath = os.path.join(cls.get_base_path, ttry.strftime('%Y%m'),
ttry.strftime('%Y%m%d%H0000.mat'))
return subpath
logger = logging.getLogger('ImaCountFile')
def __init__(self, filename):
''' Load the matlab count rate file
:param filename: The filename of the file.
URI is no more supported (no more maintained.)
'''
if filename.startswith('http://') or filename.startswith('https://') or filename.startswith('file:///'):
warnings.warn('No URI support (%s). Use local file.' %
filename)
# try:
# self.file, self.uriheader = urllib.request.urlretrieve(filename)
# except IOError as e:
# emsg = 'Failed to open the matlab file at %s' % filename
# self.logger.error(emsg)
# raise DatafileNotFoundError(emsg)
self.uri = self.file = filename
# self.mat contains all the matlab data.
try:
self.mat = scipy.io.loadmat(self.file)
''' Matlab formatte data'''
except (ValueError) as e:
""" This exception confirms the existence of data file, but unreadable due to matfile corruption.
"""
_logger.warning('The data file {} may be corrupted.'.format(self.file))
_logger.warning('emesg: {}'.format(e))
_logger.warning('No data has been read from the above file.')
raise _ex.IrfpyException('Data file {} corrupted'.format(self.file))
except (IOError, FileNotFoundError, TypeError) as e:
"""This exception is for checking the data file exists or not.
"""
emsg = 'Failed to open the matlab file: %s' % self.file
raise DatafileNotFoundError(emsg)
[docs] def getobstime(self):
''' Return the observation time list.
:returns: Observation time list (component is in
``datetime.datetime`` instance) is returned.
>>> imacnt = ImaCountFile(ImaCountFile.get_sample_filename())
>>> print(len(imacnt.getobstime()))
304
'''
obstime_mat = list(set(self.mat.get('iontime')[0]))
obstime_mat.sort()
obstime_mpl = [num2date(i - 366) for i in obstime_mat]
obstime_mpl = [i.replace(tzinfo=None) for i in obstime_mpl]
return obstime_mpl
[docs] def gettimerange(self, start_offset=0, stop_offset=0):
''' Return the time range of the observation.
:returns: Range of data file time as a list of
``datetime.datetime`` instances.
[*t0* + ``start_offset``, *t1* + ``stop_offset``] is returned, where
*t0* is the start time (the time of the first packet) and *t1* is
the end time (the time of the last packet).
This means that the end time *t1* does not agree with the end of
the observation time.
The end of the observation time is the end time plus data
acquisition time (12 sec or 24 sec, typically).
Specifying ``stop_offset`` will compensate this issue.
:keyword start_offset: Offset in seconds to be added to the start time.
:keyword stop_offset: Offset in seconds to be added to the stop time.
>>> imacnt = ImaCountFile(ImaCountFile.get_sample_filename())
>>> print(imacnt.gettimerange())
[datetime.datetime(2010, 11, 9, 14, 0, 40, 323800), datetime.datetime(2010, 11, 9, 15, 1, 16, 823860)]
>>> print(imacnt.gettimerange(stop_offset=12))
[datetime.datetime(2010, 11, 9, 14, 0, 40, 323800), datetime.datetime(2010, 11, 9, 15, 1, 28, 823860)]
'''
tlist = self.getobstime()
off0 = datetime.timedelta(seconds=start_offset)
off1 = datetime.timedelta(seconds=stop_offset)
return [min(tlist) + off0, max(tlist) + off1]
def __len__(self):
return len(self.getobstime())
[docs]def isdb():
try:
dbpath = ImaCountFile.get_base_path()
except DbPathNotFoundError:
return False
return os.path.isdir(dbpath)
[docs]def get_ima_count_file_reduced(time):
''' Return the :class:`ImaCountFile` instance covering the specified time, from reduced-HK embedded data.
:param time: Time of the data
:type time: ``datetime.datetime`` instance.
Note: This is not yet fully debugged.
'''
ttry = time
path = os.path.join(
ImaCountFile.get_base_path(),
ttry.strftime('%Y%m'),
ttry.strftime('%Y%m%d%H0000.mat'))
try:
imacnt = ImaCountFile(path)
t0, t1 = imacnt.gettimerange()
if time < t0:
ttry = time - datetime.timedelta(seconds=3600)
except DatafileNotFoundError as e:
ttry = time - datetime.timedelta(seconds=3600) # Rewind 1 hour
if ttry != time:
ttry = time - datetime.timedelta(seconds=3600)
# Data 1 hour ago should be fetched.
path = os.path.join(
ImaCountFile.get_base_path(),
ttry.strftime('%Y%m'),
ttry.strftime('%Y%m%d%H0000.mat'))
try:
imacnt = ImaCountFile(path)
except DatafileNotFoundError as e:
return None
return imacnt
[docs]def get_ima_count_file_full(time):
''' Return the :class:`ImaCountFile` instance covering the specified time, from full-HK embedded data.
:param time: Time of the data
:type time: ``datetime.datetime`` instance.
Note: This is not yet fully debugged.
'''
ttry = time
path = os.path.join(
ImaCountFile.get_base_path(),
ttry.strftime('%Y'),
ttry.strftime('%m'),
ttry.strftime('%d'),
ttry.strftime('DDS_IMA_%Y%m%dT%H0000.mat'))
try:
imacnt = ImaCountFile(path)
t0, t1 = imacnt.gettimerange()
if time < t0:
ttry = time - datetime.timedelta(seconds=3600)
except DatafileNotFoundError as e:
ttry = time - datetime.timedelta(seconds=3600) # Rewind 1 hour
if ttry != time:
ttry = time - datetime.timedelta(seconds=3600)
# Data 1 hour ago should be fetched.
path = os.path.join(
ImaCountFile.get_base_path(),
ttry.strftime('%Y'),
ttry.strftime('%m'),
ttry.strftime('%d'),
ttry.strftime('DDS_IMA_%Y%m%dT%H0000.mat'))
try:
imacnt = ImaCountFile(path)
except DatafileNotFoundError as e:
return None
return imacnt
[docs]def get_ima_count_file(time):
''' Return the :class:`ImaCountFile` instance covering the specified time.
The dataset is selected either from
- Full-HK embedded data (aka Leif-matfile, or aurora-matfile). :func:`get_ima_count_file_full`
- Reduced-HK embedded data (aka Herbert-matfile, or spaghetti-matfile). :func:`get_ima_count_file_reduced`
Because of slightly different name convention, this function judges which dataset user is
referring automatically. Full-HK is prioritized, and definitively recommended to be used.
:param time: Time of the data
:type time: ``datetime.datetime`` instance.
:returns: Ima count file.
:rtype: :class:`ImaCountFile` instance.
'''
icf = get_ima_count_file_full(time)
if icf is None:
icf = get_ima_count_file_reduced(time)
return icf
[docs]def reshape_mat3d(mat, npol=16):
''' Reshape the matlab contents into MAEPT order
See also :func:`irfpy.vima.scidata.reshape_mat3d`.
'''
return imascidata.reshape_mat3d(mat, npol=npol)
[docs]def reshape_mat(mat):
''' Reshape the matlab contents into MAET order.
MEX has, in nominal, operatied with full mode (MODE24).
It has (32, 16, 96, 16) array.
: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.
>>> f = get_ima_count_file(datetime.datetime(2006, 12, 5, 8, 15))
>>> if f is not None:
... s, m, a, e, t = reshape_mat_24(f.mat)
... else:
... s = m = a = e = t = None
>>> print(s.shape) # doctest: +SKIP
(32, 16, 96, 288)
>>> print(m.shape) # doctest: +SKIP
(2, 32, 16, 288)
>>> print(t.shape) # doctest: +SKIP
(32, 16, 288)
## The above results are in 2012-10-11.
'''
return imascidata.reshape_mat(mat)
[docs]def reshape_mat_fast24s(mat):
''' A specific reshape tool for 24s mode. MAEPT order.
:param mat: Input of the matlab formatted data
:returns: A tuple of (spectra, masslist, azimlist, elevlist, timelist).
Here ``spectra`` is the observed counts with a shape of
(M=32, A=16, E=32, P=6, T=<arb.>)
``masslist``, ``azimlist``, ``elevlist`` is not yet supported.
``timelist`` is the observation time list with a shape of (T=<arb.>).
.. todo..
Merge to generic routine.
>>> t0 = datetime.datetime(2013, 12, 29, 7, 50)
>>> f = get_ima_count_file(t0)
>>> s, m, a, e, t = reshape_mat_fast24s(f.mat)
>>> print(s.shape)
(32, 16, 32, 6, 144)
This is the MAEPT order (mass-azimth, energy, polar, time.
>>> print(t.shape)
(144,)
'''
sp = mat['ionspectra'] # E96xElse:
ndat = sp.shape[1]
if ndat % (32 * 16 * 2) != 0:
warnings.warn('Data block may be corrupted. Continue processing, but use the data carefully.')
ndat2 = ndat // 32 // 16
sp = sp.flatten().reshape(96, ndat2, 32, 16)
### Sp is now ordered, E, T, M, A
### Here, Time is really a time, 12s, including polar information.
### Energy to be 3x32.
sp = sp.reshape(3, 32, ndat2, 32, 16)
### Sp is now T0, E, T, M, A
sp = sp.transpose([1, 2, 0, 3, 4])
### Sp is now E, T, T0, M, A
sp = sp.reshape(32, ndat2 * 3, 32, 16)
### Sp is now E, T, M, A
sp = sp.reshape(32, ndat2 * 3 // 6, 6, 32, 16)
### Sp is now E, T, P, M, A
sp = np.transpose(sp, [3, 4, 0, 2, 1])
ta = mat['iontime']
ta = ta.flatten().reshape(ndat2, 32, 16)
ta = ta[:, 0, 0]
### ta is each 12 s data
ta2 = ta[::2] ### ta is each 24s data
return sp, None, None, None, ta2
[docs]def reshape_mat_24(mat):
''' (Deprecated!) Reshape the matlab contents into MAET order.
Use :func:`reshape_mat`.
'''
import warnings
warnings.warn('Use reshape_mat function.', category=DeprecationWarning)
return reshape_mat(mat)