In [1]:
%run basics
%matplotlib
In [2]:
def residuals_LloydTaylor(params,ER,T,weighted,weight_type):
rb = params[0]
E0 = params[1]
resid = ER - ER_LloydTaylor(T,rb,E0)
if weighted:
weights = CalculateWeights(resid,weight_type=weight_type)
return resid*weights
else:
return resid
def ER_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 [3]:
def residuals_RHLRC_D(params,NEE,Fsd,D,T,D0,E0,weighted,weight_type):
alpha = params[0]
beta = params[1]
k = params[2]
rb = params[3]
resid = NEE - NEE_RHLRC_D(Fsd,D,T,alpha,beta,k,D0,rb,E0)
if weighted:
weights = CalculateWeights(resid,weight_type=weight_type)
return resid*weights
else:
return resid
def NEE_RHLRC_D(Fsd,D,T,alpha,beta,k,D0,rb,E0):
NEE = -1*GPP_RHLRC_D(Fsd,D,alpha,beta,k,D0) + ER_LloydTaylor(T,rb,E0)
return NEE
def GPP_RHLRC_D(Fsd,D,alpha,beta,k,D0):
beta = beta*SHD_func_Lasslop(D,k,D0)
GPP = alpha*beta*Fsd/(alpha*Fsd+beta)
return GPP
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 [4]:
def irls_leastsq(func,p0,args=(),maxfev=3,weight_type='Huber',mode='verbose'):
weighted = False
popt,pcov = scipy.optimize.leastsq(func,p0,args=args+(weighted,weight_type))
if mode=='verbose':
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))
if mode=='verbose':
print 'Weighted call '+str(n)+' to leastsq: ',popt
ER_fit = ER_LloydTaylor(args[1],popt[0],popt[1])
resid = args[0] - ER_fit
weights = CalculateWeights(resid,weight_type=weight_type)
n = n + 1
return popt,pcov
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
def TukeyBiweight(x):
return ((numpy.abs(x)<1.0)*(1.0 - x**2)**2)
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
def OLSWeight(x):
weights = numpy.ones(len(x))
return weights
def MAD_scale(resid):
scale = numpy.median(numpy.abs(resid - numpy.median(resid))) / 0.6745
return resid/scale
In [5]:
ncname ="../Sites/Whroo/Data/Portal/Whroo_L3.nc"
ds = qcio.nc_read_series(ncname)
ldt = ds.series["DateTime"]["Data"]
ts = int(ds.globalattributes["time_step"])
site_name = ds.globalattributes["site_name"]
nrecs = int(ds.globalattributes["nc_nrecs"])
In [6]:
# get the data and synchronise the gaps
indicator = numpy.ones(nrecs,dtype=numpy.int)
Fsd,f,a = qcutils.GetSeriesasMA(ds,"Fsd")
idx = numpy.where(f!=0)[0]
indicator[idx] = numpy.int(0)
D,f,a = qcutils.GetSeriesasMA(ds,"VPD")
idx = numpy.where(f!=0)[0]
indicator[idx] = numpy.int(0)
T,f,a = qcutils.GetSeriesasMA(ds,"Ta")
idx = numpy.where(f!=0)[0]
indicator[idx] = numpy.int(0)
ustar,f,a = qcutils.GetSeriesasMA(ds,"ustar")
idx = numpy.where(f!=0)[0]
indicator[idx] = numpy.int(0)
Fc,f,a = qcutils.GetSeriesasMA(ds,"Fc")
idx = numpy.where(f!=0)[0]
indicator[idx] = numpy.int(0)
In [7]:
indicator_night = numpy.copy(indicator)
# apply a day/night filter
idx = numpy.where(Fsd>=10)[0]
indicator_night[idx] = numpy.int(0)
# apply a simple ustar filter
idx = numpy.where(ustar<=0.4)[0]
indicator_night[idx] = numpy.int(0)
In [8]:
# synchronise the gaps and apply the ustar filter
T_night = numpy.ma.masked_where(indicator_night==0,T)
ER = numpy.ma.masked_where(indicator_night==0,Fc)
In [9]:
# loop over the windows and get E0
LT_results = {"start_date":[],"mid_date":[],"end_date":[],"rb":[],"E0":[]}
LT_prior = {"rb":1.0,"E0":100}
window_length = 30
window_offset = 5
start_date = ldt[0]
last_date = ldt[-1]
end_date = start_date+datetime.timedelta(days=window_length)
last_E0_OK = False
while end_date<=last_date:
si = qcutils.GetDateIndex(ldt,str(start_date),ts=ts)
ei = qcutils.GetDateIndex(ldt,str(end_date),ts=ts)
Tsub = numpy.ma.compressed(T_night[si:ei+1])
ERsub = numpy.ma.compressed(ER[si:ei+1])
LT_results["start_date"].append(start_date)
LT_results["mid_date"].append(start_date+(start_date-end_date)/2)
LT_results["end_date"].append(end_date)
if len(ERsub)>=10:
LT_prior["rb"] = numpy.mean(ERsub)
p0 = [LT_prior["rb"],LT_prior["E0"]]
popt,pcov = irls_leastsq(residuals_LloydTaylor,p0,
args=(ERsub,Tsub),maxfev=3,
weight_type='Huber',mode='quiet')
# QC results
if popt[1]<50 or popt[1]>400:
if last_E0_OK:
popt[1] = LT_results["E0"][-1]
last_E0_OK = False
else:
if popt[1]<50: popt[1] = float(50)
if popt[1]>400: popt[1] = float(400)
last_E0_OK = False
else:
last_E0_OK = True
else:
popt[0] = LT_prior["rb"]
popt[1] = LT_prior["E0"]
LT_results["rb"].append(popt[0])
LT_results["E0"].append(popt[1])
start_date = start_date+datetime.timedelta(days=window_offset)
end_date = start_date+datetime.timedelta(days=window_length)
In [10]:
fig = plt.figure(figsize=(24,6))
ax1 = plt.subplot(211)
plt.plot(LT_results["end_date"],LT_results["rb"],'b.')
plt.ylabel("rb")
ax2 = plt.subplot(212)
plt.plot(LT_results["end_date"],LT_results["E0"],'b.')
plt.ylabel("E0")
plt.tight_layout()
plt.show()
In [10]:
indicator_day = numpy.copy(indicator)
# apply a day/night filter
idx = numpy.where(Fsd<10)[0]
indicator_day[idx] = numpy.int(0)
In [11]:
# synchronise the gaps and apply the day/night filter
Fsd_day = numpy.ma.masked_where(indicator_day==0,Fsd)
D_day = numpy.ma.masked_where(indicator_day==0,D)
T_day = numpy.ma.masked_where(indicator_day==0,T)
NEE_day = numpy.ma.masked_where(indicator_day==0,Fc)
In [53]:
LL_results = {"start_date":[],"mid_date":[],"end_date":[],
"alpha":[],"beta":[],"k":[],"rb":[]}
LL_prior = {"rb":1.0,"alpha":0.01,"beta":10,"k":0}
LL_fixed = {"E0":100,"D0":1}
window_length = 30
window_offset = 5
start_date = ldt[0]
last_date = ldt[-1]
end_date = start_date+datetime.timedelta(days=window_length)
while end_date<=last_date:
si = qcutils.GetDateIndex(ldt,str(start_date),ts=ts)
ei = qcutils.GetDateIndex(ldt,str(end_date),ts=ts)
Fsdsub = numpy.ma.compressed(Fsd_day[si:ei+1])
Dsub = numpy.ma.compressed(D_day[si:ei+1])
Tsub = numpy.ma.compressed(T_day[si:ei+1])
NEEsub = numpy.ma.compressed(NEE_day[si:ei+1])
ERsub = numpy.ma.compressed(ER[si:ei+1])
LL_results["start_date"].append(start_date)
mid_date = start_date+(start_date-end_date)/2
LL_results["mid_date"].append(mid_date)
LL_results["end_date"].append(end_date)
if len(NEEsub)>=10:
LL_prior["rb"] = numpy.mean(ERsub)
LL_prior["beta"] = numpy.abs(numpy.percentile(NEEsub,3)-numpy.percentile(NEEsub,97))
diffs = [abs(dt-mid_date) for dt in LT_results["mid_date"]]
val,idx = min((val,idx) for (idx,val) in enumerate(diffs))
LL_fixed["E0"] = LT_results["E0"][idx]
p0 = [LL_prior["alpha"],LL_prior["beta"],LL_prior["k"],LL_prior["rb"]]
popt,pcov = irls_leastsq(residuals_RHLRC_D,p0,
args=(NEEsub,Fsdsub,Dsub,Tsub,LL_fixed["D0"],LL_fixed["E0"]),
maxfev=3,weight_type='Huber',mode="quiet")
# need to put parameter QC here
# if rb is negative, discard whole parameter set
if popt[3]<=0: popt = [numpy.nan,numpy.nan,numpy.nan,numpy.nan]
# if beta is negative or greater than 250, discard whole parameter set
if popt[1]<=0:
popt = [numpy.nan,numpy.nan,numpy.nan,numpy.nan]
if popt[1]>=250:
popt = [numpy.nan,numpy.nan,numpy.nan,numpy.nan]
# if k is negative, discard whole parameter set
if popt[2]<=0: popt = [numpy.nan,numpy.nan,numpy.nan,numpy.nan]
else:
#popt = [LL_prior["alpha"],LL_prior["beta"],LL_prior["k"],LL_prior["rb"]]
popt = [numpy.nan,numpy.nan,numpy.nan,numpy.nan]
LL_results["alpha"].append(popt[0])
LL_results["beta"].append(popt[1])
LL_results["k"].append(popt[2])
LL_results["rb"].append(popt[3])
start_date = start_date+datetime.timedelta(days=window_offset)
end_date = start_date+datetime.timedelta(days=window_length)
In [54]:
fig = plt.figure(figsize=(24,6))
ax1 = plt.subplot(411)
plt.plot(LT_results["end_date"],LT_results["rb"],'b.')
plt.plot(LL_results["end_date"],LL_results["rb"],'r+')
plt.ylabel("rb")
ax2 = plt.subplot(412,sharex=ax1)
plt.plot(LL_results["end_date"],LL_results["alpha"],'b.')
plt.ylabel("alpha")
ax3 = plt.subplot(413,sharex=ax1)
plt.plot(LL_results["end_date"],LL_results["beta"],'b.')
plt.ylabel("beta")
ax4 = plt.subplot(414,sharex=ax1)
plt.plot(LL_results["end_date"],LL_results["k"],'b.')
plt.ylabel("k")
plt.tight_layout()
plt.show()
In [ ]: