'''This module is for ACE data to compare with SWIM.
ACE data, both non-shifted and shifted to the Moon position, is supported
by this module. The non-shifted data was obtained from the web, and the class Ace is for it.
The shifted data can be generated by ace2moon.py program in the script folder,
and the program dumps it to ace2moon.dat file.
The class AceMoon is for the shifted data.
The simplest way of retrieving the solar wind condition at the Moon is as follows.
>>> acemoon = AceMoon()
>>> dat = acemoon.getAverage(datetime.datetime(2009,1,25,13,0,0), datetime.datetime(2009,1,25,14,0,0))
'''
import os
import sys
import math
import urllib.request
import urllib.parse
import urllib.error
import logging
import re
import numpy
import datetime
import dateutil.parser
import calendar
import gzip
from irfpy.util.irfpyrc import *
from irfpy.util.utc import *
from irfpy.util.julday import *
from irfpy.cy1orb.MoonPos import MoonGse
__msg_isleapyear_obsolete=False
def _isleapyear(year):
''' (OBSOLETE) Return True if the specified year is a leap year.
This function is obsolete, because the equivalent function is prepared by
the python API. Use the equivalent calendar.isleap() function.
'''
global __msg_isleapyear_obsolete
if not __msg_isleapyear_obsolete:
logging.warn('irfpy.swim.ace._isleapyear() obsolete: Use calendar.isleap()')
__msg_isleapyear_obsolete=True
if year % 4 != 0:
return False
if year % 100 !=0 :
return True
if year % 400 == 0:
return True
return False
def _doi2md(year, doi):
''' DOI converts to (month, date)
@TODO the function should go to util.utc?
'''
if calendar.isleap(year):
month=[0,1,32,61,92,122,153,183,214,245,275,306,336]
if doi <=0 or doi >=367:
raise ValueError('Specified DOY %d does not exist for year %d.'%(doi, year))
else:
month=[0,1,32,60,91,121,152,182,213,244,274,305,335]
if doi <=0 or doi >=366:
raise ValueError('Specified DOY %d does not exist for year %d.'%(doi, year))
for mon in range(12,0,-1):
if month[mon] <= doi:
m=mon
d=doi-month[mon]+1
return m,d
def _md2doi(year, mon, day):
''' Convert month-day to doi.
@TODO the function should go to util.utc?
'''
if calendar.isleap(year):
month=[0,1,32,61,92,122,153,183,214,245,275,306,336]
else:
month=[0,1,32,60,91,121,152,182,213,244,274,305,335]
return month[mon]+day-1
[docs]class AceElem:
''' Class for ACE data at a time.
'''
def __init__(self, line):
data=line.split()
if len(data)!=44:
raise ValueError('Line %s does not contain data'%line)
self.__data=data
[docs] def getJulday(self):
yr=int(self.__data[0])
mo,dy=_doi2md(yr, int(self.__data[1]))
hr=int(self.__data[2])
mn=int(self.__data[3])
sec=float(self.__data[4])
return Julday(yr,mo,dy,hr,mn,sec)
def __str__(self):
return '(ACEdata@%s)'%(self.getJulday())
[docs] @classmethod
def getContents(self):
'''
'''
cont='ear day hr min sec fp_year fp_doy ACEepoch Np Tp Alpha_ratio Vp V_rtn_r V_rtn_t V_rtn_n V_gse_x V_gse_y V_gse_z V_gsm_x V_gsm_y V_gsm_z B_rtn_r B_rtn_t B_rtn_n B_gse_x B_gse_y B_gse_z B_gsm_x B_gsm_y B_gsm_z Bmag Lambda Delta dBrms pos_gse_x pos_gse_y pos_gse_z pos_gsm_x pos_gsm_y pos_gsm_z pos_hs_x pos_hs_y pos_hs_z MAG_pts'.split()
return cont
[docs] def getNp(self):
return float(self.__data[8])
[docs] def getTp(self):
return float(self.__data[9])
[docs] def getAlpha_ratio(self):
return float(self.__data[10])
[docs] def getVp(self):
return float(self.__data[11])
[docs] def getV_gse_x(self):
return float(self.__data[15])
[docs] def getV_gse_y(self):
return float(self.__data[16])
[docs] def getV_gse_z(self):
return float(self.__data[17])
[docs] def getB_gse_x(self):
return float(self.__data[24])
[docs] def getB_gse_y(self):
return float(self.__data[25])
[docs] def getB_gse_z(self):
return float(self.__data[26])
[docs] def getBmag(self):
return float(self.__data[30])
[docs] def getpos_gse_x(self):
return float(self.__data[34])
[docs] def getpos_gse_y(self):
return float(self.__data[35])
[docs] def getpos_gse_z(self):
return float(self.__data[36])
[docs]class AceElemAtMoon:
''' Ace data at Moon at an instance. Using AceElem class, the time of the Moon is calculated.
'''
__moonpos=None
def __init__(self, aceElem):
''' Initialize the data. ValueError is thrown if the data is not available.
'''
if AceElemAtMoon.__moonpos is None:
logging.info('%%%%% Moon data is not there.')
AceElemAtMoon.__moonpos=MoonGse()
AceElemAtMoon.__moonpos.setFromDefaultUri()
else:
logging.info('%%%%% Moon data is there.')
moonpos=AceElemAtMoon.__moonpos
t=aceElem.getJulday() # The observation time
logging.debug( t)
xa=aceElem.getpos_gse_x() # Position of ACE in km
logging.debug( xa)
xm=moonpos.getGse(t)[0]
logging.debug( xm)
### X-distance
xdist = xm - xa
logging.debug( 'Distance is %f km'%xdist)
vx = aceElem.getV_gse_x()
if vx == -9999.9:
raise ValueError('Velocity data by ACE is not available at %s'%t)
logging.debug( 'Velocity is %f km/s'%vx)
tx = xdist/vx
logging.debug( 'Time difference is %f s'%tx)
self.jd_at_moon = t.dayAfter(tx*Julday.SECOND)
logging.debug( 'Time at ACE is %s / Time at Moon is %s'%(t,self.jd_at_moon))
self.aceElem=aceElem
[docs] def getJulday(self):
return self.jd_at_moon
[docs] def getNp(self):
return self.aceElem.getNp()
[docs] def getTp(self):
return self.aceElem.getTp()
[docs] def getAlpha_ratio(self):
return self.aceElem.getAlpha_ratio()
[docs] def getVp(self):
return self.aceElem.getVp()
[docs] def getV_gse_x(self):
return self.aceElem.getV_gse_x()
[docs] def getV_gse_y(self):
return self.aceElem.getV_gse_y()
[docs] def getV_gse_z(self):
return self.aceElem.getV_gse_z()
[docs] def getB_gse_x(self):
return self.aceElem.getB_gse_x()
[docs] def getB_gse_y(self):
return self.aceElem.getB_gse_y()
[docs] def getB_gse_z(self):
return self.aceElem.getB_gse_z()
[docs] def getBmag(self):
return self.aceElem.getBmg()
[docs] def getpos_gse_x(self):
return self.aceElem.getpos_gse_x()
[docs] def getpos_gse_y(self):
return self.aceElem.getpos_gse_y()
[docs] def getpos_gse_z(self):
return self.aceElem.getpos_gse_z()
[docs]class AceMoonData:
def __init__(self, data_tuple):
self.__data=data_tuple
[docs] def getJulday(self):
t=dateutil.parser.parse(self.__data[1])
return convert(t, Julday)
[docs] def getNp(self):
return float(self.__data[4])
[docs] def getTp(self):
return float(self.__data[5])
[docs] def getAlpha_ratio(self):
return float(self.__data[6])
[docs] def getVp(self):
return float(self.__data[7])
[docs] def getV_gse_x(self):
return float(self.__data[8])
[docs] def getV_gse_y(self):
return float(self.__data[9])
[docs] def getV_gse_z(self):
return float(self.__data[10])
[docs] def getB_gse_x(self):
return float(self.__data[11])
[docs] def getB_gse_y(self):
return float(self.__data[12])
[docs] def getB_gse_z(self):
return float(self.__data[13])
[docs] def getBmag(self):
return float(self.__data[14])
[docs] def getpos_gse_x(self):
return float(self.__data[15])
[docs] def getpos_gse_y(self):
return float(self.__data[16])
[docs] def getpos_gse_z(self):
return float(self.__data[17])
def __str__(self):
return 'AceMoonData@%s [N=%s,V=%s,T=%s,R=%s,B=%s]'%(self.__data[1],self.__data[4], self.__data[7], self.__data[5], self.__data[6], self.__data[14])
[docs]class AceMoon:
''' Class for ACE Moon data.
The shifted ACE data is generated by ace2moon.py in script folder from the non-shifted ACE data.
The ACE data can be infered to the Moon position considering the ACE and Moon position,
and ace2moon.py will dump the result to shifted ace2moon.dat file.
The file location can be specified by .irfpyrc, where [swim] section acemoon property.
'''
def __init__(self, rcfile=None):
self.__rc = Rc()
if rcfile:
self.__rc.append(rcfile)
self.__acemoonuri = self.__rc.get('swim', 'acemoon')
logging.debug('ACE MOON data will be retrieved from %s'%self.__acemoonuri)
self.__acemoon_local,status = urllib.request.urlretrieve(self.__acemoonuri)
logging.debug('Retrieved to %s, with status of "%s"'%(self.__acemoon_local, status))
self.__cache=JdSeries()
self.__cached_day=[]
[docs] def clearCache(self):
self.__cache=JdSeries()
self.__cached_day=[]
[docs] def read(self, year, month, day):
''' Read a data from cache for the specified day.
'''
f=open(self.__acemoon_local)
daystring='%04d-%02d-%02d'%(year, month,day)
logging.debug('Day string=%s'%daystring)
if daystring in self.__cached_day:
logging.debug('Data of AceMoon for %s already there.'%daystring)
return
for line in f:
sline=line.split()
if len(sline) != 18:
continue
if sline[1].startswith(daystring):
data=AceMoonData(sline)
jd=data.getJulday()
if not self.__cache.hasElement(jd):
self.__cache.add(data.getJulday(),data)
self.__cached_day.append(daystring)
logging.debug('Cached data number=%d'%self.__cache.size())
[docs] def get(self, t):
t_dt = convert(t, datetime.datetime)
self.read(t_dt.year, t_dt.month, t_dt.day)
return self.__cache.getNeighbor(convert(t_dt, Julday))
[docs] def getList(self, start_time, stop_time):
''' Get ACE data at Moon between start_time and stop_time. Data format is JdSeries with AceMoonData class.
'''
# Read all the date specified.
t0=convert(start_time, datetime.datetime)
t1=convert(stop_time, datetime.datetime)
t=t0
while t <= t1:
y=t.year
m=t.month
d=t.day
doi=_md2doi(y,m,d)
self.read(y,m,d)
t=t+datetime.timedelta(0,43200) # Data file segment is 1 days. Half value (0.5 day) is fine if you include redunduncy.
y=t1.year # The last day should be explictly loaded.
m=t1.month
d=t1.day
doi=_md2doi(y,m,d)
self.read(y,m,d)
### Get the time of the data.
tlist = self.__cache.getJuldayList()
t0_jd=convert(start_time, julday.Julday)
t1_jd=convert(stop_time, julday.Julday)
series=JdSeries()
for t in tlist:
if t>=t0_jd and t<=t1_jd:
series.add(self.__cache.getNeighbor(t))
return series
def __s(self, sumX, sumX2, nX):
''' Calculate the standard deviation.
From given sum(X), sum(X^2), and number_of(X), standard deviation,
sqrt( sum(X^2)/nX - (sum(X)/nX)^2 ) is returned.
'''
sbx2=sumX2/nX-sumX/nX*sumX/nX
if sbx2 < 0:
logging.warn('Variation is minus (%g). Probably small numerical error.'%sbx2)
sbx2 = 0
return math.sqrt(sbx2)
[docs] def getAverage(self, t0, t1):
''' Returns the average and the standard deviation of the solar wind parameter.
The returned value is a tuple ( ave_np, std_np, ave_vx, std_vx, ave_vy, std_vy, ave_vz, std_vz, ave_Tp, std_Tp, ave_rHe, std_rHe, ave_bx, std_bx, ave_by, std_by, ave_bz, std_bz )
where rHe is the helium ratio.
'''
fulldata=self.getList(t0,t1)
np=0.
np2=0.
n_np=0
vx=0.
vx2=0.
n_vx=0
vy=0.
vy2=0.
n_vy=0
vz=0.
vz2=0.
n_vz=0
bx=0.
bx2=0.
n_bx=0
by=0.
by2=0.
n_by=0
bz=0.
bz2=0.
n_bz=0
tp=0.
tp2=0.
n_tp=0
pa=0.
pa2=0.
n_pa=0
timefulldata = fulldata.getJuldayList()
for jd in timefulldata:
data=fulldata.getNeighbor(jd)
prm=data.getData()
_np=prm.getNp()
if _np != -9999.9:
np=np+_np
np2=np2+_np*_np
n_np=n_np+1
_vx=prm.getV_gse_x()
if _vx != -9999.9:
vx=vx+_vx
vx2=vx2+_vx*_vx
n_vx=n_vx+1
_vy=prm.getV_gse_y()
if _vy != -9999.9:
vy=vy+_vy
vy2=vy2+_vy*_vy
n_vy=n_vy+1
_vz=prm.getV_gse_z()
if _vz != -9999.9:
vz=vz+_vz
vz2=vz2+_vz*_vz
n_vz=n_vz+1
_tp=prm.getTp()
if _tp != -9999.9:
tp=tp+_tp
tp2=tp2+_tp*_tp
n_tp=n_tp+1
_pa=prm.getAlpha_ratio()
if _pa != -9999.9:
pa=pa+_pa
pa2=pa2+_pa*_pa
n_pa=n_pa+1
_bx=prm.getB_gse_x()
if _bx != -9999.9:
bx=bx+_bx
bx2=bx2+_bx*_bx
n_bx=n_bx+1
_by=prm.getB_gse_y()
if _by != -9999.9:
by=by+_by
by2=by2+_by*_by
n_by=n_by+1
_bz=prm.getB_gse_z()
if _bz != -9999.9:
bz=bz+_bz
bz2=bz2+_bz*_bz
n_bz=n_bz+1
if n_np!=0:
anp=np/n_np
snp=self.__s(np, np2, n_np)
else:
anp=float("nan")
snp=float("nan")
if n_vx!=0:
avx=vx/n_vx
svx=self.__s(vx, vx2, n_vx)
else:
avx=float('nan')
svx=float('nan')
if n_vy!=0:
avy=vy/n_vy
svy=self.__s(vy, vy2, n_vy)
else:
avy=float('nan')
svy=float('nan')
if n_vz!=0:
avz=vz/n_vz
svz=self.__s(vz, vz2, n_vz)
else:
avz=float('nan')
svz=float('nan')
if n_tp!=0:
atp=tp/n_tp
stp=self.__s(tp, tp2, n_tp)
else:
atp=float('nan')
stp=float('nan')
if n_pa!=0:
apa=pa/n_pa
spa=self.__s(pa, pa2, n_pa)
else:
apa=float('nan')
spa=float('nan')
if n_bx!=0:
abx=bx/n_bx
sbx=self.__s(bx, bx2, n_bx)
else:
abx=float('nan')
sbx=float('nan')
if n_by!=0:
aby=by/n_by
sby=self.__s(by, by2, n_by)
else:
aby=float('nan')
sby=float('nan')
if n_bz!=0:
abz=bz/n_bz
sbz=self.__s(bz, bz2, n_bz)
else:
abz=float('nan')
sbz=float('nan')
return (anp,snp,avx,svx,avy,svy,avz,svz,atp,stp,apa,spa,abx,sbx,aby,sby,abz,sbz)
[docs]class Ace:
''' Class for ACE L2 MERGE data.
Data is set by irfpyrc setup file, from the entry [swim] acemerge.
'''
def __init__(self, rcfile=None):
self.__rc = Rc()
if rcfile:
self.__rc.append(rcfile)
self.__aceurl = self.__rc.get('swim', 'acemerge')
logging.debug('ACE data will be retrieved from %s'% self.__aceurl)
self.__cache=JdSeries()
self.__cache_contents=[]
# def readAll(self):
# ''' Read a data from 'AceMerge.txt'. The file can be set through .irfpyrc setting.
# '''
# self.__acefile, status=urllib.urlretrieve(self.__aceurl)
# logging.info('ACE data retrieved. %s with status "%s"'%(self.__acefile, status))
# f=open(self.__acefile)
# self.__cache=JdSeries()
# for l in f:
# try:
# a=AceElem(l)
# jd=a.getJulday()
# if not self.__cache.hasElement(jd):
# self.__cache.add(a.getJulday(), a)
# except ValueError,e:
# continue
[docs] def read(self, year, month, day):
''' Read a data from segmented 'AceMerge.txt_yyyy_doi.gz' file.
'''
doi = _md2doi(year, month, day)
suffix='_%04d_%02dx.gz'%(year, doi/10)
if suffix in self.__cache_contents:
logging.info('Data already loaded for %d-%d-%d (%s)'%(year, month,day, suffix))
return
uri=self.__aceurl+suffix
fname,status=urllib.request.urlretrieve(uri)
logging.info('ACE data retrieved. %s with status "%s"'%(fname, status))
f=gzip.GzipFile(fname, mode='rb')
for l in f:
try:
a=AceElem(l)
jd=a.getJulday()
if not self.__cache.hasElement(jd):
self.__cache.add(a.getJulday(), a)
except ValueError as e:
continue
f.close()
self.__cache_contents.append(suffix)
[docs] def get(self, t):
''' Get ACE data at the time of t. Data format is a JdObject with AceElem class.
'''
t_dt=convert(t, datetime.datetime)
y=t_dt.year
m=t_dt.month
d=t_dt.day
doi = _md2doi(y,m,d)
self.read(y,m,d)
t_jd=convert(t, Julday)
jdobj = self.__cache.getNeighbor(t_jd)
logging.info('Data retrieved: %s'%jdobj)
return jdobj
[docs] def getList(self, start_time, stop_time):
''' Get ACE data between start_time and stop_time. Data format is JdSeries with AceElem class.
'''
# Read all the date specified.
t0=convert(start_time, datetime.datetime)
t1=convert(stop_time, datetime.datetime)
t=t0
while t <= t1:
y=t.year
m=t.month
d=t.day
doi=_md2doi(y,m,d)
self.read(y,m,d)
t=t+datetime.timedelta(5) # Data file length is 10 days. Half value (5 day) is fine if you include redunduncy.
y=t1.year # The last day should be explictly loaded.
m=t1.month
d=t1.day
doi=_md2doi(y,m,d)
self.read(y,m,d)
### Get the time of the data.
tlist = self.__cache.getJuldayList()
t0_jd=convert(start_time, julday.Julday)
t1_jd=convert(stop_time, julday.Julday)
series=JdSeries()
for t in tlist:
if t>=t0_jd and t<=t1_jd:
series.add(self.__cache.getNeighbor(t))
return series
[docs] def getAverage(self, t0, t1):
''' Returns the average and the standard deviation of the solar wind parameter.
The returned value is a tuple ( ave_np, std_np, ave_vx, std_vx, ave_vy, std_vy, ave_vz, std_vz, ave_Tp, std_Tp, ave_rHe, std_rHe, ave_bx, std_bx, ave_by, std_by, ave_bz, std_bz )
where rHe is the helium ratio.
'''
fulldata=self.getList(t0,t1)
np=0.
np2=0.
n_np=0
vx=0.
vx2=0.
n_vx=0
vy=0.
vy2=0.
n_vy=0
vz=0.
vz2=0.
n_vz=0
bx=0.
bx2=0.
n_bx=0
by=0.
by2=0.
n_by=0
bz=0.
bz2=0.
n_bz=0
tp=0.
tp2=0.
n_tp=0
pa=0.
pa2=0.
n_pa=0
for data in fulldata:
jd=data.getJd()
prm=data.getData()
_np=prm.getNp()
if _np != -9999.9:
np=np+_np
np2=np2+_np*_np
n_np=n_np+1
_vx=prm.getV_gse_x()
if _vx != -9999.9:
vx=vx+_vx
vx2=vx2+_vx*_vx
n_vx=n_vx+1
_vy=prm.getV_gse_y()
if _vy != -9999.9:
vy=vy+_vy
vy2=vy2+_vy*_vy
n_vy=n_vy+1
_vz=prm.getV_gse_z()
if _vz != -9999.9:
vz=vz+_vz
vz2=vz2+_vz*_vz
n_vz=n_vz+1
_tp=prm.getTp()
if _tp != -9999.9:
tp=tp+_tp
tp2=tp2+_tp*_tp
n_tp=n_tp+1
_pa=prm.getAlpha_ratio()
if _pa != -9999.9:
pa=pa+_pa
pa2=pa2+_pa*_pa
n_pa=n_pa+1
_bx=prm.getB_gse_x()
if _bx != -9999.9:
bx=bx+_bx
bx2=bx2+_bx*_bx
n_bx=n_bx+1
_by=prm.getB_gse_y()
if _by != -9999.9:
by=by+_by
by2=by2+_by*_by
n_by=n_by+1
_bz=prm.getB_gse_z()
if _bz != -9999.9:
bz=bz+_bz
bz2=bz2+_bz*_bz
n_bz=n_bz+1
if n_np!=0:
anp=np/n_np
snp=math.sqrt(np2/n_np - np/n_np*np/n_np)
else:
anp=float("nan")
snp=float("nan")
if n_vx!=0:
avx=vx/n_vx
svx=math.sqrt(vx2/n_vx-vx/n_vx*vx/n_vx)
else:
avx=float('nan')
svx=float('nan')
if n_vy!=0:
avy=vy/n_vy
svy=math.sqrt(vy2/n_vy - vy/n_vy*vy/n_vy)
else:
avy=float('nan')
svy=float('nan')
if n_vz!=0:
avz=vz/n_vz
svz=math.sqrt(vz2/n_vz - vz/n_vz*vz/n_vz)
else:
avz=float('nan')
svz=float('nan')
if n_tp!=0:
atp=tp/n_tp
stp=math.sqrt(tp2/n_tp - tp/n_tp*tp/n_tp)
else:
atp=float('nan')
stp=float('nan')
if n_pa!=0:
apa=pa/n_pa
spa=math.sqrt(pa2/n_pa - pa/n_pa*pa/n_pa)
else:
apa=float('nan')
spa=float('nan')
if n_bx!=0:
abx=bx/n_bx
sbx=math.sqrt(bx2/n_bx-bx/n_bx*bx/n_bx)
else:
abx=float('nan')
sbx=float('nan')
if n_by!=0:
aby=by/n_by
sby=math.sqrt(by2/n_by - by/n_by*by/n_by)
else:
aby=float('nan')
sby=float('nan')
if n_bz!=0:
abz=bz/n_bz
sbz=math.sqrt(bz2/n_bz - bz/n_bz*bz/n_bz)
else:
abz=float('nan')
sbz=float('nan')
return (anp,snp,avx,svx,avy,svy,avz,svz,atp,stp,apa,spa,abx,sbx,aby,sby,abz,sbz)
[docs]class AceData:
''' Class for ACE Browse data. For L2 data, use Ace class.
You can refer to the data index by class varibales of (YEAR_F, DOY, BX, BY, BZ, BT, BP, B, NPROTON, HEPRATIO, VSW, TSW).
'''
## Default URI of ACE browse data. Original data is in spaghetti:/moon/support/ACE/Browse/PARAM/ACE_Browse_Data.txt
DEFAULT_BROWSE_URI='/Volumes/scidata/data/moon/support/ACE/Browse/PARAM/ACE_Browse_Data.txt'
# DEFAULT_BROWSE_URI='/Volumes/moon/support/ACE/Browse/PARAM/ACE_Browse_Data.txt'
(YEAR_F, DOY, BX, BY, BZ, BT, BP, B, NPROTON, HEPRATIO, VSW, TSW)=list(range(12))
LABEL=('YEAR', 'DOY', 'Bx', 'By', 'Bz', 'Bt', 'Bp', '|B|', 'NH+', 'NHe++/NH+', 'Vsw', 'Tsw')
def __init__(self):
self.__datatable=None
[docs] def setFromUri(self, uri=None):
if uri is None:
uri = self.DEFAULT_BROWSE_URI
logging.debug('URI=%s'%uri)
f=urllib.request.urlretrieve(uri)
logging.debug('Local file=%s'%f[0])
fp=open(f[0], 'r')
buf=''
endhead=re.compile('^BEGIN DATA.*$')
while not endhead.match(buf):
buf=fp.readline()
# logging.debug('buf=%s'%buf)
tbl=numpy.loadtxt(fp)
logging.debug('Data read to matrix: shape=%s'%str(tbl.shape))
# Read data to Nx12 array.
self.__datatable=tbl.copy()
self.__jdser=JdSeries()
for i in range(self.__datatable.shape[0]):
yrf = self.__datatable[i,self.YEAR_F]
yr = int(yrf)
val = self.__datatable[i,self.DOY]
doy = int(val)
val = (val-doy)*24
hr = int(val)
val = (val-hr)*60
mn = int(val)
val = (val-mn)*60
jd = Julday(yr, 1, doy, hr, mn, val)
# logging.debug('JD=%s'%str(jd))
self.__jdser.add(jd, self.__datatable[i, :])
logging.debug('Loaded data size=%d'% self.__jdser.size())
[docs] def getData(self, dt):
jd=dt2jd(dt)
return self.__jdser.getNeighbor(jd)
[docs] def getTimeseries(self, datetime0, datetime1, value=(BX,BY,BZ,NPROTON,VSW)):
'''Returns the JdObject of the dataset (subdataset) for the specified time interval.
@param[in] value Array-like to specify which values are retrieved. Default is (BX,BY,BZ,NPROTON,VSW).
'''
jdlist=sorted(self.__jdser.getJuldayList())
timeser = JdSeries()
tt0 = dt2tt(datetime0)
tt1 = dt2tt(datetime1)
for jd in jdlist:
# dt=jd2dt(jd)
tt = dt2tt(jd.getDatetime())
# logging.debug(dt)
if tt>=tt0 and tt<=tt1:
vals = self.__jdser.getNeighbor(jd).getData()
data = [ vals[i] for i in value ]
# print data
timeser.add(jd, data)
return timeser