In [1]:
%run basics
%matplotlib
In [2]:
def NEE_Lasslop(Fsd,D,T,D0,E0,alpha,beta0,k,rb):
NEE = -1*GPP_Lasslop(Fsd,D,alpha,beta0,k,D0) + Reco_LloydTaylor(T,rb,E0)
return NEE
In [3]:
def NEE_Isaac(Fsd,D,T,D0,E0,alpha,beta0,rb):
#def NEE_Isaac(args,params):
NEE = -1*GPP_Isaac(Fsd,D,alpha,beta0,D0) + Reco_LloydTaylor(T,rb,E0)
return NEE
In [4]:
def NEE_NRH_Isaac(Fsd,D,T,D0,E0,alpha,beta0,eta,rb):
NEE = -1*GPP_NRH_Isaac(Fsd,D,alpha,beta0,eta,D0) + Reco_LloydTaylor(T,rb,E0)
return NEE
In [5]:
def GPP_Lasslop(Fsd,D,alpha,beta0,k,D0):
beta = beta0*SHD_func_Lasslop(D,k,D0)
GPP = alpha*beta*Fsd/(alpha*Fsd+beta)
return GPP
In [6]:
def GPP_Isaac(Fsd,D,alpha,beta0,D0):
beta = beta0*SHD_func_Isaac(D,D0)
GPP = alpha*beta*Fsd/(alpha*Fsd+beta)
return GPP
In [7]:
def GPP_NRH_Isaac(Fsd,D,alpha,beta0,eta,D0):
beta = beta0*SHD_func_Isaac(D,D0)
GPP = (1/(2*eta))*(alpha*Fsd+beta-numpy.sqrt((alpha*Fsd+beta)**2-4*alpha*beta*eta*Fsd))
return GPP
In [8]:
def SHD_func_Lasslop(D,k,D0):
SHD_func = numpy.ones(len(D))
idx = numpy.where(D>D0)[0]
SHD_func[idx] = numpy.exp(-k*(D[idx]-D0))
return SHD_func
In [9]:
def SHD_func_Isaac(D,D0):
SHD_func = numpy.ones(len(D))
idx = numpy.where(D>D0)[0]
SHD_func[idx] = 1/(1+(D[idx]-D0)/D0)
return SHD_func
In [10]:
def Reco_LloydTaylor(T,rb,E0):
# 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 [11]:
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 [12]:
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
NEE_plt = NEE_Isaac(args[1],args[2],args[3],args[4],args[5],popt[0],popt[1],popt[2])
# NEE_plt = NEE_Isaac(args[1:],popt)
resid = args[0] - NEE_plt
weights = CalculateWeights(resid,weight_type=weight_type)
fig = plt.figure()
plt.scatter(args[1],args[0],c=weights)
plt.plot(args[1],NEE_plt,'r+')
plt.show()
n = n + 1
return popt,pcov
In [13]:
def residuals_LloydTaylor(params,Reco,T,weighted,weight_type):
rb = params[0]
E0 = params[1]
resid = Reco - Reco_LloydTaylor(T,rb,E0)
if weighted:
weights = CalculateWeights(resid,weight_type=weight_type)
return resid*weights
else:
return resid
In [14]:
def residuals_Lasslop(params,NEE,Fsd,D,T,D0,E0,weighted,weight_type):
alpha = params[0]
beta0 = params[1]
k = params[2]
rb = params[3]
resid = NEE - NEE_Lasslop(Fsd,D,T,D0,E0,alpha,beta0,k,rb)
if weighted:
weights = CalculateWeights(resid,weight_type=weight_type)
return resid*weights
else:
return resid
In [15]:
def residuals_Isaac(params,NEE,Fsd,D,T,D0,E0,weighted,weight_type):
alpha = params[0]
beta0 = params[1]
eta = params[2]
rb = params[3]
resid = NEE - NEE_NRH_Isaac(Fsd,D,T,D0,E0,alpha,beta0,eta,rb)
if weighted:
weights = CalculateWeights(resid,weight_type=weight_type)
return resid*weights
else:
return resid
In [16]:
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 [17]:
def TukeyBiweight(x):
return ((numpy.abs(x)<1.0)*(1.0 - x**2)**2)
In [18]:
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 [19]:
def OLSWeight(x):
weights = numpy.ones(len(x))
return weights
In [20]:
def MAD_scale(resid):
scale = numpy.median(numpy.abs(resid - numpy.median(resid))) / 0.6745
return resid/scale
In [21]:
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 [22]:
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")
In [23]:
fig = plt.figure(1)
ax1=plt.subplot(411)
plt.plot(ldt,Fsd)
ax2=plt.subplot(412,sharex=ax1)
plt.plot(ldt,D)
ax3=plt.subplot(413,sharex=ax1)
plt.plot(ldt,T)
ax4=plt.subplot(414,sharex=ax1)
plt.plot(ldt,Fc)
plt.show()
In [24]:
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)
NEE_day = 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 [25]:
fig = plt.figure(1)
ax1=plt.subplot(411)
plt.plot(ldt,Fsd)
ax2=plt.subplot(412,sharex=ax1)
plt.plot(ldt,D)
ax3=plt.subplot(413,sharex=ax1)
plt.plot(ldt,T)
ax4=plt.subplot(414,sharex=ax1)
plt.plot(ldt,NEE_day)
plt.show()
In [26]:
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)
ldt_plt = ldt[si:ei+1]
NEE_plt = NEE_day[si:ei+1]
Fsd_plt = Fsd[si:ei+1]
D_plt = D[si:ei+1]
T_plt = T[si:ei+1]
idx = numpy.ma.where(NEE_plt.mask==False)[0]
ldt_plt = [ldt_plt[i] for i in idx]
NEE_plt = numpy.ma.compressed(NEE_plt)
Fsd_plt = numpy.ma.compressed(Fsd_plt)
D_plt = numpy.ma.compressed(D_plt)
T_plt = numpy.ma.compressed(T_plt)
In [27]:
fig=plt.figure(1)
ax1=plt.subplot(411)
plt.plot(ldt_plt,Fsd_plt,'b.')
ax2=plt.subplot(412,sharex=ax1)
plt.plot(ldt_plt,D_plt,'b.')
ax3=plt.subplot(413,sharex=ax1)
plt.plot(ldt_plt,T_plt,'b.')
ax4=plt.subplot(414,sharex=ax1)
plt.plot(ldt_plt,NEE_plt,'b.')
plt.show()
In [26]:
fig=plt.figure(1)
plt.plot(Fsd_plt,NEE_plt,'b.')
plt.show()
In [31]:
#popt,pcov = irls_leastsq(residuals_Lasslop,[0.1,100,0,1],args=(NEE_plt,Fsd_plt,D_plt,T_plt,0.005,210),maxfev=3,weight_type='Huber')
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)
In [44]:
fig=plt.figure()
plt.plot(D_plt,SHD_func_Isaac(D_plt,0.005),'b.')
plt.show()
In [ ]: