update api, add median line

Examples/SimpleSNR-2013-04-14T0854.py

@@ -13,21 +13,24 @@
 'zlim_pl': [None,None],
 'vlim_pl': [72,90],
 'flim_pl': [3.5,5.5],
-'int': True,
+'int': False,
 'samples': True,
 'odir': 'out/2013-04-14T0854',
 'vlim': [30, 60],
-'zlim': (90, None),
-'tlim': ['2013-04-14T08:54:00Z',
-         '2013-04-14T08:54:50Z'],
+'zlim': (200, 350),
+'tlim': ['2013-04-14T08:53:00Z',
+         '2013-04-14T08:56:00Z'],
+#'tlim': ['2013-04-14T08:54:00Z',
+#         '2013-04-14T08:54:50Z'],
+'medthres': 2.,
 'verbose': True,
 #   'tmark': [(datetime(2013,4,14,8,54,30,tzinfo=UTC),300.,'onset',-1),
 #             (datetime(2013,4,14,8,54,41,tzinfo=UTC),300.,'quiescence',1)]
   }
 #%% iterate over list. Files are ID'd by file extension (See README.rst)
 flist = (
-#'d0346834.dt3.h5', # 480 us long pulse
-'d0346834.dt1.h5',
+'d0346834.dt3.h5', # 480 us long pulse
+#'d0346834.dt1.h5',
 #'d0346834.dt0.h5', #alt code
 
 #'20130413.001_ac_30sec.h5',

isrutils/looper.py

@@ -11,7 +11,7 @@
 
 from . import Path,str2dt
 from .switchyard import isrselect
-from .plots import plotsnr,plotplasmaline
+from .plots import plotsnr,plotplasmaline,plotsumionline
 
 def simpleloop(flist,P):
     flist=sorted(flist) #in case glob
@@ -38,7 +38,9 @@ def simpleloop(flist,P):
     for f in flist:
       #  ft = ftype(f)
        # ax[ft] = {}
-        specdown,specup,snrsamp,azel,isrlla,snrint,snr30int = isrselect(Path(P['path'])/f, P)
+        specdown,specup,snrsamp,azel,isrlla,snrint,snr30int,ionsum = isrselect(Path(P['path'])/f, P)
+        # summed ion line over altitude range
+        plotsumionline(ionsum,None,f,P)
         # 15 sec integration
         plotsnr(snrint,f,Pint,azel,ctxt='int_')
         # 200 ms integration

isrutils/overlayISRopt.py

@@ -5,7 +5,7 @@
 matplotlib.use('agg') # NOTE comment out this line to enable visible plots
 #
 from .plasmaline import readplasmaline#,plotplasmaline
-from .summed import sumlongpulse,dojointplot
+from .summed import sumionline,dojointplot
 #
 from GeoData.utilityfuncs import readNeoCMOS
 #
@@ -33,7 +33,7 @@ def overlayisrhist(P):
     spec,freq = readplasmaline(P['isrfn'],P)
     #plotplasmaline(spec,freq,isrfn,P)
 #%% (2-3) read ISR long pulse
-    lpsum,beamazel,isrlla = sumlongpulse(P)
+    lpsum,beamazel,isrlla = sumionline(P['isrfn'],P)
 #%% (4) load optical data
     if optfn is not None:
         #hst = []; hstazel=[]; hstlla=[]; hstut=[]

isrutils/plots.py

@@ -6,12 +6,12 @@
 import h5py
 from datetime import datetime,timedelta
 from pytz import UTC
-from numpy import log10,absolute, meshgrid, sin, radians,unique,atleast_1d
+from numpy import log10,absolute, meshgrid, sin, radians,unique,atleast_1d, median
 from numpy.ma import masked_invalid
 from pandas import DataFrame
 from xarray import DataArray
 #
-from matplotlib.pyplot import figure,subplots
+from matplotlib.pyplot import figure,subplots,gcf
 from matplotlib.dates import MinuteLocator,SecondLocator
 from mpl_toolkits.mplot3d import Axes3D
 from matplotlib.dates import DateFormatter
@@ -401,4 +401,51 @@ def timeticks(tdiff):
         return SecondLocator(bysecond=range(0,60,5)),  SecondLocator(bysecond=range(0,60,2))
 
     else:
-        return SecondLocator(bysecond=range(0,60,2)),  SecondLocator(bysecond=range(0,60,1))
+        return SecondLocator(bysecond=range(0,60,2)),  SecondLocator(bysecond=range(0,60,1))
+
+
+def plotsumionline(dsum,ax,fn,P):
+    if dsum is None:
+        return
+
+    assert isinstance(dsum,DataArray) and dsum.ndim==1,'incorrect input type'
+    assert dsum.size > 1,'must have at least two data points to plot'
+
+    if not ax:
+        fg = figure()
+        ax = fg.gca()
+    else:
+        fg = gcf()
+
+    ax.plot(dsum.time,dsum.values,label='$\sum_{range} |P_{rx}|$')
+
+    if P['verbose']:
+        med = median(dsum.values)
+        ax.axhline(med,color='gold',linestyle='--',label='median')
+        if 'medthres' in P:
+            medthres = P['medthres']*med
+            ax.axhline(medthres,color='red',linestyle='--',label='threshold')
+
+    ax.set_ylabel('summed power')
+    ax.set_xlabel('time [UTC]')
+    ax.set_title('{} summed over ranges ({}..{})km'.format(expfn(fn),P['zlim'][0],P['zlim'][1]))
+
+    ax.set_yscale('log')
+    ax.grid(True)
+    ax.legend()
+
+    fg.autofmt_xdate()
+
+    writeplots(fg, dsum.time[0].item(), P['odir'],'summedAlt')
+
+def plotsumplasmaline(plsum):
+    assert isinstance(plsum,DataArray)
+
+    fg = figure()
+    ax = fg.gca()
+    plsum.plot(ax=ax)
+    ax.set_ylabel('summed power')
+    ax.set_xlabel('time [UTC]')
+    ax.set_title('plasma line summed over altitude (200..350)km and frequency (3.5..5.5)MHz')
+
+    #ax.xaxis.set_major_locator(timeticks(plsum.time[-1]-plsum.time[0])[0])

isrutils/summed.py

@@ -15,8 +15,9 @@
 import matplotlib.animation as anim
 #
 from .plasmaline import readplasmaline
-from .common import timeticks,findindex2Dsphere,timesync,projectisrhist
+from .common import findindex2Dsphere,timesync,projectisrhist
 from .snrpower import readpower_samples
+from .plots import plotsumionline
 #
 from GeoData.plotting import plotazelscale
 
@@ -37,7 +38,7 @@ def dojointplot(ds,spec,freq,beamazel,optical,optazel,optlla,isrlla,heightkm,uto
     gs = gridspec.GridSpec(2, 1, height_ratios=[3,1])
 #%% setup radar plot(s)
     a1 = fg.add_subplot(gs[1])
-    plotsumlongpulse(ds,a1,expfn(P['isrfn']),P['zlim'])
+    plotsumionline(ds,a1,expfn(P['isrfn']),P['zlim'])
 
     h1 = a1.axvline(nan,color='k',linestyle='--')
     t1 = a1.text(0.05,0.95,'time=',transform=a1.transAxes,va='top',ha='left')
@@ -124,34 +125,11 @@ def compclim(imgs,lower=0.5,upper=99.9,Nsamples=50):
     return clim
 
 #%% dt3
-def sumlongpulse(P):
-    snrsamp,azel,lla = readpower_samples(P['isrfn'],P)
+def sumionline(fn,P):
+    snrsamp,azel,lla = readpower_samples(fn,P)
     assert isinstance(snrsamp,DataArray)
 
-    return snrsamp.sum(axis=0),azel,lla
-
-def plotsumlongpulse(dsum,ax,rmode,zlim):
-    assert isinstance(dsum,DataArray) and dsum.ndim==1,'incorrect input type'
-    assert dsum.size > 1,'must have at least two data points to plot'
-
-    if not ax:
-        fg = figure()
-        ax = fg.gca()
-
-    try:
-        dsum.plot(ax=ax)
-    except ValueError as e:
-        print('Windows bug?ValueError: ordinal must be >= 1     your error is {}'.format(e) )
-
-    ax.set_ylabel('summed power')
-    ax.set_xlabel('time [UTC]')
-    ax.set_title('{} summed over altitude ({}..{})km'.format(rmode,zlim[0],zlim[1]))
-
-    ax.set_yscale('log')
-    ax.grid(True)
-
-    ax.xaxis.set_major_locator(timeticks(dsum.time[-1] - dsum.time[0])[0])
-    return ax
+    return snrsamp.sum(dim='srng'),azel,lla
 
 #%% plasma line
 def sumplasmaline(fn,P):
@@ -166,18 +144,6 @@ def sumplasmaline(fn,P):
 
     for s in spec:
         find = (P['flim'][0] <= absolute(freq[s]/1.e6)) & (absolute(freq[s]/1.e6) < P['flim'][1])
-        specsum.loc[:,s] = spec.loc[:,:,zind,find].sum(axis=3).sum(axis=2)
+        specsum.loc[:,s] = spec.loc[:,:,zind,find].sum(axis=3).sum(axis=2) #FIXME .sum(dim=)
 
     return specsum
-
-def plotsumplasmaline(plsum):
-    assert isinstance(plsum,DataArray)
-
-    fg = figure()
-    ax = fg.gca()
-    plsum.plot(ax=ax)
-    ax.set_ylabel('summed power')
-    ax.set_xlabel('time [UTC]')
-    ax.set_title('plasma line summed over altitude (200..350)km and frequency (3.5..5.5)MHz')
-
-    ax.xaxis.set_major_locator(timeticks(plsum.time[-1]-plsum.time[0])[0])

isrutils/switchyard.py

@@ -6,10 +6,10 @@
 from .rawacf import readACF
 from .plasmaline import readplasmaline
 from .snrpower import readpower_samples,readsnr_int,snrvtime_fit
+from .summed import sumionline
 from .plots import plotsnr,plotsnr1d,plotplasmaline
 
 
-
 def isrstacker(flist,P):
 
     for fn in flist:
@@ -54,11 +54,11 @@ def isrselect(fn,P):
 #%% handle path, detect file type
     ft = ftype(fn)
 #%% plasma line
-    specdown=specup=None
+    specdown=None; specup=None
     if ft in ('dt1','dt2'):
         specdown,specup,azel = readplasmaline(fn,P)
 #%% ~ 200 millisecond raw altcode and longpulse
-    snrsamp=isrlla=None
+    snrsamp=None; isrlla=None
     if ft in ('dt0','dt3'):
 #        tic = time()
         snrsamp,azel,isrlla = readpower_samples(fn,P)
@@ -71,13 +71,14 @@ def isrselect(fn,P):
         if P['verbose']:
             print('ACF/PSD read & plot took {:.1f} sec.'.format(time()-tic))
 #%% multi-second integration (numerous integrated pulses)
-    snrint=None
+    snrint=None; ionsum=None
     if ft in ('dt0','dt3'):
         snrint = readsnr_int(fn,P['beamid'])
+        ionsum,azel,isrlla = sumionline(fn,P)
 #%% 30 second integration plots
     if fn.stem.rsplit('_',1)[-1] == '30sec':
         snr30int = snrvtime_fit(fn,P['beamid'])
     else:
         snr30int=None
 
-    return specdown,specup,snrsamp,azel,isrlla,snrint,snr30int
+    return specdown,specup,snrsamp,azel,isrlla,snrint,snr30int,ionsum