In [1]:
%run basics

In [2]:
def match_masks(series_list):
    if len(series_list)<2:
        print "match_masks: only 1 series, nothing to do"
        return
    mask=numpy.ma.mask_or(series_list[0].mask,series_list[1].mask)
    if len(series_list)>2:
        for item in series_list[2:]:
            mask=numpy.ma.mask_or(mask,item.mask)
    for item in series_list:
        item.mask=mask
    return

In [3]:
access_name=qcio.get_filename_dialog()

In [4]:
aws_name=qcio.get_filename_dialog()

In [5]:
tower_name=qcio.get_filename_dialog()

In [6]:
ds={}
ds["access"]=qcio.nc_read_series(access_name)
ds["aws"]=qcio.nc_read_series(aws_name)
ds["tower"]=qcio.nc_read_series(tower_name)


ERROR:qc.utils:get_UTCfromlocaltime: time_zone not in global attributes, checking elsewhere ...
ERROR:qc.utils:get_UTCfromlocaltime: time_zone not in global attributes, checking elsewhere ...

In [11]:
ts = int(ds["tower"].globalattributes["time_step"])
ldt={"tower":{},"access":{},"aws":{}}
for item in ["tower","access","aws"]:
    ldt[item]["30minute"]=ds[item].series["DateTime"]["Data"]
ldt_start=ldt["tower"]["30minute"][0]
ldt_end=ldt["tower"]["30minute"][-1]
for item in ["access","aws"]:
    ldt_start=max(ldt_start,ldt[item]["30minute"][0])
    ldt_end=min(ldt_end,ldt[item]["30minute"][-1])
print ldt_start,ldt_end
si={}
ei={}
for item in ["tower","access","aws"]:
    si[item]=qcutils.GetDateIndex(ldt[item]["30minute"],str(ldt_start),match="startnextday",ts=ts)
    ei[item]=qcutils.GetDateIndex(ldt[item]["30minute"],str(ldt_end),match="endpreviousday",ts=ts)
    ldt[item]["30minute"]=ldt[item]["30minute"][si[item]:ei[item]+1]
    print ldt[item]["30minute"][0],ldt[item]["30minute"][-1]


2011-01-01 08:00:00 2013-11-30 23:00:00
2011-01-02 00:30:00 2013-11-30 00:00:00
2011-01-02 00:30:00 2013-11-30 00:00:00
2011-01-02 00:30:00 2013-11-30 00:00:00

In [12]:
nPerHr = int(float(60)/ts+0.5)
nPerDay = int(float(24)*nPerHr+0.5)
nDays = len(ldt["tower"]["30minute"])/nPerDay

In [19]:
RH={"tower":{},"access":{},"aws":{}}
Ta={"tower":{},"access":{},"aws":{}}
ps={"tower":{},"access":{},"aws":{}}
q={"tower":{},"access":{},"aws":{}}
for item,label in zip(["tower","access","aws"],["RH","RH_11","RH_0"]):
    RH[item]["30minute"],f,a=qcutils.GetSeriesasMA(ds[item],label,si=si[item],ei=ei[item])
for item,label in zip(["tower","access","aws"],["Ta","Ta_11","Ta_0"]):
    Ta[item]["30minute"],f,a=qcutils.GetSeriesasMA(ds[item],label,si=si[item],ei=ei[item])
for item,label in zip(["tower","access","aws"],["ps","ps_11","ps_0"]):
    ps[item]["30minute"],f,a=qcutils.GetSeriesasMA(ds[item],label,si=si[item],ei=ei[item])
for item,label in zip(["tower","access"],["q","q_11"]):
    q[item]["30minute"],f,a=qcutils.GetSeriesasMA(ds[item],label,si=si[item],ei=ei[item])
q["aws"]["30minute"]=mf.specifichumidityfromRH(RH["aws"]["30minute"],Ta["aws"]["30minute"],ps["aws"]["30minute"])
ldt["tower"]["day"]=ldt["tower"]["30minute"][0::nPerDay]
match_masks([q["tower"]["30minute"],q["access"]["30minute"],q["aws"]["30minute"]])
for item in ["tower","access","aws"]:
    q[item]["day"]=q[item]["30minute"].reshape(nDays,nPerDay)
    q[item]["day_avg"]=numpy.ma.average(q[item]["day"],axis=1)
for item in ["tower","access","aws"]:
    RH[item]["day"]=RH[item]["30minute"].reshape(nDays,nPerDay)
    RH[item]["day_avg"]=numpy.ma.average(RH[item]["day"],axis=1)
for item in ["tower","access","aws"]:
    Ta[item]["day"]=Ta[item]["30minute"].reshape(nDays,nPerDay)
    Ta[item]["day_avg"]=numpy.ma.average(Ta[item]["day"],axis=1)

In [15]:
fig=plt.figure()
plt.plot(ldt["tower"]["day"],q["tower"]["day_avg"],'b')
plt.plot(ldt["tower"]["day"],q["access"]["day_avg"],'r+')
plt.plot(ldt["tower"]["day"],q["aws"]["day_avg"],'go')
plt.show()

In [18]:
fig=plt.figure()
plt.plot(ldt["tower"]["day"],RH["tower"]["day_avg"],'b')
plt.plot(ldt["tower"]["day"],RH["access"]["day_avg"],'r+')
plt.plot(ldt["tower"]["day"],RH["aws"]["day_avg"],'go')
plt.show()

In [20]:
fig=plt.figure()
plt.plot(ldt["tower"]["day"],Ta["tower"]["day_avg"],'b')
plt.plot(ldt["tower"]["day"],Ta["access"]["day_avg"],'r+')
plt.plot(ldt["tower"]["day"],Ta["aws"]["day_avg"],'go')
plt.show()

In [ ]: