In [1]:
%run basics
%matplotlib
In [6]:
name_list = ["/mnt/OzFlux/Sites/Calperum/Data/Portal/Calperum_L3.nc"]
In [2]:
name_list = ["/mnt/OzFlux/Sites/AdelaideRiver/Data/Portal/AdelaideRiver_L3.nc",
"/mnt/OzFlux/Sites/AliceSpringsMulga/Data/Portal/AliceSpringsMulga_L3.nc",
"/mnt/OzFlux/Sites/Calperum/Data/Portal/Calperum_L3.nc",
"/mnt/OzFlux/Sites/CapeTribulation/Data/Portal/CapeTribulation_L3.nc",
"/mnt/OzFlux/Sites/CowBay/Data/Portal/CowBay_L3.nc",
"/mnt/OzFlux/Sites/CumberlandPlains/Data/Portal/CumberlandPlains_L3.nc",
"/mnt/OzFlux/Sites/DalyPasture/Data/Portal/DalyPasture_L3.nc",
"/mnt/OzFlux/Sites/DalyUncleared/Data/Portal/DalyUncleared_L3.nc",
"/mnt/OzFlux/Sites/DryRiver/Data/Portal/DryRiver_L3.nc",
"/mnt/OzFlux/Sites/Emerald/Data/Portal/Emerald_L3.nc",
"/mnt/OzFlux/Sites/FoggDam/Data/Portal/FoggDam_L3.nc",
"/mnt/OzFlux/Sites/Gingin/Data/Portal/Gingin_L3.nc",
"/mnt/OzFlux/Sites/GreatWesternWoodlands/Data/Portal/GreatWesternWoodlands_L3.nc",
"/mnt/OzFlux/Sites/HowardSprings/Data/Portal/HowardSprings_L3.nc",
"/mnt/OzFlux/Sites/Otway/Data/Portal/Otway_L3.nc",
"/mnt/OzFlux/Sites/RedDirtMelonFarm/Data/Portal/RedDirtMelonFarm_L3.nc",
"/mnt/OzFlux/Sites/RiggsCreek/Data/Portal/RiggsCreek_L3.nc",
"/mnt/OzFlux/Sites/RobsonCreek/Data/Portal/RobsonCreek_L3.nc",
"/mnt/OzFlux/Sites/Samford/Data/Portal/Samford_L3.nc",
"/mnt/OzFlux/Sites/SturtPlains/Data/Portal/SturtPlains_L3.nc",
"/mnt/OzFlux/Sites/Tumbarumba/Data/Portal/Tumbarumba_L3.nc",
"/mnt/OzFlux/Sites/Whroo/Data/Portal/Whroo_L3.nc",
"/mnt/OzFlux/Sites/WombatStateForest/Data/Portal/WombatStateForest_L3.nc",
"/mnt/OzFlux/Sites/Yanco/Data/Portal/Yanco_L3.nc"]
In [3]:
rad = ['Fsd','Fsu','Fld','Flu','Fn']
met = ['Ah','Cc','Precip','ps','Ta','Ws','Wd']
soil = ['Fg','Ts','Sws']
flux = ['Fm','ustar','Fh','Fe','Fc']
rad_gap_len = numpy.array([])
met_gap_len = numpy.array([])
soil_gap_len = numpy.array([])
flux_gap_len = numpy.array([])
for nc_name in name_list:
print "Processing ",nc_name
ds = qcio.nc_read_series(nc_name)
ts = int(ds.globalattributes["time_step"])
# number of time steps per hour
n = 60/ts
bins = [0,2*n+1,24*n+1,5*24*n+1,30*24*n+1,365*24*n]
bin_labels = ["<2hrs","<1day","<5days","<30days",">30days","365 days"]
# radiation
for item in rad:
data,flag,attr = qcutils.GetSeriesasMA(ds,item)
mask = numpy.ma.getmaskarray(data)
gap_startend = qcutils.contiguous_regions(mask)
for i in range(len(gap_startend)):
length = gap_startend[i][1]-gap_startend[i][0]
rad_gap_len = numpy.append(rad_gap_len,length)
# meteorology
for item in met:
data,flag,attr = qcutils.GetSeriesasMA(ds,item)
mask = numpy.ma.getmaskarray(data)
gap_startend = qcutils.contiguous_regions(mask)
for i in range(len(gap_startend)):
length = gap_startend[i][1]-gap_startend[i][0]
met_gap_len = numpy.append(met_gap_len,length)
# soil
for item in soil:
data,flag,attr = qcutils.GetSeriesasMA(ds,item)
mask = numpy.ma.getmaskarray(data)
gap_startend = qcutils.contiguous_regions(mask)
for i in range(len(gap_startend)):
length = gap_startend[i][1]-gap_startend[i][0]
soil_gap_len = numpy.append(soil_gap_len,length)
# fluxes
for item in flux:
data,flag,attr = qcutils.GetSeriesasMA(ds,item)
mask = numpy.ma.getmaskarray(data)
gap_startend = qcutils.contiguous_regions(mask)
for i in range(len(gap_startend)):
length = gap_startend[i][1]-gap_startend[i][0]
flux_gap_len = numpy.append(flux_gap_len,length)
In [7]:
rad_hist,rad_edges = numpy.histogram(rad_gap_len,bins=bins)
rad_hist = rad_hist.astype(float)/float(numpy.sum(rad_hist))
met_hist,met_edges = numpy.histogram(met_gap_len,bins=bins)
met_hist = met_hist.astype(float)/float(numpy.sum(met_hist))
soil_hist,soil_edges = numpy.histogram(soil_gap_len,bins=bins)
soil_hist = soil_hist.astype(float)/float(numpy.sum(soil_hist))
flux_hist,flux_edges = numpy.histogram(flux_gap_len,bins=bins)
flux_hist = flux_hist.astype(float)/float(numpy.sum(flux_hist))
bar_width = 0.35
fig,axs = plt.subplots(2,2,figsize=(6,6))
index = numpy.arange(len(rad_hist))
axs[0,0].bar(index,rad_hist)
axs[0,0].set_title("Radiation")
axs[0,0].set_ylabel("Frequency",fontsize=16)
#axs[0,0].set_ylabel("Occurrences",fontsize=16)
axs[0,0].set_xticks(index+bar_width)
axs[0,0].set_xticklabels(bin_labels)
axs[0,1].bar(index,met_hist)
axs[0,1].set_title("Meteorology")
axs[0,1].set_xticks(index+bar_width)
axs[0,1].set_xticklabels(bin_labels)
axs[1,0].bar(index,soil_hist)
axs[1,0].set_title("Soil")
axs[1,0].set_xlabel("Gap length",fontsize=16)
axs[1,0].set_ylabel("Frequency",fontsize=16)
#axs[1,0].set_ylabel("Occurrences",fontsize=16)
axs[1,0].set_xticks(index+bar_width)
axs[1,0].set_xticklabels(bin_labels)
axs[1,1].bar(index,flux_hist)
axs[1,1].set_title("Fluxes")
axs[1,1].set_xlabel("Gap length",fontsize=16)
axs[1,1].set_xticks(index+bar_width)
axs[1,1].set_xticklabels(bin_labels)
plt.tight_layout()
plt.show()
In [44]:
from scipy import stats
rad_sum,rad_edges,rad_num = stats.binned_statistic(rad_gap_len,rad_gap_len,bins=bins,statistic="sum")
rad_sum = rad_sum.astype(float)/float(numpy.sum(rad_sum))
met_sum,met_edges,met_num = stats.binned_statistic(met_gap_len,met_gap_len,bins=bins,statistic="sum")
met_sum = met_sum.astype(float)/float(numpy.sum(met_sum))
soil_sum,soil_edges,soil_num = stats.binned_statistic(soil_gap_len,soil_gap_len,bins=bins,statistic="sum")
soil_sum = soil_sum.astype(float)/float(numpy.sum(soil_sum))
flux_sum,flux_edges,flux_num = stats.binned_statistic(flux_gap_len,flux_gap_len,bins=bins,statistic="sum")
flux_sum = flux_sum.astype(float)/float(numpy.sum(flux_sum))
bar_width = 0.35
fig,axs = plt.subplots(2,2,figsize=(10,10))
index = numpy.arange(len(rad_hist[0]))
axs[0,0].bar(index,rad_sum)
axs[0,0].set_title("Radiation")
#axs[0,0].set_ylabel("Frequency",fontsize=16)
axs[0,0].set_ylabel("Fraction of all gaps",fontsize=16)
axs[0,0].set_xticks(index+bar_width)
axs[0,0].set_xticklabels(bin_labels)
axs[0,1].bar(index,met_sum)
axs[0,1].set_title("Meteorology")
axs[0,1].set_xticks(index+bar_width)
axs[0,1].set_xticklabels(bin_labels)
axs[1,0].bar(index,soil_sum)
axs[1,0].set_title("Soil")
axs[1,0].set_xlabel("Gap length",fontsize=16)
#axs[1,0].set_ylabel("Frequency",fontsize=16)
axs[1,0].set_ylabel("Fraction of all gaps",fontsize=16)
axs[1,0].set_xticks(index+bar_width)
axs[1,0].set_xticklabels(bin_labels)
axs[1,1].bar(index,flux_sum)
axs[1,1].set_title("Fluxes")
axs[1,1].set_xlabel("Gap length",fontsize=16)
axs[1,1].set_xticks(index+bar_width)
axs[1,1].set_xticklabels(bin_labels)
plt.tight_layout()
plt.show()
In [ ]: