In [1]:
%run basics
%matplotlib
In [2]:
def Fre_LloydTaylor(T,rb,E0):
"""
Lloyd-Taylor expression for respiration as a function of temperature.
E0 is the activation energy (sensitivity of respiration to temperature)
rb is the baseline respiration
"""
# return rb*numpy.exp(E0*(1/(c.Tref-c.T0)-1/(T-c.T0)))
return rb*numpy.exp(E0*(1/(288.13-227.13)-1/((T+273.15)-227.13)))
In [3]:
def plot_Reco(fig,i,Reco,T,Sws,weight_type='Huber'):
ax = fig.add_subplot(2,2,i)
popt,pcov = irls_leastsq(residuals_LT,[2,200],args=(Reco,T),maxfev=3,weight_type=weight_type)
resid = Reco - Reco_LloydTaylor(T,popt[0],popt[1])
weights = CalculateWeights(resid,weight_type=weight_type)
plt.scatter(T,Reco,c=weights)
#mean,edges,number = scipy.stats.binned_statistic(T,Reco,statistic='mean')
#mid = [(edges[i]+edges[i+1])/2 for i in range(0,len(edges)-1)]
#plt.plot(mid,mean,'ro')
Xfit = numpy.linspace(numpy.min(T),numpy.max(T),100)
Yfit = Reco_LloydTaylor(Xfit,popt[0],popt[1])
plt.plot(Xfit,Yfit,'r-')
text_str = ('%.2f < Sws < %.2f'%(Sws[0],Sws[-1]))
plt.text(0.5,0.9,text_str,horizontalalignment='center',transform=ax.transAxes)
text_str = ('rb=%.2f,E0=%.2f'%(popt[0],popt[-1]))
plt.text(0.5,0.8,text_str,horizontalalignment='center',transform=ax.transAxes)
In [4]:
def irls_leastsq(func,p0,args=(),maxfev=3,weight_type='Huber'):
weighted = False
popt,pcov = scipy.optimize.leastsq(func,p0,args=args+(weighted,weight_type))
print 'After non-weighted call to leastsq: ',popt
n=1
weighted = True
while n<=maxfev:
popt,pcov = scipy.optimize.leastsq(func,popt,args=args+(weighted,weight_type))
print 'Weighted call '+str(n)+' to leastsq: ',popt
n = n + 1
return popt,pcov
In [5]:
def residuals_LloydTaylor(params,Fre,T,weighted,weight_type):
rb = params[0]
E0 = params[1]
resid = Fre - Fre_LloydTaylor(T,rb,E0)
if weighted:
weights = CalculateWeights(resid,weight_type=weight_type)
return resid*weights
else:
return resid
In [6]:
def CalculateWeights(resid,weight_type='Huber'):
if weight_type.lower()=='tukey':
weights = TukeyBiweight(MAD_scale(resid))
elif weight_type.lower()=='huber':
weights = HuberTWeight(MAD_scale(resid))
elif weight_type.lower()=='ols':
weights = OLSWeight(MAD_scale(resid))
else:
print "CalculateWeights: Unknown weight type, used Huber's t weight"
weights = HuberTWeight(MAD_scale(resid))
return weights
In [7]:
def TukeyBiweight(x):
return ((numpy.abs(x)<1.0)*(1.0 - x**2)**2)
In [8]:
def HuberTWeight(x):
weights = numpy.ones(len(x))
idx = numpy.where(numpy.abs(x)>1.345)[0]
weights[idx] = 1.345/numpy.abs(x[idx])
return weights
In [9]:
def OLSWeight(x):
weights = numpy.ones(len(x))
return weights
In [10]:
def MAD_scale(resid):
scale = numpy.median(numpy.abs(resid - numpy.median(resid))) / 0.6745
return resid/scale
In [11]:
ncname = qcio.get_filename_dialog(path="../../Sites/HowardSprings/Data/Processed/all")
print ncname
ds = qcio.nc_read_series(ncname)
ldt = ds.series["DateTime"]["Data"]
ts = int(ds.globalattributes["time_step"])
site_name = ds.globalattributes["site_name"]
In [12]:
Fsd,Fsd_flag,Fsd_attr = qcutils.GetSeriesasMA(ds,"Fsd")
Fc,Fc_flag,Fc_attr = qcutils.GetSeriesasMA(ds,"Fc")
T,T_flag,T_attr = qcutils.GetSeriesasMA(ds,"Ts")
ustar,ustar_flag,ustar_attr = qcutils.GetSeriesasMA(ds,"ustar")
In [13]:
fig = plt.figure(1)
ax1=plt.subplot(411)
plt.plot(ldt,Fsd)
ax2=plt.subplot(412,sharex=ax1)
plt.plot(ldt,Fc)
ax3=plt.subplot(413,sharex=ax1)
plt.plot(ldt,ustar)
ax4=plt.subplot(414,sharex=ax1)
plt.plot(ldt,T)
plt.show()
In [16]:
# only use Fc when Fc QC flag=0
Fc_best = numpy.ma.masked_where(Fc_flag!=0,Fc)
Fc_best = numpy.ma.masked_where(Fc_best<0,Fc_best)
# mask day time values
Fc_nocturnal = numpy.ma.masked_where(Fsd>=10,Fc_best)
# mask low ustar values
# NOTE: ustar threshold value to be read from control file
Fre = numpy.ma.masked_where(ustar<0.25,Fc_nocturnal)
# synchronise masks across drivers and targets
mask = Fre.mask
for item in [T]:
mask = numpy.ma.mask_or(mask,item.mask)
# apply synchronised mask across drivers and targets
Fre = numpy.ma.array(Fre,mask=mask)
T = numpy.ma.array(T,mask=mask)
In [17]:
fig = plt.figure(1)
ax1=plt.subplot(211)
plt.plot(ldt,Fre)
ax2=plt.subplot(212,sharex=ax1)
plt.plot(ldt,T)
plt.show()
In [35]:
# get the yearly E0 values
first_date = ldt[0]
last_date = ldt[-1]
years = [yr for yr in range(first_date.year,last_date.year+1)]
for year in [yr for yr in range(first_date.year,last_date.year+1)]:
idx = numpy.where(ds.series["Year"]["Data"]==year)[0]
Fre_year = numpy.ma.compressed(Fre[idx])
T_year = numpy.ma.compressed(T[idx])
popt,pcov = irls_leastsq(residuals_LloydTaylor,[1,100],args=(Fre_year,T_year),maxfev=3,weight_type='Huber')
Fre_fit = Fre_LloydTaylor(T_year,popt[0],popt[1])
fig = plt.figure()
plt.plot(T_year,Fre_year,'b.')
plt.plot(T_year,Fre_fit,'r-')
plt.show()
In [36]:
def update_end_date(end_date,window_units,window_length,window_overlap):
if window_units.lower()=="days":
end_date = start_date + dateutil.relativedelta.relativedelta(days=window_length)
elif window_units.lower()=="months":
end_date = start_date + dateutil.relativedelta.relativedelta(months=window_length)
elif window_units.lower()=="years":
end_date = start_date + dateutil.relativedelta.relativedelta(years=window_length)
else:
print "Unrecognised option for window length units"
return end_date
def update_start_date(start_date,window_units,window_length,window_overlap):
window_jump = window_length - window_overlap
if window_units.lower()=="days":
start_date = start_date + dateutil.relativedelta.relativedelta(days=window_jump)
elif window_units.lower()=="months":
start_date = start_date + dateutil.relativedelta.relativedelta(months=window_jump)
elif window_units.lower()=="years":
start_date = start_date + dateutil.relativedelta.relativedelta(years=window_jump)
else:
print "Unrecognised option for window length units"
return start_date
In [45]:
# Get the long term E0
# set the window units, length and overlap
# these will be read from the control file
window_units = "years"
window_length = 1
window_overlap = 0
# get the first and last dates in the data set
first_date = ldt[0]
last_date = ldt[-1]
print first_date,last_date
# get the initial start date
start_date = ldt[0]
# get the initial end date
end_date = update_end_date(start_date,window_units,window_length,window_overlap)
# create empty lists to hold the start and end dates of the data window
start_dates = []
end_dates = []
# loop until the last date to create the lists of start and end dates
while start_date<last_date:
start_dates.append(start_date)
end_dates.append(end_date)
start_date = update_start_date(start_date,window_units,window_length,window_overlap)
end_date = update_end_date(end_date,window_units,window_length,window_overlap)
for sd,ed in zip(start_dates,end_dates):
si = qcutils.GetDateIndex(ldt,str(sd),ts=ts,default=0)
ei = qcutils.GetDateIndex(ldt,str(ed),ts=ts,default=len(ldt))
print sd,si,ed,ei
Fre_year = numpy.ma.compressed(Fre[si:ei+1])
T_year = numpy.ma.compressed(T[si:ei+1])
popt,pcov = irls_leastsq(residuals_LloydTaylor,[1,100],args=(Fre_year,T_year),maxfev=3,weight_type='Huber')
Fre_fit = Fre_LloydTaylor(T_year,popt[0],popt[1])
fig = plt.figure()
plt.plot(T_year,Fre_year,'b.')
plt.plot(T_year,Fre_fit,'r-')
plt.show()
# one year
#end_date = datetime.datetime(start_date.year+1,1,1)+datetime.timedelta(minutes=ts)
#end_date = datetime.datetime(start_date.year+1,1,1)+datetime.timedelta(minutes=ts)
#this_date = start_date + dateutil.relativedelta.relativedelta(months=1)
#while this_date<end_date:
#si = qcutils.GetDateIndex(ldt,str(start_date),ts=ts)
#ei = qcutils.GetDateIndex(ldt,str(this_date),ts=ts)
#ldt_plt = ldt[si:ei+1]
#Fre_plt = Fc_nocturnal_filtered[si:ei+1]
#T_plt = T[si:ei+1]
#idx = numpy.ma.where(Fre_plt.mask==False)[0]
#ldt_plt = [ldt_plt[i] for i in idx]
#Fre_plt = numpy.ma.compressed(Fre_plt)
#T_plt = numpy.ma.compressed(T_plt)
In [27]:
print [year for year in range(first_date.year,last_date.year+1)]
In [17]:
fig=plt.figure(1)
ax1=plt.subplot(211)
plt.plot(ldt_plt,Fre_plt,'b.')
ax2=plt.subplot(212,sharex=ax1)
plt.plot(ldt_plt,T_plt,'b.')
plt.show()
In [ ]:
fig=plt.figure(1)
plt.plot(T_plt,Fre_plt,'b.')
plt.show()
In [ ]:
popt,pcov = irls_leastsq(residuals_Isaac,[0.1,100,1,1],args=(NEE_plt,Fsd_plt,D_plt,T_plt,0.005,50),maxfev=3,weight_type='Huber')
print start_date,this_date,popt
start_date = start_date + datetime.timedelta(days=5)
this_date = start_date + datetime.timedelta(days=10)