|
#!/usr/bin/env ipython |
|
# -*- coding: utf-8 -*- |
|
import numpy as np |
|
from numpy import ( |
|
nanmean, nanstd, |
|
) |
|
import argparse, os |
|
from h5py import File as h5 |
|
from datetime import datetime, timedelta |
|
from shared_funcs.funcs import ( |
|
My2DArray, date2utc, utc2date, rebine, ncep_gh, nans |
|
) |
|
import itertools |
|
#--- we need X to build the figure |
|
if 'DISPLAY' in os.environ: |
|
from pylab import figure, show, close |
|
else: |
|
raise SystemExit( |
|
" >>> ERROR: the DISPLAY environment variable doesn't exist!" |
|
) |
|
|
|
#--- retrieve args |
|
parser = argparse.ArgumentParser( |
|
formatter_class=argparse.ArgumentDefaultsHelpFormatter |
|
) |
|
parser.add_argument( |
|
'-y', '--years', |
|
type=int, |
|
nargs=2, |
|
default=[2006,2013], |
|
help='years for period of data-processing. So that the profile\ |
|
will be built from YearIni/Jan/01 until YearEnd/Dec/31)', |
|
) |
|
parser.add_argument( |
|
'-fid', '--fname_inp', |
|
type=str, |
|
default='./test.h5', |
|
help='filename of .h5 input (contains geopotential-height-corrected Scalers)', |
|
) |
|
parser.add_argument( |
|
'-obs', '--obs', |
|
type=str, |
|
nargs='+', |
|
default='wAoP_wPrs_wGh', |
|
help='observables to process; can be more than one, and must be one\ |
|
of the keys inside the .h5 file. They will be plotted in the\ |
|
same figure in superposed fashion.', |
|
) |
|
parser.add_argument( |
|
'-nb', '--nbin', |
|
type=int, |
|
default=100, |
|
help='total number of bins for time domain in the long-term profile' |
|
) |
|
parser.add_argument( |
|
'-fig', '--fname_fig', |
|
type=str, |
|
default='test.png', |
|
help='filename for output figure. Or "screen" for interactive plot.' |
|
) |
|
|
|
pa = parser.parse_args() |
|
|
|
#--- |
|
finp_d = h5(pa.fname_inp,'r') |
|
date_ini = datetime(pa.years[0], 1, 1) |
|
date_end = datetime(pa.years[1], 12, 31) |
|
assert date_end>date_ini, '`date_end` must be greater than `date_ini`!!' |
|
|
|
#--- |
|
nbin = pa.nbin |
|
|
|
class d(object): |
|
""" |
|
data buffer for each bin of the profile |
|
""" |
|
def __init__(self): |
|
self.n = 0 # number of original data points contributing |
|
self.rates = My2DArray(1, dtype=np.float32) |
|
self.mean = np.nan |
|
self.std = np.nan |
|
|
|
#obs = pa.obs #'wAoP_wPrs_wGh' |
|
utc_ini = date2utc(datetime(1970,1,1,0,0,0)) |
|
ini,end=datetime(pa.years[0],1,1),datetime(pa.years[1],12,31) |
|
Dt = (end-ini).total_seconds() |
|
dt = Dt/nbin |
|
|
|
# list of all pairs (year,month) of interest |
|
YYYY_MM = list(itertools.product( |
|
range(pa.years[0],pa.years[1]+1,1), |
|
range(1,12+1,1), |
|
)) |
|
|
|
dd_ = {} # several dd's |
|
for obs in pa.obs: |
|
# buffer for each observable |
|
dd = [ d() for i in range(nbin) ] |
|
for yyyy, mm in YYYY_MM: |
|
rpath='%04d/%02d' % (yyyy,mm) |
|
tutc = finp_d['t_utc/'+rpath][...] # [utc sec] |
|
r = finp_d[obs+'/'+rpath][...] # rate |
|
ind = np.array((tutc-date2utc(ini))/dt, dtype=int) # (*) |
|
ind_list = np.unique(ind) # (**) |
|
print yyyy, mm, ind_list |
|
if ind_list.size>1: |
|
assert ind_list[0]<ind_list[1],"wrong."# consistency |
|
for i in ind_list: |
|
cc = ind==i # collect all with similar "time bins" |
|
ri = r[cc] # rates of this i-block |
|
n = dd[i].n |
|
dd[i].n += ri.size |
|
dd[i].rates.resize_rows(n+ri.size) |
|
dd[i].rates[n:n+ri.size] = ri |
|
|
|
if len(ind_list)==1: # [likely] |
|
continue |
|
else: # [unlikely] |
|
for j in ind_list[:-1]: |
|
n = dd[j].n |
|
dd[j].rates = dd[j].rates[:n] |
|
dd[j].mean = nanmean(dd[j].rates) |
|
dd[j].std = nanstd(dd[j].rates) |
|
dd[j].rates = None # delete data |
|
dd_[obs] = dd |
|
""" |
|
(*): list of integers corresponding to the nearest integer value |
|
of the utc-sec discretized in units of `dt` time-windows. |
|
For example, taking utc-sec="equivalent to 09/nov/2008 05:40", |
|
`ini=datetime(2008,nov,1,0,0)`, and `dt=1day` results in |
|
an `ind` value of 9. |
|
(**): to collect all with similar `ind` values |
|
""" |
|
|
|
# build the figure |
|
do = {} |
|
fig = figure(1, figsize=(8,4)) |
|
ax = fig.add_subplot(111) |
|
for obs in pa.obs: |
|
prof = np.array([dd_[obs][n].mean for n in range(nbin)]) |
|
t, r = np.arange(prof.size)*dt/(86400.*365), 100.*prof |
|
do[obs] = np.array([t,r]).T |
|
ax.plot(t, r, '-o', ms=2, alpha=.7, label=obs) |
|
np.savetxt('test_%s.txt'%obs,do[obs],fmt='%12.2f') |
|
|
|
ax.legend(loc='best') |
|
ax.grid(True) |
|
ax.set_xlabel('years since '+date_ini.strftime('%b/%Y')) |
|
ax.set_ylabel('count rate') |
|
if pa.fname_fig=='screen': |
|
print " >>> showing interactively..." |
|
show() |
|
else: |
|
fig.savefig(pa.fname_fig, dpi=135, bbox_inches='tight') |
|
print "\n ---> we generated: "+pa.fname_fig |
|
|
|
close(fig) |
|
|
|
|
|
#EOF |
