Figure Set 1: FATES_NPP (for C-only at that site and RD/ECA)
Figure Set 2: FATES_NPP
Figure Set 3: FATES_STOREC
Figure Set 3: FATES_STOREC_TF
Figure Set 4: FATES_STOREN
Figure Set 4: FATES_STOREN_TF
Figure Set 5: FATES_STOREC_TF / FATES_STOREN_TF
Figure Set 3: FATES_L2FR
Figure Set 6: FATES_FROOTC_ALLOC / FATES_LEAFC_ALLOC
Figure Set 7: FATES_FROOTC / FATES_LEAFC
Figure Set 7: FATES_FROOTC
Figure Set 7: FATES_LEAFC
Figure Set 8: FATES_NO3UPTAKE
Figure Set 8: FATES_NH4UPTAKE
Figure Set 9: NET_NMIN
Figure Set 9: GROSS_NMIN
Figure Set10: SMIN_NO3_LEACHED
Figure Set10: DENIT
for a cell, you can choose which variables you are interested in for a single simulation.
some plots from postprocessor
In [100]:
importos,globimportxarrayasxrimportmatplotlib.pyplotaspltimportnumpyasnpimportpandasaspdimportpyreadr# to read .rds filesimportdatetimeasdtnp.set_printoptions(suppress=True)
Select the timeperiod you are interested to plot for daily timeseries¶
In [107]:
# Logging happened on Jan 1, 1855# Start Date Index idx_start_date=4*365+1idx_logging=5*365+1# For multiple variables on a same plot upto 5 years after logging Yrs_Next=5# upto 5 years after logging
print("\nAll Plots >>>>>>")fig,axs=plt.subplots(9,1,figsize=(15,30),dpi=400)# Subplot 1 >>>idx_ax=0var="FATES_NPP"# C-onlytry:sim_conly=sims.replace("RD","Conly")ifsims=="DUK_PIDE_RD_AgBgW_testing":sim_conly="DUK_PIDE_Conly_AgBgW"except:sim_conly="DUK_PIDE_Conly_AgBgW"key=sim_conlyts_data=ds[key][var][idx_start_date:idx_start_date+Yrs_Next*365]# Daily data range from logging date to next five yearsaxs[idx_ax].plot(ts_data.time,ts_data,label=var+" C-only")# RDkey=simsts_data=ds[key][var][idx_start_date:idx_start_date+Yrs_Next*365]# Daily data range from logging date to next five yearsaxs[idx_ax].plot(ts_data.time,ts_data,label=var)axs[idx_ax].set_ylabel(f"{ds[key][var].units}")axs[idx_ax].set_title(f"{key} Daily Logging + 5 years\n",fontsize=16)axs[idx_ax].axhline(y=0,color='r',lw=1,alpha=.1)print(key)# Subplot 1 <<<# Subplot 2 >>>idx_ax=1vars_plot=("""FATES_STOREC_TFFATES_STOREC""").split('\n')vars_plot=vars_plot[1:-1]var=vars_plot[0]ts_data=ds[key][var][idx_start_date:idx_start_date+Yrs_Next*365]# Daily data range from logging date to next five yearsaxs[idx_ax].plot(ts_data.time,ts_data,label=var,color='b')axs[idx_ax].set_ylabel(f"{var}({ds[key][var].units})",color='b')ax_tmp=axs[idx_ax].twinx()var=vars_plot[1]ts_data=ds[key][var][idx_start_date:idx_start_date+Yrs_Next*365]# Daily data range from logging date to next five yearsax_tmp.plot(ts_data.time,ts_data,'k.',label=var,alpha=.3)ax_tmp.legend(loc=5)ax_tmp.set_ylabel(f"{var}({ds[key][var].units})")ax_tmp.axhline(y=0,color='r',lw=1,alpha=.1)# Subplot 2 <<<# Subplot 3 >>>idx_ax=2vars_plot=("""FATES_STOREN_TFFATES_STOREN""").split('\n')vars_plot=vars_plot[1:-1]var=vars_plot[0]ts_data=ds[key][var][idx_start_date:idx_start_date+Yrs_Next*365]# Daily data range from logging date to next five yearsaxs[idx_ax].plot(ts_data.time,ts_data,label=var,color='b')axs[idx_ax].set_ylabel(f"{var}({ds[key][var].units})",color='b')ax_tmp=axs[idx_ax].twinx()var=vars_plot[1]ts_data=ds[key][var][idx_start_date:idx_start_date+Yrs_Next*365]# Daily data range from logging date to next five yearsax_tmp.plot(ts_data.time,ts_data,'k.',label=var,alpha=.3)ax_tmp.set_ylabel(f"{var}({ds[key][var].units})")ax_tmp.legend(loc=5)ax_tmp.axhline(y=0,color='r',lw=1,alpha=.1)# Subplot 3 <<<# Subplot 4 >>>idx_ax=3ts_data=(ds[key]["FATES_STOREC_TF"]/ds[key]["FATES_STOREN_TF"])[idx_start_date:idx_start_date+Yrs_Next*365]# Daily data range from logging date to next five yearsaxs[idx_ax].plot(ts_data.time,ts_data,label="FATES_STOREC_TF/FATES_STOREN_TF",color='b')axs[idx_ax].set_ylabel(f"FATES_STOREC_TF/FATES_STOREN_TF",color='b')ax_tmp=axs[idx_ax].twinx()ts_data=(ds[key]["FATES_STOREC"]/ds[key]["FATES_STOREN"])[idx_start_date:idx_start_date+Yrs_Next*365]# Daily data range from logging date to next five yearsax_tmp.plot(ts_data.time,ts_data,'k.',label="FATES_STOREC/FATES_STOREN")ax_tmp.legend(loc=5)ax_tmp.set_ylabel(f"{var}({ds[key][var].units})")# Subplot 4 <<<# Subplot 5 >>>idx_ax=4ts_data=(ds[key]["FATES_FROOT_ALLOC"]/ds[key]["FATES_LEAF_ALLOC"])[idx_start_date:idx_start_date+Yrs_Next*365]# Daily data range from logging date to next five yearsaxs[idx_ax].plot(ts_data.time,ts_data,label="FATES_FROOT_ALLOC/FATES_LEAF_ALLOC",color='b')axs[idx_ax].set_ylabel(f"FATES_FROOT_ALLOC/FATES_LEAF_ALLOC",color='b')ax_tmp=axs[idx_ax].twinx()var="FATES_L2FR"ts_data=ds[key][var][idx_start_date:idx_start_date+Yrs_Next*365]# Daily data range from logging date to next five yearsax_tmp.plot(ts_data.time,ts_data,'k.',label=var)ax_tmp.set_ylabel(f"{var}({ds[key][var].units})")ax_tmp.legend(loc=5)# Subplot 5 <<<# Subplot 6 >>>idx_ax=5vars_plot=("""FATES_LEAFCFATES_FROOTC""").split('\n')vars_plot=vars_plot[1:-1]fori_var,varinenumerate(vars_plot):ts_data=ds[key][var][idx_start_date:idx_start_date+Yrs_Next*365]# Daily data range from logging date to next five yearsaxs[idx_ax].plot(ts_data.time,ts_data,label=var)axs[idx_ax].set_ylabel(f"{var}({ds[key][var].units})",color='orange')ax_tmp=axs[idx_ax].twinx()ts_data=(ds[key]["FATES_FROOTC"]/ds[key]["FATES_LEAFC"])[idx_start_date:idx_start_date+Yrs_Next*365]# Daily data range from logging date to next five yearsax_tmp.plot(ts_data.time,ts_data,'k.',label="FATES_FROOTC/FATES_LEAFC")ax_tmp.set_ylabel(f"FATES_FROOTC/FATES_LEAFC")ax_tmp.legend(loc=5)ax_tmp.set_ylabel(f"Unitless")# Subplot 6 <<<# Subplot 7 >>>idx_ax=6vars_plot=("""FATES_NO3UPTAKEFATES_NH4UPTAKE""").split('\n')vars_plot=vars_plot[1:-1]var=vars_plot[0]ts_data=ds[key][var][idx_start_date:idx_start_date+Yrs_Next*365]# Daily data range from logging date to next five yearsaxs[idx_ax].plot(ts_data.time,ts_data,label=var,color='b')axs[idx_ax].set_ylabel(f"{var}({ds[key][var].units})",color='b')axs[idx_ax].axhline(y=0,color='r',lw=1,alpha=.1)ax_tmp=axs[idx_ax].twinx()var=vars_plot[1]ts_data=ds[key][var][idx_start_date:idx_start_date+Yrs_Next*365]# Daily data range from logging date to next five yearsax_tmp.plot(ts_data.time,ts_data,'k.',label=var)ax_tmp.set_ylabel(f"{var}({ds[key][var].units})")ax_tmp.legend(loc=5)# Subplot 7 <<<# Subplot 8 >>>idx_ax=7vars_plot=("""GROSS_NMINNET_NMIN""").split('\n')vars_plot=vars_plot[1:-1]fori_var,varinenumerate(vars_plot):ts_data=ds[key][var][idx_start_date:idx_start_date+Yrs_Next*365]# Daily data range from logging date to next five yearsaxs[idx_ax].plot(ts_data.time,ts_data,label=var)axs[idx_ax].set_ylabel(f"{ds[key][var].units}")axs[idx_ax].axhline(y=0,color='r',lw=1,alpha=.1)# Subplot 8 <<<# Subplot 9 >>>idx_ax=8vars_plot=("""SMIN_NO3_LEACHEDDENIT""").split('\n')vars_plot=vars_plot[1:-1]fori_var,varinenumerate(vars_plot):ts_data=ds[key][var][idx_start_date:idx_start_date+Yrs_Next*365]# Daily data range from logging date to next five yearsaxs[idx_ax].plot(ts_data.time,ts_data,label=var)axs[idx_ax].set_ylabel(f"{ds[key][var].units}")axs[idx_ax].axhline(y=0,color='r',lw=1,alpha=.1)# Subplot 9 <<<print(key)foriinrange(len(axs)):axs[i].legend()axs[i].axvline(x=dt.date(1855,1,1),color='g',lw=5,alpha=.1)#fig.savefig('./Plots/' + f'{key}_daily.pdf',bbox_inches='tight')#fig.savefig('./Plots/' + f'{key}_daily.png',bbox_inches='tight')
All Plots >>>>>>
DUK_PIDC_RD_AgBgW_CN53
/Users/ud4/opt/anaconda3/envs/pyces/lib/python3.9/site-packages/dask/core.py:121: RuntimeWarning: invalid value encountered in true_divide
return func(*(_execute_task(a, cache) for a in args))
cx_int = cx_int + cx_logratio
! Reset the integrator if its sign changes
if( abs(cx_logratio)>nearzero .and. abs(cx0)>nearzero) then
if( abs(cx_logratio/abs(cx_logratio) - cx0/abs(cx0)) > nearzero ) then
cx_int = cx_logratio
end if
end if
dcxdt_ratio = cx_logratio-cx0
ema_dcxdt = pid_drv_wgt*dcxdt_ratio + (1._r8-pid_drv_wgt)*ema_dcxdt
cx0 = cx_logratio
defSafeLog(val):# As defined in PRTAllometricCNPModsafelog_min=0.001safelog_max=1000clipped_val=np.clip(val,safelog_min,safelog_max)logval=np.log(clipped_val)returnlogval
In [124]:
ifc2x=="CN"andc2x!="CP":print("case1")cn_ratio=(ds[key]["FATES_STOREC_TF"]/ds[key]["FATES_STOREN_TF"])[idx_start_date:idx_start_date+Yrs_Next*365]cx_logratio_ar=SafeLog(cn_ratio)ifc2x!="CN"andc2x=="CP":print("case2")cp_ratio=(ds[key]["FATES_STOREC_TF"]/ds[key]["FATES_STOREP_TF"])[idx_start_date:idx_start_date+Yrs_Next*365]cx_logratio_ar=SafeLog(cp_ratio)ifc2x=="CNP":print("case3")cn_ratio=(ds[key]["FATES_STOREC_TF"]/ds[key]["FATES_STOREN_TF"])[idx_start_date:idx_start_date+Yrs_Next*365]cp_ratio=(ds[key]["FATES_STOREC_TF"]/ds[key]["FATES_STOREP_TF"])[idx_start_date:idx_start_date+Yrs_Next*365]#cx_logratio_ar = SafeLog(cn_ratio)#keeping priority if CN Ratio only#cx_logratio_ar = SafeLog(cn_ratio) # fcncx_logratio_ar=SafeLog(cp_ratio)# fcncx_logratio_ar=np.ravel(cx_logratio_ar)#cx_logratio_ar = np.insert(cx_logratio_ar, 0, 0) # adding first value as 0#derivativedcxdt_ratio_mat=cx_logratio_ar[1:]-cx_logratio_ar[:-1]
#integrationpid_drv_wgt=1/20.# n-day smoothing of the derivative of the process function in the PID controllercx_int=cx_logratio_ar[0]#0cx0=cx_logratio_ar[0]#0.ema_dcxdt=ema_dcxdt_init#cx_logratio_ar[0]*(pid_drv_wgt)#0. #cx_logratio_ar = cx_logratio_ar[1:] # removing the zero that we added earliernearzero=1e-30cx_int_ar=[]dcxdt_ratio_ar=[]ema_dcxdt_ar=[]foridx,cx_logratioinenumerate(cx_logratio_ar[1:]):cx_int=cx_int+cx_logratio#Reset the integrator if its sign changesifabs(cx_logratio)>nearzeroandabs(cx0)>nearzero:ifabs(cx_logratio/abs(cx_logratio)-cx0/abs(cx0))>nearzero:cx_int=cx_logratiodcxdt_ratio=cx_logratio-cx0ema_dcxdt=pid_drv_wgt*dcxdt_ratio+(1.-pid_drv_wgt)*ema_dcxdtcx0=cx_logratioifidx==1:print(dcxdt_ratio,ema_dcxdt)cx_int_ar.append(cx_int)dcxdt_ratio_ar.append(dcxdt_ratio)ema_dcxdt_ar.append(ema_dcxdt)cx_int_ar=np.array(cx_int_ar)dcxdt_ratio_ar=np.array(dcxdt_ratio_ar)ema_dcxdt_ar=np.array(ema_dcxdt_ar)# Repeating the first value twice in the beginingcx_int_ar=np.insert(cx_int_ar,0,cx_int_ar[0])dcxdt_ratio_ar=np.insert(dcxdt_ratio_ar,0,dcxdt_ratio_ar[0])ema_dcxdt_ar=np.insert(ema_dcxdt_ar,0,ema_dcxdt_ar[0])
l2fr_delta_ar=[]foridxinrange(len(cx_logratio_ar[1:])):cx_logratio=cx_logratio_ar[idx]cx_int=cx_int_ar[idx]ema_dcxdt=ema_dcxdt_ar[idx]l2fr_delta=pid_kp*cx_logratio+pid_ki*cx_int+pid_kd*ema_dcxdtl2fr_delta_ar.append(l2fr_delta)l2fr_delta_ar=np.array(l2fr_delta_ar)# Repeating the first value twice in the beginingl2fr_delta_ar=np.insert(l2fr_delta_ar,0,l2fr_delta_ar[0])
Plots with L2FR calculated based on C/C' and N/N' from the model outputs¶
The variables that are calcuated based on the C/C' and N/N'
cx_logratio_ar : CN ratio
ema_dcxdt_ar : Smoothened derivative term PID
cx_int_ar : Integral term PID
l2fr_delta_ar : Delta l2fr
Other variables are directly ploted from model variables (except l2fr_delta_model, which is just the difference of FATES_L2FR one month apart
Leaf C
C Store
C/C'
CN Ratio
derivative/integral
proportion of there terms
multiple terms
delta l2fr
l2fr
In [130]:
fig,axs=plt.subplots(10,1,figsize=(15,40),dpi=400)idx_ax=0#FATES_LEAFCvars_plot=("""FATES_LEAFC""").split('\n')vars_plot=vars_plot[1:-1]fori_var,varinenumerate(vars_plot):ts_data=ds[key][var][idx_start_date:idx_start_date+Yrs_Next*365]# Daily data range from logging date to next five yearsaxs[idx_ax].plot(ts_data.time,ts_data,label=var)axs[idx_ax].set_ylabel(f"{ds[key][var].units}")axs[idx_ax].set_title(f"{key} Daily Logging + 5 years (PID control parameters) \n",fontsize=16)axs[idx_ax].axhline(y=0,color='k',ls='--',lw=1,alpha=.3)#ax_tmp = axs[idx_ax].twinx()#ts_data = (ds[key]["FATES_FROOTC"]/ds[key]["FATES_LEAFC"])[idx_start_date:idx_start_date+Yrs_Next*365] # Daily data range from logging date to next five years#ax_tmp.plot(ts_data.time, ts_data, 'k+', label = "FATES_FROOTC/FATES_LEAFC", markersize = 4, alpha =.1)#ax_tmp.legend(loc = 5)#ax_tmp.set_ylabel (f"Unitless")# Subplot >>>idx_ax=1+idx_axts_data=(ds[key]["FATES_STOREC"]/ds[key]["FATES_STOREC_TF"])[idx_start_date:idx_start_date+Yrs_Next*365]# Daily data range from logging date to next five yearsaxs[idx_ax].plot(ts_data.time,ts_data,label="FATES_STOREC_TARGET",alpha=.3)axs[idx_ax].set_ylabel(f"FATES_STOREC({ds[key]['FATES_STOREC'].units})")ax_tmp=axs[idx_ax].twinx()ts_data=(ds[key]["FATES_STOREN"]/ds[key]["FATES_STOREN_TF"])[idx_start_date:idx_start_date+Yrs_Next*365]# Daily data range from logging date to next five yearsax_tmp.plot(ts_data.time,ts_data,'r+',label="FATES_STOREN_TARGET",alpha=.1,markersize=6)ax_tmp.legend(loc=5)ax_tmp.set_ylabel(f"{'FATES_STOREN'}({ds[key]['FATES_STOREN'].units})",color='r')ax_tmp.axhline(y=0,color='k',ls='--',lw=1,alpha=.3)# Subplot 2 <<<# Subplot >>>idx_ax=1+idx_axts_data=(ds[key]["FATES_STOREC_TF"])[idx_start_date:idx_start_date+Yrs_Next*365]# Daily data range from logging date to next five yearsaxs[idx_ax].plot(ts_data.time,ts_data,label="FATES_STOREC_TF",alpha=.3)axs[idx_ax].set_ylabel(f"FATES_STOREC_TF({ds[key]['FATES_STOREC_TF'].units})")ax_tmp=axs[idx_ax].twinx()ts_data=(ds[key]["FATES_STOREN_TF"])[idx_start_date:idx_start_date+Yrs_Next*365]# Daily data range from logging date to next five yearsax_tmp.plot(ts_data.time,ts_data,'r+',label="FATES_STOREN_TF",alpha=.1,markersize=6)ax_tmp.legend(loc=5)ax_tmp.set_ylabel(f"{'FATES_STOREN_TF'}({ds[key]['FATES_STOREN_TF'].units})",color='r')ax_tmp.axhline(y=0,color='k',ls='--',lw=1,alpha=.3)# Subplot <<<# Subplot >>>idx_ax=1+idx_axkey=simsts_data=ds[key][var][idx_start_date:idx_start_date+Yrs_Next*365]# Daily data range from logging date to next five yearsaxs[idx_ax].plot(ts_data.time,cx_logratio_ar,label="SafeLog(FATES Store TF C:N* or $F_{CN}$*")axs[idx_ax].axhline(y=0,color='k',ls='--',lw=1,alpha=.3)# Subplot <<<# Subplot >>>idx_ax=1+idx_axaxs[idx_ax].plot(ts_data.time,ema_dcxdt_ar,label="Smoothened Derivative of Log(C:N) or $dF_{CN}$*")axs[idx_ax].axhline(y=0,color='r',lw=1,alpha=.1)#ax_tmp = axs[idx_ax].twinx()#ax_tmp.plot(ts_data.time, pid_kd * ema_dcxdt_ar, 'r+', label = "pid_kd * ema_dcxdt_ar", markersize = 4, alpha =.1 )#ax_tmp.legend(loc = 5)# Subplot <<<# Subplot >>>idx_ax=1+idx_axaxs[idx_ax].plot(ts_data.time,cx_int_ar,label="Integral of Log(C:N)* or $\int$$F_{CN}$*")axs[idx_ax].axhline(y=0,color='k',ls='--',lw=1,alpha=.3)# Subplot <<<# Subplot >>>idx_ax=1+idx_axKpxFcn=pid_kp*cx_logratio_arKdxDerivative=pid_kd*ema_dcxdt_arKixIntegral=pid_ki*cx_int_arprint(f"pid_kp : {pid_kp}\n pid_kd : {pid_kd}\n pid_ki : {pid_ki} ")axs[idx_ax].plot(ts_data.time,KpxFcn,label=r"Kp$\times$ $F_{CN}$*")axs[idx_ax].plot(ts_data.time,KdxDerivative,label=r"Kd $\times$ $dF_{CN}$*")axs[idx_ax].plot(ts_data.time,KixIntegral,label=r"Ki $\times$ $\int$$F_{CN}$*")axs[idx_ax].axhline(y=0,color='k',ls='--',lw=1,alpha=.3)# Subplot <<<# Subplot >>>idx_ax=1+idx_axKpxFcn_prop=np.abs(pid_kp*cx_logratio_ar)/np.sum(np.abs((KpxFcn,KdxDerivative,KixIntegral)),axis=0)KdxDerivative_prop=pid_kd*ema_dcxdt_ar/np.sum(np.abs((KpxFcn,KdxDerivative,KixIntegral)),axis=0)KixIntegral_prop=pid_ki*cx_int_ar/np.sum(np.abs((KpxFcn,KdxDerivative,KixIntegral)),axis=0)print(f"pid_kp : {pid_kp}\n pid_kd : {pid_kd}\n pid_ki : {pid_ki} ")axs[idx_ax].plot(ts_data.time,KpxFcn_prop,label=r"(Kp $\times$ $F_{CN}$)/$\Sigma$|Kp$\times$ $F_{CN}$ + Kd $\times$ $dF_{CN}$ + Ki $\times$ $\int$$F_{CN}$|*")axs[idx_ax].plot(ts_data.time,KdxDerivative_prop,label=r"(Kd $\times$ $dF_{CN}$/$\Sigma$|Kp$\times$ $F_{CN}$ + Kd $\times$ $dF_{CN}$ + Ki $\times$ $\int$$F_{CN}$|*")axs[idx_ax].plot(ts_data.time,KixIntegral_prop,label=r"(Ki $\times$ $\int$$F_{CN}$/$\Sigma$|Kp$\times$ $F_{CN}$ + Kd $\times$ $dF_{CN}$ + Ki $\times$ $\int$$F_{CN}$|*")axs[idx_ax].axhline(y=0,color='k',ls='--',lw=1,alpha=.3)# Subplot <<<# Subplot >>>idx_ax=1+idx_axaxs[idx_ax].plot(ts_data.time,l2fr_delta_ar,label=r"$\Delta$ L2FR*")axs[idx_ax].axhline(y=0,color='k',ls='--',lw=1,alpha=.3)axs[idx_ax].set_ylabel(r"$\Delta$ L2FR*")ax_tmp=axs[idx_ax].twinx()ts_data_l2fr=ds[key]['FATES_L2FR'][idx_start_date:idx_start_date+Yrs_Next*365]ts=ts_data_l2fr.valuesdelta_l2fr_model=ts[1:]-ts[:-1]ax_tmp.plot(ts_data_l2fr.time[:-1],delta_l2fr_model,'r+',label="$\Delta$ FATES_L2FR",alpha=.1,markersize=6)ax_tmp.set_ylabel("$\Delta$ FATES_L2FR",color='r')ax_tmp.legend(loc=5)# Subplot <<<# Subplot >>>idx_ax=1+idx_axL2FR_ar=l2fr_delta_ar.cumsum()axs[idx_ax].plot(ts_data.time,L2FR_ar,label="L2FR*")axs[idx_ax].axhline(y=0,color='k',ls='--',lw=1,alpha=.3)axs[idx_ax].set_ylabel(r"L2FR*")ax_tmp=axs[idx_ax].twinx()ts_data_l2fr=ds[key]['FATES_L2FR'][idx_start_date:idx_start_date+Yrs_Next*365]ax_tmp.plot(ts_data_l2fr.time,ts_data_l2fr,'r+',label="FATES_L2FR",alpha=.1,markersize=6)ax_tmp.set_ylabel(f"{'FATES_L2FR'}({ts_data_l2fr.units})",color='r')foriinrange(len(axs)):axs[i].legend()axs[i].axvline(x=dt.date(1855,1,1),color='g',lw=5,alpha=.1)
df_vals=pd.read_csv(filenames_sim_out[0],header=None,delim_whitespace=True)df_vals.columns=cols# removing the leading string from first columndf_vals['timex']=df_vals['timex'].str.replace('sinkhole','')df_vals['timex']=df_vals['timex'].str.replace('_01:00:00','')df_vals['timex']=pd.to_datetime(df_vals['timex'])# since there are many cohorts per time step, categorizing every cohort with a unique valuedf_vals['timex_code']=pd.factorize(df_vals['timex'])[0]print(f"Sanity Check: {((2007-1850)+1)*365} should be equal to {df_vals['timex_code'].values[-1]+1}")df_vals.head(30)
plt.figure(figsize=(15,6))plt.plot(np.unique(df_vals['timex'],return_counts=1)[0],np.unique(df_vals['timex_code'],return_counts=1)[1])plt.ylabel("Number of cohorts",fontsize=14)plt.axvline(x=dt.date(1855,1,1),color='r',lw=5,alpha=.2)plt.title("Cohorts Vs Time (daily)\n",fontsize=16)
# Set the date range you want to subset betweenstart_date='1854-01-01'end_date='1859-01-01'# Create a boolean mask based on the date rangemask=(df_vals['timex']>=start_date)&(df_vals['timex']<=end_date)# Subset the DataFrame based on the boolean masksubset_df=df_vals[mask]print(f"Maximum Cohorts: {max(np.unique(subset_df['timex_code'],return_counts=1)[1])}")
Maximum Cohorts: 14
In [141]:
# Plotting based on the cohort index# initializing the cohort index to -1subset_df['cohort_id']=-1# for every time step assign a unique index in the cohort idfor_,grpinsubset_df.groupby('timex'):subset_df.loc[grp.index,"cohort_id"]=range(len(grp))
/var/folders/f1/01gxw8vn74q_x_rf_p5ztryjr405zq/T/ipykernel_28297/3966830711.py:4: SettingWithCopyWarning:
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead
See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
/Users/ud4/opt/anaconda3/envs/pyces/lib/python3.9/site-packages/pandas/core/indexing.py:1773: SettingWithCopyWarning:
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead
See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
# Set the date range you want to subset betweenstart_date='1854-01-01'end_date='1894-01-01'# Create a boolean mask based on the date rangemask=(df_vals['timex']>=start_date)&(df_vals['timex']<=end_date)# Subset the DataFrame based on the boolean masksubset_df=df_vals[mask]print(f"Maximum Cohorts: {max(np.unique(subset_df['timex_code'],return_counts=1)[1])}")
# Set the date range you want to subset betweenstart_date='1855-01-02'end_date='1894-01-01'# Create a boolean mask based on the date rangemask=(df_vals['timex']>=start_date)&(df_vals['timex']<=end_date)# Subset the DataFrame based on the boolean masksubset_df=df_vals[mask]print(f"Maximum Cohorts: {max(np.unique(subset_df['timex_code'],return_counts=1)[1])}")
# Checking if the l2fr values (sum over all cohorts) from the print are same as l2fr variable outputmax_vals=5000sum_l2fr=df_vals.groupby('timex_code')['l2fr'].sum().reset_index()[:max_vals]mean_l2fr=df_vals.groupby('timex_code')['l2fr'].mean().reset_index()[:max_vals]min_l2fr=df_vals.groupby('timex_code')['l2fr'].min().reset_index()[:max_vals]max_l2fr=df_vals.groupby('timex_code')['l2fr'].max().reset_index()[:max_vals]
print("\nAll Plots >>>>>>")idx_start_date=5*365# Jan 1, 1855# For multiple variables on a same plot upto 5 years after logging Yrs_Next=5# upto 5 years after logging fig,axs=plt.subplots(9,1,figsize=(15,30),dpi=400)# Subplot 1 >>>idx_ax=0var="FATES_NPP"# C-onlysim_conly=sims.replace("RD","Conly")ifsims=="DUK_PIDE_RD_AgBgW_testing":sim_conly="DUK_PIDE_Conly_AgBgW"key=sim_conlyts_data=ds[key][var][idx_start_date:idx_start_date+Yrs_Next*365]# Daily data range from logging date to next five yearsaxs[idx_ax].scatter(ts_data.time,ts_data,marker=".",label=var+" C-only",alpha=.2)# RDkey=simsts_data=ds[key][var][idx_start_date:idx_start_date+Yrs_Next*365]# Daily data range from logging date to next five yearsaxs[idx_ax].scatter(ts_data.time,ts_data,marker=".",label=var,alpha=.2)axs[idx_ax].set_ylabel(f"{ds[key][var].units}")axs[idx_ax].set_title(f"{key} Daily Logging + 5 years (Scatter) \n",fontsize=16)axs[idx_ax].axhline(y=0,color='r',lw=1,alpha=.1)print(key)# Subplot 1 <<<# Subplot 2 >>>idx_ax=1vars_plot=("""FATES_STOREC_TFFATES_STOREC""").split('\n')vars_plot=vars_plot[1:-1]var=vars_plot[0]ts_data=ds[key][var][idx_start_date:idx_start_date+Yrs_Next*365]# Daily data range from logging date to next five yearsaxs[idx_ax].scatter(ts_data.time,ts_data,marker=".",label=var,alpha=.2)axs[idx_ax].set_ylabel(f"{ds[key][var].units}")ax_tmp=axs[idx_ax].twinx()var=vars_plot[1]ts_data=ds[key][var][idx_start_date:idx_start_date+Yrs_Next*365]# Daily data range from logging date to next five yearsax_tmp.scatter(ts_data.time,ts_data,marker=".",color='k',label=var,alpha=.3)ax_tmp.legend(loc=5)ax_tmp.set_ylabel(f"{var}({ds[key][var].units})")ax_tmp.axhline(y=0,color='r',lw=1,alpha=.1)# Subplot 2 <<<# Subplot 3 >>>idx_ax=2vars_plot=("""FATES_STOREN_TFFATES_STOREN""").split('\n')vars_plot=vars_plot[1:-1]var=vars_plot[0]ts_data=ds[key][var][idx_start_date:idx_start_date+Yrs_Next*365]# Daily data range from logging date to next five yearsaxs[idx_ax].scatter(ts_data.time,ts_data,marker=".",label=var,alpha=.2)axs[idx_ax].set_ylabel(f"{ds[key][var].units}")ax_tmp=axs[idx_ax].twinx()var=vars_plot[1]ts_data=ds[key][var][idx_start_date:idx_start_date+Yrs_Next*365]# Daily data range from logging date to next five yearsax_tmp.scatter(ts_data.time,ts_data,marker=".",color='k',label=var,alpha=.3)ax_tmp.set_ylabel(f"{var}({ds[key][var].units})")ax_tmp.legend(loc=5)ax_tmp.axhline(y=0,color='r',lw=1,alpha=.1)# Subplot 3 <<<# Subplot 4 >>>idx_ax=3ts_data=(ds[key]["FATES_STOREC_TF"]/ds[key]["FATES_STOREN_TF"])[idx_start_date:idx_start_date+Yrs_Next*365]# Daily data range from logging date to next five yearsaxs[idx_ax].scatter(ts_data.time,ts_data,label="FATES_STOREC_TF/FATES_STOREN_TF",alpha=.2)axs[idx_ax].set_ylabel(f"FATES_STOREC_TF/FATES_STOREN_TF")ax_tmp=axs[idx_ax].twinx()ts_data=(ds[key]["FATES_STOREC"]/ds[key]["FATES_STOREN"])[idx_start_date:idx_start_date+Yrs_Next*365]# Daily data range from logging date to next five yearsax_tmp.scatter(ts_data.time,ts_data,marker=".",color='k',label="FATES_STOREC/FATES_STOREN",alpha=.2)ax_tmp.legend(loc=5)ax_tmp.set_ylabel(f"FATES_STOREC/FATES_STOREN")# Subplot 4 <<<# Subplot 5 >>>idx_ax=4ts_data=(ds[key]["FATES_FROOT_ALLOC"]/ds[key]["FATES_LEAF_ALLOC"])[idx_start_date:idx_start_date+Yrs_Next*365]# Daily data range from logging date to next five yearsaxs[idx_ax].scatter(ts_data.time,ts_data,label="FATES_FROOT_ALLOC/FATES_LEAF_ALLOC",alpha=.2)axs[idx_ax].set_ylabel(f"FATES_FROOT_ALLOC/FATES_LEAF_ALLOC")ax_tmp=axs[idx_ax].twinx()ts_data=ds[key]["FATES_L2FR"][idx_start_date:idx_start_date+Yrs_Next*365]# Daily data range from logging date to next five yearsax_tmp.scatter(ts_data.time,ts_data,marker=".",color='k',label="FATES_L2FR",alpha=.2)ax_tmp.set_ylabel(f"{ds[key]['FATES_L2FR'].units}")ax_tmp.legend(loc=5)# Subplot 5 <<<# Subplot 6 >>>idx_ax=5vars_plot=("""FATES_LEAFCFATES_FROOTC""").split('\n')vars_plot=vars_plot[1:-1]fori_var,varinenumerate(vars_plot):ts_data=ds[key][var][idx_start_date:idx_start_date+Yrs_Next*365]# Daily data range from logging date to next five yearsaxs[idx_ax].scatter(ts_data.time,ts_data,label=var,alpha=.2)axs[idx_ax].set_ylabel(f"{ds[key][var].units}")ax_tmp=axs[idx_ax].twinx()ts_data=(ds[key]["FATES_FROOTC"]/ds[key]["FATES_LEAFC"])[idx_start_date:idx_start_date+Yrs_Next*365]# Daily data range from logging date to next five yearsax_tmp.scatter(ts_data.time,ts_data,marker=".",color='k',label="FATES_FROOTC/FATES_LEAFC",alpha=.2)ax_tmp.legend(loc=5)ax_tmp.set_ylabel(f"Unitless")# Subplot 6 <<<# Subplot 7 >>>idx_ax=6vars_plot=("""FATES_NO3UPTAKEFATES_NH4UPTAKE""").split('\n')vars_plot=vars_plot[1:-1]var=vars_plot[0]ts_data=ds[key][var][idx_start_date:idx_start_date+Yrs_Next*365]# Daily data range from logging date to next five yearsaxs[idx_ax].scatter(ts_data.time,ts_data,label=var,alpha=.2)axs[idx_ax].set_ylabel(f"{ds[key][var].units}")axs[idx_ax].axhline(y=0,color='r',lw=1,alpha=.1)ax_tmp=axs[idx_ax].twinx()var=vars_plot[1]ts_data=ds[key][var][idx_start_date:idx_start_date+Yrs_Next*365]# Daily data range from logging date to next five yearsax_tmp.scatter(ts_data.time,ts_data,marker=".",color='k',label=var,alpha=.2)ax_tmp.set_ylabel(f"{var}({ds[key][var].units})")ax_tmp.legend(loc=5)# Subplot 7 <<<# Subplot 8 >>>idx_ax=7vars_plot=("""GROSS_NMINNET_NMIN""").split('\n')vars_plot=vars_plot[1:-1]fori_var,varinenumerate(vars_plot):ts_data=ds[key][var][idx_start_date:idx_start_date+Yrs_Next*365]# Daily data range from logging date to next five yearsaxs[idx_ax].scatter(ts_data.time,ts_data,label=var,alpha=.2)axs[idx_ax].set_ylabel(f"{ds[key][var].units}")axs[idx_ax].axhline(y=0,color='r',lw=1,alpha=.1)# Subplot 8 <<<# Subplot 9 >>>idx_ax=8vars_plot=("""SMIN_NO3_LEACHEDDENIT""").split('\n')vars_plot=vars_plot[1:-1]fori_var,varinenumerate(vars_plot):ts_data=ds[key][var][idx_start_date:idx_start_date+Yrs_Next*365]# Daily data range from logging date to next five yearsaxs[idx_ax].scatter(ts_data.time,ts_data,label=var,alpha=.2)axs[idx_ax].set_ylabel(f"{ds[key][var].units}")axs[idx_ax].axhline(y=0,color='r',lw=1,alpha=.1)# Subplot 9 <<<print(key)foriinrange(len(axs)):axs[i].legend()axs[i].axvline(x=dt.date(1855,1,1),color='g',lw=5,alpha=.1)#fig.savefig('./Plots/' + f'{key}_Scatter_daily.pdf',bbox_inches='tight')#fig.savefig('./Plots/' + f'{key}_Scatter_daily.png',bbox_inches='tight')
All Plots >>>>>>
DUK_PIDC_RD_AgBgW_CN53
/Users/ud4/opt/anaconda3/envs/pyces/lib/python3.9/site-packages/dask/core.py:121: RuntimeWarning:
invalid value encountered in true_divide
print("\nFigure Set 2 >>>>>>")# For multiple variables on a same plot#sims = "ORN_PIDE_RD_AgBgW" # Difined at the topkey=simsvars_plot=("""FATES_GPPFATES_AUTORESPFATES_NPP""").split('\n')vars_plot=vars_plot[1:-1]ymin=9e20ymax=-9e20sum_ts=0fori_var,varinenumerate(vars_plot):ts_data=ds[key][var].groupby("time.year").mean('time')sum_ts=sum_ts+ts_datats_data.plot(figsize=(20,3))ifnp.min(ts_data.values)<ymin:ymin=np.min(ts_data.values)ifnp.max(ts_data.values)>ymax:ymax=np.max(ts_data.values)plt.title(f"{ds[key][var].long_name} ({var}) - {key} - AnnualSUM | Units: {ds[key][var].units}")plt.axvline(x=logging_year,color='r',lw=5,label='logging year',alpha=.2)#plt.axvline(x = 1996, color = 'g',lw=5, label = 'logging year', alpha=.2)#plt.ylim(.2e-6,1.5e-6)ifi_var!=len(vars_plot)-1:plt.xticks([])plt.xlabel(None)fig=plt.figure(figsize=(20,9))plt.title(f"Common Plot for simulation {key}",fontsize=15)fori_var,varinenumerate(vars_plot):ts_data=ds[key][var].groupby("time.year").mean('time')plt.plot(ts_data.year,ts_data,label=var)plt.ylim(ymin*.95,ymax*1.05)plt.legend(fontsize=14)plt.axvline(x=logging_year,color='r',lw=5,label='logging year',alpha=.2)fig1=plt.figure(figsize=(20,9))plt.title(f"Fractional Plot for simulation {key}",fontsize=15)fori_var,varinenumerate(vars_plot):ts_data=ds[key][var].groupby("time.year").mean('time')plt.plot(ts_data.year,ts_data/sum_ts,label=var)#plt.ylim(ymin*.95,ymax*1.05)plt.legend(fontsize=14)plt.axvline(x=logging_year,color='r',lw=5,label='logging year',alpha=.2)
print("\nFigure Set 3 >>>>>>")# For multiple variables on a same plot#sims = "ORN_PIDE_RD_AgBgW" # Difined at the topkey=simsvars_plot=("""FATES_STORECFATES_STOREC_TF""").split('\n')vars_plot=vars_plot[1:-1]ymin=9e20ymax=-9e20sum_ts=0fori_var,varinenumerate(vars_plot):ts_data=ds[key][var].groupby("time.year").mean('time')sum_ts=sum_ts+ts_datats_data.plot(figsize=(20,3))ifnp.min(ts_data.values)<ymin:ymin=np.min(ts_data.values)ifnp.max(ts_data.values)>ymax:ymax=np.max(ts_data.values)plt.title(f"{ds[key][var].long_name} ({var}) - {key} - AnnualSUM | Units: {ds[key][var].units}")plt.axvline(x=logging_year,color='r',lw=5,label='logging year',alpha=.2)#plt.axvline(x = 1996, color = 'g',lw=5, label = 'logging year', alpha=.2)#plt.ylim(.2e-6,1.5e-6)ifi_var!=len(vars_plot)-1:plt.xticks([])plt.xlabel(None)fig=plt.figure(figsize=(20,9))plt.title(f"Common Plot for simulation {key}",fontsize=15)fori_var,varinenumerate(vars_plot):ifvar=="FATES_L2FR":continuets_data=ds[key][var].groupby("time.year").mean('time')plt.plot(ts_data.year,ts_data,label=var)plt.ylim(ymin*.95,ymax*1.05)plt.legend(fontsize=14)plt.axvline(x=logging_year,color='r',lw=5,label='logging year',alpha=.2)fig1=plt.figure(figsize=(20,9))plt.title(f"Fractional Plot for simulation {key}",fontsize=15)fori_var,varinenumerate(vars_plot):ts_data=ds[key][var].groupby("time.year").mean('time')plt.plot(ts_data.year,ts_data/sum_ts,label=var)#plt.ylim(ymin*.95,ymax*1.05)plt.legend(fontsize=14)plt.axvline(x=logging_year,color='r',lw=5,label='logging year',alpha=.2)fig2=plt.figure(figsize=(20,9))plt.title(f"Ratio Plot for simulation {key}",fontsize=15)#for i_var,var in enumerate(vars_plot):ifTrue:ts_data=(ds[key]["FATES_STOREC"]/ds[key]["FATES_STOREC_TF"]).groupby("time.year").mean('time')plt.plot(ts_data.year,ts_data,label="FATES_STOREC/FATES_STOREC_TF")#plt.ylim(ymin*.95,ymax*1.05)plt.legend(fontsize=14)plt.title(f"{key}: FATES_STOREC/FATES_STOREC_TF")plt.axvline(x=logging_year,color='r',lw=5,label='logging year',alpha=.2)
print("\nFigure Set 4 >>>>>>")# For multiple variables on a same plot#sims = "ORN_PIDE_RD_AgBgW" # Difined at the topkey=simsvars_plot=("""FATES_STORENFATES_STOREN_TF""").split('\n')vars_plot=vars_plot[1:-1]ymin=9e20ymax=-9e20sum_ts=0fori_var,varinenumerate(vars_plot):ts_data=ds[key][var].groupby("time.year").mean('time')sum_ts=sum_ts+ts_datats_data.plot(figsize=(20,3))ifnp.min(ts_data.values)<ymin:ymin=np.min(ts_data.values)ifnp.max(ts_data.values)>ymax:ymax=np.max(ts_data.values)plt.title(f"{ds[key][var].long_name} ({var}) - {key} - AnnualSUM | Units: {ds[key][var].units}")plt.axvline(x=logging_year,color='r',lw=5,label='logging year',alpha=.2)#plt.axvline(x = 1996, color = 'g',lw=5, label = 'logging year', alpha=.2)#plt.ylim(.2e-6,1.5e-6)ifi_var!=len(vars_plot)-1:plt.xticks([])plt.xlabel(None)fig=plt.figure(figsize=(20,9))plt.title(f"Common Plot for simulation {key}",fontsize=15)fori_var,varinenumerate(vars_plot):ifvar=="FATES_L2FR":continuets_data=ds[key][var].groupby("time.year").mean('time')plt.plot(ts_data.year,ts_data,label=var)plt.ylim(ymin*.95,ymax*1.05)plt.legend(fontsize=14)plt.axvline(x=logging_year,color='r',lw=5,label='logging year',alpha=.2)fig1=plt.figure(figsize=(20,9))plt.title(f"Fractional Plot for simulation {key}",fontsize=15)fori_var,varinenumerate(vars_plot):ts_data=ds[key][var].groupby("time.year").mean('time')plt.plot(ts_data.year,ts_data/sum_ts,label=var)#plt.ylim(ymin*.95,ymax*1.05)plt.legend(fontsize=14)plt.axvline(x=logging_year,color='r',lw=5,label='logging year',alpha=.2)fig2=plt.figure(figsize=(20,9))plt.title(f"Ratio Plot for simulation {key}",fontsize=15)#for i_var,var in enumerate(vars_plot):ifTrue:ts_data=(ds[key]["FATES_STOREN"]/ds[key]["FATES_STOREN_TF"]).groupby("time.year").mean('time')plt.plot(ts_data.year,ts_data,label="FATES_STOREN/FATES_STOREN_TF")plt.axvline(x=logging_year,color='r',lw=5,label='logging year',alpha=.2)#plt.ylim(ymin*.95,ymax*1.05)plt.legend(fontsize=14)plt.title(f"{key}: FATES_STOREN/FATES_STOREN_TF")
print("\nFigure Set 5 >>>>>>")# For multiple variables on a same plot#sims = "ORN_PIDE_RD_AgBgW" # Difined at the topkey=simsvars_plot=("""FATES_STOREC_TFFATES_STOREN_TF""").split('\n')vars_plot=vars_plot[1:-1]ymin=9e20ymax=-9e20sum_ts=0fori_var,varinenumerate(vars_plot):ts_data=ds[key][var].groupby("time.year").mean('time')sum_ts=sum_ts+ts_datats_data.plot(figsize=(20,3))ifnp.min(ts_data.values)<ymin:ymin=np.min(ts_data.values)ifnp.max(ts_data.values)>ymax:ymax=np.max(ts_data.values)plt.title(f"{ds[key][var].long_name} ({var}) - {key} - AnnualSUM | Units: {ds[key][var].units}")plt.axvline(x=logging_year,color='r',lw=5,label='logging year',alpha=.2)#plt.axvline(x = 1996, color = 'g',lw=5, label = 'logging year', alpha=.2)#plt.ylim(.2e-6,1.5e-6)ifi_var!=len(vars_plot)-1:plt.xticks([])plt.xlabel(None)fig=plt.figure(figsize=(20,9))plt.title(f"Common Plot for simulation {key}",fontsize=15)fori_var,varinenumerate(vars_plot):ifvar=="FATES_L2FR":continuets_data=ds[key][var].groupby("time.year").mean('time')plt.plot(ts_data.year,ts_data,label=var)plt.ylim(ymin*.95,ymax*1.05)plt.legend(fontsize=14)plt.axvline(x=logging_year,color='r',lw=5,label='logging year',alpha=.2)fig1=plt.figure(figsize=(20,9))plt.title(f"Fractional Plot for simulation {key}",fontsize=15)fori_var,varinenumerate(vars_plot):ts_data=ds[key][var].groupby("time.year").mean('time')plt.plot(ts_data.year,ts_data/sum_ts,label=var)#plt.ylim(ymin*.95,ymax*1.05)plt.legend(fontsize=14)plt.axvline(x=logging_year,color='r',lw=5,label='logging year',alpha=.2)fig2=plt.figure(figsize=(20,9))plt.title(f"Ratio Plot for simulation {key}",fontsize=15)#for i_var,var in enumerate(vars_plot):ifTrue:ts_data=(ds[key]["FATES_STOREC_TF"]/ds[key]["FATES_STOREN_TF"]).groupby("time.year").mean('time')plt.plot(ts_data.year,ts_data,label="FATES_STOREC_TF/FATES_STOREN_TF")plt.axvline(x=logging_year,color='r',lw=5,label='logging year',alpha=.2)#plt.ylim(ymin*.95,ymax*1.05)plt.legend(fontsize=14)plt.title(f"{key}: FATES_STOREC_TF/FATES_STOREN_TF")fig3=plt.figure(figsize=(20,9))plt.title(f"Ratio Plot for simulation {key}",fontsize=15)#for i_var,var in enumerate(vars_plot):ifTrue:ts_data=(ds[key]["FATES_STOREC"]/ds[key]["FATES_STOREN"]).groupby("time.year").mean('time')plt.plot(ts_data.year,ts_data,label="FATES_STOREC/FATES_STOREN")plt.axvline(x=logging_year,color='r',lw=5,label='logging year',alpha=.2)#plt.ylim(ymin*.95,ymax*1.05)plt.legend(fontsize=14)plt.title(f"{key}: FATES_STOREC/FATES_STOREN")
print("\nFigure Set 6 >>>>>>")# For multiple variables on a same plot#sims = "ORN_PIDE_RD_AgBgW" # Difined at the topkey=sims'''vars_plot = ["FATES_CROOT_ALLOC","FATES_FROOT_ALLOC","FATES_FROOT_ALLOC","FATES_SEED_ALLOC","FATES_STEM_ALLOC","FATES_STORE_ALLOC"]'''vars_plot=("""FATES_CROOT_ALLOCFATES_FROOT_ALLOCFATES_FROOT_ALLOCFATES_SEED_ALLOCFATES_STEM_ALLOCFATES_STORE_ALLOC""").split('\n')vars_plot=vars_plot[1:-1]ymin=9e20ymax=-9e20sum_ts=0fori_var,varinenumerate(vars_plot):ts_data=ds[key][var].groupby("time.year").mean('time')sum_ts=sum_ts+ts_datats_data.plot(figsize=(20,3))ifnp.min(ts_data.values)<ymin:ymin=np.min(ts_data.values)ifnp.max(ts_data.values)>ymax:ymax=np.max(ts_data.values)plt.title(f"{ds[key][var].long_name} ({var}) - {key} - AnnualSUM | Units: {ds[key][var].units}")plt.axvline(x=logging_year,color='r',lw=5,label='logging year',alpha=.2)#plt.axvline(x = 1996, color = 'g',lw=5, label = 'logging year', alpha=.2)#plt.ylim(.2e-6,1.5e-6)ifi_var!=len(vars_plot)-1:plt.xticks([])plt.xlabel(None)fig=plt.figure(figsize=(20,9))plt.title(f"Common Plot for simulation {key}",fontsize=15)fori_var,varinenumerate(vars_plot):ts_data=ds[key][var].groupby("time.year").mean('time')plt.plot(ts_data.year,ts_data,label=var)plt.ylim(ymin*.95,ymax*1.05)plt.legend(fontsize=14)plt.axvline(x=logging_year,color='r',lw=5,label='logging year',alpha=.2)fig1=plt.figure(figsize=(20,9))plt.title(f"Fractional Plot for simulation {key}",fontsize=15)fori_var,varinenumerate(vars_plot):ts_data=ds[key][var].groupby("time.year").mean('time')plt.plot(ts_data.year,ts_data/sum_ts,label=var)#plt.ylim(ymin*.95,ymax*1.05)plt.legend(fontsize=14)plt.axvline(x=logging_year,color='r',lw=5,label='logging year',alpha=.2)fig2=plt.figure(figsize=(20,9))plt.title(f"Ratio Plot for simulation {key}",fontsize=15)#for i_var,var in enumerate(vars_plot):ifTrue:ts_data=(ds[key]["FATES_FROOT_ALLOC"]/ds[key]["FATES_LEAF_ALLOC"]).groupby("time.year").mean('time')plt.plot(ts_data.year,ts_data,label="FATES_FROOT_ALLOC/FATES_LEAF_ALLOC")ts_data=ds[key]["FATES_L2FR"].groupby("time.year").mean('time')plt.plot(ts_data.year,ts_data,label="FATES_L2FR")#plt.ylim(ymin*.95,ymax*1.05)plt.legend(fontsize=14)plt.title(f"{key}: FATES_FROOT_ALLOC/FATES_LEAF_ALLOC")plt.axvline(x=logging_year,color='r',lw=5,label='logging year',alpha=.2)
Figure Set 6 >>>>>>
/Users/ud4/opt/anaconda3/envs/pyces/lib/python3.9/site-packages/dask/core.py:121: RuntimeWarning:
invalid value encountered in true_divide
print("\nFigure Set 7 >>>>>>")# For multiple variables on a same plot#sims = "ORN_PIDE_RD_AgBgW" # Difined at the topkey=simsvars_plot=("""FATES_LEAFCFATES_FROOTC""").split('\n')vars_plot=vars_plot[1:-1]ymin=9e20ymax=-9e20sum_ts=0fori_var,varinenumerate(vars_plot):ts_data=ds[key][var].groupby("time.year").mean('time')sum_ts=sum_ts+ts_datats_data.plot(figsize=(20,3))ifnp.min(ts_data.values)<ymin:ymin=np.min(ts_data.values)ifnp.max(ts_data.values)>ymax:ymax=np.max(ts_data.values)plt.title(f"{ds[key][var].long_name} ({var}) - {key} - AnnualSUM | Units: {ds[key][var].units}")plt.axvline(x=logging_year,color='r',lw=5,label='logging year',alpha=.2)#plt.axvline(x = 1996, color = 'g',lw=5, label = 'logging year', alpha=.2)#plt.ylim(.2e-6,1.5e-6)ifi_var!=len(vars_plot)-1:plt.xticks([])plt.xlabel(None)fig=plt.figure(figsize=(20,9))plt.title(f"Common Plot for simulation {key}",fontsize=15)fori_var,varinenumerate(vars_plot):ifvar=="FATES_L2FR":continuets_data=ds[key][var].groupby("time.year").mean('time')plt.plot(ts_data.year,ts_data,label=var)plt.ylim(ymin*.95,ymax*1.05)plt.legend(fontsize=14)plt.axvline(x=logging_year,color='r',lw=5,label='logging year',alpha=.2)fig1=plt.figure(figsize=(20,9))plt.title(f"Fractional Plot for simulation {key}",fontsize=15)fori_var,varinenumerate(vars_plot):ts_data=ds[key][var].groupby("time.year").mean('time')plt.plot(ts_data.year,ts_data/sum_ts,label=var)#plt.ylim(ymin*.95,ymax*1.05)plt.legend(fontsize=14)plt.axvline(x=logging_year,color='r',lw=5,label='logging year',alpha=.2)fig2=plt.figure(figsize=(20,9))plt.title(f"Ratio Plot for simulation {key}",fontsize=15)#for i_var,var in enumerate(vars_plot):ifTrue:ts_data=(ds[key]["FATES_FROOTC"]/ds[key]["FATES_LEAFC"]).groupby("time.year").mean('time')plt.plot(ts_data.year,ts_data,label="FATES_FROOTC/FATES_LEAFC")ts_data=ds[key]["FATES_L2FR"].groupby("time.year").mean('time')plt.plot(ts_data.year,ts_data,label="FATES_L2FR")#plt.ylim(ymin*.95,ymax*1.05)plt.legend(fontsize=14)plt.title(f"{key}: FATES_FROOTC/FATES_LEAFC")plt.axvline(x=logging_year,color='r',lw=5,label='logging year',alpha=.2)
print("\nFigure Set 8 >>>>>>")# For multiple variables on a same plot#sims = "ORN_PIDE_RD_AgBgW" # Difined at the topkey=simsvars_plot=("""FATES_NO3UPTAKEFATES_NH4UPTAKE""").split('\n')vars_plot=vars_plot[1:-1]ymin=9e20ymax=-9e20sum_ts=0fori_var,varinenumerate(vars_plot):ts_data=ds[key][var].groupby("time.year").mean('time')sum_ts=sum_ts+ts_datats_data.plot(figsize=(20,3))ifnp.min(ts_data.values)<ymin:ymin=np.min(ts_data.values)ifnp.max(ts_data.values)>ymax:ymax=np.max(ts_data.values)plt.title(f"{ds[key][var].long_name} ({var}) - {key} - AnnualSUM | Units: {ds[key][var].units}")plt.axvline(x=logging_year,color='r',lw=5,label='logging year',alpha=.2)#plt.axvline(x = 1996, color = 'g',lw=5, label = 'logging year', alpha=.2)#plt.ylim(.2e-6,1.5e-6)ifi_var!=len(vars_plot)-1:plt.xticks([])plt.xlabel(None)fig=plt.figure(figsize=(20,9))plt.title(f"Common Plot for simulation {key}",fontsize=15)fori_var,varinenumerate(vars_plot):ifvar=="FATES_L2FR":continuets_data=ds[key][var].groupby("time.year").mean('time')plt.plot(ts_data.year,ts_data,label=var)plt.ylim(ymin*.95,ymax*1.05)plt.legend(fontsize=14)plt.axvline(x=logging_year,color='r',lw=5,label='logging year',alpha=.2)fig1=plt.figure(figsize=(20,9))plt.title(f"Fractional Plot for simulation {key}",fontsize=15)fori_var,varinenumerate(vars_plot):ts_data=ds[key][var].groupby("time.year").mean('time')plt.plot(ts_data.year,ts_data/sum_ts,label=var)#plt.ylim(ymin*.95,ymax*1.05)plt.legend(fontsize=14)plt.axvline(x=logging_year,color='r',lw=5,label='logging year',alpha=.2)
Figure Set 8 >>>>>>
Out[169]:
<matplotlib.lines.Line2D at 0x188d94190>
In [170]:
### *Figure Set 9*
In [171]:
print("\nFigure Set 9 >>>>>>")# For multiple variables on a same plot#sims = "ORN_PIDE_RD_AgBgW" # Difined at the topkey=simsvars_plot=("""GROSS_NMINNET_NMIN""").split('\n')vars_plot=vars_plot[1:-1]ymin=9e20ymax=-9e20sum_ts=0fori_var,varinenumerate(vars_plot):ts_data=ds[key][var].groupby("time.year").mean('time')sum_ts=sum_ts+ts_datats_data.plot(figsize=(20,3))ifnp.min(ts_data.values)<ymin:ymin=np.min(ts_data.values)ifnp.max(ts_data.values)>ymax:ymax=np.max(ts_data.values)plt.title(f"{ds[key][var].long_name} ({var}) - {key} - AnnualSUM | Units: {ds[key][var].units}")plt.axvline(x=logging_year,color='r',lw=5,label='logging year',alpha=.2)#plt.axvline(x = 1996, color = 'g',lw=5, label = 'logging year', alpha=.2)#plt.ylim(.2e-6,1.5e-6)ifi_var!=len(vars_plot)-1:plt.xticks([])plt.xlabel(None)fig=plt.figure(figsize=(20,9))plt.title(f"Common Plot for simulation {key}",fontsize=15)fori_var,varinenumerate(vars_plot):ifvar=="FATES_L2FR":continuets_data=ds[key][var].groupby("time.year").mean('time')plt.plot(ts_data.year,ts_data,label=var)plt.ylim(ymin*.95,ymax*1.05)plt.legend(fontsize=14)plt.axvline(x=logging_year,color='r',lw=5,label='logging year',alpha=.2)fig1=plt.figure(figsize=(20,9))plt.title(f"Fractional Plot for simulation {key}",fontsize=15)fori_var,varinenumerate(vars_plot):ts_data=ds[key][var].groupby("time.year").mean('time')plt.plot(ts_data.year,ts_data/sum_ts,label=var)#plt.ylim(ymin*.95,ymax*1.05)plt.legend(fontsize=14)plt.axvline(x=logging_year,color='r',lw=5,label='logging year',alpha=.2)
print("\nFigure Set 10 >>>>>>")# For multiple variables on a same plot#sims = "ORN_PIDE_RD_AgBgW" # Difined at the topkey=simsvars_plot=("""SMIN_NO3_LEACHEDDENIT""").split('\n')vars_plot=vars_plot[1:-1]ymin=9e20ymax=-9e20sum_ts=0fori_var,varinenumerate(vars_plot):ts_data=ds[key][var].groupby("time.year").mean('time')sum_ts=sum_ts+ts_datats_data.plot(figsize=(20,3))ifnp.min(ts_data.values)<ymin:ymin=np.min(ts_data.values)ifnp.max(ts_data.values)>ymax:ymax=np.max(ts_data.values)plt.title(f"{ds[key][var].long_name} ({var}) - {key} - AnnualSUM | Units: {ds[key][var].units}")plt.axvline(x=logging_year,color='r',lw=5,label='logging year',alpha=.2)#plt.axvline(x = 1996, color = 'g',lw=5, label = 'logging year', alpha=.2)#plt.ylim(.2e-6,1.5e-6)ifi_var!=len(vars_plot)-1:plt.xticks([])plt.xlabel(None)fig=plt.figure(figsize=(20,9))plt.title(f"Common Plot for simulation {key}",fontsize=15)fori_var,varinenumerate(vars_plot):ifvar=="FATES_L2FR":continuets_data=ds[key][var].groupby("time.year").mean('time')plt.plot(ts_data.year,ts_data,label=var)plt.ylim(ymin*.95,ymax*1.05)plt.legend(fontsize=14)plt.axvline(x=logging_year,color='r',lw=5,label='logging year',alpha=.2)fig1=plt.figure(figsize=(20,9))plt.title(f"Fractional Plot for simulation {key}",fontsize=15)fori_var,varinenumerate(vars_plot):ts_data=ds[key][var].groupby("time.year").mean('time')plt.plot(ts_data.year,ts_data/sum_ts,label=var)#plt.ylim(ymin*.95,ymax*1.05)plt.legend(fontsize=14)plt.axvline(x=logging_year,color='r',lw=5,label='logging year',alpha=.2)
Figure Set 10 >>>>>>
Out[172]:
<matplotlib.lines.Line2D at 0x195238970>
In [ ]:
In [ ]:
In [173]:
# For multiple variables on a same plot#sims = "ORN_PIDE_RD_AgBgW" # Difined at the topkey=simsvars_plot=("""FATES_FROOTCFATES_LEAFCFATES_NONSTRUCTCFATES_REPROCFATES_SAPWOODCFATES_STORECFATES_STRUCTCFATES_VEGC""").split('\n')vars_plot=vars_plot[1:-1]ymin=9e20ymax=-9e20sum_ts=0fori_var,varinenumerate(vars_plot):print(var)ts_data=ds[key][var].groupby("time.year").mean('time')sum_ts=sum_ts+ts_datats_data.plot(figsize=(20,3))ifnp.min(ts_data.values)<ymin:ymin=np.min(ts_data.values)ifnp.max(ts_data.values)>ymax:ymax=np.max(ts_data.values)plt.title(f"{ds[key][var].long_name} ({var}) - {key} - AnnualSUM | Units: {ds[key][var].units}")plt.axvline(x=logging_year,color='r',lw=5,label='logging year',alpha=.2)#plt.axvline(x = 1996, color = 'g',lw=5, label = 'logging year', alpha=.2)#plt.ylim(.2e-6,1.5e-6)ifi_var!=len(vars_plot)-1:plt.xticks([])plt.xlabel(None)fig=plt.figure(figsize=(20,9))plt.title(f"Common Plot for simulation {key}",fontsize=15)fori_var,varinenumerate(vars_plot):ts_data=ds[key][var].groupby("time.year").mean('time')plt.plot(ts_data.year,ts_data,label=var)plt.ylim(ymin*.95,ymax*1.05)plt.legend(fontsize=14)plt.axvline(x=logging_year,color='r',lw=5,label='logging year',alpha=.2)fig1=plt.figure(figsize=(20,9))plt.title(f"Fractional Plot for simulation {key}",fontsize=15)fori_var,varinenumerate(vars_plot):ts_data=ds[key][var].groupby("time.year").mean('time')plt.plot(ts_data.year,ts_data/sum_ts,label=var)#plt.ylim(ymin*.95,ymax*1.05)plt.legend(fontsize=14)plt.axvline(x=logging_year,color='r',lw=5,label='logging year',alpha=.2)
FATES_NPP (for C-only at that site and RD/ECA) [done]
FATES_NPP [done]
FATES_STOREC
FATES_STOREC_TF
FATES_STOREN
FATES_STOREN_TF
FATES_STOREC_TF / FATES_STOREN_TF
FATES_L2FR
FATES_FROOTC_ALLOC / FATES_LEAFC_ALLOC
FATES_FROOTC / FATES_LEAFC
FATES_FROOTC
FATES_LEAFC
FATES_NO3UPTAKE
FATES_NH4UPTAKE
NET_NMIN
GROSS_NMIN
It's also be good to look at N leaching (NLEACH?) and denitrification (DENIT?). Not sure those are the variable names but the variables (however named) should be in the output files I hope.
And let's take a look at PID-D for both Duke and ORNL, and PID-B for Duke.
