In [1]:
%run basics
%matplotlib
In [2]:
def NEE_Lasslop(Fsd,D,T,alpha,beta0,k,rb,D0,E0):
NEE = -1*GPP_Lasslop(Fsd,D,alpha,beta0,k,D0) + Reco_LloydTaylor(T,rb,E0)
return NEE
In [3]:
def GPP_Lasslop(Fsd,D,alpha,beta0,k,D0):
beta = beta0*VPD_Lasslop(D,k,D0)
GPP = alpha*beta*Fsd/(alpha*Fsd+beta)
return GPP
In [4]:
def VPD_Lasslop(D,k,D0):
VPD = numpy.ones(len(D))
idx = numpy.where(D>D0)[0]
VPD[idx] = numpy.exp(-k*(D[idx]-D0))
return VPD
In [5]:
def Reco_LloydTaylor(T,rb,E0):
return rb*numpy.exp(E0*(1/(c.Tref-c.T0)-1/(T-c.T0)))
In [6]:
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 [7]:
def irls_leastsq(func,p0,args=(),maxfev=3,weight_type='Huber'):
weighted = False
p,cov = scipy.optimize.leastsq(func,p0,args=args+(weighted,weight_type))
print 'After non-weighted call to leastsq: ',p[0],p[1]
n=1
weighted = True
while n<=maxfev:
p,cov = scipy.optimize.leastsq(func,[p[0],p[1]],args=args+(weighted,weight_type))
print 'Weighted call '+str(n)+' to leastsq: ',p[0],p[1]
n = n + 1
return p,cov
In [8]:
def residuals_LT(p,Reco,T,weighted,weight_type):
rb = p[0]
E0 = p[1]
resid = Reco - Reco_LloydTaylor(T,rb,E0)
if weighted:
weights = CalculateWeights(resid,weight_type=weight_type)
return resid*weights
else:
return resid
In [9]:
def residuals_LL(p,NEE,Fsd,D,T,D0,E0,weighted,weight_type):
alpha = p[0]
beta0 = p[1]
k = p[2]
rb = p[3]
resid = NEE - NEE_Lasslop(Fsd,D,T,alpha,beta0,k,rb,D0,E0)
if weighted:
weights = CalculateWeights(resid,weight_type=weight_type)
return resid*weights
else:
return resid
In [10]:
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 [11]:
def TukeyBiweight(x):
return ((numpy.abs(x)<1.0)*(1.0 - x**2)**2)
In [12]:
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 [13]:
def OLSWeight(x):
weights = numpy.ones(len(x))
return weights
In [14]:
def MAD_scale(resid):
scale = numpy.median(numpy.abs(resid - numpy.median(resid))) / 0.6745
return resid/scale
In [38]:
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 [39]:
Fsd,f,a = qcutils.GetSeriesasMA(ds,"Fsd")
D,f,a = qcutils.GetSeriesasMA(ds,"SHD")
T,f,a = qcutils.GetSeriesasMA(ds,"Ts")
Fc,f,a = qcutils.GetSeriesasMA(ds,"Fc")
fig1,axs = plt.subplots(nrows=4,ncols=1,sharex=True)
axs[0].plot(Fsd)
axs[1].plot(D)
axs[2].plot(T)
axs[3].plot(Fc)
plt.show
Fc = numpy.ma.masked_where(f!=0,Fc)
Fc_day = numpy.ma.masked_where(Fsd<=10,Fc)
mask = Fc_day.mask
for item in [D,T]:
mask = numpy.ma.mask_or(mask,item.mask)
NEP_day = -1*numpy.ma.array(Fc_day,mask=mask)
Fsd = numpy.ma.array(Fsd,mask=mask)
D = numpy.ma.array(D,mask=mask)
T = numpy.ma.array(T,mask=mask)
In [30]:
plt.plot(T)
plt.show()
In [33]:
start_date = ldt[0]
end_date = datetime.datetime(start_date.year+1,1,1)+datetime.timedelta(minutes=ts)
this_date = start_date + datetime.timedelta(days=10)
while this_date<end_date:
si = qcutils.GetDateIndex(ldt,str(start_date),ts=ts)
ei = qcutils.GetDateIndex(ldt,str(this_date),ts=ts)
NEP_plt = numpy.ma.compressed(NEP_day[si:ei+1])
Fsd_plt = numpy.ma.compressed(Fsd[si:ei+1])
D_plt = numpy.ma.compressed(D[si:ei+1])
T_plt = numpy.ma.compressed(T[si:ei+1])
popt,pcov = irls_leastsq(residuals_LL,[0.1,100,0,2],args=(NEP_plt,Fsd_plt,D_plt,T_plt,0.01,210),maxfev=3,weight_type='Huber')
print start_date,this_date,popt[0],popt[1]
start_date = start_date + datetime.timedelta(days=5)
this_date = start_date + datetime.timedelta(days=10)
In [19]:
print len(NEP_day)
In [26]:
plt.plot(Fsd)
plt.show()
In [ ]: