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)
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]
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 [ ]: