In [1]:
%run basics
%matplotlib
from scipy.optimize import curve_fit
In [2]:
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 NEE_RHLRC_D(data,alpha,beta,k,D0,rb,E0):
Fsd = data["Fsd"]
D = data["D"]
T = data["T"]
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]
if isinstance(k,numpy.ndarray):
SHD_func[idx] = numpy.exp(-k[idx]*(D[idx]-D0))
else:
SHD_func[idx] = numpy.exp(-k*(D[idx]-D0))
return SHD_func
In [4]:
def interp_params(param_rslt_array):
def do_interp(array_1D):
xp = numpy.arange(len(arr))
fp = array_1D[:]
nan_index = numpy.isnan(fp)
fp[nan_index] = numpy.interp(xp[nan_index], xp[~nan_index], fp[~nan_index])
return fp
arr = param_rslt_array.copy()
num_vars = numpy.shape(arr)
if len(num_vars) == 1:
arr = do_interp(arr)
else:
num_vars = num_vars[1]
for i in range(num_vars):
arr[:, i] = do_interp(arr[:, i])
return arr
In [5]:
def ApplyQuantileFilter(data,quantile):
lower = float(quantile)
upper = float(100)-lower
q = numpy.percentile(numpy.ma.compressed(data),[lower,upper])
return numpy.ma.masked_where((data<q[0])|(data>q[1]),data)
In [6]:
def get_LT_params(ldt,ER,T,info):
"""
Purpose:
Returns rb and E0 for the Lloyd & Taylor respiration function.
Usage:
Author: PRI
Date: April 2016
"""
mta = numpy.array([])
LT_results = {"start_date":mta,"mid_date":mta,"end_date":mta,
"rb":mta,"E0":mta}
LT_prior = {"rb":1.0,"E0":100}
start_date = ldt[0]
last_date = ldt[-1]
end_date = start_date+datetime.timedelta(days=info["window_length"])
last_E0_OK = False
while end_date<=last_date:
LT_results["start_date"] = numpy.append(LT_results["start_date"],start_date)
LT_results["mid_date"] = numpy.append(LT_results["mid_date"],start_date+(end_date-start_date)/2)
LT_results["end_date"] = numpy.append(LT_results["end_date"],end_date)
si = qcutils.GetDateIndex(ldt,str(start_date),ts=ts)
ei = qcutils.GetDateIndex(ldt,str(end_date),ts=ts)
Tsub = numpy.ma.compressed(T[si:ei+1])
ERsub = numpy.ma.compressed(ER[si:ei+1])
if len(ERsub)>=10:
LT_prior["rb"] = numpy.mean(ERsub)
p0 = [LT_prior["rb"],LT_prior["E0"]]
popt,pcov = curve_fit(ER_LloydTaylor,Tsub,ERsub,p0=p0)
# 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
LT_results["rb"] = numpy.append(LT_results["rb"],popt[0])
LT_results["E0"] = numpy.append(LT_results["E0"],popt[1])
else:
LT_results["rb"] = numpy.append(LT_results["rb"],numpy.nan)
LT_results["E0"] = numpy.append(LT_results["E0"],numpy.nan)
start_date = start_date+datetime.timedelta(days=info["window_offset"])
end_date = start_date+datetime.timedelta(days=info["window_length"])
# start_date = end_date
# end_date = start_date+dateutil.relativedelta.relativedelta(years=1)
return LT_results
In [7]:
def get_LL_params(ldt,Fsd,D,T,NEE,ER,LT_results,info):
# Lasslop as it was written in Lasslop et al (2010), mostly ...
# Actually, the only intended difference is the window length and offset
# Lasslop et al used window_length=4, window_offset=2
mta = numpy.array([])
LL_results = {"start_date":mta,"mid_date":mta,"end_date":mta,
"alpha":mta,"beta":mta,"k":mta,"rb":mta,
"alpha_low":mta,"rb_low":mta,"rb_prior":mta,"E0":mta}
LL_prior = {"rb":1.0,"alpha":0.01,"beta":10,"k":0}
LL_fixed = {"D0":1}
D0 = LL_fixed["D0"]
drivers = {}
start_date = ldt[0]
last_date = ldt[-1]
end_date = start_date+datetime.timedelta(days=info["window_length"])
while end_date<=last_date:
#print start_date,end_date
sub_results = {"RMSE":[],"alpha":[],"beta":[],"k":[],"rb":[]}
si = qcutils.GetDateIndex(ldt,str(start_date),ts=ts)
ei = qcutils.GetDateIndex(ldt,str(end_date),ts=ts)
drivers["Fsd"] = numpy.ma.compressed(Fsd[si:ei+1])
drivers["D"] = numpy.ma.compressed(D[si:ei+1])
drivers["T"] = numpy.ma.compressed(T[si:ei+1])
NEEsub = numpy.ma.compressed(NEE[si:ei+1])
ERsub = numpy.ma.compressed(ER[si:ei+1])
mid_date = start_date+(end_date-start_date)/2
# get the value of E0 for the period closest to the mid-point of this period
diffs = [abs(dt-mid_date) for dt in LT_results["mid_date"]]
val,idx = min((val,idx) for (idx,val) in enumerate(diffs))
LL_results["E0"] = numpy.append(LL_results["E0"],LT_results["E0_int"][idx])
LL_results["start_date"] = numpy.append(LL_results["start_date"],start_date)
LL_results["mid_date"] = numpy.append(LL_results["mid_date"],mid_date)
LL_results["end_date"] = numpy.append(LL_results["end_date"],end_date)
if len(NEEsub)>=10:
# alpha and rb from linear fit between NEE and Fsd at low light levels
idx = numpy.where(drivers["Fsd"]<100)[0]
alpha_low,rb_low = numpy.polyfit(drivers["Fsd"][idx],NEEsub[idx],1)
if len(ERsub)>=10: LL_prior["rb"] = numpy.mean(ERsub)
for bm in [0.5,1,2]:
#print "Doing beta multiplier: ",bm
LL_prior["beta"] = numpy.abs(numpy.percentile(NEEsub,3)-numpy.percentile(NEEsub,97))
LL_prior["beta"] = bm*LL_prior["beta"]
E0 = LL_results["E0"][-1]
p0 = [LL_prior["alpha"],LL_prior["beta"],LL_prior["k"],LL_prior["rb"]]
try:
fopt = lambda x,alpha,beta,k,rb:NEE_RHLRC_D(x,alpha,beta,k,D0,rb,E0)
popt,pcov = curve_fit(fopt,drivers,NEEsub,p0=p0)
alpha,beta,k,rb = popt[0],popt[1],popt[2],popt[3]
last_alpha_OK = True
except RuntimeError:
#print " Setting all parameters to NaN: 1"
alpha,beta,k,rb = numpy.nan,numpy.nan,numpy.nan,numpy.nan
last_alpha_OK = False
# QC the parameters
# k first
if numpy.isnan(k) or k<0 or k>2:
k = 0
try:
p0 = [LL_prior["alpha"],LL_prior["beta"],LL_prior["rb"]]
fopt = lambda x,alpha,beta,rb:NEE_RHLRC_D(x,alpha,beta,k,D0,rb,E0)
popt,pcov = curve_fit(fopt,drivers,NEEsub,p0=p0)
alpha,beta,rb = popt[0],popt[1],popt[2]
last_alpha_OK = True
except RuntimeError:
#print " Setting all parameters to NaN: 2"
alpha,beta,k,rb = numpy.nan,numpy.nan,numpy.nan,numpy.nan
last_alpha_OK = False
# then alpha
if numpy.isnan(alpha) or alpha<0 or alpha>0.22:
if last_alpha_OK==True:
alpha = LL_results["alpha"][-1]
else:
alpha = 0
try:
p0 = [LL_prior["beta"],LL_prior["k"],LL_prior["rb"]]
fopt = lambda x,beta,k,rb:NEE_RHLRC_D(x,alpha,beta,k,D0,rb,E0)
popt,pcov = curve_fit(fopt,drivers,NEEsub,p0=p0)
beta,k,rb = popt[0],popt[1],popt[2]
except RuntimeError:
#print " Setting all parameters to NaN: 3"
alpha,beta,k,rb = numpy.nan,numpy.nan,numpy.nan,numpy.nan
# then beta
if beta<0:
beta = 0
try:
p0 = [LL_prior["alpha"],LL_prior["k"],LL_prior["rb"]]
fopt = lambda x,alpha,k,rb:NEE_RHLRC_D(x,alpha,beta,k,D0,rb,E0)
popt,pcov = curve_fit(fopt,drivers,NEEsub,p0=p0)
alpha,k,rb = popt[0],popt[1],popt[2]
except RuntimeError:
#print " Setting all parameters to NaN: 4"
alpha,beta,k,rb = numpy.nan,numpy.nan,numpy.nan,numpy.nan
elif beta>250:
#print " Setting all parameters to NaN: 5"
alpha,beta,k,rb = numpy.nan,numpy.nan,numpy.nan,numpy.nan
# and finally rb
if rb<0:
#print " Setting all parameters to NaN: 6"
alpha,beta,k,rb = numpy.nan,numpy.nan,numpy.nan,numpy.nan
# now get the RMSE for this set of parameters
if not numpy.isnan(alpha) and not numpy.isnan(beta) and not numpy.isnan(k) and not numpy.isnan(rb):
NEEest = NEE_RHLRC_D(drivers,alpha,beta,k,D0,rb,E0)
sub_results["RMSE"].append(numpy.sqrt(numpy.mean((NEEsub-NEEest)**2)))
sub_results["alpha"].append(alpha)
sub_results["beta"].append(beta)
sub_results["k"].append(k)
sub_results["rb"].append(rb)
# now find the minimum RMSE and the set of parameters for the minimum
if len(sub_results["RMSE"])!=0:
min_RMSE = min(sub_results["RMSE"])
idx = sub_results["RMSE"].index(min_RMSE)
LL_results["alpha"] = numpy.append(LL_results["alpha"],sub_results["alpha"][idx])
LL_results["alpha_low"] = numpy.append(LL_results["alpha_low"],float(-1)*alpha_low)
LL_results["rb"] = numpy.append(LL_results["rb"],sub_results["rb"][idx])
LL_results["rb_low"] = numpy.append(LL_results["rb_low"],rb_low)
LL_results["rb_prior"] = numpy.append(LL_results["rb_prior"],LL_prior["rb"])
LL_results["beta"] = numpy.append(LL_results["beta"],sub_results["beta"][idx])
LL_results["k"] = numpy.append(LL_results["k"],sub_results["k"][idx])
else:
LL_results["alpha"] = numpy.append(LL_results["alpha"],numpy.nan)
LL_results["alpha_low"] = numpy.append(LL_results["alpha_low"],float(-1)*alpha_low)
LL_results["rb"] = numpy.append(LL_results["rb"],numpy.nan)
LL_results["rb_low"] = numpy.append(LL_results["rb_low"],rb_low)
LL_results["rb_prior"] = numpy.append(LL_results["rb_prior"],LL_prior["rb"])
LL_results["beta"] = numpy.append(LL_results["beta"],numpy.nan)
LL_results["k"] = numpy.append(LL_results["k"],numpy.nan)
else:
LL_results["alpha"] = numpy.append(LL_results["alpha"],numpy.nan)
LL_results["alpha_low"] = numpy.append(LL_results["alpha_low"],float(-1)*alpha_low)
LL_results["rb"] = numpy.append(LL_results["rb"],numpy.nan)
LL_results["rb_low"] = numpy.append(LL_results["rb_low"],rb_low)
LL_results["rb_prior"] = numpy.append(LL_results["rb_prior"],LL_prior["rb"])
LL_results["beta"] = numpy.append(LL_results["beta"],numpy.nan)
LL_results["k"] = numpy.append(LL_results["k"],numpy.nan)
# update the start and end datetimes
start_date = start_date+datetime.timedelta(days=info["window_offset"])
end_date = start_date+datetime.timedelta(days=info["window_length"])
LL_results["D0"] = D0
return LL_results
In [8]:
def plot_LTparams_ER(ldt,ER,ER_LT,LT_results):
fig, axs = plt.subplots(3,1,sharex=True,figsize=(24,6))
axs[0].plot(LT_results["mid_date"],LT_results["rb"],'bo')
axs[0].set_ylabel("rb (umol/m2/s)")
axs[1].plot(LT_results["mid_date"],LT_results["E0"],'bo')
axs[1].set_ylabel("E0 (C)")
axs[2].plot(ldt,ER,'bo')
axs[2].plot(ldt,ER_LT,'r--')
axs[2].axhline(y=0,linewidth=4,color="r")
axs[2].set_ylabel("ER (umol/m2/s)")
plt.tight_layout()
plt.draw()
In [ ]:
def plot_LLparams(LT_results,LL_results):
fig, axs = plt.subplots(4,1,sharex=True,figsize=(24,6))
axs[0].plot(LT_results["mid_date"],LT_results["rb"],'bo')
axs[0].plot(LL_results["mid_date"],LL_results["rb"],'ro')
axs[0].plot(LL_results["mid_date"],LL_results["rb_low"],'go')
axs[0].plot(LL_results["mid_date"],LL_results["rb_prior"],'yo')
axs[0].set_ylabel("rb")
axs[1].plot(LL_results["mid_date"],LL_results["alpha"],'bo')
axs[1].plot(LL_results["mid_date"],LL_results["alpha_low"],'ro')
axs[1].set_ylabel("alpha")
axs[2].plot(LL_results["mid_date"],LL_results["beta"],'bo')
axs[2].set_ylabel("beta")
axs[3].plot(LL_results["mid_date"],LL_results["k"],'bo')
axs[3].set_ylabel("k")
plt.tight_layout()
plt.show()
In [9]:
#ncname ="../Sites/Whroo/Data/Portal/Whroo_L5.nc"
ncname ="../Sites/Calperum/Data/Portal/Calperum_L5.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 [10]:
# 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 [11]:
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.24)[0]
indicator_night[idx] = numpy.int(0)
In [12]:
# synchronise the gaps and apply the ustar filter
T_night = numpy.ma.masked_where(indicator_night==0,T)
ustar_night = numpy.ma.masked_where(indicator_night==0,ustar)
ER = numpy.ma.masked_where(indicator_night==0,Fc)
In [13]:
fig, axs = plt.subplots(3,1,sharex=True,figsize=(24,6))
axs[0].plot(ldt,ustar_night,'bo')
axs[0].set_ylabel("ustar (m/s)")
axs[1].plot(ldt,T_night,'bo')
axs[1].set_ylabel("T (C)")
axs[2].plot(ldt,ER,'bo')
axs[2].axhline(y=0,linewidth=4,color="r")
axs[2].set_ylabel("ER (umol/m2/s)")
plt.tight_layout()
plt.draw()
In [14]:
info = {"window_length":30,"window_offset":5}
LT_results = get_LT_params(ldt,ER,T_night,info)
In [15]:
LT_results["rb_int"] = interp_params(LT_results["rb"])
LT_results["E0_int"] = interp_params(LT_results["E0"])
In [16]:
fig,axs = plt.subplots(2,sharex=True,figsize=(24,6))
axs[0].plot(LT_results["mid_date"],LT_results["rb"],'b.')
axs[0].plot(LT_results["mid_date"],LT_results["rb_int"],'r+')
axs[1].plot(LT_results["mid_date"],LT_results["E0"],'b.')
axs[1].plot(LT_results["mid_date"],LT_results["E0_int"],'r+')
plt.tight_layout()
plt.draw()
In [17]:
ntsperday = float(24)*float(60)/float(ts)
days_at_beginning = float(info["window_length"])/2 - float(info["window_offset"])/2
rb_beginning = numpy.ones(days_at_beginning*ntsperday)*LT_results["rb_int"][0]
rb_middle = numpy.repeat(LT_results["rb_int"],info["window_offset"]*ntsperday)
nend = len(ldt) - (len(rb_beginning)+len(rb_middle))
rb_end = numpy.ones(nend)*LT_results["rb_int"][-1]
rb_tts = numpy.concatenate((rb_beginning,rb_middle,rb_end))
E0_beginning = numpy.ones(days_at_beginning*ntsperday)*LT_results["E0_int"][0]
E0_middle = numpy.repeat(LT_results["E0_int"],info["window_offset"]*ntsperday)
nend = len(ldt) - (len(E0_beginning)+len(E0_middle))
E0_end = numpy.ones(nend)*LT_results["E0_int"][-1]
E0_tts = numpy.concatenate((E0_beginning,E0_middle,E0_end))
ER_LT = ER_LloydTaylor(T,rb_tts,E0_tts)
In [18]:
plot_LTparams_ER(ldt,ER,ER_LT,LT_results)
In [19]:
indicator_day = numpy.copy(indicator)
# apply a day/night filter
idx = numpy.where(Fsd<=10)[0]
indicator_day[idx] = numpy.int(0)
In [20]:
# 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 [21]:
LL_results = get_LL_params(ldt,Fsd_day,D_day,T_day,NEE_day,ER,LT_results,info)
In [22]:
LL_results["alpha_int"] = interp_params(LL_results["alpha"])
LL_results["beta_int"] = interp_params(LL_results["beta"])
LL_results["k_int"] = interp_params(LL_results["k"])
LL_results["rb_int"] = interp_params(LL_results["rb"])
LL_results["E0_int"] = interp_params(LL_results["E0"])
In [23]:
plot_LLparams(LT_results,LL_results)
In [24]:
ntsperday = float(24)*float(60)/float(ts)
days_at_beginning = float(info["window_length"])/2 - float(info["window_offset"])/2
int_list = ["alpha_int","beta_int","k_int","rb_int","E0_int"]
tts_list = ["alpha_tts","beta_tts","k_tts","rb_tts","E0_tts"]
for tts_item,int_item in zip(tts_list,int_list):
beginning = numpy.ones(days_at_beginning*ntsperday)*LL_results[int_item][0]
middle = numpy.repeat(LL_results[int_item],info["window_offset"]*ntsperday)
nend = len(ldt) - (len(rb_beginning)+len(rb_middle))
end = numpy.ones(nend)*LL_results[int_item][-1]
LL_results[tts_item] = numpy.concatenate((beginning,middle,end))
In [25]:
D0 = LL_results["D0"]
rb = LL_results["rb_tts"]
E0 = LL_results["E0_tts"]
ER_LL = ER_LloydTaylor(T,rb,E0)
alpha = LL_results["alpha_tts"]
beta = LL_results["beta_tts"]
k = LL_results["k_tts"]
GPP_LL = GPP_RHLRC_D(Fsd,D,alpha,beta,k,D0)
data = {"Fsd":Fsd,"T":T,"D":D}
NEE_LL = NEE_RHLRC_D(data,alpha,beta,k,D0,rb,E0)
In [26]:
fig, axs = plt.subplots(2,1,sharex=True,figsize=(24,6))
axs[0].plot(ldt,ER,'bo')
axs[0].plot(ldt,ER_LT,'r--')
axs[0].plot(ldt,ER_LL,'g-')
axs[0].set_ylabel("ER (umol/m2/s)")
axs[0].axhline(y=0,linewidth=4,color="r")
axs[1].plot(ldt,Fc,'bo')
axs[1].plot(ldt,NEE_LL,'r-',linewidth=4)
axs[1].set_ylabel("NEE (umol/m2/s)")
axs[1].axhline(y=0,linewidth=3,color="g")
plt.tight_layout()
plt.show()
In [ ]:
# Lasslop with k fixed at 0.2
# Actually, the only intended difference is the window length and offset
# Lasslop et al used window_length=4, window_offset=2
LL_results = {"start_date":[],"mid_date":[],"end_date":[],
"alpha":[],"beta":[],"k":[],"rb":[],
"alpha_low":[],"rb_low":[]}
LL_prior = {"rb":1.0,"alpha":0.01,"beta":10}
LL_fixed = {"E0":100,"D0":1,"k":0.2}
drivers = {}
window_length = 30
window_offset = 5
start_date = ldt[0]
last_date = ldt[-1]
#last_date = start_date+datetime.timedelta(days=5*window_length)
end_date = start_date+datetime.timedelta(days=window_length)
while end_date<=last_date:
#print start_date,end_date
sub_results = {"RMSE":[],"alpha":[],"beta":[],"k":[],"rb":[]}
si = qcutils.GetDateIndex(ldt,str(start_date),ts=ts)
ei = qcutils.GetDateIndex(ldt,str(end_date),ts=ts)
drivers["Fsd"] = numpy.ma.compressed(Fsd_day[si:ei+1])
drivers["D"] = numpy.ma.compressed(D_day[si:ei+1])
drivers["T"] = 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])
mid_date = start_date+(start_date-end_date)/2
if len(NEEsub)>=10:
# alpha and rb from linear fit between NEE and Fsd at low light levels
idx = numpy.where(drivers["Fsd"]<100)[0]
alpha_low,rb_low = numpy.polyfit(drivers["Fsd"][idx],NEEsub[idx],1)
LL_prior["rb"] = numpy.mean(ERsub)
for bm in [0.5,1,2]:
#print "Doing beta multiplier: ",bm
LL_prior["beta"] = numpy.abs(numpy.percentile(NEEsub,3)-numpy.percentile(NEEsub,97))
LL_prior["beta"] = bm*LL_prior["beta"]
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]
E0 = LL_fixed["E0"]
D0 = LL_fixed["D0"]
k = LL_fixed["k"]
p0 = [LL_prior["alpha"],LL_prior["beta"],LL_prior["rb"]]
try:
fopt = lambda x,alpha,beta,rb:NEE_RHLRC_D(x,alpha,beta,k,D0,rb,E0)
popt,pcov = curve_fit(fopt,drivers,NEEsub,p0=p0)
alpha,beta,rb = popt[0],popt[1],popt[2]
last_alpha_OK = True
except RuntimeError:
#print " Setting all parameters to NaN: 1"
alpha,beta,rb = numpy.nan,numpy.nan,numpy.nan
last_alpha_OK = False
# QC the parameters
# then alpha
if numpy.isnan(alpha) or alpha<0 or alpha>0.22:
if last_alpha_OK==True:
alpha = LL_results["alpha"][-1]
else:
alpha = 0
try:
p0 = [LL_prior["beta"],LL_prior["rb"]]
fopt = lambda x,beta,rb:NEE_RHLRC_D(x,alpha,beta,k,D0,rb,E0)
popt,pcov = curve_fit(fopt,drivers,NEEsub,p0=p0)
beta,rb = popt[0],popt[1]
except RuntimeError:
#print " Setting all parameters to NaN: 3"
alpha,beta,rb = numpy.nan,numpy.nan,numpy.nan
# then beta
if beta<0:
beta = 0
try:
p0 = [LL_prior["alpha"],LL_prior["rb"]]
fopt = lambda x,alpha,rb:NEE_RHLRC_D(x,alpha,beta,k,D0,rb,E0)
popt,pcov = curve_fit(fopt,drivers,NEEsub,p0=p0)
alpha,rb = popt[0],popt[1]
except RuntimeError:
#print " Setting all parameters to NaN: 4"
alpha,beta,rb = numpy.nan,numpy.nan,numpy.nan
elif beta>250:
#print " Setting all parameters to NaN: 5"
alpha,beta,rb = numpy.nan,numpy.nan,numpy.nan
# and finally rb
if rb<0:
#print " Setting all parameters to NaN: 6"
alpha,beta,rb = numpy.nan,numpy.nan,numpy.nan
# now get the RMSE for this set of parameters
if not numpy.isnan(alpha) and not numpy.isnan(beta) and not numpy.isnan(rb):
NEEest = NEE_RHLRC_D(drivers,alpha,beta,k,D0,rb,E0)
sub_results["RMSE"].append(numpy.sqrt(numpy.mean((NEEsub-NEEest)**2)))
sub_results["alpha"].append(alpha)
sub_results["beta"].append(beta)
sub_results["rb"].append(rb)
# now find the minimum RMSE and the set of parameters for the minimum
if len(sub_results["RMSE"])!=0:
LL_results["start_date"].append(start_date)
LL_results["mid_date"].append(mid_date)
LL_results["end_date"].append(end_date)
min_RMSE = min(sub_results["RMSE"])
idx = sub_results["RMSE"].index(min_RMSE)
LL_results["alpha"].append(sub_results["alpha"][idx])
alpha_low = float(-1)*alpha_low
LL_results["alpha_low"].append(alpha_low)
LL_results["rb_low"].append(rb_low)
LL_results["beta"].append(sub_results["beta"][idx])
LL_results["rb"].append(sub_results["rb"][idx])
# update the start and end datetimes
start_date = start_date+datetime.timedelta(days=window_offset)
end_date = start_date+datetime.timedelta(days=window_length)
In [ ]:
fig = plt.figure(figsize=(24,6))
ax1 = plt.subplot(311)
plt.plot(LT_results["end_date"],LT_results["rb"],'b.')
plt.plot(LL_results["end_date"],LL_results["rb"],'r+')
plt.plot(LL_results["end_date"],LL_results["rb_low"],'g^')
plt.ylabel("rb")
ax2 = plt.subplot(312,sharex=ax1)
plt.plot(LL_results["end_date"],LL_results["alpha"],'b.')
plt.plot(LL_results["end_date"],LL_results["alpha_low"],'r+')
plt.ylabel("alpha")
ax3 = plt.subplot(313,sharex=ax1)
plt.plot(LL_results["end_date"],LL_results["beta"],'b.')
plt.ylabel("beta")
plt.tight_layout()
plt.show()
In [ ]: