'''IMA extra data file access.
See also a tutorial :ref:`ima_extra`.
A demo script to plot E-t diagram in ``vima_demo/vima_dd_simpple.py``.
IMA extra is a calibrated database that Andrei Fedorov at IRAP prepared.
This module provides the access to the VEX/IMA extra dataset.
The database location should be specified by ``.irfpyrc`` configuration system.
::
[vima]
imaddncbase = /path/to/the/imaextra/database
IMA extra consists of four kinds of dataset.
- IMAPARAM (supported by irfpy)
- IMAEXTRA (supported by irfpy)
- IMAANG (not supported by irfpy yet)
- IMAEXTRALIGHT (not supported by irfpy yet)
The supported datasets are explained in the following.
*IMAPARAM* dataset
Parameter (moment) values that IRAP produced. It is *not* the official version.
.. warning::
There is no responsibility at IRF and irfpy develop team,
as well as ASPERA-4 PI group, for the accuracy of the daset.
Dataset should be used by the user's responsibility.
It is strongly recommended to contact to the ASPERA-4 PI team
for the use of the dataset.
There are eleven data in the IMA parameter database.
- ScanFlag: ?
- Quality_P: ?
- Density_P: Density H+
- Velocity_P: Velocity H+
- Temperature_P: Temperature H+
- EnergyThermal_P: Thermal Energy H+
- Quality_O: ?
- Density_O: Density O+
- Velocity_O: Velocity O+
- Temperature_O: Temperature O+
- EnergyThermal_O: Thermal Energy O+
The simplest way of getting data is using :func:`getImaParam` function.
>>> import datetime
>>> t0 = datetime.datetime(2008, 1, 27, 3, 0, 0)
>>> t1 = datetime.datetime(2008, 1, 27, 3, 10, 0)
>>> ima_parameters = getImaParam(t0, t1)
:attr:`obstime` is a list of the observation time.
>>> print(ima_parameters.obstime)
[datetime.datetime(2008, 1, 27, 3, 0, 18, 779000), datetime.datetime(2008, 1, 27, 3, 3, 30, 875000), datetime.datetime(2008, 1, 27, 3, 6, 42, 875000), datetime.datetime(2008, 1, 27, 3, 9, 54, 812000)]
:attr:`data` is a dictionary of the dataset.
>>> print(ima_parameters.data['Density_P']) # doctest: +NORMALIZE_WHITESPACE
[0.00756063 0.00653219 0.00567522 0.01819916]
*IMAEXTRA* dataset
IMA extra dataset is a separated counts to Proton and Oxygen.
Example of data access follows.
>>> from irfpy.vima import imaextra
>>> import datetime
>>> t0 = datetime.datetime(2011, 11, 5, 6)
>>> t1 = datetime.datetime(2011, 11, 5, 7)
>>> extdata = imaextra.getImaExtra(t0, t1)
You can know the observation start time as
>>> obstime = extdata.getobstime() # Returns an array
>>> print(obstime[0])
2011-11-05 06:01:48.605000
You can get the proton count spectrum as
>>> proton = extdata.getHpSpec()
>>> print(proton.shape) # The shape should be (A16, E96, P16, Time)
(16, 96, 16, 19)
>>> print(proton.sum()) # Total counts that owes to Proton is
11609479.0
Also, the oxygen (heavy) count spectrum as
>>> heavy = extdata.getHeavySpec()
>>> print(heavy.shape) # The shape should be (A16, E96, P16, Time)
(16, 96, 16, 19)
>>> print(heavy.max()) # Maximum counts in the bin that owes to Oxygen (heavy) is
3348.076
The count that could not be assigned either of proton or oxygen (rest matrix) is
>>> rest = extdata.getRestRm()
>>> print(rest.shape) # Rest matrix is ordered as (M32, A16, E96, P16, Time)
(32, 16, 96, 16, 19)
>>> print(rest.sum()) # Total counts of rest matrix
799203.06
*Direct access to the IMA extra data*
You can access the data directly via the attribute :attr:`data`;
while this is not recommended.
Rather, use the access functions as :meth:`getHpSpec` or :meth:`getHeavySpec`.
Nevertheless, the ``data`` is a dictironary, with key and entries as follows
- 'AngTableN': Angular table index. (N,)
- 'ETableN': Energy table index. 0 (ver-1) or 1 (ver-3). (N,)
- 'FracO2': O2+ fraction (TBC). (N, P16, E96)
- 'HeavySpec': Heavy ion spectrum. (N, P16, E96, A16)
- 'HpSpec': Proton ion spectrum. (N, P16, E96, A16)
- 'Pacc': Post acceleration (N,)
- 'RestRm': The rest spectrum. (N, P16, E96, A16, M32)
- 'Time': Time. (N,)
Here *N* is the number of data.
The order in the raw matrix is the "inverted way" as MAEPT order
(MAEPT is the standard order of array for irfpy).
.. warning::
Again, be careful using the "RAW" data, since the order of the axis of the array is
different from the irfpy-standard.
Therefore, again and again, it is recommended to use the accessing methods, like
:meth:`ImaExtra.getHpSpec` or similar, since these have the standard MAEPT order.
*The following text is for developers*
>>> vimaddncpath = irfpy.util.irfpyrc.Rc().get('vima', 'imaddncbase')
>>> vimaextrapath = os.path.join(vimaddncpath, 'IMAEXTRA')
>>> vimaparampath = os.path.join(vimaddncpath, 'IMAPARAM')
:func:`getImaExtraFile` and :func:`getImaParamFile` will open the data file.
>>> fe = getImaExtraFile('imaextra20123640214.nc.gz')
>>> fp = getImaParamFile('imaparam20100630434.nc.gz')
You may also use the time directly.
>>> fe = getImaExtraFile(t=datetime.datetime(2007, 9, 10, 3, 0))
>>> fp = getImaParamFile(t=datetime.datetime(2008, 1, 27, 3, 0))
Each IMAEXTRA file corresponds to :class:`ImaExtraFile`.
This can be instanced manually, but use
:class:`ImaExtraFileFactory`, which provides accessor and cache to
:class:`ImaExtraFile` object.
>>> fact = ImaExtraFileFactory.createFactory()
>>> iex = fact.getImaExtraFile(os.path.join(vimaextrapath, 'imaextra20123640214.nc.gz'))
>>> print(len(iex.var['Time']))
39
>>> fact2 = ImaParamFileFactory.createFactory()
>>> ipr = fact2.getImaParamFile(os.path.join(vimaparampath, 'imaparam20100630434.nc.gz'))
>>> print(len(ipr.var['Time']))
38
The dataset include many files (thus, :class:`ImaExtraFile` objects).
The :class:`ImaExtraDatabase` provides a simple database of the dataset.
This class has functionality converting the time to filename
(:meth:`get_filename`).
>>> db = ImaExtraDatabase.createDatabase()
>>> fn = db.get_filename(datetime.datetime(2007, 9, 10, 3, 0))
>>> print(os.path.basename(fn)) # The filename
imaextra20072520101.nc.gz
>>> db2 = ImaParamDatabase.createDatabase()
>>> fn = db2.get_filename(datetime.datetime(2008, 1, 27, 3, 0))
>>> print(os.path.basename(fn))
imaparam20080260247.nc.gz
'''
import os
import datetime
import logging
_logger = logging.getLogger(__name__)
import numpy as np
import irfpy.imacommon.imaextra2
from irfpy.imacommon.imaextra2 import ImaDdNcFileFactory, ImaDdNcDatabase
import irfpy.util.irfpyrc
import irfpy.util.timeseriesdb
import irfpy.util.ringcache
from irfpy.util import exception as _ex
[docs]def isdb():
rc = irfpy.util.irfpyrc.Rc()
database = rc.get('vima', 'imaddncbase')
if database is None:
return False
if os.path.isdir(database):
return True
else:
return False
[docs]class ImaParam(irfpy.imacommon.imaextra2.ImaParamCommon):
''' The object stores data of IMA parameter
>>> t0 = datetime.datetime(2010, 7, 10, 4, 30)
>>> file0 = getImaParamFile(t=t0)
>>> ip0 = file0.getImaParam()
>>> print(ip0.ndat)
32
>>> t1 = datetime.datetime(2010, 7, 13, 8, 40)
>>> file1 = getImaParamFile(t=t1)
>>> ip1 = file1.getImaParam()
>>> print(ip1.ndat)
37
>>> ip01 = ip0 + ip1
>>> print(ip01.data['Velocity_P'].shape)
(69, 3)
>>> ip10 = ip1 + ip0
>>> print(ip01.data['Velocity_O'].shape)
(69, 3)
>>> print((ip01.data['Temperature_P'] == ip10.data['Temperature_P']).all())
True
>>> print(ip01.obstime == ip10.obstime)
True
>>> # specifying the same dataset for adding, duplication check works.
>>> ip11 = ip1 + ip1 # If you add two identical dataset
>>> print(ip11.data['ScanFlag'].shape) # the result is the same as before.
(37,)
>>> print(ip1.obstime == ip11.obstime)
True
'''
data_keys = ['ScanFlag',
'Quality_P', 'Density_P', 'Velocity_P', 'Temperature_P', 'EnergyThermal_P',
'Quality_O', 'Density_O', 'Velocity_O', 'Temperature_O', 'EnergyThermal_O',
]
def __init__(self):
irfpy.imacommon.imaextra2.ImaParamCommon.__init__(self, self.data_keys)
def __add__(self, other):
new = ImaParam()
new.append(self)
new.append(other)
return new
[docs]class ImaParamFile(irfpy.imacommon.imaextra2.ImaParamFileCommon):
''' An IMA param file for VEX.
IMA param file class. Usually, you want to instance via
:class:`ImaDdNcFileFactory` then caching will be effective.
Member of ``dim`` contains the dimension information and ``var`` for variables.
'''
def __init__(self, filename, gunzip=True):
''' Open the file and read the data.
'''
irfpy.imacommon.imaextra2.ImaParamFileCommon.__init__(self, filename, gunzip=gunzip)
[docs] def getImaParam(self):
''' Obtain the :class:`ImaParam` object
'''
ip = ImaParam()
ip.obstime = self.getobstime()
ip.ndat = len(ip.obstime)
for k in list(ip.data.keys()):
ip.data[k] = np.array(self.var[k], copy=True)
return ip
[docs]class ImaParamFileFactory(ImaDdNcFileFactory):
__singleton_instance = None
def __init__(self, *args, **kwds):
ImaDdNcFileFactory.__init__(self, ImaParamFile, *args, **kwds)
[docs] @classmethod
def createFactory(cls):
''' Return the "master" factory.
'''
if cls.__singleton_instance is None:
cls.__singleton_instance = ImaParamFileFactory()
return cls.__singleton_instance
[docs] def getImaParamFile(self, filename, gunzip=True):
return self.getDdNcFile(filename, gunzip=gunzip)
[docs]class ImaParamDatabase(ImaDdNcDatabase):
__singleton_instance = None
def __init__(self, database=None):
if database is None:
rc = irfpy.util.irfpyrc.Rc()
database = rc.get('vima', 'imaddncbase')
database = os.path.join(database, 'IMAPARAM')
ImaDdNcDatabase.__init__(self, database, 'imaparam', ImaParamFileFactory())
[docs] @classmethod
def createDatabase(cls, refresh=False):
''' Create the database.
'''
if cls.__singleton_instance is None or refresh:
cls.__singleton_instance = ImaParamDatabase()
return cls.__singleton_instance
[docs]def getImaParamFile(filename=None, t=None, gunzip=True):
''' Get the IMA parameter in a form of :class:`ImaParam` object.
It is a user function, but slightly low-level. High level functions are to be made
.. todo
Create a high-level user function that returns IMA parameter in a simple way
:param t0: Time to start
:param t1: Time to end
:returns: IMA parameter object
:rtype: :class:`ImaParam`
>>> imaParam = getImaParam(datetime.datetime(2011, 11, 5, 6), datetime.datetime(2011, 11, 5, 8))
'''
if filename is not None:
try:
return ImaParamFileFactory.createFactory().getImaParamFile(filename, gunzip=gunzip)
except IOError as e:
logger = logging.getLogger('getImaParamFile')
logger.warning('No file found %s.' % filename)
rc = irfpy.util.irfpyrc.Rc()
dbpath = rc.get('vima', 'imaddncbase')
filename = os.path.join(dbpath, 'IMAPARAM', filename)
logger.warning('Try %s now.' % filename)
return ImaParamFileFactory.createFactory().getImaParamFile(filename, gunzip=gunzip)
elif t is not None:
db = ImaParamDatabase.createDatabase()
filename = db.get_filename(t)
return ImaParamFileFactory.createFactory().getImaParamFile(filename, gunzip=gunzip)
else:
raise ValueError('Either ``filename`` or ``t`` should be specified')
[docs]def getImaParam(t0, t1):
''' Return the :class:`ImaParam` object.
:param t0: Time to start
:param t1: Time to end
:returns: IMA parameter object
:rtype: :class:`ImaParam`
'''
db = ImaParamDatabase.createDatabase() # Data base for IMA parameter files
if len(db) == 0:
raise _ex.IrfpyException('Data file not found. Check the setup (e.g., .irfpyrc, or the data file)')
fns = db.get_filenames(t0, t1) # Return the file name covering t0 to t1.
_logger.debug('Number of file: {}'.format(len(fns)))
ip = ImaParam()
imaparam_reader = lambda filename: ImaParamFileFactory.createFactory().getImaParamFile(fn).getImaParam()
for fn in fns:
_logger.debug('Start reading file ({})'.format(fn))
from irfpy.util import filepairv
import irfpy.vima
iprm = filepairv.filepair(fn, fn + '.filepair.gz', converter=imaparam_reader, version=irfpy.vima.__version__)
# iprm = ImaParamFileFactory.createFactory().getImaParamFile(fn).getImaParam()
_logger.debug('Finish reading file ({}). Data length={}'.format(fn, iprm.ndat))
ip.append(iprm)
ip.trim(t0, t1)
return ip
[docs]class iterparam:
""" Iterate the IMA parameter file.
Iterated is a list of ``(t, dat)`` were ``t`` is the time of the observation start,
and ``dat`` is a dictionary that contains the data.
The key of ``dat`` is :attr:`ImaParam.data_keys`. The ``dat`` is *not* a :class:`ImaParam` object, but
rather, :attr:`ImaParam.data` data.
:param t0: Time to start. Iterated data's start time is always later than or equal to t0.
:param t1: Time to stop. Iterated data's stop time is always earlier than or equal to t1.
:returns: A list, ``(t, dat)``. ``t`` is time, and ``dat`` is a dictionary containing the IMAPARAM data.
Its key is :attr:`ImaParam.data_keys`.
:rtype: ``list`` of (``datetime.datetime``, ``dict``)
If you want to iterate data from 2011-11-05T05:15 for 10 minutes,
you can do as follows:
>>> t0 = datetime.datetime(2011, 11, 5, 5, 15)
>>> t1 = datetime.datetime(2011, 11, 5, 5, 25)
>>> for t, mom in iterparam(t0, t1):
... print(t, mom['Density_P'])
2011-11-05 05:17:00.200000 9.732157
2011-11-05 05:20:12.262000 11.949523
2011-11-05 05:23:24.262000 10.680418
"""
data_keys = ImaParam.data_keys
""" Key of the data
"""
_logger = logging.getLogger(__name__ + '.iterparam')
def __init__(self, t0, t1):
self.t0 = t0
self.t1 = t1
db = ImaParamDatabase.createDatabase() # Create a database for IMA extra files.
if len(db) == 0:
raise _ex.IrfpyException('Data file not found. Check the setup (e.g., .irfpyrc, or the data file)')
self.filename_list = db.get_filenames(self.t0, self.t1) # Return the file names to read
self.current_file_contents = None
self.current_file_position = -1
self.current_obstime = self.t0
def __iter__(self):
return self
def __next__(self):
while len(self.filename_list) > 0:
self._logger.debug('Remained files = {}'.format(len(self.filename_list)))
self._logger.debug('The target file = {}'.format(self.filename_list[0]))
# If current file contents are not loaded
if self.current_file_contents is None:
self._logger.debug('No file contents have been loaded. Start to load')
self.current_file_contents = getImaParamFile(filename=self.filename_list[0]).getImaParam()
self.current_file_position = -1
## HERE YOU SHOULD PUT if FAILED!
self._logger.debug('File contents read finished.')
# You proceed 1 data, and return.
self._logger.debug('Progress 1 data')
self.current_file_position += 1
try:
self._logger.debug('Get the time the loaded contents')
t = self.current_file_contents.obstime[self.current_file_position]
self._logger.debug('Finish time list: Result {}'.format(t))
if t <= self.current_obstime:
self._logger.debug('The obtainted time is not appropriate. Skip this')
continue
if t > self.t1:
self._logger.debug('The obtainted time is past the required range. Stop iteration')
raise StopIteration()
self._logger.debug('Start loading the IMA param file.')
dat = self.current_file_contents
self._logger.debug('Finished.')
self._logger.debug('Preparation of the returned data')
dat2 = {}
for k in self.data_keys:
dat2[k] = dat.data[k][self.current_file_position]
self._logger.debug('Finish preparation')
self.current_obstime = t
return t, dat2
except IndexError:
self._logger.debug('Whole the data is returned. Go next file')
self.filename_list.pop(0) # Remove the current file
self.current_file_contents = None
self.current_file_position = -1
raise StopIteration()
class _iterextra:
data_keys = ImaExtra.data_keys
""" Key of the data
"""
_logger = logging.getLogger(__name__ + '.iterextra')
def __init__(self, t0, t1):
self.t0 = t0
self.t1 = t1
db = ImaExtraDatabase.createDatabase() # Create a database for IMA extra files.
self.filename_list = db.get_filenames(self.t0, self.t1) # Return the file names to read
self.current_file_contents = None
self.current_file_position = -1
self.current_obstime = self.t0
def __iter__(self):
return self
def __next__(self):
while len(self.filename_list) > 0:
self._logger.debug('Remained files = {}'.format(len(self.filename_list)))
self._logger.debug('The target file = {}'.format(self.filename_list[0]))
# If current file contents are not loaded
if self.current_file_contents is None:
self._logger.debug('No file contents have been loaded. Start to load')
self.current_file_contents = getImaExtraFile(filename=self.filename_list[0]).getImaExtra()
self.current_file_position = -1
## HERE YOU SHOULD PUT if FAILED!
self._logger.debug('File contents read finished.')
# You proceed 1 data, and return.
self._logger.debug('Progress 1 data')
self.current_file_position += 1
try:
self._logger.debug('Get the time the loaded contents')
t = self.current_file_contents.obstime[self.current_file_position]
self._logger.debug('Finish time list: Result {}'.format(t))
if t <= self.current_obstime:
self._logger.debug('The obtainted time is not appropriate. Skip this')
continue
if t > self.t1:
self._logger.debug('The obtainted time is past the required range. Stop iteration')
raise StopIteration()
self._logger.debug('Start loading the IMA extra file.')
dat = self.current_file_contents
self._logger.debug('Finished.')
self._logger.debug('Preparation of the returned data')
dat2 = {}
for k in self.data_keys:
dat2[k] = dat.data[k][self.current_file_position]
self._logger.debug('Finish preparation')
self.current_obstime = t
return t, dat2
except IndexError:
self._logger.debug('Whole the data is returned. Go next file')
self.filename_list.pop(0) # Remove the current file
self.current_file_contents = None
self.current_file_position = -1
raise StopIteration()