diff --git a/Examples/SimpleSNR-2013-04-14T0854.py b/Examples/SimpleSNR-2013-04-14T0854.py index 47b3eb3..fa48a6f 100755 --- a/Examples/SimpleSNR-2013-04-14T0854.py +++ b/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', diff --git a/isrutils/looper.py b/isrutils/looper.py index 9f313aa..a91af56 100644 --- a/isrutils/looper.py +++ b/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 diff --git a/isrutils/overlayISRopt.py b/isrutils/overlayISRopt.py index c9cb790..fb1b2ec 100644 --- a/isrutils/overlayISRopt.py +++ b/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=[] diff --git a/isrutils/plots.py b/isrutils/plots.py index 531cc4f..ef0d801 100644 --- a/isrutils/plots.py +++ b/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)) \ No newline at end of file + 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]) \ No newline at end of file diff --git a/isrutils/summed.py b/isrutils/summed.py index 676469a..a345070 100644 --- a/isrutils/summed.py +++ b/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]) diff --git a/isrutils/switchyard.py b/isrutils/switchyard.py index bc973da..8432697 100644 --- a/isrutils/switchyard.py +++ b/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