In [1]:
%run basics
%matplotlib
import collections
import qcrpLL
from matplotlib.patches import Polygon
import xlrd
In [2]:
def normal_dist(bins,mu,sigma):
return 1/(sigma*numpy.sqrt(2*numpy.pi))*numpy.exp(-(bins-mu)**2 /(2*sigma**2))
In [3]:
def colour_boxes(ax,bp,data,nsites,nalts):
boxColors = {1:['darkkhaki'],2:['darkkhaki','royalblue'],
3:['red','darkkhaki','royalblue']}
numBoxes = nsites*nalts
medians = list(range(numBoxes))
for i in range(numBoxes):
box = bp['boxes'][i]
boxX = []
boxY = []
for j in range(5):
boxX.append(box.get_xdata()[j])
boxY.append(box.get_ydata()[j])
boxCoords = list(zip(boxX, boxY))
# Alternate between Dark Khaki and Royal Blue
k = i % nalts
boxPolygon = Polygon(boxCoords, facecolor=boxColors[nalts][k])
ax.add_patch(boxPolygon)
# Now draw the median lines back over what we just filled in
med = bp['medians'][i]
medianX = []
medianY = []
for j in range(2):
medianX.append(med.get_xdata()[j])
medianY.append(med.get_ydata()[j])
ax.plot(medianX, medianY, 'k')
medians[i] = medianY[0]
In [4]:
def get_data_xlsheet(sheet,frow,label):
header = sheet.row_values(frow-1)
col = header.index(label)
lrow = int(sheet.nrows)
a = numpy.array(active_sheet.col_values(col)[frow:lrow])
idx = numpy.where(~numpy.isnan(a)&(a!=-9999))[0]
return a[idx]
In [5]:
mu,sigma = 0.41,0.1
s = numpy.random.normal(mu,sigma,100)
idx = numpy.where(s>0)[0]
s2 = s[idx]
fig,axs = plt.subplots(2,figsize=(6,6))
count, bins, ignored = axs[0].hist(s,10,normed=True)
axs[0].plot(bins,normal_dist(bins,mu,sigma),linewidth=2, color='r')
count, bins, ignored = axs[1].hist(s2,10,normed=True)
axs[1].plot(bins,normal_dist(bins,mu,sigma),linewidth=2, color='r')
plt.show()
In [6]:
cf_name = "controlfiles/Whroo/Portal/L6_uncertainty.txt"
cf = qcio.get_controlfilecontents(cf_name)
in_name = qcio.get_infilenamefromcf(cf)
ds = qcio.nc_read_series(in_name)
ts = int(ds.globalattributes["time_step"])
ntsperday = float(24)*float(60)/float(ts)
ldt = ds.series["DateTime"]["Data"]
qcrp.ParseL6ControlFile(cf,ds)
In [7]:
months = collections.OrderedDict()
months["ER_LT"] = collections.OrderedDict()
months["GPP_LT"] = collections.OrderedDict()
months["NEE_LT"] = collections.OrderedDict()
info = {"window_length":30,"window_offset":5,"ts":ts}
days_at_beginning = float(info["window_length"])/2 - float(info["window_offset"])/2
int_list = ["rb_int","E0_int"]
tts_list = ["rb_tts","E0_tts"]
si = qcutils.GetDateIndex(ldt,str(ldt[0]),ts=ts,default=0,match="startnextmonth")
ei = qcutils.GetDateIndex(ldt,str(ldt[-1]),ts=ts,default=len(ldt)-1,match="endpreviousmonth")
ldt = ldt[si:ei+1]
Fsd,Fsd_flag,Fsd_attr = qcutils.GetSeriesasMA(ds,"Fsd",si=si,ei=ei)
Ta,Ta_flag,Ta_attr = qcutils.GetSeriesasMA(ds,"Ta",si=si,ei=ei)
if "solar_altitude" not in ds.series.keys(): qcts.get_synthetic_fsd(ds)
Fsd_syn,f,a = qcutils.GetSeriesasMA(ds,"Fsd_syn",si=si,ei=ei)
sa,f,a = qcutils.GetSeriesasMA(ds,"solar_altitude",si=si,ei=ei)
ustar,ustar_flag,ustar_attr = qcutils.GetSeriesasMA(ds,"ustar",si=si,ei=ei)
Fc,Fc_flag,Fc_attr = qcutils.GetSeriesasMA(ds,"Fc",si=si,ei=ei)
Fc_ok = numpy.ma.masked_where((Fc_flag!=0)|(ustar_flag!=0)|(Ta_flag!=0),Fc,copy=True)
Ta_ok = numpy.ma.masked_where((Fc_flag!=0)|(ustar_flag!=0)|(Ta_flag!=0),Ta,copy=True)
In [8]:
nthresholds = len(s2)
for i,ustar_threshold in enumerate(s2):
data = {}
if (i%10)==0: print datetime.datetime.now(),i
#print datetime.datetime.now(),i
ustar_dict = qcrp.get_ustar_thresholds_annual(ldt,ustar_threshold)
# get the ustar indicator (1==>turbulent)
turb_ind = qcrp.get_turbulence_indicator_ustar(ldt,ustar,ustar_dict,ts)
Fc_filtered = numpy.ma.masked_where(turb_ind["values"]==0,Fc_ok,copy=True)
Ta_filtered = numpy.ma.masked_where(turb_ind["values"]==0,Ta_ok,copy=True)
# get the daytime indicator (1==>daytime)
day_ind = qcrp.get_day_indicator(cf,Fsd,Fsd_syn,sa)
ER_night = numpy.ma.masked_where(day_ind["values"]==1,Fc_filtered,copy=True)
Ta_night = numpy.ma.masked_where(day_ind["values"]==1,Ta_filtered,copy=True)
# get rb and E0 for the windows
LT_results = qcrpLL.get_LT_params(ldt,ER_night,Ta_night,info)
# interpolate rb and E0 across missing windows
LT_results["rb_int"] = qcrpLL.interp_params(LT_results["rb"])
LT_results["E0_int"] = qcrpLL.interp_params(LT_results["E0"])
# repeat windowed rb and Eo across the tower time step
for tts_item,int_item in zip(tts_list,int_list):
beginning = numpy.ones(days_at_beginning*ntsperday)*LT_results[int_item][0]
middle = numpy.repeat(LT_results[int_item],info["window_offset"]*ntsperday)
nend = len(ldt) - (len(beginning)+len(middle))
end = numpy.ones(nend)*LT_results[int_item][-1]
LT_results[tts_item] = numpy.concatenate((beginning,middle,end))
# get the Lloyd-Taylor ecosystem respiration
ER_LT_all = qcrpLL.ER_LloydTaylor(Ta,LT_results["rb_tts"],LT_results["E0_tts"])
# merge with observations
data["ER_LT"] = numpy.ma.copy(ER_night)
idx = numpy.where(numpy.ma.getmaskarray(data["ER_LT"])==True)[0]
data["ER_LT"][idx] = ER_LT_all[idx]
# get the NEE
data["NEE_LT"] = numpy.ma.copy(Fc)
idx = numpy.where(day_ind["values"]==0)[0]
data["NEE_LT"][idx] = data["ER_LT"][idx]
# partition NEE to get GPP
data["GPP_LT"] = float(-1)*data["NEE_LT"]+data["ER_LT"]
# get the monthly totals
series_list = data.keys()
series_list.sort()
# loop over the months in the data file
start_date = ldt[0]
end_date = start_date+dateutil.relativedelta.relativedelta(months=1)
end_date = end_date-dateutil.relativedelta.relativedelta(minutes=ts)
last_date = ldt[-1]
while start_date<=last_date:
si = qcutils.GetDateIndex(ldt,str(start_date),ts=ts,default=0)
ei = qcutils.GetDateIndex(ldt,str(end_date),ts=ts,default=len(ldt)-1)
td = ldt[si].strftime("%Y-%m-%d")
for item in ["ER_LT","GPP_LT","NEE_LT"]:
if not td in months[item].keys(): months[item][td] = numpy.array([])
data_1d = data[item][si:ei+1]*12.01*ts*60/1E6
months[item][td] = numpy.append(months[item][td],numpy.ma.sum(data_1d))
start_date = end_date+dateutil.relativedelta.relativedelta(minutes=ts)
end_date = start_date+dateutil.relativedelta.relativedelta(months=1)
end_date = end_date-dateutil.relativedelta.relativedelta(minutes=ts)
print datetime.datetime.now(),i
In [10]:
import cPickle as pickle
with open("../Sites/Whroo/Data/months.pickle", 'wb') as handle:
pickle.dump(months, handle, protocol=pickle.HIGHEST_PROTOCOL)
In [ ]:
with open("../Sites/Whroo/Data/months.pickle", 'wb') as handle:
months = pickle.load(handle)
In [37]:
xl_path = "../Sites/Whroo/Data/Portal/"
xl_name = xl_path+"Whroo_L6_filteredL5_Summary.xls"
xl_book = xlrd.open_workbook(xl_name)
datemode = xl_book.datemode
In [38]:
active_sheet = xl_book.sheet_by_name("Monthly")
xldt_l6summary = get_data_xlsheet(active_sheet,2,"Months")
nee_lt_reddyproc = get_data_xlsheet(active_sheet,2,"NEE_REddyProc")
nee_lt_fluxnet = get_data_xlsheet(active_sheet,2,"NEE_LT_FluxNet")
idx = numpy.where(nee_lt_fluxnet==0)[0]
nee_lt_fluxnet[idx] = nee_lt_reddyproc[idx]
nee_ll_fluxnet = get_data_xlsheet(active_sheet,2,"NEE_LL_FluxNet")
idx = numpy.where(nee_ll_fluxnet==0)[0]
nee_ll_fluxnet[idx] = nee_lt_reddyproc[idx]
nee_lt_ozfluxqc = get_data_xlsheet(active_sheet,2,"NEE_LT")
nee_ll_ozfluxqc = get_data_xlsheet(active_sheet,2,"NEE_LL")
nee_solo_ozfluxqc = get_data_xlsheet(active_sheet,2,"NEE_SOLO")
nee_dingo = get_data_xlsheet(active_sheet,2,"NEE_DINGO")
nee_alt_max = numpy.copy(nee_lt_reddyproc)
nee_alt_min = numpy.copy(nee_lt_reddyproc)
for item in [nee_lt_fluxnet,nee_ll_fluxnet,nee_lt_ozfluxqc,nee_ll_ozfluxqc,nee_solo_ozfluxqc,nee_dingo]:
nee_alt_max = numpy.maximum(nee_alt_max,item)
nee_alt_min = numpy.minimum(nee_alt_min,item)
x_alt_nee = range(1,len(nee_alt_max)+1)
In [42]:
gpp_lt_reddyproc = get_data_xlsheet(active_sheet,2,"GPP_REddyProc")
gpp_lt_fluxnet = get_data_xlsheet(active_sheet,2,"GPP_LT_FluxNet")
idx = numpy.where(gpp_lt_fluxnet==0)[0]
gpp_lt_fluxnet[idx] = gpp_lt_reddyproc[idx]
gpp_ll_fluxnet = get_data_xlsheet(active_sheet,2,"GPP_LL_FluxNet")
idx = numpy.where(gpp_ll_fluxnet==0)[0]
gpp_ll_fluxnet[idx] = gpp_lt_reddyproc[idx]
gpp_lt_ozfluxqc = get_data_xlsheet(active_sheet,2,"GPP_LT")
gpp_ll_ozfluxqc = get_data_xlsheet(active_sheet,2,"GPP_LL")
gpp_solo_ozfluxqc = get_data_xlsheet(active_sheet,2,"GPP_SOLO")
gpp_dingo = get_data_xlsheet(active_sheet,2,"GPP_DINGO")
gpp_dingo = float(-1)*gpp_dingo
gpp_alt_max = numpy.copy(gpp_lt_reddyproc)
gpp_alt_min = numpy.copy(gpp_lt_reddyproc)
for item in [gpp_lt_fluxnet,gpp_ll_fluxnet,gpp_lt_ozfluxqc,gpp_ll_ozfluxqc,gpp_solo_ozfluxqc,gpp_dingo]:
gpp_alt_max = numpy.maximum(gpp_alt_max,item)
gpp_alt_min = numpy.minimum(gpp_alt_min,item)
x_alt_gpp = range(1,len(gpp_alt_max)+1)
In [40]:
er_lt_reddyproc = get_data_xlsheet(active_sheet,2,"ER_REddyProc")
er_lt_fluxnet = get_data_xlsheet(active_sheet,2,"ER_LT_FluxNet")
idx = numpy.where(er_lt_fluxnet==0)[0]
er_lt_fluxnet[idx] = er_lt_reddyproc[idx]
er_ll_fluxnet = get_data_xlsheet(active_sheet,2,"ER_LL_FluxNet")
idx = numpy.where(er_ll_fluxnet==0)[0]
er_ll_fluxnet[idx] = er_lt_reddyproc[idx]
er_lt_ozfluxqc = get_data_xlsheet(active_sheet,2,"ER_LT")
er_ll_ozfluxqc = get_data_xlsheet(active_sheet,2,"ER_LL")
er_solo_ozfluxqc = get_data_xlsheet(active_sheet,2,"ER_SOLO")
er_dingo = get_data_xlsheet(active_sheet,2,"ER_DINGO")
er_alt_max = numpy.copy(er_lt_reddyproc)
er_alt_min = numpy.copy(er_lt_reddyproc)
for item in [er_lt_fluxnet,er_ll_fluxnet,er_lt_ozfluxqc,er_ll_ozfluxqc,er_solo_ozfluxqc,er_dingo]:
er_alt_max = numpy.maximum(er_alt_max,item)
er_alt_min = numpy.minimum(er_alt_min,item)
x_alt_er = range(1,len(er_alt_max)+1)
In [20]:
basedate = datetime.datetime(1899, 12, 30)
dt_l6summary = [None]*len(xldt_l6summary)
for i in range(len(xldt_l6summary)):
dt_l6summary[i] = basedate + datetime.timedelta(days=xldt_l6summary[i] + 1462 * datemode)
In [45]:
nee = []; nee_max = []; nee_min = []
er = []; er_max = []; er_min = []
gpp = []; gpp_max = []; gpp_min = []
xlabels_months = []; xlabels_years = []
for item in months["GPP_LT"].keys():
gpp.append(months["GPP_LT"][item])
gpp_max.append(numpy.max(months["GPP_LT"][item]))
gpp_min.append(numpy.min(months["GPP_LT"][item]))
str_parts = item.split("-")
if str_parts[1]=="06": str_parts[1] = str_parts[1]+"\n"+str_parts[0]
xlabels_months.append(str_parts[1])
xlabels_years.append(str_parts[0])
for item in months["ER_LT"].keys():
er.append(months["ER_LT"][item])
er_max.append(numpy.max(months["ER_LT"][item]))
er_min.append(numpy.min(months["ER_LT"][item]))
for item in months["NEE_LT"].keys():
nee.append(months["NEE_LT"][item])
nee_max.append(numpy.max(months["NEE_LT"][item]))
nee_min.append(numpy.min(months["NEE_LT"][item]))
In [103]:
# plot NEE, ER and GPP on a single axis
fig,axs = plt.subplots(1,figsize=(8,4),sharex=True)
axs.set_title("Whroo: Monthly NEE, GPP and ER")
axs.plot([],[],color='red',linewidth=5,label="NEE")
axs.fill_between(x_alt_nee[0:-1],nee_alt_max[0:-1],nee_alt_min[0:-1],facecolor='lightsalmon', alpha=0.5)
axs.fill_between(x,nee_max,nee_min,facecolor='red')
axs.plot([],[],color='green',linewidth=5,label="GPP")
axs.fill_between(x_alt_gpp[0:-1],gpp_alt_max[0:-1],gpp_alt_min[0:-1],facecolor='lightgreen', alpha=0.5)
axs.fill_between(x,gpp_max,gpp_min,facecolor='green')
axs.plot([],[],color='blue',linewidth=5,label="ER")
axs.fill_between(x_alt_er[0:-1],er_alt_max[0:-1],er_alt_min[0:-1],facecolor='lightblue', alpha=0.5)
axs.fill_between(x,er_max,er_min,facecolor='blue')
legend = axs.legend(loc='upper left',ncol=3,frameon=False)
axs.set_ylabel("$gCm^{-2}month^{-1}$",fontsize=16,fontweight='bold')
axs.set_xlim([1,len(x)])
plt.xticks(numpy.arange(1,len(x)+1,1))
xtickNames = plt.setp(axs, xticklabels=xlabels_months)
plt.setp(xtickNames, rotation=0, fontsize=12)
plt.setp(axs.get_xticklabels()[::2], visible=False)
fig.savefig("../Sites/Whroo/Plots/uncertainty_versus_range.svg",format="svg")
plt.show()
In [72]:
nalts = 1
nsites = len(months["NEE_LT"].keys())
minor_ticks = [x for x in range(nalts,nalts*nsites,nalts)]
x = range(1,len(nee_lt_ozfluxqc))
fig,axs = plt.subplots(3,figsize=(8,6),sharex=True)
#fig.canvas.set_window_title(var)
plt.subplots_adjust(hspace=0.15)
axs[0].set_title("Whroo: Monthly NEE, GPP and ER")
axs[0].fill_between(x,nee_max,nee_min)
axs[0].fill_between(x_alt_nee[0:-1],nee_alt_max[0:-1],nee_alt_min[0:-1],facecolor='lightblue', alpha=0.5)
axs[0].set_ylabel("NEE")
axs[0].set_xlim([1,44])
axs[1].fill_between(x,gpp_max,gpp_min)
axs[1].fill_between(x_alt_gpp[0:-1],gpp_alt_max[0:-1],gpp_alt_min[0:-1],facecolor='lightblue', alpha=0.5)
axs[1].set_ylabel("GPP")
axs[2].fill_between(x,er_max,er_min)
axs[2].fill_between(x_alt_er[0:-1],er_alt_max[0:-1],er_alt_min[0:-1],facecolor='lightblue', alpha=0.5)
axs[2].set_ylabel("ER")
xticklabels = axs[0].get_xticklabels()+axs[1].get_xticklabels()
plt.setp(xticklabels, visible=False)
plt.xticks(numpy.arange(1,45,1))
xtickNames = plt.setp(axs[2], xticklabels=xlabels_months)
plt.setp(xtickNames, rotation=0, fontsize=10)
plt.setp(axs[2].get_xticklabels()[::2], visible=False)
plt.show()
In [84]:
fig, ax1 = plt.subplots()
ax1.fill_between(x,nee_max,nee_min)
ax1.fill_between(x,er_max,er_min)
ax1.set_ylim(-100,250)
ax2 = ax1.twinx()
ax2.fill_between(x,gpp_max,gpp_min)
ax2.set_ylim(-150,200)
plt.show()
In [85]:
fig,(ax1,ax2) = plt.subplots(2, 1, sharex=True)
ax1.fill_between(x,gpp_max,gpp_min)
ax2.fill_between(x,nee_max,nee_min)
ax2.fill_between(x,er_max,er_min)
ax1.set_ylim(0,200)
ax2.set_ylim(-100,150)
ax1.spines['bottom'].set_visible(False)
ax2.spines['top'].set_visible(False)
ax1.xaxis.tick_top()
ax1.tick_params(labeltop='off') # don't put tick labels at the top
ax2.xaxis.tick_bottom()
plt.show()
In [ ]:
fig,axs = plt.subplots(3,figsize=(12,6),sharex=True)
axs[0].plot(ldt,data["NEE_LT"],'bo')
axs[0].set_ylabel("NEE")
axs[1].plot(ldt,data["GPP_LT"],'bo')
axs[1].set_ylabel("GPP")
axs[2].plot(ldt,data["ER_LT"],'bo')
axs[2].set_ylabel("ER")
plt.tight_layout()
plt.show()
In [ ]: