''' cena_scicnt module is an implementaion for bridging Sadppac.
An implementation for bridging Sadppac.
Benefit is syntax sugar, caching, and using irfpy.util.irfpyrc module for path setting.
Entry of [cena] bmuobase in .irfpyrc file is loaded.
'''
import logging as _logging
import urllib.request
import urllib.parse
import urllib.error
import numpy
from Sadppac import *
from irfpy.util.irfpyrc import Rc
from irfpy.util.julday import Julday, JdSeries
from irfpy.util.utc import convert
_logger = _logging.getLogger(__name__)
[docs]class CenaLoader:
''' A CENA data loader.
CenaLoader is a class that treats CENA data.
CENA data is saved by the orbit sorted files in BMU format.
The reading routine is implemented in C++ library called libsadppac.
Reading routine in the libsadppac is based on the file-based, therefore
it is not convenient for daily analysis.
This class may bridge the libsadppac to usual analysis, which is based on observation time.
>>> cl=CenaLoader()
>>> print(cl.load(50))
None
>>> arr=cl.load(1945) # This will load cena data for orbit 1945.
>>> int(arr.size()) # Now v3.2, 1737. The original database returned 1769, but new data base on 100710 (nfs) returned 1768.
1737
'''
def __init__(self, rcfile=None):
''' Initializer
@param rcfile Configure file name.
'''
self.__rc = Rc()
if rcfile:
self.__rc.append(rcfile)
self.__baseuri=self.__rc.get('cena', 'bmuobase')
self.__cache = {}
_logger.debug("BMU OBASE=%s"%self.__baseuri)
if self.__baseuri is None:
_logger.error('BMU base not found. %s'%self.__baseuri)
raise RuntimeError('BMU base not found. %s'%self.__baseuri)
[docs] def clearCache(self):
self.__cache={}
[docs] def printCache(self):
print('Cache size=%d'%len(self.__cache))
for i in sorted(self.__cache):
print(i, self.__cache[i])
[docs] def load(self, orbitnr):
''' Load a data of the specified orbit number. If no data exists, None is returned.
'''
if orbitnr in self.__cache:
_logger.debug('Orbit number %d found'%orbitnr)
return self.__cache[orbitnr]
else:
uri='%s/orb_%04d/orb_%04d.cena'%(self.__baseuri, orbitnr, orbitnr)
_logger.debug('Reading data from %s'%uri)
try:
fnam, status = urllib.request.urlretrieve(uri)
_logger.debug('URL retrieved: %s / %s'%(fnam, status))
reader=CenaBmuReader()
reader.readFile(fnam)
carr=reader.getArray()
self.__cache[orbitnr]=carr
return carr
except IOError as e:
_logger.error('IOError happened! %s'%e)
self.__cache[orbitnr]=None
return None
[docs] def getObservationTime(self, orbitnr):
''' Obtain a tuple of observation time in the format of Julday.
'''
arr = self.load(orbitnr)
if arr is not None:
_logger.debug('Array size=%d'%arr.size())
else:
_logger.warning('No data loaded for %d' % orbitnr)
return []
t=[]
for ipac in range(arr.size()):
try:
ut=convert(arr.getPacket(ipac).getUtc(), Julday)
t.append(ut)
except Exception as e:
_logger.warning('Expetion caught during the convertion of observation time.\n***%s'%e)
continue
return t
[docs] def getTMmodes(self, orbitnr):
''' Obtain a JdSeries of TM mode for the specified orbit.
Sensor mode is defined in Sadppac.CenaPacket.TMMODE_XXX.
TMMODE_MASS_ACC, TMMODE_TOF_ACC, TMMODE_COUNT_ACC and TMMODE_NON_PROCESS is defined.
>>> c=CenaLoader()
>>> tm=c.getTMmodes(1945)
>>> isinstance(tm, JdSeries)
True
>>> print(tm.size()) # Now v3.2, 1737. The original value was 1769, but after the database update, 1768 is returned via nfs.
1737
>>> print(tm.first().getData()) # MASS_ACC=0x02
2
>>> tm=c.getTMmodes(1948)
>>> print(tm.first().getData()) # CNTR=0x03
3
>>> tm=c.getTMmodes(2040) # In fact, this is no data orbit for CENA.
>>> print(tm.size())
0
'''
arr = self.load(orbitnr)
if arr is None: # Being nodata, return 0-size JdSeries()
_logger.debug('No data found for orbit %d'%orbitnr)
return JdSeries()
_logger.debug('Array size=%d'%arr.size())
tmm=JdSeries()
for ipac in range(arr.size()):
try:
pac=arr.getPacket(ipac)
ut=convert(pac.getUtc(), Julday)
except Exception as e:
_logger.warning('Expetion caught during the convertion of observation time.\n***%s'%e)
continue
mod=pac.getTMmode()
_logger.debug('### %s %d'%(ut, mod))
tmm.add(ut, mod)
return tmm
[docs] def getFullMassMatrix(self, orbitnr):
''' Get mass matrix of the specified orbit.
Return the mass matrix for full model. Only data obtained in mass accumulation mode
with binning parameters nE=8, nD=7, nM=128, nP=1 is returned.
@retval JdSeries of mass matrix. Mass matrix is E=8xD=7xM=128 numpy array.
>>> c=CenaLoader()
>>> mm=c.getFullMassMatrix(1945) # mm is JdSeries with element of numpy.array
>>> isinstance(mm, JdSeries)
True
>>> print(mm.size()) ### v3.2, 1735. The original value was 1767, but after the database update, 1766 is returned.
1735
>>> mass_arr = mm.first().getData()
>>> print(mass_arr.shape)
(8, 7, 128)
>>> mm2=c.getFullMassMatrix(1948) # mm2 is JdSeries, but 1948 is in cntr mode.
>>> print(mm2.size())
0
'''
mm=JdSeries()
arr=self.load(orbitnr)
for ipac in range(arr.size()):
pac_org=arr.getPacket(ipac)
try:
jd=convert(pac_org.getUtc(), Julday)
except Exception as e:
_logger.warning('Conversion to Julday made exception. Skipped.\n*** %s'%e)
continue
if mm.hasElement(jd):
_logger.warning('Data at %s already exists. Skipped.'%jd)
continue
if pac_org.getTMmode() != CenaPacket.TMMODE_MASS_ACC:
_logger.info('Not in MASS mode. Skipping for %d at %s.'%(pac_org.get4sCycleCntr(), str(jd)))
continue
pac = pac_org.decompress() # Not clear if decompression failed. I think None will be returned.
if pac is None:
_logger.warning('Decompression failed. Skipping for %d at %s.'%(pac_org.get4sCycleCntr(), str(jd)))
continue
### Evaluate CenaPacket
nE = pac.getNrEne()
nC = pac.getNrCh()
nP = pac.getNrPhase()
nM = pac.getNrMass()
if nE!=8 or nC!=7 or nP!=1 or nM!=128:
_logger.warning('Binning parameter not supported for observation of %s. nE=%d(!=8), nC=%d(!=7), nP=%d(!=1), nM=%d(!=128)\n'%(jd,nE,nC,nP,nM))
continue
matx = pac.getCenaSciMassAcc()
if matx is None:
_logger.warning('Failed to get CenaSciMassAcc by unknown reason for %d at %s.'%(pac.get4sCycleCntr(), str(jd)))
continue
nC2 = matx.getNrChIntr() # Internal expression of channel number. In fact, this is 8.
data = matx.getCountAsArray()
data = numpy.array(data)
data=data.reshape(nE, nP, nC2, nM)
data = data[:, :, :nC, :].sum(1) # TODO: Inefficient way ??? Just to ignore the nP=1 axis.
mm.add(jd, data)
return mm
if __name__=='__main__':
loader=CenaLoader()
print(loader.load(1948))
print(loader.load(1948))
print(loader.load(1948))
print(loader.load(1948))
print(loader.load(1949))
print(loader.load(1950))
print(loader.load(1951))
print(loader.load(1952))
print(loader.load(1953))
print('--')
loader.printCache()
print('--')
loader.clearCache()
loader.printCache()
for jd in loader.getObservationTime(1950):
print(jd)