In [35]:
import qcio
import qcutils
import numpy
import statsmodels.api as sm
import scipy.interpolate as inter
inname=qcio.get_filename_dialog()
ds=qcio.nc_read_series(inname)
ts = int(ds.globalattributes['time_step'])
DateTime = ds.series['DateTime']['Data']
si = qcutils.GetDateIndex(DateTime,str(DateTime[0]),ts=ts,default=0,match='startnextday')
ei = qcutils.GetDateIndex(DateTime,str(DateTime[-1]),ts=ts,default=-1,match='endpreviousday')
DateTime = DateTime[si:ei+1]
nrecs = len(DateTime)
nperhr = int(float(60)/ts+0.5)
nperday = int(float(24)*nperhr+0.5)
ndays = nrecs/nperday
nrecs=ndays*nperday
Ah_HMP,f=qcutils.GetSeriesasMA(ds,"Ah_HMP_3m",si=si,ei=ei)
Ah_7500,f=qcutils.GetSeriesasMA(ds,"Ah_7500_Av",si=si,ei=ei)
mask = numpy.ma.mask_or(Ah_7500.mask,Ah_HMP.mask)
Ah_7500 = numpy.ma.array(Ah_7500,mask=mask)
Ah_HMP = numpy.ma.array(Ah_HMP,mask=mask)
Ah_7500_2d = numpy.ma.reshape(Ah_7500,[ndays,nperday])
Ah_HMP_2d = numpy.ma.reshape(Ah_HMP,[ndays,nperday])
Ah_7500_daily_std = numpy.ma.std(Ah_7500_2d,axis=1)
Ah_HMP_daily_std = numpy.ma.std(Ah_HMP_2d,axis=1)
Ah_stdratio_daily = Ah_HMP_daily_std/Ah_7500_daily_std
DT_daily = DateTime[0:nrecs:nperday]
slope = numpy.ones(ndays)
offset = numpy.zeros(ndays)
correl = numpy.ones(ndays)
number = numpy.zeros(ndays)
for i in range(0,ndays-1):
x = Ah_7500_2d[i,:]
y = Ah_HMP_2d[i,:]
x_nm = numpy.ma.compressed(x)
x_nm = sm.add_constant(x_nm,prepend=False)
y_nm = numpy.ma.compressed(y)
if len(y_nm)>1:
resrlm = sm.RLM(y_nm,x_nm,M=sm.robust.norms.TukeyBiweight()).fit()
coefs = resrlm.params
r = numpy.ma.corrcoef(x,y)
number[i] = numpy.ma.count(x)
slope[i] = coefs[0]
offset[i] = coefs[1]
correl[i] = r[0][1]
correl2 = numpy.ma.masked_where((correl<0.95)|(number<15),correl)
number2 = numpy.ma.masked_where((correl<0.95)|(number<15),number)
slope2 = numpy.ma.masked_where((correl<0.95)|(number<15),slope)
offset2 = numpy.ma.masked_where((correl<0.95)|(number<15),offset)
sdratio2 = numpy.ma.masked_where((correl<0.95)|(number<15),Ah_stdratio_daily)
x=numpy.arange(0,len(DT_daily))
x2=numpy.ma.array(x,mask=sdratio2.mask)
x2_nm=numpy.ma.compressed(x2)
y2_nm=numpy.ma.compressed(sdratio2)
s2 = inter.UnivariateSpline (x2_nm, y2_nm, s=1)
fig=figure()
plot(DT_daily,sdratio2,'b.')
plot(DT_daily,s2(x),'r--')
Out[35]:
In [ ]: