swim_plot_monthly_swΒΆ
This script plots monthly solar wind spectrogram into one picture.
#!/usr/bin/env python
r''' This script plots monthly solar wind spectrogram into one picture.
'''
import os
import sys
import math
import logging
logging.basicConfig()
logger = logging.getLogger(os.path.basename(sys.argv[0]))
logger.setLevel(logging.DEBUG)
import datetime
import pickle as pickle
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import numpy as np
import irfpy.swim.ace
import irfpy.cy1orb.MoonPos
import irfpy.cy1orb.Cy1Orbit
import irfpy.cy1orb.Cy1OrbitNr
from irfpy.util.julday import Julday
import irfpy.swim.swim_start
import irfpy.swim.SwimEnergy
datafilename = 'swim_plot_monthly_sw.dat'
def get_deflux(tlist):
''' return the defeferntial flux (energy-wise) for the given tlist. The returned is 1-D 16 element array
'''
maxflx = np.zeros([16])
for t in tlist:
flx = np.array(irfpy.swim.swim_start.getdeflux(t, filter=[irfpy.swim.swim_start.filter.remove_background_5numsum, irfpy.swim.swim_start.filter.remove_one_count]).getData())
flx = flx.max(axis=1) # Along D. Now E=16 dependence.
for ie in range(16):
if flx[ie] > maxflx[ie]:
maxflx[ie] = flx[ie]
return maxflx
def getace(datafilename):
logger = logging.getLogger('getace')
logger.setLevel(logging.DEBUG)
# ACE data
if not os.path.exists(datafilename):
a = irfpy.swim.ace.AceMoon()
ace_ene = {}
c = irfpy.cy1orb.Cy1OrbitNr.Cy1OrbitNr()
c.setFromDefaultUri()
for inr in range(700, 3300):
t0 = c.getStartTime(inr)
t1 = c.getStopTime(inr)
vel = a.getAverage(t0, t1)[2]
if math.isnan(vel):
logger.warn('No data for %d' % inr)
continue
ace_ene[inr] = (vel / 13.8) ** 2
logger.debug('%d %g' % (inr, vel))
fff = open(datafilename, 'w')
pickle.dump(ace_ene, fff)
else:
fff = open(datafilename)
ace_ene = pickle.load(fff)
ace_lene = np.zeros([3500])
for inr in range(3500):
ace_lene[inr] = ace_ene.get(inr, 1e50)
ace_lene = np.ma.log10(ace_lene)
ace_lene = np.ma.masked_greater(ace_lene, 40)
return ace_lene
def main(*args, **kwds):
'''
'''
espec = getespec(datafilename + '_espec')
moonphase = getmoonphase(datafilename + '_moonphase')
ace_lene = getace(datafilename + '_ace')
xarr = list(range(0, 3500))
# Plotting part
especarr = []
for inr in xarr:
if inr in espec:
es = espec[inr][1] + 1e-5 # 1e-5 to avoid zero.
else:
es = np.zeros([16]) - 1.
especarr.append(list(es))
especarr = np.array(especarr)
# yarr = range(16) # Simple index.
etbl = irfpy.swim.SwimEnergy.EnergyTable.table()
yarr = etbl.getTable() # Log of energy
yarr = np.log10(yarr)
X, Y = np.meshgrid(xarr, yarr)
Zm = np.ma.masked_less_equal(especarr, 0)
Zm = np.log10(Zm)
print(Zm.max(), Zm.min())
print(moonphase)
mp = np.zeros([max(moonphase.keys())])
for idx in range(max(moonphase.keys())):
mp[idx] = moonphase.get(idx, -1e30)
# Calculate the fullmoon orbit.
fullmoon = []
for idx in range(len(mp) - 1):
if mp[idx] != None and mp[idx+1] != None:
if mp[idx] > mp[idx+1]:
print(idx, mp[idx], idx+1, mp[idx+1])
fullmoon.append(idx+1)
fig = plt.figure(figsize=(9, 14))
nplt = len(fullmoon) - 1
for imon in range(nplt):
orb_start = fullmoon[imon]
orb_end = fullmoon[imon+1]
ax = fig.add_subplot(nplt, 1, imon+1)
X, Y = np.meshgrid(xarr[orb_start:orb_end], yarr)
Z = Zm[orb_start:orb_end, :]
print(X.shape, X.max(), X.min())
print(Y.shape, Y.max(), Y.min())
print(Z.shape, Z.max(), Z.min(), Z.mask.all())
if Z.mask.all() == True: # No data included in Z, then, skip
continue
pc = ax.pcolor(X, Y, Z.T, vmax=10, vmin=2)
ax.plot(xarr[orb_start:orb_end], ace_lene[orb_start:orb_end], 'r--')
ax.set_ylabel('Log E[eV]')
fig.colorbar(pc)
ax.set_xlabel('Orbit number')
fig.savefig('a.png')
def readdata():
logger.info('Loading from %s' % datafilename)
f = open(datafilename)
espec = pickle.load(f)
moonphase = pickle.load(f)
return espec, moonphase
def getespec(datafilename):
logger.info('Making data from scratch to save to %s' % datafilename)
if os.path.exists(datafilename):
f = open(datafilename)
espec = pickle.load(f)
f.close()
else:
# Dt is 10 min (pls and minus)
dt = datetime.timedelta(seconds=600)
# Orbit number information
onr = irfpy.cy1orb.Cy1OrbitNr.Cy1OrbitNr()
onr.setFromDefaultUri()
orb = irfpy.cy1orb.Cy1Orbit.Cy1Orbit()
# A cache of the energy spectra. [orbnr] = np.array([16, 16])
espec_cache = {}
for inr in range(onr.getFirstOrbitNr(), onr.getLastOrbitNr()+1):
# for inr in range(3020, 3023):
try:
tc = Julday(orb.getDayEqOfOrb(inr)[0]).getDatetime()
except IOError:
continue
t0 = tc - dt
t1 = tc + dt
# Getting the SWIM data availability.
tlist = irfpy.swim.swim_start.getobstime(timerange=[t0, t1])
logger.debug(' Length of data = %d (%s-%s)' % (len(tlist), t0, t1))
if len(tlist) == 0:
continue
espec = get_deflux(tlist) # Energy spectra for the "orbit".
espec_cache[inr] = (tc, espec) # inr -> (time, spectra[16])
f = open(datafilename, 'w')
pickle.dump(espec_cache, f)
f.close()
return espec
def getmoonphase(datafilename):
if os.path.exists(datafilename):
f = open(datafilename)
moonphase = pickle.load(f)
f.close()
else:
orb = irfpy.cy1orb.Cy1Orbit.Cy1Orbit()
pos = irfpy.cy1orb.MoonPos.MoonGse()
pos.setFromDefaultUri()
moonphase = {}
for inr in range(onr.getFirstOrbitNr(), onr.getLastOrbitNr()+1):
# for inr in range(3020, 3050):
try:
teq = orb.getDayEqOfOrb(inr)
teq = Julday(teq[0])
teq = teq.getDatetime()
except IOError:
continue
moongse = pos.getGse(teq)
moonph = math.atan2(moongse[1], moongse[0]) * 180 / math.pi
moonphase[inr] = moonph
print(inr, moonph)
f = open(datafilename, 'w')
pickle.dump(moonphase, f)
f.close()
return moonphase
if __name__ == '__main__':
logger.debug('Start of the program %s' % sys.argv[0])
main(sys.argv[1:])