pwd

'/Users/ud4/repos/GitHub/FATESFACE'

from IPython.display import Image, display
display(Image(filename='/Users/ud4/Downloads/PIDs.png',width=800))
display(Image(filename='/Users/ud4/Downloads/BiomassVsPIDs.png',width=800))
print ("Source : Knox 2023 (under review; https://d197for5662m48.cloudfront.net/documents/publicationstatus/129549/preprint_pdf/4892791c8be5504aff47afd6065937b7.pdf)")

Source : Knox 2023 (under review; https://d197for5662m48.cloudfront.net/documents/publicationstatus/129549/preprint_pdf/4892791c8be5504aff47afd6065937b7.pdf)

ALL PID cases

PID	Kp	Kd
A	0.0005	0.1
B	0.0005	0.005
C	0.0001	0.01
D	0.001	0.1
E	0.001	0.5
F	0.005	0.1
G	0.005	0.5
H	0.001	1.0

Index

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

Features of this Jupyter Notebook

• for a cell, you can choose which variables you are interested in for a single simulation.
• some plots from postprocessor

import os,glob
import xarray as xr
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import pyreadr # to read .rds files
import datetime as dt

np.set_printoptions(suppress=True)

Define FilePaths

path_in = "/Users/ud4/FATESMDS_analysis/outputs/runs/tests_alp/230309/"

# RD only

fnames={}
fnames["DUK_PIDA_Conly_AgBgW"] = f"{path_in}Bharat_AW_Nalloc_api25e3sm_mf0df100_r0615_Base_PIDA_AgBgW_processed/Bharat_AW_Nalloc_api25e3sm_mf0df100_r0615_Base_PIDA_AgBgW_US-DUK_trans.nc"
fnames["ORN_PIDA_Conly_AgBgW"] = f"{path_in}Bharat_AW_Nalloc_api25e3sm_mf0df100_r0615_Base_PIDA_AgBgW_processed/Bharat_AW_Nalloc_api25e3sm_mf0df100_r0615_Base_PIDA_AgBgW_US-ORN_trans.nc"

fnames["DUK_PIDB_Conly"] = f"{path_in}Bharat_AW_Nalloc_api25e3sm_mf0df100_r0614_Base_PIDB_AgBgW_processed/Bharat_AW_Nalloc_api25e3sm_mf0df100_r0614_Base_PIDB_AgBgW_US-DUK_trans.nc"
fnames["ORN_PIDB_Conly"] = f"{path_in}Bharat_AW_Nalloc_api25e3sm_mf0df100_r0614_Base_PIDB_AgBgW_processed/Bharat_AW_Nalloc_api25e3sm_mf0df100_r0614_Base_PIDB_AgBgW_US-ORN_trans.nc"

fnames["DUK_PIDC_Conly"] = f"{path_in}Bharat_AW_Nalloc_api25e3sm_mf0df100_r0614_Base_PIDC_AgBgW_processed/Bharat_AW_Nalloc_api25e3sm_mf0df100_r0614_Base_PIDC_AgBgW_US-DUK_trans.nc"
fnames["ORN_PIDC_Conly"] = f"{path_in}Bharat_AW_Nalloc_api25e3sm_mf0df100_r0614_Base_PIDC_AgBgW_processed/Bharat_AW_Nalloc_api25e3sm_mf0df100_r0614_Base_PIDC_AgBgW_US-ORN_trans.nc"

fnames["DUK_PIDD_Conly"] = f"{path_in}Bharat_AW_Nalloc_api25e3sm_mf0df100_r0614_Base_PIDD_AgBgW_processed/Bharat_AW_Nalloc_api25e3sm_mf0df100_r0614_Base_PIDD_AgBgW_US-DUK_trans.nc"
fnames["ORN_PIDD_Conly"] = f"{path_in}Bharat_AW_Nalloc_api25e3sm_mf0df100_r0614_Base_PIDD_AgBgW_processed/Bharat_AW_Nalloc_api25e3sm_mf0df100_r0614_Base_PIDD_AgBgW_US-ORN_trans.nc"

fnames["DUK_PIDA_RD_AgBgW"] = f"{path_in}Bharat_AW_Nalloc_api25e3sm_mf0df100_r0615_Base_PIDA_AgBgW_RD_processed/Bharat_AW_Nalloc_api25e3sm_mf0df100_r0615_Base_PIDA_AgBgW_RD_US-DUK_trans.nc"
fnames["ORN_PIDA_RD_AgBgW"] = f"{path_in}Bharat_AW_Nalloc_api25e3sm_mf0df100_r0615_Base_PIDA_AgBgW_RD_processed/Bharat_AW_Nalloc_api25e3sm_mf0df100_r0615_Base_PIDA_AgBgW_RD_US-ORN_trans.nc"

fnames["DUK_PIDB_RD_AgBgW"] = f"{path_in}Bharat_AW_Nalloc_api25e3sm_mf0df100_r0614_Base_PIDB_AgBgW_RD_processed/Bharat_AW_Nalloc_api25e3sm_mf0df100_r0614_Base_PIDB_AgBgW_RD_US-DUK_trans.nc"
fnames["ORN_PIDB_RD_AgBgW"] = f"{path_in}Bharat_AW_Nalloc_api25e3sm_mf0df100_r0614_Base_PIDB_AgBgW_RD_processed/Bharat_AW_Nalloc_api25e3sm_mf0df100_r0614_Base_PIDB_AgBgW_RD_US-ORN_trans.nc"

#"ORN_PIDC_RD" : Simulation missing *Running
fnames["DUK_PIDC_Conly_AgBgW"] = f"{path_in}Bharat_AW_Nalloc_api25e3sm_mf0df100_r0614_Base_PIDC_AgBgW_processed/Bharat_AW_Nalloc_api25e3sm_mf0df100_r0614_Base_PIDC_AgBgW_US-DUK_trans.nc"
#fnames["ORN_PIDC_RD_AgBgW"] = f"{path_in}Bharat_AW_Nalloc_api25e3sm_mf0df100_r0614_Base_PIDC_AgBgW_RD_processed/Bharat_AW_Nalloc_api25e3sm_mf0df100_r0614_Base_PIDC_AgBgW_RD_US-DUK_trans.nc"

fnames["DUK_PIDC_RD_AgBgW"] = f"{path_in}Bharat_AW_Nalloc_api25e3sm_mf0df100_r0614_Base_PIDC_AgBgW_RD_processed/Bharat_AW_Nalloc_api25e3sm_mf0df100_r0614_Base_PIDC_AgBgW_RD_US-DUK_trans.nc"
#fnames["ORN_PIDC_RD_AgBgW"] = f"{path_in}Bharat_AW_Nalloc_api25e3sm_mf0df100_r0614_Base_PIDC_AgBgW_RD_processed/Bharat_AW_Nalloc_api25e3sm_mf0df100_r0614_Base_PIDC_AgBgW_RD_US-DUK_trans.nc"

fnames["DUK_PIDD_RD_AgBgW"] = f"{path_in}Bharat_AW_Nalloc_api25e3sm_mf0df100_r0614_Base_PIDD_AgBgW_RD_processed/Bharat_AW_Nalloc_api25e3sm_mf0df100_r0614_Base_PIDD_AgBgW_RD_US-DUK_trans.nc"
fnames["ORN_PIDD_RD_AgBgW"] = f"{path_in}Bharat_AW_Nalloc_api25e3sm_mf0df100_r0614_Base_PIDD_AgBgW_RD_processed/Bharat_AW_Nalloc_api25e3sm_mf0df100_r0614_Base_PIDD_AgBgW_RD_US-ORN_trans.nc"

#fnames["DUK_PIDE_RD_AgBgW"] = f"{path_in}Bharat_AW_Nalloc_api25e3sm_mf0df100_r0621_Base_PIDE_AgBgW_RD_processed/Bharat_AW_Nalloc_api25e3sm_mf0df100_r0621_Base_PIDE_AgBgW_RD_US-DUK_trans.nc"
fnames["ORN_PIDE_RD_AgBgW"] = f"{path_in}Bharat_AW_Nalloc_api25e3sm_mf0df100_r0621_Base_PIDE_AgBgW_RD_processed/Bharat_AW_Nalloc_api25e3sm_mf0df100_r0621_Base_PIDE_AgBgW_RD_US-ORN_trans.nc"

fnames["DUK_PIDF_Conly_AgBgW"] = f"{path_in}Bharat_AW_Nalloc_api25e3sm_mf0df100_r0621_Base_PIDF_AgBgW_processed/Bharat_AW_Nalloc_api25e3sm_mf0df100_r0621_Base_PIDF_AgBgW_US-DUK_trans.nc"
fnames["ORN_PIDF_Conly_AgBgW"] = f"{path_in}Bharat_AW_Nalloc_api25e3sm_mf0df100_r0621_Base_PIDF_AgBgW_processed/Bharat_AW_Nalloc_api25e3sm_mf0df100_r0621_Base_PIDF_AgBgW_US-ORN_trans.nc"

fnames["DUK_PIDF_RD_AgBgW"] = f"{path_in}Bharat_AW_Nalloc_api25e3sm_mf0df100_r0621_Base_PIDF_AgBgW_RD_processed/Bharat_AW_Nalloc_api25e3sm_mf0df100_r0621_Base_PIDF_AgBgW_RD_US-DUK_trans.nc"
fnames["ORN_PIDF_RD_AgBgW"] = f"{path_in}Bharat_AW_Nalloc_api25e3sm_mf0df100_r0621_Base_PIDF_AgBgW_RD_processed/Bharat_AW_Nalloc_api25e3sm_mf0df100_r0621_Base_PIDF_AgBgW_RD_US-ORN_trans.nc"

fnames["DUK_PIDG_Conly_AgBgW"] = f"{path_in}Bharat_AW_Nalloc_api25e3sm_mf0df100_r0621_Base_PIDG_AgBgW_processed/Bharat_AW_Nalloc_api25e3sm_mf0df100_r0621_Base_PIDG_AgBgW_US-DUK_trans.nc"
fnames["ORN_PIDG_Conly_AgBgW"] = f"{path_in}Bharat_AW_Nalloc_api25e3sm_mf0df100_r0621_Base_PIDG_AgBgW_processed/Bharat_AW_Nalloc_api25e3sm_mf0df100_r0621_Base_PIDG_AgBgW_US-ORN_trans.nc"

fnames["DUK_PIDG_RD_AgBgW"] = f"{path_in}Bharat_AW_Nalloc_api25e3sm_mf0df100_r0621_Base_PIDG_AgBgW_RD_processed/Bharat_AW_Nalloc_api25e3sm_mf0df100_r0621_Base_PIDG_AgBgW_RD_US-DUK_trans.nc"
fnames["ORN_PIDG_RD_AgBgW"] = f"{path_in}Bharat_AW_Nalloc_api25e3sm_mf0df100_r0621_Base_PIDG_AgBgW_RD_processed/Bharat_AW_Nalloc_api25e3sm_mf0df100_r0621_Base_PIDG_AgBgW_RD_US-ORN_trans.nc"

fnames["DUK_PIDH_Conly_AgBgW"] = f"{path_in}Bharat_AW_Nalloc_api25e3sm_mf0df100_r0621_Base_PIDH_AgBgW_processed/Bharat_AW_Nalloc_api25e3sm_mf0df100_r0621_Base_PIDH_AgBgW_US-DUK_trans.nc"
fnames["ORN_PIDH_Conly_AgBgW"] = f"{path_in}Bharat_AW_Nalloc_api25e3sm_mf0df100_r0621_Base_PIDH_AgBgW_processed/Bharat_AW_Nalloc_api25e3sm_mf0df100_r0621_Base_PIDH_AgBgW_US-ORN_trans.nc"

fnames["DUK_PIDH_RD_AgBgW"] = f"{path_in}Bharat_AW_Nalloc_api25e3sm_mf0df100_r0621_Base_PIDH_AgBgW_RD_processed/Bharat_AW_Nalloc_api25e3sm_mf0df100_r0621_Base_PIDH_AgBgW_RD_US-DUK_trans.nc"
fnames["ORN_PIDH_RD_AgBgW"] = f"{path_in}Bharat_AW_Nalloc_api25e3sm_mf0df100_r0621_Base_PIDH_AgBgW_RD_processed/Bharat_AW_Nalloc_api25e3sm_mf0df100_r0621_Base_PIDH_AgBgW_RD_US-ORN_trans.nc"


# AgBgW

fnames["DUK_PIDE_Conly_AgBgW"] = f"{path_in}Bharat_AW_Nalloc_mf0df100_r0707_Base_PIDE_AgBgW_processed/Bharat_AW_Nalloc_mf0df100_r0707_Base_PIDE_AgBgW_US-DUK_trans.nc"
fnames["DUK_PIDE_RD_AgBgW"] = f"{path_in}Bharat_AW_Nalloc_mf0df100_r0707_Base_PIDE_AgBgW_RD_processed/Bharat_AW_Nalloc_mf0df100_r0707_Base_PIDE_AgBgW_RD_US-DUK_trans.nc"



# BgW only

fnames["DUK_PIDE_Conly_BgW"] = f"{path_in}Bharat_AW_Nalloc_mf0df100_r0629_PIDE_BgW_processed/Bharat_AW_Nalloc_mf0df100_r0629_PIDE_BgW_US-DUK_trans.nc"
fnames["ORN_PIDE_Conly_BgW"] = f"{path_in}Bharat_AW_Nalloc_mf0df100_r0629_PIDE_BgW_processed/Bharat_AW_Nalloc_mf0df100_r0629_PIDE_BgW_US-ORN_trans.nc"

fnames["DUK_PIDE_RD_BgW"] = f"{path_in}Bharat_AW_Nalloc_mf0df100_r0629_PIDE_BgW_RD_processed/Bharat_AW_Nalloc_mf0df100_r0629_PIDE_BgW_RD_US-DUK_trans.nc"
fnames["ORN_PIDE_RD_BgW"] = f"{path_in}Bharat_AW_Nalloc_mf0df100_r0629_PIDE_BgW_RD_processed/Bharat_AW_Nalloc_mf0df100_r0629_PIDE_BgW_RD_US-ORN_trans.nc"

# RootW only

fnames["DUK_PIDE_Conly_RootW"] = f"{path_in}Bharat_AW_Nalloc_mf0df100_r0629_PIDE_RootW_processed/Bharat_AW_Nalloc_mf0df100_r0629_PIDE_RootW_US-DUK_trans.nc"
fnames["ORN_PIDE_Conly_RootW"] = f"{path_in}Bharat_AW_Nalloc_mf0df100_r0629_PIDE_RootW_processed/Bharat_AW_Nalloc_mf0df100_r0629_PIDE_RootW_US-ORN_trans.nc"

fnames["DUK_PIDE_RD_RootW"] = f"{path_in}Bharat_AW_Nalloc_mf0df100_r0629_PIDE_RootW_RD_processed/Bharat_AW_Nalloc_mf0df100_r0629_PIDE_RootW_RD_US-DUK_trans.nc"
fnames["ORN_PIDE_RD_RootW"] = f"{path_in}Bharat_AW_Nalloc_mf0df100_r0629_PIDE_RootW_RD_processed/Bharat_AW_Nalloc_mf0df100_r0629_PIDE_RootW_RD_US-ORN_trans.nc"

# default is CN 53
fnames["DUK_PIDE_Conly_AgBgW_CN55"] = f"{path_in}Bharat_AW_Nalloc_mf0df100_r0628_Base_PIDE_AgBgW_CN55_processed/Bharat_AW_Nalloc_mf0df100_r0628_Base_PIDE_AgBgW_CN55_US-DUK_trans.nc"
fnames["ORN_PIDE_Conly_AgBgW_CN55"] = f"{path_in}Bharat_AW_Nalloc_mf0df100_r0628_Base_PIDE_AgBgW_CN55_processed/Bharat_AW_Nalloc_mf0df100_r0628_Base_PIDE_AgBgW_CN55_US-ORN_trans.nc"

fnames["DUK_PIDE_Conly_AgBgW_CN35"] = f"{path_in}Bharat_AW_Nalloc_mf0df100_r0628_Base_PIDE_AgBgW_CN35_processed/Bharat_AW_Nalloc_mf0df100_r0628_Base_PIDE_AgBgW_CN35_US-DUK_trans.nc"
fnames["ORN_PIDE_Conly_AgBgW_CN35"] = f"{path_in}Bharat_AW_Nalloc_mf0df100_r0628_Base_PIDE_AgBgW_CN35_processed/Bharat_AW_Nalloc_mf0df100_r0628_Base_PIDE_AgBgW_CN35_US-ORN_trans.nc"

# default is CN 53
fnames["DUK_PIDE_Conly_AgBgW_CN55"] = f"{path_in}Bharat_AW_Nalloc_mf0df100_r0628_Base_PIDE_AgBgW_CN55_processed/Bharat_AW_Nalloc_mf0df100_r0628_Base_PIDE_AgBgW_CN55_US-DUK_trans.nc"
fnames["ORN_PIDE_Conly_AgBgW_CN55"] = f"{path_in}Bharat_AW_Nalloc_mf0df100_r0628_Base_PIDE_AgBgW_CN55_processed/Bharat_AW_Nalloc_mf0df100_r0628_Base_PIDE_AgBgW_CN55_US-ORN_trans.nc"

fnames["DUK_PIDE_RD_AgBgW_CN55"] = f"{path_in}Bharat_AW_Nalloc_mf0df100_r0628_Base_PIDE_AgBgW_CN55_RD_processed/Bharat_AW_Nalloc_mf0df100_r0628_Base_PIDE_AgBgW_CN55_RD_US-DUK_trans.nc"
#fnames["ORN_PIDE_RD_AgBgW_CN55"] = f"{path_in}Bharat_AW_Nalloc_mf0df100_r0628_Base_PIDE_AgBgW_CN55_RD_processed/Bharat_AW_Nalloc_mf0df100_r0628_Base_PIDE_AgBgW_CN55_RD_US-ORN_trans.nc"

fnames["DUK_PIDE_Conly_AgBgW_CN35"] = f"{path_in}Bharat_AW_Nalloc_mf0df100_r0628_Base_PIDE_AgBgW_CN35_processed/Bharat_AW_Nalloc_mf0df100_r0628_Base_PIDE_AgBgW_CN35_US-DUK_trans.nc"
fnames["ORN_PIDE_Conly_AgBgW_CN35"] = f"{path_in}Bharat_AW_Nalloc_mf0df100_r0628_Base_PIDE_AgBgW_CN35_processed/Bharat_AW_Nalloc_mf0df100_r0628_Base_PIDE_AgBgW_CN35_US-ORN_trans.nc"

fnames["DUK_PIDE_RD_AgBgW_CN35"] = f"{path_in}Bharat_AW_Nalloc_mf0df100_r0628_Base_PIDE_AgBgW_CN35_RD_processed/Bharat_AW_Nalloc_mf0df100_r0628_Base_PIDE_AgBgW_CN35_RD_US-DUK_trans.nc"
fnames["ORN_PIDE_RD_AgBgW_CN35"] = f"{path_in}Bharat_AW_Nalloc_mf0df100_r0628_Base_PIDE_AgBgW_CN35_RD_processed/Bharat_AW_Nalloc_mf0df100_r0628_Base_PIDE_AgBgW_CN35_RD_US-ORN_trans.nc"

#fnames["DUK_PIDzero_RD"] = : Simulation missing *Running
#fnames["ORN_PIDzero_RD"] = f"{path_in}Bharat_AW_Nalloc_api25e3sm_mf0df100_r06122_Base_PIDzero_RD_processed/Bharat_AW_Nalloc_api25e3sm_mf0df100_r06122_Base_PIDzero_RD_US-ORN_trans.nc"

## This is for testing_vals.txt

fnames["DUK_PIDE_RD_AgBgW_testing"] = f"{path_in}Bharat_AW_Nalloc_mf0df100_r0726_l2fr_PIDE_CN53_AgBgW_RD_processed/Bharat_AW_Nalloc_mf0df100_r0726_l2fr_PIDE_CN53_AgBgW_RD_US-DUK_trans.nc"

# Other global attributes for plots

logging_year= 1855

ds = {}

for idx, key in enumerate(fnames.keys()):
    print (key)
    ds[key] = xr.open_mfdataset(fnames[key],decode_times=True)
for idx, key in enumerate(ds.keys()):
    ds[key]['time'] =  pd.to_datetime(ds[key].time.values.astype(str))

DUK_PIDA_Conly_AgBgW
ORN_PIDA_Conly_AgBgW
DUK_PIDB_Conly
ORN_PIDB_Conly
DUK_PIDC_Conly
ORN_PIDC_Conly
DUK_PIDD_Conly
ORN_PIDD_Conly
DUK_PIDA_RD_AgBgW
ORN_PIDA_RD_AgBgW
DUK_PIDB_RD_AgBgW
ORN_PIDB_RD_AgBgW
DUK_PIDC_Conly_AgBgW
DUK_PIDC_RD_AgBgW
DUK_PIDD_RD_AgBgW
ORN_PIDD_RD_AgBgW
ORN_PIDE_RD_AgBgW
DUK_PIDF_Conly_AgBgW
ORN_PIDF_Conly_AgBgW
DUK_PIDF_RD_AgBgW
ORN_PIDF_RD_AgBgW
DUK_PIDG_Conly_AgBgW
ORN_PIDG_Conly_AgBgW
DUK_PIDG_RD_AgBgW
ORN_PIDG_RD_AgBgW
DUK_PIDH_Conly_AgBgW
ORN_PIDH_Conly_AgBgW
DUK_PIDH_RD_AgBgW
ORN_PIDH_RD_AgBgW
DUK_PIDE_Conly_AgBgW
DUK_PIDE_RD_AgBgW
DUK_PIDE_Conly_BgW
ORN_PIDE_Conly_BgW
DUK_PIDE_RD_BgW
ORN_PIDE_RD_BgW
DUK_PIDE_Conly_RootW
ORN_PIDE_Conly_RootW
DUK_PIDE_RD_RootW
ORN_PIDE_RD_RootW
DUK_PIDE_Conly_AgBgW_CN55
ORN_PIDE_Conly_AgBgW_CN55
DUK_PIDE_Conly_AgBgW_CN35
ORN_PIDE_Conly_AgBgW_CN35
DUK_PIDE_RD_AgBgW_CN55
DUK_PIDE_RD_AgBgW_CN35
ORN_PIDE_RD_AgBgW_CN35
DUK_PIDE_RD_AgBgW_testing

Enter name of the single simulation that you want to investigate:

sims = "ORN_PIDB_RD_AgBgW"
sims = "DUK_PIDB_RD_AgBgW"
sims = "DUK_PIDC_RD_AgBgW"
#sims = "DUK_PIDD_RD_AgBgW"
#sims = "ORN_PIDD_RD_AgBgW"
#sims = "DUK_PIDE_RD_AgBgW"
#sims = "ORN_PIDE_RD_AgBgW"
#sims = "DUK_PIDF_RD_AgBgW"

#sims = "ORN_PIDF_RD_AgBgW"
#sims = "DUK_PIDG_RD_AgBgW"
#sims = "ORN_PIDH_RD_AgBgW"
#sims = "DUK_PIDH_RD_AgBgW"

#sims = "DUK_PIDE_Conly_BgW"
#sims = "ORN_PIDE_Conly_BgW"

#sims = "DUK_PIDE_RD_BgW"
#sims = "ORN_PIDE_RD_BgW"

#sims = "DUK_PIDE_Conly_RootW"
#sims = "ORN_PIDE_Conly_RootW"

#sims = "DUK_PIDE_RD_RootW"
#sims = "ORN_PIDE_RD_RootW" 

#sims = "DUK_PIDE_RD_AgBgW_CN55"
#sims = "ORN_PIDE_RD_AgBgW_CN55" # Not available

#sims = "DUK_PIDE_RD_AgBgW_CN35"
#sims = "ORN_PIDE_RD_AgBgW_CN35"

sims = "DUK_PIDE_RD_AgBgW_testing"

Select the timeperiod you are interested to plot for daily timeseries

# Logging happened on Jan 1, 1855

# Start Date Index  
idx_start_date = 4*365+1 

idx_logging = 5*365+1 

# For multiple variables on a same plot upto 5 years after logging 
Yrs_Next = 5 # upto 5 years after logging 

Daily Plots

print ("\nAll Plots >>>>>>")

fig, axs = plt.subplots(9,1 , figsize=(15,30), dpi=400)

# Subplot 1 >>>
idx_ax = 0
var = "FATES_NPP"
# C-only
try:
    sim_conly = sims.replace("RD", "Conly")
    if sims == "DUK_PIDE_RD_AgBgW_testing" : 
        sim_conly = "DUK_PIDE_Conly_AgBgW"
except:
    sim_conly = "DUK_PIDE_Conly_AgBgW"
    
key=sim_conly
ts_data = ds[key][var][idx_start_date:idx_start_date+Yrs_Next*365] # Daily data range from logging date to next five years
axs[idx_ax].plot(ts_data.time, ts_data, label = var + " C-only")
# RD
key=sims
ts_data = ds[key][var][idx_start_date:idx_start_date+Yrs_Next*365] # Daily data range from logging date to next five years
axs[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 = 1
vars_plot = (
"""
FATES_STOREC_TF
FATES_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 years
axs[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 years
ax_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 = 2
vars_plot = (
"""
FATES_STOREN_TF
FATES_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 years
axs[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 years
ax_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 = 3
ts_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 years
axs[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 years
ax_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 = 4
ts_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 years
axs[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()
ts_data = ds[key]["FATES_L2FR"][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_L2FR")
ax_tmp.set_ylabel (f"{var}({ds[key][var].units})")
ax_tmp.legend(loc = 5)
# Subplot 5 <<<

# Subplot 6 >>>
idx_ax = 5
vars_plot = (
"""
FATES_LEAFC
FATES_FROOTC
"""
).split('\n')
vars_plot = vars_plot[1:-1]
for i_var,var in enumerate(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 years
    axs[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 years
ax_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 = 6
vars_plot = (
"""
FATES_NO3UPTAKE
FATES_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 years
axs[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 years
ax_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 = 7
vars_plot = (
"""
GROSS_NMIN
NET_NMIN
"""
).split('\n')
vars_plot = vars_plot[1:-1]
for i_var,var in enumerate(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 years
    axs[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 = 8
vars_plot = (
"""
SMIN_NO3_LEACHED
DENIT
"""
).split('\n')
vars_plot = vars_plot[1:-1]
for i_var,var in enumerate(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 years
    axs[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)
for i in range(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_PIDE_RD_AgBgW_testing

/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))

DUK_PIDE_RD_AgBgW_testing

Working on the PID Variables

• Kp, Kd, Ki : inputs Calculate:
• fcn = SafeLog(MAXnp((C/C')/(M/M'))
• integrate:
• derivative:
• multiples

Hint:

'''

   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

'''

# PID Params

#E|0.001|0.5
print (sims.split('_')[1])
if sims.split('_')[1] == 'PIDE':
    pid_kp = 0.001
    pid_kd = 0.5
    pid_ki = 0.0
if sims.split('_')[1] == 'PIDH':
    pid_kp = 0.001
    pid_kd = 1.0
    pid_ki = 0.0
    
if sims.split('_')[1] == 'PIDC':
    pid_kp = 0.0001
    pid_kd = 0.01
    pid_ki = 0.0
    
if sims.split('_')[1] == 'PIDF':
    pid_kp = 0.005
    pid_kd = 0.1
    pid_ki = 0.0
    
if sims.split('_')[1] == 'PIDG':
    pid_kp = 0.005
    pid_kd = 0.5
    pid_ki = 0.0

carbon_nitrogen = True        
carbon_phosphorus = True    

if carbon_nitrogen == True:
    c2x = "CN"
    
if carbon_phosphorus == True:
    c2x = "CP"

if carbon_nitrogen ==True and carbon_phosphorus == True:
    c2x = "CNP"

PIDE

def SafeLog(val):
    # As defined in PRTAllometricCNPMod
    safelog_min = 0.001
    safelog_max = 1000
    clipped_val = np.clip(val, safelog_min, safelog_max)
    logval = np.log(clipped_val)
    return logval

if c2x == "CN" and c2x != "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)
if c2x != "CN" and c2x == "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)
if c2x == "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) # fcn
cx_logratio_ar = SafeLog(cp_ratio) # fcn
cx_logratio_ar = np.ravel(cx_logratio_ar)
#cx_logratio_ar = np.insert(cx_logratio_ar, 0, 0) # adding first value as 0

#derivative
dcxdt_ratio_mat =  cx_logratio_ar[1:] - cx_logratio_ar[:-1]

case3

ts_data_l2fr = ds[key]['FATES_L2FR'][idx_start_date:idx_start_date+Yrs_Next*365]
ts = ts_data_l2fr.values
delta_l2fr_model = ts[1:] - ts[:-1]

delta_l2fr_model[:,0][0]
ema_dcxdt_init = (delta_l2fr_model[:,0][0] - pid_kp*dcxdt_ratio_mat[0])/pid_kd
ema_dcxdt_init

-0.000216938316822052

#integration
pid_drv_wgt = 1/20.   # n-day smoothing of the derivative of the process function in the PID controller
cx_int = cx_logratio_ar[0]#0
cx0 = 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 earlier
nearzero = 1e-30
cx_int_ar = []
dcxdt_ratio_ar = []
ema_dcxdt_ar = []
for idx, cx_logratio in enumerate (cx_logratio_ar[1:]):
    cx_int = cx_int + cx_logratio
    #Reset the integrator if its sign changes
    if abs(cx_logratio)> nearzero and abs(cx0)> nearzero :
        if abs(cx_logratio/abs(cx_logratio) - cx0/abs(cx0)) > nearzero :
             cx_int = cx_logratio
    dcxdt_ratio = cx_logratio-cx0         
    ema_dcxdt = pid_drv_wgt*dcxdt_ratio + (1.-pid_drv_wgt)*ema_dcxdt
    cx0 = cx_logratio
    if idx == 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 begining
cx_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])

5.8472157e-05 -0.0001849060775339603

# Sanity check
sum(dcxdt_ratio_ar[1:] == dcxdt_ratio_mat)

1824

l2fr_delta_ar = []

for idx in range(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_dcxdt
    l2fr_delta_ar.append(l2fr_delta)
l2fr_delta_ar = np.array(l2fr_delta_ar)    

# Repeating the first value twice in the begining
l2fr_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

fig, axs = plt.subplots(10,1 , figsize=(15,40), dpi=400)


idx_ax = 0
#FATES_LEAFC
vars_plot = (
"""
FATES_LEAFC
"""
).split('\n')
vars_plot = vars_plot[1:-1]
for i_var,var in enumerate(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 years
    axs[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_ax
ts_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 years
axs[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 years
ax_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_ax
ts_data = (ds[key]["FATES_STOREC_TF"])[idx_start_date:idx_start_date+Yrs_Next*365] # Daily data range from logging date to next five years
axs[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 years
ax_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_ax
key=sims
ts_data = ds[key][var][idx_start_date:idx_start_date+Yrs_Next*365] # Daily data range from logging date to next five years
axs[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_ax
axs[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_ax
axs[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_ax
KpxFcn = pid_kp * cx_logratio_ar    
KdxDerivative = pid_kd * ema_dcxdt_ar    
KixIntegral = pid_ki * cx_int_ar    
print (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_ax
KpxFcn_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_ax
axs[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.values
delta_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_ax
L2FR_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')



for i in range(len(axs)): 
    axs[i].legend()
    axs[i].axvline(x = dt.date(1855,1,1), color = 'g',lw=5, alpha=.1)

pid_kp : 0.001
 pid_kd : 0.5
 pid_ki : 0.0 
pid_kp : 0.001
 pid_kd : 0.5
 pid_ki : 0.0 

Printing Variables that go in L2FR calculations

Printing the following lines of the routine CNPAdjustFRootTargets in the file parteh/PRTAllometricCNPMod.F90

852     print*, 'sinkhole',trim(dateTimeString) ,cx_logratio, cp_ratio, cn_ratio, cx_int, &
853     dcxdt_ratio, ema_dcxdt, cx0, l2fr_delta, l2fr, store_c_max, store_c_act, &
854     store_nut_max, store_nut_act, l2fr_min, store_N_max, store_N_act, &
855     store_P_max, store_P_act, dbh

filenames_sim_out = glob.glob(f"{path_in}*{sims.split('_')[1]}*r0811_l2fr.txt")

f"{path_in}*{sims.split('_')[1]}*l2fr.txt"

'/Users/ud4/FATESMDS_analysis/outputs/runs/tests_alp/230309/*PIDE*l2fr.txt'

filenames_sim_out

['/Users/ud4/FATESMDS_analysis/outputs/runs/tests_alp/230309/PIDE_AgBgW_r0811_l2fr.txt']

cols = ['timex', 'cx_logratio', 'cp_ratio', 'cn_ratio', 'cx_int', 
        'dcxdt_ratio', 'ema_dcxdt', 'cx0', 'l2fr_delta', 'l2fr',
       'store_c_max', 'store_c_act', 'store_nut_max', 'store_nut_act', 'l2fr_min',
       'store_N_max', 'store_N_act', 'store_P_max', 'store_P_act','dbh' ]

df_vals = pd.read_csv(filenames_sim_out[0], header=None, delim_whitespace=True)
df_vals.columns = cols
# removing the leading string from first column
df_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 value
df_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)

Sanity Check: 57670 should be equal to 57670

	timex	cx_logratio	cp_ratio	cn_ratio	cx_int	dcxdt_ratio	ema_dcxdt	cx0	l2fr_delta	l2fr	...	store_c_act	store_nut_max	store_nut_act	l2fr_min	store_N_max	store_N_act	store_P_max	store_P_act	dbh	timex_code
0	1850-01-01	0.001595	0.458194	1.001596	0.025248	-0.000403	8.545468e-06	0.001595	0.000006	0.212009	...	7.851091	0.003764	0.008225	0.01	0.058807	0.058789	0.003764	0.008225	11.981116	0
1	1850-01-01	0.001515	0.457973	1.001516	0.024135	-0.000409	8.235517e-06	0.001515	0.000006	0.213019	...	9.326194	0.004471	0.009775	0.01	0.069861	0.069841	0.004471	0.009775	13.705120	0
2	1850-01-01	0.001391	0.454351	1.001392	0.008578	-0.000423	7.622402e-06	0.001391	0.000005	0.232612	...	11.930909	0.005721	0.012604	0.01	0.089384	0.089357	0.005721	0.012604	16.611559	0
3	1850-01-01	0.001267	0.452071	1.001267	0.007863	-0.000434	7.100165e-06	0.001267	0.000005	0.244996	...	15.056028	0.007220	0.015986	0.01	0.112811	0.112777	0.007220	0.015986	19.920853	0
4	1850-01-01	0.001173	0.462080	1.001173	0.007309	-0.000441	6.804564e-06	0.001173	0.000005	0.190013	...	18.201289	0.008729	0.018907	0.01	0.136390	0.136350	0.008729	0.018907	23.101532	0
5	1850-01-01	0.001043	0.463048	1.001044	0.006572	-0.000455	6.167436e-06	0.001043	0.000004	0.184538	...	22.708153	0.010892	0.023539	0.01	0.170184	0.170134	0.010892	0.023539	27.458106	0
6	1850-01-01	-0.285429	0.432666	0.751692	-241.552691	-0.000502	-1.326177e-04	-0.285429	-0.000352	0.357088	...	26.556579	0.012739	0.029462	0.01	0.199052	0.264968	0.012739	0.029462	31.029836	0
7	1850-01-01	0.003750	0.379747	1.003757	0.047627	-0.000075	4.534833e-05	0.003750	0.000026	0.733943	...	28.300767	0.013538	0.035772	0.01	0.211538	0.211461	0.013538	0.035772	32.538742	0
8	1850-01-02	0.004025	0.064169	1.004034	0.004025	0.004025	2.012715e-04	0.004025	0.000105	1.086434	...	0.076287	0.000037	0.000571	0.01	0.000570	0.000570	0.000037	0.000571	0.390945	1
9	1850-01-02	0.001722	0.236868	1.001723	0.026970	0.000127	1.445319e-05	0.001722	0.000009	0.212018	...	7.852568	0.003764	0.015913	0.01	0.058807	0.058793	0.003764	0.015913	11.981116	1
10	1850-01-02	0.001643	0.236257	1.001644	0.025778	0.000128	1.423871e-05	0.001643	0.000009	0.213028	...	9.327979	0.004471	0.018952	0.01	0.069861	0.069845	0.004471	0.018952	13.705120	1
11	1850-01-02	0.001493	0.225283	1.001494	0.010071	0.000103	1.237038e-05	0.001493	0.000008	0.232619	...	11.933253	0.005721	0.025426	0.01	0.089384	0.089366	0.005721	0.025426	16.611559	1
12	1850-01-02	0.001355	0.218846	1.001356	0.009218	0.000088	1.115221e-05	0.001355	0.000007	0.245003	...	15.059062	0.007220	0.033029	0.01	0.112811	0.112790	0.007220	0.033029	19.920853	1
13	1850-01-02	0.001350	0.250470	1.001351	0.008659	0.000178	1.535049e-05	0.001350	0.000009	0.190022	...	18.205020	0.008729	0.034888	0.01	0.136390	0.136353	0.008729	0.034888	23.101532	1
14	1850-01-02	0.001234	0.254098	1.001235	0.007806	0.000191	1.542809e-05	0.001234	0.000009	0.184547	...	22.712927	0.010892	0.042906	0.01	0.170184	0.170137	0.010892	0.042906	27.458106	1
15	1850-01-02	-0.285365	0.173968	0.751740	-241.838056	0.000065	-1.227611e-04	-0.285365	-0.000347	0.356741	...	26.562306	0.012739	0.073289	0.01	0.199052	0.265008	0.012739	0.073289	31.029836	1
16	1850-01-02	0.003434	0.103309	1.003440	0.041916	-0.000305	2.781351e-05	0.003434	0.000017	0.733946	...	28.302045	0.013539	0.131503	0.01	0.211544	0.211544	0.013539	0.131503	32.539486	1
17	1850-01-03	-0.001947	0.063845	0.998054	-0.001947	-0.001947	-9.737284e-05	-0.001947	-0.000051	1.086279	...	0.075991	0.000036	0.000571	0.01	0.000570	0.000571	0.000036	0.000571	0.390608	2
18	1850-01-03	-0.001588	0.060027	0.998413	-0.001588	-0.005614	-8.948259e-05	-0.001588	-0.000046	1.086387	...	0.076090	0.000037	0.000608	0.01	0.000570	0.000572	0.000037	0.000608	0.390945	2
19	1850-01-03	-0.000326	0.236395	0.999674	-0.000326	-0.002048	-8.865855e-05	-0.000326	-0.000045	0.211973	...	7.836984	0.003764	0.015913	0.01	0.058807	0.058797	0.003764	0.015913	11.981116	2
20	1850-01-03	-0.000405	0.235785	0.999596	-0.000405	-0.002047	-8.883953e-05	-0.000405	-0.000045	0.212983	...	9.309487	0.004471	0.018952	0.01	0.069861	0.069849	0.004471	0.018952	13.705120	2
21	1850-01-03	-0.000582	0.224834	0.999418	-0.000582	-0.002075	-9.200611e-05	-0.000582	-0.000047	0.232573	...	11.909632	0.005721	0.025426	0.01	0.089384	0.089374	0.005721	0.025426	16.611559	2
22	1850-01-03	-0.000737	0.218411	0.999263	-0.000737	-0.002092	-9.398426e-05	-0.000737	-0.000048	0.244955	...	15.029302	0.007220	0.033030	0.01	0.112811	0.112803	0.007220	0.033030	19.920853	2
23	1850-01-03	-0.000652	0.249973	0.999348	-0.000652	-0.002002	-8.552831e-05	-0.000652	-0.000043	0.189978	...	18.169092	0.008729	0.034888	0.01	0.136390	0.136357	0.008729	0.034888	23.101532	2
24	1850-01-03	-0.000756	0.253594	0.999244	-0.000756	-0.001991	-8.488657e-05	-0.000756	-0.000043	0.184504	...	22.668165	0.010892	0.042906	0.01	0.170184	0.170140	0.010892	0.042906	27.458106	2
25	1850-01-03	-0.289091	0.173476	0.748944	-242.127147	-0.003726	-3.029400e-04	-0.289091	-0.000441	0.356301	...	26.471357	0.012739	0.073245	0.01	0.199052	0.265087	0.012739	0.073245	31.029836	2
26	1850-01-03	0.000355	0.102997	1.000355	0.041221	-0.003078	-1.274662e-04	0.000355	-0.000063	0.733881	...	28.215919	0.013540	0.131505	0.01	0.211559	0.211559	0.013540	0.131505	32.541223	2
27	1850-01-04	0.003733	0.064150	1.003740	0.003733	0.003733	1.866365e-04	0.003733	0.000097	1.086426	...	0.076265	0.000037	0.000571	0.01	0.000570	0.000570	0.000037	0.000571	0.390945	3
28	1850-01-04	0.001533	0.060297	1.001535	0.001533	0.003356	7.789724e-05	0.001533	0.000040	1.086357	...	0.076287	0.000037	0.000608	0.01	0.000571	0.000572	0.000037	0.000608	0.391031	3
29	1850-01-04	0.001350	0.236827	1.001351	0.001350	0.001676	-4.253656e-07	0.001350	0.000001	0.211974	...	7.850623	0.003764	0.015912	0.01	0.058807	0.058800	0.003764	0.015912	11.981116	3

30 rows × 21 columns

Number of Cohorts over time

print (f"Maximum Cohorts: {max(np.unique(df_vals['timex_code'], return_counts=1)[1])}")

Maximum Cohorts: 32

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)

Text(0.5, 1.0, 'Cohorts Vs Time (daily)\n')

FATES_NPP just for reference

ts_data = ds["DUK_PIDE_RD_AgBgW_testing"]["FATES_NPP"]
#ts_data
plt.figure(figsize=(15,6))
plt.scatter(x=ts_data.time, y=ts_data, s=2)
plt.title( "FATES_NPP_Daily\n", fontsize  = 16)

plt.figure(figsize=(15,6))
ts_data = ts_data.groupby("time.year").mean('time')
plt.scatter(x=ts_data.year, y=ts_data)
plt.title( "FATES_NPP Annual Mean\n", fontsize  = 16)

Text(0.5, 1.0, 'FATES_NPP Annual Mean\n')

Plot of L2FR from print file

# Set the date range you want to subset between
start_date = '1854-01-01'
end_date = '1859-01-01'

# Create a boolean mask based on the date range
mask = (df_vals['timex'] >= start_date) & (df_vals['timex'] <= end_date)

# Subset the DataFrame based on the boolean mask
subset_df = df_vals[mask]
print (f"Maximum Cohorts: {max(np.unique(subset_df['timex_code'], return_counts=1)[1])}")

Maximum Cohorts: 17

# Plotting based on the cohort index

# initializing the cohort index to -1
subset_df['cohort_id'] = -1

# for every time step assign a unique index in the cohort id
for _, grp in subset_df.groupby('timex'):
    subset_df.loc[grp.index, "cohort_id"] = range(len(grp))

/var/folders/f1/01gxw8vn74q_x_rf_p5ztryjr405zq/T/ipykernel_67082/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
  subset_df['cohort_id'] = -1
/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
  self._setitem_single_column(ilocs[0], value, pi)

# Which CMAP coloring to chose
cmap_select = 'Wistia'


fig, axs = plt.subplots(13,1 , figsize=(20,50), dpi=100)

idx_ax = 0

axs[idx_ax].scatter(x=subset_df['timex'], y=subset_df['dbh'],  
                    marker='.', label = 'dbh', c=subset_df['dbh'], 
                    cmap=cmap_select, s=5)
axs[idx_ax].set_title(f"{key} \n",fontsize=16)
'''
axs[idx_ax].scatter(x=subset_df['timex'], y=subset_df['store_c_act'],  
                    marker='.', label = 'store_c_act', c=subset_df['cohort_id'], 
                    cmap=cmap_select, s=5)

idx_ax = idx_ax +1
axs[idx_ax].scatter(x=subset_df['timex'], y=subset_df['store_c_max'],  
                    marker='.', label = 'store_c_max', c=subset_df['cohort_id'], 
                    cmap=cmap_select, s=5)
'''
idx_ax = idx_ax +1
axs[idx_ax].scatter(x=subset_df['timex'], y=subset_df['store_c_act']/subset_df['store_c_max'],  
                    marker='.', label = 'store_c_TF', c=subset_df['dbh'], 
                    cmap=cmap_select, s=5)

idx_ax = idx_ax +1


'''
axs[idx_ax].scatter(x=subset_df['timex'], y=subset_df['store_nut_act'],  
                    marker='.', label = 'store_nut_act', c=subset_df['cohort_id'], 
                    cmap=cmap_select, s=5)

idx_ax = idx_ax +1
axs[idx_ax].scatter(x=subset_df['timex'], y=subset_df['store_nut_max'],  
                    marker='.', label = 'store_nut_max', c=subset_df['cohort_id'], 
                    cmap=cmap_select, s=5)
'''
axs[idx_ax].scatter(x=subset_df['timex'], y=subset_df['store_N_act']/subset_df['store_N_max'],  
                    marker='.', label = 'store_N_TF', c=subset_df['dbh'], 
                    cmap=cmap_select, s=5)

idx_ax = idx_ax +1

axs[idx_ax].scatter(x=subset_df['timex'], y=subset_df['store_P_act']/subset_df['store_P_max'],  
                    marker='.', label = 'store_P_TF', c=subset_df['dbh'], 
                    cmap=cmap_select, s=5)

idx_ax = idx_ax +1

axs[idx_ax].scatter(x=subset_df['timex'], y=subset_df['cn_ratio'],  
                    marker='.', label = 'cn_ratio', c=subset_df['dbh'], 
                    cmap=cmap_select, s=5)
idx_ax = idx_ax +1
axs[idx_ax].scatter(x=subset_df['timex'], y=subset_df['cp_ratio'],  
                    marker='.', label = 'cp_ratio', c=subset_df['dbh'], 
                    cmap=cmap_select, s=5)


idx_ax = idx_ax +1
axs[idx_ax].scatter(x=subset_df['timex'], y=subset_df['cx_logratio'],  
                    marker='.', label = 'cx_logratio', 
                    c=subset_df['dbh'],  cmap=cmap_select, s=5)

idx_ax = idx_ax +1
axs[idx_ax].scatter(x=subset_df['timex'], y=subset_df['dcxdt_ratio'],  
                    marker='o', label = 'dcxdt_ratio', c=subset_df['dbh'], 
                    cmap=cmap_select, s=15)

idx_ax = idx_ax +1
axs[idx_ax].scatter(x=subset_df['timex'], y=subset_df['ema_dcxdt'],  
                    marker='.', label = 'ema_dcxdt', c=subset_df['dbh'], 
                    cmap=cmap_select, s=5)

idx_ax = idx_ax +1
axs[idx_ax].scatter(x=subset_df['timex'], y=subset_df['cx0'],  
                    marker='.', label = 'cx0', c=subset_df['dbh'], 
                    cmap=cmap_select, s=5)

idx_ax = idx_ax +1
axs[idx_ax].scatter(x=subset_df['timex'], y=subset_df['cx_int'],  
                    marker='.', label = 'cx_int', c=subset_df['dbh'], 
                    cmap=cmap_select, s=5)

idx_ax = idx_ax +1
axs[idx_ax].scatter(x=subset_df['timex'], y=subset_df['l2fr_delta'],  
                    marker='.', label = 'l2fr_delta', 
                    c=subset_df['dbh'],  cmap=cmap_select, s=5)

idx_ax = idx_ax +1
axs[idx_ax].scatter(x=subset_df['timex'], y=subset_df['l2fr'],  
                    marker='.', label = 'l2fr', c=subset_df['dbh'], 
                    cmap=cmap_select, s=5)





for i in range(len(axs)): 
    axs[i].legend(loc=5, fontsize=14)
    axs[i].axvline(x = dt.date(1855,1,1), color = 'r',lw=5, alpha=.1)

# Which CMAP coloring to chose
cmap_select = 'tab20b'


fig, axs = plt.subplots(13,1 , figsize=(20,50), dpi=100)

idx_ax = 0

axs[idx_ax].scatter(x=subset_df['timex'], y=subset_df['dbh'],  
                    marker='.', label = 'dbh', c=subset_df['dbh'], 
                    cmap=cmap_select, s=5)
axs[idx_ax].set_title(f"{key}  \n",fontsize=16)
'''
axs[idx_ax].scatter(x=subset_df['timex'], y=subset_df['store_c_act'],  
                    marker='.', label = 'store_c_act', c=subset_df['cohort_id'], 
                    cmap=cmap_select, s=5)

idx_ax = idx_ax +1
axs[idx_ax].scatter(x=subset_df['timex'], y=subset_df['store_c_max'],  
                    marker='.', label = 'store_c_max', c=subset_df['cohort_id'], 
                    cmap=cmap_select, s=5)
'''
idx_ax = idx_ax +1
axs[idx_ax].scatter(x=subset_df['timex'], y=subset_df['store_c_act']/subset_df['store_c_max'],  
                    marker='.', label = 'store_c_TF', c=subset_df['dbh'], 
                    cmap=cmap_select, s=5)

idx_ax = idx_ax +1


'''
axs[idx_ax].scatter(x=subset_df['timex'], y=subset_df['store_nut_act'],  
                    marker='.', label = 'store_nut_act', c=subset_df['cohort_id'], 
                    cmap=cmap_select, s=5)

idx_ax = idx_ax +1
axs[idx_ax].scatter(x=subset_df['timex'], y=subset_df['store_nut_max'],  
                    marker='.', label = 'store_nut_max', c=subset_df['cohort_id'], 
                    cmap=cmap_select, s=5)
'''
axs[idx_ax].scatter(x=subset_df['timex'], y=subset_df['store_N_act']/subset_df['store_N_max'],  
                    marker='.', label = 'store_N_TF', c=subset_df['dbh'], 
                    cmap=cmap_select, s=5)

idx_ax = idx_ax +1

axs[idx_ax].scatter(x=subset_df['timex'], y=subset_df['store_P_act']/subset_df['store_P_max'],  
                    marker='.', label = 'store_P_TF', c=subset_df['dbh'], 
                    cmap=cmap_select, s=5)

idx_ax = idx_ax +1

axs[idx_ax].scatter(x=subset_df['timex'], y=subset_df['cn_ratio'],  
                    marker='.', label = 'cn_ratio', c=subset_df['dbh'], 
                    cmap=cmap_select, s=5)
idx_ax = idx_ax +1
axs[idx_ax].scatter(x=subset_df['timex'], y=subset_df['cp_ratio'],  
                    marker='.', label = 'cp_ratio', c=subset_df['dbh'], 
                    cmap=cmap_select, s=5)


idx_ax = idx_ax +1
axs[idx_ax].scatter(x=subset_df['timex'], y=subset_df['cx_logratio'],  
                    marker='.', label = 'cx_logratio', 
                    c=subset_df['dbh'],  cmap=cmap_select, s=5)

idx_ax = idx_ax +1
axs[idx_ax].scatter(x=subset_df['timex'], y=subset_df['dcxdt_ratio'],  
                    marker='o', label = 'dcxdt_ratio', c=subset_df['dbh'], 
                    cmap=cmap_select, s=15)

idx_ax = idx_ax +1
axs[idx_ax].scatter(x=subset_df['timex'], y=subset_df['ema_dcxdt'],  
                    marker='.', label = 'ema_dcxdt', c=subset_df['dbh'], 
                    cmap=cmap_select, s=5)

idx_ax = idx_ax +1
axs[idx_ax].scatter(x=subset_df['timex'], y=subset_df['cx0'],  
                    marker='.', label = 'cx0', c=subset_df['dbh'], 
                    cmap=cmap_select, s=5)

idx_ax = idx_ax +1
axs[idx_ax].scatter(x=subset_df['timex'], y=subset_df['cx_int'],  
                    marker='.', label = 'cx_int', c=subset_df['dbh'], 
                    cmap=cmap_select, s=5)

idx_ax = idx_ax +1
axs[idx_ax].scatter(x=subset_df['timex'], y=subset_df['l2fr_delta'],  
                    marker='.', label = 'l2fr_delta', 
                    c=subset_df['dbh'],  cmap=cmap_select, s=5)

idx_ax = idx_ax +1
axs[idx_ax].scatter(x=subset_df['timex'], y=subset_df['l2fr'],  
                    marker='.', label = 'l2fr', c=subset_df['dbh'], 
                    cmap=cmap_select, s=5)





for i in range(len(axs)): 
    axs[i].legend(loc=5, fontsize=14)
    axs[i].axvline(x = dt.date(1855,1,1), color = 'r',lw=5, alpha=.1)

Plotting for multiple decades

# Set the date range you want to subset between
start_date = '1854-01-01'
end_date = '1894-01-01'

# Create a boolean mask based on the date range
mask = (df_vals['timex'] >= start_date) & (df_vals['timex'] <= end_date)

# Subset the DataFrame based on the boolean mask
subset_df = df_vals[mask]
print (f"Maximum Cohorts: {max(np.unique(subset_df['timex_code'], return_counts=1)[1])}")

Maximum Cohorts: 29

# Which CMAP coloring to chose
cmap_select = 'viridis'


fig, axs = plt.subplots(13,1 , figsize=(20,50), dpi=100)

idx_ax = 0

axs[idx_ax].scatter(x=subset_df['timex'], y=subset_df['dbh'],  
                    marker='.', label = 'dbh', c=subset_df['dbh'], 
                    cmap=cmap_select, s=5)
axs[idx_ax].set_title(f"{key} {start_date} - {end_date} \n",fontsize=16)
'''
axs[idx_ax].scatter(x=subset_df['timex'], y=subset_df['store_c_act'],  
                    marker='.', label = 'store_c_act', c=subset_df['cohort_id'], 
                    cmap=cmap_select, s=5)

idx_ax = idx_ax +1
axs[idx_ax].scatter(x=subset_df['timex'], y=subset_df['store_c_max'],  
                    marker='.', label = 'store_c_max', c=subset_df['cohort_id'], 
                    cmap=cmap_select, s=5)
'''
idx_ax = idx_ax +1
axs[idx_ax].scatter(x=subset_df['timex'], y=subset_df['store_c_act']/subset_df['store_c_max'],  
                    marker='.', label = 'store_c_TF', c=subset_df['dbh'], 
                    cmap=cmap_select, s=5)

idx_ax = idx_ax +1


'''
axs[idx_ax].scatter(x=subset_df['timex'], y=subset_df['store_nut_act'],  
                    marker='.', label = 'store_nut_act', c=subset_df['cohort_id'], 
                    cmap=cmap_select, s=5)

idx_ax = idx_ax +1
axs[idx_ax].scatter(x=subset_df['timex'], y=subset_df['store_nut_max'],  
                    marker='.', label = 'store_nut_max', c=subset_df['cohort_id'], 
                    cmap=cmap_select, s=5)
'''
axs[idx_ax].scatter(x=subset_df['timex'], y=subset_df['store_N_act']/subset_df['store_N_max'],  
                    marker='.', label = 'store_N_TF', c=subset_df['dbh'], 
                    cmap=cmap_select, s=5)

idx_ax = idx_ax +1

axs[idx_ax].scatter(x=subset_df['timex'], y=subset_df['store_P_act']/subset_df['store_P_max'],  
                    marker='.', label = 'store_P_TF', c=subset_df['dbh'], 
                    cmap=cmap_select, s=5)

idx_ax = idx_ax +1

axs[idx_ax].scatter(x=subset_df['timex'], y=subset_df['cn_ratio'],  
                    marker='.', label = 'cn_ratio', c=subset_df['dbh'], 
                    cmap=cmap_select, s=5)
idx_ax = idx_ax +1
axs[idx_ax].scatter(x=subset_df['timex'], y=subset_df['cp_ratio'],  
                    marker='.', label = 'cp_ratio', c=subset_df['dbh'], 
                    cmap=cmap_select, s=5)


idx_ax = idx_ax +1
axs[idx_ax].scatter(x=subset_df['timex'], y=subset_df['cx_logratio'],  
                    marker='.', label = 'cx_logratio', 
                    c=subset_df['dbh'],  cmap=cmap_select, s=5)

idx_ax = idx_ax +1
axs[idx_ax].scatter(x=subset_df['timex'], y=subset_df['dcxdt_ratio'],  
                    marker='o', label = 'dcxdt_ratio', c=subset_df['dbh'], 
                    cmap=cmap_select, s=15)

idx_ax = idx_ax +1
axs[idx_ax].scatter(x=subset_df['timex'], y=subset_df['ema_dcxdt'],  
                    marker='.', label = 'ema_dcxdt', c=subset_df['dbh'], 
                    cmap=cmap_select, s=5)

idx_ax = idx_ax +1
axs[idx_ax].scatter(x=subset_df['timex'], y=subset_df['cx0'],  
                    marker='.', label = 'cx0', c=subset_df['dbh'], 
                    cmap=cmap_select, s=5)

idx_ax = idx_ax +1
axs[idx_ax].scatter(x=subset_df['timex'], y=subset_df['cx_int'],  
                    marker='.', label = 'cx_int', c=subset_df['dbh'], 
                    cmap=cmap_select, s=5)

idx_ax = idx_ax +1
axs[idx_ax].scatter(x=subset_df['timex'], y=subset_df['l2fr_delta'],  
                    marker='.', label = 'l2fr_delta', 
                    c=subset_df['dbh'],  cmap=cmap_select, s=5)

idx_ax = idx_ax +1
axs[idx_ax].scatter(x=subset_df['timex'], y=subset_df['l2fr'],  
                    marker='.', label = 'l2fr', c=subset_df['dbh'], 
                    cmap=cmap_select, s=5)





for i in range(len(axs)): 
    axs[i].legend(loc=5, fontsize=14)
    axs[i].axvline(x = dt.date(1855,1,1), color = 'r',lw=5, alpha=.1)

# Checking if the l2fr values (sum over all cohorts) from the print are same as l2fr variable output
max_vals = 5000
sum_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]

# Reading daily l2fr files FATES output
#list_of_d_files= glob.glob(f"/Users/ud4/FATESMDS_analysis/outputs/runs/tests_alp/230309/run/Bharat_AW_Nalloc_mf0df100_r0725_Test_PIDE_AgBgW_RD_US-DUK_I20TRELMFATES.elm.h0.*.nc")
#ds_FATES_daily = xr.open_mfdataset(f"/Users/ud4/FATESMDS_analysis/outputs/runs/tests_alp/230309/run/Bharat_AW_Nalloc_mf0df100_r0725_Test_PIDE_AgBgW_RD_US-DUK_I20TRELMFATES.elm.h0.*.nc")

#plt.plot(ds_FATES_daily['FATES_L2FR']) #sum_l2fr['l2fr'].plot() fig1, ax = plt.subplots(figsize=(30, 6)) scatter = ax.scatter(x=range(len(ds_FATES_daily['FATES_L2FR'])), y = ds_FATES_daily['FATES_L2FR'], c=ds_FATES_daily['FATES_L2FR'], cmap='viridis', marker='o', label = 'ds_FATES_daily') plt.legend()

#sum_l2fr['l2fr'].plot()

fig1, ax = plt.subplots(figsize=(30, 6))
scatter = ax.scatter(x=range(len(sum_l2fr['l2fr'])), y = sum_l2fr['l2fr'], c=sum_l2fr['l2fr'], cmap='viridis', marker='o', label = 'sum_l2fr')
plt.legend()

<matplotlib.legend.Legend at 0x17f99fc40>

#mean_l2fr['l2fr'].plot()

fig1, ax = plt.subplots(figsize=(30, 6))
scatter = ax.scatter(x=range(len(mean_l2fr['l2fr'])), y = mean_l2fr['l2fr'], c=mean_l2fr['l2fr'], cmap='viridis', marker='o', label = 'mean_l2fr')
plt.legend()

<matplotlib.legend.Legend at 0x17de1d310>

#min_l2fr['l2fr'].plot()


fig, ax = plt.subplots(figsize=(30, 6))
scatter = ax.scatter(x=range(len(min_l2fr['l2fr'])), y = min_l2fr['l2fr'], c=min_l2fr['l2fr'], cmap='viridis', marker='o', label = 'min_l2fr')
plt.legend()

<matplotlib.legend.Legend at 0x17cd30c70>

#max_l2fr['l2fr'].plot()


fig, ax = plt.subplots(figsize=(30, 6))
scatter = ax.scatter(x=range(len(max_l2fr['l2fr'])), y = max_l2fr['l2fr'], c=max_l2fr['l2fr'], cmap='viridis', marker='o', label = 'max_l2fr')
plt.legend()

<matplotlib.legend.Legend at 0x17c548430>

fig, ax = plt.subplots(figsize=(30, 6))
scatter = ax.scatter(ds['DUK_PIDE_RD_AgBgW_testing']['FATES_L2FR'].time, ds['DUK_PIDE_RD_AgBgW_testing']['FATES_L2FR'], c=ds['DUK_PIDE_RD_AgBgW_testing']['FATES_L2FR'], cmap='viridis', marker='o', label = 'FATES_L2FR_yrly')
plt.legend()

<matplotlib.legend.Legend at 0x17cd42160>

ts_data_l2fr = ds['DUK_PIDE_RD_AgBgW_testing']['FATES_L2FR'][idx_start_date:idx_start_date+Yrs_Next*365]

fig, ax = plt.subplots(figsize=(30, 6))
scatter = ax.scatter(ts_data_l2fr.time, ts_data_l2fr, c=ts_data_l2fr, cmap='viridis', marker='o', label = 'FATES_L2FR_Daily')
plt.legend()
ax.axvline(x = dt.date(1855,1,1), color = 'r',lw=5, alpha=.1)

<matplotlib.lines.Line2D at 0x18038e250>

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 = 0
var = "FATES_NPP"
# C-only
sim_conly = sims.replace("RD", "Conly")
if sims == "DUK_PIDE_RD_AgBgW_testing" : 
    sim_conly = "DUK_PIDE_Conly_AgBgW"
key=sim_conly
ts_data = ds[key][var][idx_start_date:idx_start_date+Yrs_Next*365] # Daily data range from logging date to next five years
axs[idx_ax].scatter(ts_data.time, ts_data, marker=".", label = var + " C-only",alpha=.2)
# RD
key=sims
ts_data = ds[key][var][idx_start_date:idx_start_date+Yrs_Next*365] # Daily data range from logging date to next five years
axs[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 = 1
vars_plot = (
"""
FATES_STOREC_TF
FATES_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 years
axs[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 years
ax_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 = 2
vars_plot = (
"""
FATES_STOREN_TF
FATES_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 years
axs[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 years
ax_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 = 3
ts_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 years
axs[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 years
ax_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 = 4
ts_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 years
axs[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 years
ax_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 = 5
vars_plot = (
"""
FATES_LEAFC
FATES_FROOTC
"""
).split('\n')
vars_plot = vars_plot[1:-1]
for i_var,var in enumerate(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 years
    axs[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 years
ax_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 = 6
vars_plot = (
"""
FATES_NO3UPTAKE
FATES_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 years
axs[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 years
ax_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 = 7
vars_plot = (
"""
GROSS_NMIN
NET_NMIN
"""
).split('\n')
vars_plot = vars_plot[1:-1]
for i_var,var in enumerate(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 years
    axs[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 = 8
vars_plot = (
"""
SMIN_NO3_LEACHED
DENIT
"""
).split('\n')
vars_plot = vars_plot[1:-1]
for i_var,var in enumerate(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 years
    axs[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)
for i in range(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_PIDE_RD_AgBgW_testing

/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))

DUK_PIDE_RD_AgBgW_testing

print ("\nAll Plots >>>>>>")
# For multiple variables on a same plot
fig, axs = plt.subplots(9,1 , figsize=(15,30), dpi=400)

# Subplot 1 >>>
idx_ax = 0
var = "FATES_NPP"
# C-only
sim_conly = sims.replace("RD", "Conly")
if sims == "DUK_PIDE_RD_AgBgW_testing" : 
    sim_conly = "DUK_PIDE_Conly_AgBgW"
key=sim_conly
ts_data = ds[key][var].groupby("time.year").mean('time')
axs[idx_ax].plot(ts_data.year, ts_data, label = var + " C-only")
# RD
key=sims
ts_data = ds[key][var].groupby("time.year").mean('time')
axs[idx_ax].plot(ts_data.year, ts_data, label = var )
axs[idx_ax].set_ylabel (f"{ds[key][var].units}")
axs[idx_ax].set_title(f"{key} Annual Mean\n",fontsize=16)

# Subplot 1 <<<

# Subplot 2 >>>
idx_ax = 1
vars_plot = (
"""
FATES_STOREC_TF
FATES_STOREC
"""
).split('\n')
vars_plot = vars_plot[1:-1]
var = vars_plot[0]
ts_data = ds[key][var].groupby("time.year").mean('time')
axs[idx_ax].plot(ts_data.year, ts_data, label = var)
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].groupby("time.year").mean('time')
ax_tmp.plot(ts_data.year,ts_data,'k+',  label = var)
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 = 2
vars_plot = (
"""
FATES_STOREN_TF
FATES_STOREN
"""
).split('\n')
vars_plot = vars_plot[1:-1]
var = vars_plot[0]
ts_data = ds[key][var].groupby("time.year").mean('time')
axs[idx_ax].plot(ts_data.year, ts_data, label = var)
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].groupby("time.year").mean('time')
ax_tmp.plot(ts_data.year,ts_data,'k+',  label = var)
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 = 3
ts_data = (ds[key]["FATES_STOREC_TF"]/ds[key]["FATES_STOREN_TF"]).groupby("time.year").mean('time')
axs[idx_ax].plot(ts_data.year, ts_data, label = "FATES_STOREC_TF/FATES_STOREN_TF")
axs[idx_ax].set_ylabel (f"{ds[key][var].units}")
ax_tmp = axs[idx_ax].twinx()
ts_data = (ds[key]["FATES_STOREC"]/ds[key]["FATES_STOREN"]).groupby("time.year").mean('time')
ax_tmp.plot(ts_data.year, 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 = 4
ts_data = (ds[key]["FATES_FROOT_ALLOC"]/ds[key]["FATES_LEAF_ALLOC"]).groupby("time.year").mean('time')
axs[idx_ax].plot(ts_data.year, ts_data, label = "FATES_FROOT_ALLOC/FATES_LEAF_ALLOC")
#axs[idx_ax].set_ylabel (f"{ds[key][var].units}")
ax_tmp = axs[idx_ax].twinx()
ts_data = ds[key]["FATES_L2FR"].groupby("time.year").mean('time')
ax_tmp.plot(ts_data.year, ts_data,'k+', label = "FATES_L2FR")
ax_tmp.set_ylabel (f"{ds[key]['FATES_L2FR'].units}")
ax_tmp.legend(loc = 5)
# Subplot 5 <<<

# Subplot 6 >>>
idx_ax = 5
vars_plot = (
"""
FATES_LEAFC
FATES_FROOTC
"""
).split('\n')
vars_plot = vars_plot[1:-1]
for i_var,var in enumerate(vars_plot):
    ts_data = ds[key][var].groupby("time.year").mean('time')
    axs[idx_ax].plot(ts_data.year, ts_data, label = var)
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"]).groupby("time.year").mean('time')
ax_tmp.plot(ts_data.year, ts_data, 'k+', label = "FATES_FROOTC/FATES_LEAFC")
ax_tmp.legend(loc = 1)
ax_tmp.set_ylabel (f"Unitless")

# Subplot 6 <<<

# Subplot 7 >>>
idx_ax = 6
vars_plot = (
"""
FATES_NO3UPTAKE
FATES_NH4UPTAKE
"""
).split('\n')
vars_plot = vars_plot[1:-1]
var = vars_plot[0]
ts_data = ds[key][var].groupby("time.year").mean('time')
axs[idx_ax].plot(ts_data.year, ts_data, label = var)
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].groupby("time.year").mean('time')
ax_tmp.plot(ts_data.year,ts_data,'k+',  label = var)
ax_tmp.legend(loc = 5)
# Subplot 7 <<<

# Subplot 8 >>>
idx_ax = 7
vars_plot = (
"""
GROSS_NMIN
NET_NMIN
"""
).split('\n')
vars_plot = vars_plot[1:-1]
for i_var,var in enumerate(vars_plot):
    ts_data = ds[key][var].groupby("time.year").mean('time')
    axs[idx_ax].plot(ts_data.year, ts_data, label = var)
axs[idx_ax].set_ylabel (f"{ds[key][var].units}")
# Subplot 8 <<<

# Subplot 9 >>>
idx_ax = 8
vars_plot = (
"""
SMIN_NO3_LEACHED
DENIT
"""
).split('\n')
vars_plot = vars_plot[1:-1]
for i_var,var in enumerate(vars_plot):
    ts_data = ds[key][var].groupby("time.year").mean('time')
    axs[idx_ax].plot(ts_data.year, ts_data, label = var)
axs[idx_ax].set_ylabel (f"{ds[key][var].units}")
# Subplot 9 <<<


for i in range(len(axs)): 
    axs[i].legend()
    axs[i].axvline(x = logging_year, color = 'r',lw=5, alpha=.2)
    
fig.savefig('./Plots/' + f'{key}.pdf',bbox_inches='tight')
fig.savefig('./Plots/' + f'{key}.png',bbox_inches='tight')

All Plots >>>>>>

/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))

#print(breakit)

The next cell is for the carbon only simuations for the chosen site

Figure Set 1

print ("\nFigure Set 1 >>>>>>")
# For multiple variables on a same plot

# C-only
sim_conly = sims.replace("RD", "Conly")
if sims == "DUK_PIDE_RD_AgBgW_testing" : 
    sim_conly = "DUK_PIDE_Conly_AgBgW"
key=sim_conly
#key=sims
vars_plot = (
"""
FATES_GPP
FATES_AUTORESP
FATES_NPP
"""
).split('\n')
vars_plot = vars_plot[1:-1]



ymin = 9e20
ymax = -9e20
sum_ts = 0
for i_var,var in enumerate(vars_plot):
    ts_data = ds[key][var].groupby("time.year").mean('time')
    sum_ts= sum_ts + ts_data
    ts_data.plot(figsize=(20,3))
    if np.min(ts_data.values) < ymin:
        ymin = np.min(ts_data.values)
    if np.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)
    if i_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)
for i_var,var in enumerate(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)
#plt.close(fig)
    
fig1 = plt.figure(figsize=(20,9))
plt.title (f"Fractional Plot for simulation {key}", fontsize=15)
for i_var,var in enumerate(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)
    
#plt.close(fig1)
#plt.axvline(x = logging_year, color = 'r',lw=5, label = 'logging year', alpha=.2)

Figure Set 1 >>>>>>

<matplotlib.lines.Line2D at 0x17ec48760>

ts_data.year

<xarray.DataArray 'year' (year: 148)>
array([1850, 1851, 1852, 1853, 1854, 1855, 1856, 1857, 1858, 1859, 1860, 1861,
       1862, 1863, 1864, 1865, 1866, 1867, 1868, 1869, 1870, 1871, 1872, 1873,
       1874, 1875, 1876, 1877, 1878, 1879, 1880, 1881, 1882, 1883, 1884, 1885,
       1886, 1887, 1888, 1889, 1890, 1891, 1892, 1893, 1894, 1895, 1896, 1897,
       1898, 1899, 1900, 1901, 1902, 1903, 1904, 1905, 1906, 1907, 1908, 1909,
       1910, 1911, 1912, 1913, 1914, 1915, 1916, 1917, 1918, 1919, 1920, 1921,
       1922, 1923, 1924, 1925, 1926, 1927, 1928, 1929, 1930, 1931, 1932, 1933,
       1934, 1935, 1936, 1937, 1938, 1939, 1940, 1941, 1942, 1943, 1944, 1945,
       1946, 1947, 1948, 1949, 1950, 1951, 1952, 1953, 1954, 1955, 1956, 1957,
       1958, 1959, 1960, 1961, 1962, 1963, 1964, 1965, 1966, 1967, 1968, 1969,
       1970, 1971, 1972, 1973, 1974, 1975, 1976, 1977, 1978, 1979, 1980, 1981,
       1982, 1983, 1984, 1985, 1986, 1987, 1988, 1989, 1990, 1991, 1992, 1993,
       1994, 1995, 1996, 1997])
Coordinates:
  * year     (year) int64 1850 1851 1852 1853 1854 ... 1993 1994 1995 1996 1997

Rest of the results are for the chosen sim!

Figure Set 2

print ("\nFigure Set 2 >>>>>>")

# For multiple variables on a same plot

#sims = "ORN_PIDE_RD_AgBgW" # Difined at the top
key=sims

vars_plot = (
"""
FATES_GPP
FATES_AUTORESP
FATES_NPP
"""
).split('\n')
vars_plot = vars_plot[1:-1]



ymin = 9e20
ymax = -9e20
sum_ts = 0
for i_var,var in enumerate(vars_plot):
    ts_data = ds[key][var].groupby("time.year").mean('time')
    sum_ts= sum_ts + ts_data
    ts_data.plot(figsize=(20,3))
    if np.min(ts_data.values) < ymin:
        ymin = np.min(ts_data.values)
    if np.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)
    if i_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)
for i_var,var in enumerate(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)
for i_var,var in enumerate(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 2 >>>>>>

<matplotlib.lines.Line2D at 0x182a52610>

Figure Set 3

print ("\nFigure Set 3 >>>>>>")
# For multiple variables on a same plot

#sims = "ORN_PIDE_RD_AgBgW" # Difined at the top
key=sims

vars_plot = (
"""
FATES_STOREC
FATES_STOREC_TF
"""
).split('\n')
vars_plot = vars_plot[1:-1]



ymin = 9e20
ymax = -9e20
sum_ts = 0
for i_var,var in enumerate(vars_plot):
    ts_data = ds[key][var].groupby("time.year").mean('time')
    sum_ts= sum_ts + ts_data
    ts_data.plot(figsize=(20,3))
    if np.min(ts_data.values) < ymin:
        ymin = np.min(ts_data.values)
    if np.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)
    if i_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)
for i_var,var in enumerate(vars_plot):
    if var == "FATES_L2FR" : continue
    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)
for i_var,var in enumerate(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):
if True:
    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)
    
        

Figure Set 3 >>>>>>

<matplotlib.lines.Line2D at 0x17d6dd340>

Figure Set 4

print ("\nFigure Set 4 >>>>>>")
# For multiple variables on a same plot

#sims = "ORN_PIDE_RD_AgBgW" # Difined at the top
key=sims

vars_plot = (
"""
FATES_STOREN
FATES_STOREN_TF
"""
).split('\n')
vars_plot = vars_plot[1:-1]



ymin = 9e20
ymax = -9e20
sum_ts = 0
for i_var,var in enumerate(vars_plot):
    ts_data = ds[key][var].groupby("time.year").mean('time')
    sum_ts= sum_ts + ts_data
    ts_data.plot(figsize=(20,3))
    if np.min(ts_data.values) < ymin:
        ymin = np.min(ts_data.values)
    if np.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)
    if i_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)
for i_var,var in enumerate(vars_plot):
    if var == "FATES_L2FR" : continue
    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)
for i_var,var in enumerate(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):
if True:
    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")
    
        

Figure Set 4 >>>>>>

Figure Set 5

print ("\nFigure Set 5 >>>>>>")
# For multiple variables on a same plot

#sims = "ORN_PIDE_RD_AgBgW" # Difined at the top
key=sims

vars_plot = (
"""
FATES_STOREC_TF
FATES_STOREN_TF
"""
).split('\n')
vars_plot = vars_plot[1:-1]



ymin = 9e20
ymax = -9e20
sum_ts = 0
for i_var,var in enumerate(vars_plot):
    ts_data = ds[key][var].groupby("time.year").mean('time')
    sum_ts= sum_ts + ts_data
    ts_data.plot(figsize=(20,3))
    if np.min(ts_data.values) < ymin:
        ymin = np.min(ts_data.values)
    if np.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)
    if i_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)
for i_var,var in enumerate(vars_plot):
    if var == "FATES_L2FR" : continue
    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)
for i_var,var in enumerate(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):
if True:
    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):
if True:
    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")
    
        

Figure Set 5 >>>>>>

Figure Set 6

print ("\nFigure Set 6 >>>>>>")
# For multiple variables on a same plot

#sims = "ORN_PIDE_RD_AgBgW" # Difined at the top
key=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_ALLOC
FATES_FROOT_ALLOC
FATES_FROOT_ALLOC
FATES_SEED_ALLOC
FATES_STEM_ALLOC
FATES_STORE_ALLOC
"""
).split('\n')
vars_plot = vars_plot[1:-1]



ymin = 9e20
ymax = -9e20
sum_ts = 0
for i_var,var in enumerate(vars_plot):
    ts_data = ds[key][var].groupby("time.year").mean('time')
    sum_ts= sum_ts + ts_data
    ts_data.plot(figsize=(20,3))
    if np.min(ts_data.values) < ymin:
        ymin = np.min(ts_data.values)
    if np.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)
    if i_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)
for i_var,var in enumerate(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)
for i_var,var in enumerate(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):
if True:
    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
  return func(*(_execute_task(a, cache) for a in args))

<matplotlib.lines.Line2D at 0x17ce92910>

Figure Set 7

print ("\nFigure Set 7 >>>>>>")

# For multiple variables on a same plot

#sims = "ORN_PIDE_RD_AgBgW" # Difined at the top
key=sims

vars_plot = (
"""
FATES_LEAFC
FATES_FROOTC
"""
).split('\n')
vars_plot = vars_plot[1:-1]



ymin = 9e20
ymax = -9e20
sum_ts = 0
for i_var,var in enumerate(vars_plot):
    ts_data = ds[key][var].groupby("time.year").mean('time')
    sum_ts= sum_ts + ts_data
    ts_data.plot(figsize=(20,3))
    if np.min(ts_data.values) < ymin:
        ymin = np.min(ts_data.values)
    if np.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)
    if i_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)
for i_var,var in enumerate(vars_plot):
    if var == "FATES_L2FR" : continue
    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)
for i_var,var in enumerate(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):
if True:
    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)
    
        

Figure Set 7 >>>>>>

<matplotlib.lines.Line2D at 0x1824a0850>

Figure Set 8

print ("\nFigure Set 8 >>>>>>")

# For multiple variables on a same plot

#sims = "ORN_PIDE_RD_AgBgW" # Difined at the top
key=sims

vars_plot = (
"""
FATES_NO3UPTAKE
FATES_NH4UPTAKE
"""
).split('\n')
vars_plot = vars_plot[1:-1]



ymin = 9e20
ymax = -9e20
sum_ts = 0
for i_var,var in enumerate(vars_plot):
    ts_data = ds[key][var].groupby("time.year").mean('time')
    sum_ts= sum_ts + ts_data
    ts_data.plot(figsize=(20,3))
    if np.min(ts_data.values) < ymin:
        ymin = np.min(ts_data.values)
    if np.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)
    if i_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)
for i_var,var in enumerate(vars_plot):
    if var == "FATES_L2FR" : continue
    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)
for i_var,var in enumerate(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 >>>>>>

<matplotlib.lines.Line2D at 0x17dff1a60>

### *Figure Set 9*

print ("\nFigure Set 9 >>>>>>")

# For multiple variables on a same plot

#sims = "ORN_PIDE_RD_AgBgW" # Difined at the top
key=sims

vars_plot = (
"""
GROSS_NMIN
NET_NMIN
"""
).split('\n')
vars_plot = vars_plot[1:-1]



ymin = 9e20
ymax = -9e20
sum_ts = 0
for i_var,var in enumerate(vars_plot):
    ts_data = ds[key][var].groupby("time.year").mean('time')
    sum_ts= sum_ts + ts_data
    ts_data.plot(figsize=(20,3))
    if np.min(ts_data.values) < ymin:
        ymin = np.min(ts_data.values)
    if np.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)
    if i_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)
for i_var,var in enumerate(vars_plot):
    if var == "FATES_L2FR" : continue
    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)
for i_var,var in enumerate(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 9 >>>>>>

<matplotlib.lines.Line2D at 0x18b906190>

Figure Set 10

print ("\nFigure Set 10 >>>>>>")

# For multiple variables on a same plot

#sims = "ORN_PIDE_RD_AgBgW" # Difined at the top
key=sims

vars_plot = (
"""
SMIN_NO3_LEACHED
DENIT
"""
).split('\n')
vars_plot = vars_plot[1:-1]



ymin = 9e20
ymax = -9e20
sum_ts = 0
for i_var,var in enumerate(vars_plot):
    ts_data = ds[key][var].groupby("time.year").mean('time')
    sum_ts= sum_ts + ts_data
    ts_data.plot(figsize=(20,3))
    if np.min(ts_data.values) < ymin:
        ymin = np.min(ts_data.values)
    if np.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)
    if i_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)
for i_var,var in enumerate(vars_plot):
    if var == "FATES_L2FR" : continue
    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)
for i_var,var in enumerate(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 >>>>>>

<matplotlib.lines.Line2D at 0x189591ee0>

# For multiple variables on a same plot

#sims = "ORN_PIDE_RD_AgBgW" # Difined at the top
key=sims

vars_plot = (
"""
FATES_FROOTC
FATES_LEAFC
FATES_NONSTRUCTC
FATES_REPROC
FATES_SAPWOODC
FATES_STOREC
FATES_STRUCTC
FATES_VEGC
"""
).split('\n')
vars_plot = vars_plot[1:-1]



ymin = 9e20
ymax = -9e20
sum_ts = 0
for i_var,var in enumerate(vars_plot):
    print(var)
    ts_data = ds[key][var].groupby("time.year").mean('time')
    sum_ts= sum_ts + ts_data
    ts_data.plot(figsize=(20,3))
    if np.min(ts_data.values) < ymin:
        ymin = np.min(ts_data.values)
    if np.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)
    if i_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)
for i_var,var in enumerate(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)
for i_var,var in enumerate(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_FROOTC
FATES_LEAFC
FATES_NONSTRUCTC
FATES_REPROC
FATES_SAPWOODC
FATES_STOREC
FATES_STRUCTC
FATES_VEGC

<matplotlib.lines.Line2D at 0x17fdc1f40>

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.

Subplot 1: FATES_NPP (for C-only at that site and RD/ECA) Subplot 1: FATES_NPP Subplot 2: FATES_STOREC Subplot 2: FATES_STOREC_TF Subplot 3: FATES_STOREN Subplot 3: FATES_STOREN_TF Subplot 4: FATES_STOREC_TF / FATES_STOREN_TF Subplot 5: FATES_FROOTC_ALLOC / FATES_LEAFC_ALLOC Subplot 5: FATES_L2FR Subplot 6: FATES_FROOTC / FATES_LEAFC Subplot 6: FATES_FROOTC Subplot 6: FATES_LEAFC Subplot 7: FATES_NO3UPTAKE Subplot 7: FATES_NH4UPTAKE Subplot 8: NET_NMIN Subplot 8: GROSS_NMIN Subplot 9: SMIN_NO3_LEACHED Subplot 9: DENIT

print ("\nAll Plots >>>>>>")
# For multiple variables on a same plot
fig, axs = plt.subplots(9,1 , figsize=(15,30), dpi=400)

# Subplot 1 >>>
idx_ax = 0
var = "FATES_NPP"

# C-only
sim_conly = sims.replace("RD", "Conly")
if sims == "DUK_PIDE_RD_AgBgW_testing" : 
    sim_conly = "DUK_PIDE_Conly_AgBgW"
key=sim_conly
ts_data = ds[key][var].groupby("time.year").mean('time')
axs[idx_ax].plot(ts_data.year, ts_data, label = var + " C-only")
# RD
key=sims
ts_data = ds[key][var].groupby("time.year").mean('time')
axs[idx_ax].plot(ts_data.year, ts_data, label = var )
axs[idx_ax].set_ylabel (f"{ds[key][var].units}")
axs[idx_ax].set_title(f"{key}\n",fontsize=16)

# Subplot 1 <<<

# Subplot 2 >>>
idx_ax = 1
vars_plot = (
"""
FATES_STOREC_TF
FATES_STOREC
"""
).split('\n')
vars_plot = vars_plot[1:-1]
for i_var,var in enumerate(vars_plot):
    ts_data = ds[key][var].groupby("time.year").mean('time')
    axs[idx_ax].plot(ts_data.year, ts_data, label = var)
axs[idx_ax].set_ylabel (f"{ds[key][var].units}")
# Subplot 2 <<<

# Subplot 3 >>>
idx_ax = 2
vars_plot = (
"""
FATES_STOREN_TF
FATES_STOREN
"""
).split('\n')
vars_plot = vars_plot[1:-1]
var = vars_plot[0]
ts_data = ds[key][var].groupby("time.year").mean('time')
axs[idx_ax].plot(ts_data.year, ts_data, label = var)
ax_tmp = axs[idx_ax].twinx()
var = vars_plot[1]
ts_data = ds[key][var].groupby("time.year").mean('time')
ax_tmp.plot(ts_data.year,ts_data,'k+',  label = var)
ax_tmp.legend(loc = 5)
axs[idx_ax].set_ylabel (f"{ds[key][var].units}")
# Subplot 3 <<<

# Subplot 4 >>>
idx_ax = 3
ts_data = (ds[key]["FATES_STOREC_TF"]/ds[key]["FATES_STOREN_TF"]).groupby("time.year").mean('time')
axs[idx_ax].plot(ts_data.year, ts_data, label = "FATES_STOREC_TF/FATES_STOREN_TF")
ax_tmp = axs[idx_ax].twinx()
ts_data = (ds[key]["FATES_STOREC"]/ds[key]["FATES_STOREN"]).groupby("time.year").mean('time')
ax_tmp.plot(ts_data.year, ts_data, 'k+', label = "FATES_STOREC/FATES_STOREN")
ax_tmp.legend(loc = 5)
# Subplot 4 <<<

# Subplot 5 >>>
idx_ax = 4
ts_data = (ds[key]["FATES_FROOT_ALLOC"]/ds[key]["FATES_LEAF_ALLOC"]).groupby("time.year").mean('time')
axs[idx_ax].plot(ts_data.year, ts_data, label = "FATES_FROOT_ALLOC/FATES_LEAF_ALLOC")
ts_data = ds[key]["FATES_L2FR"].groupby("time.year").mean('time')
axs[idx_ax].plot(ts_data.year, ts_data, label = "FATES_L2FR")
# Subplot 5 <<<

# Subplot 6 >>>
idx_ax = 5
vars_plot = (
"""
FATES_LEAFC
FATES_FROOTC
"""
).split('\n')
vars_plot = vars_plot[1:-1]
for i_var,var in enumerate(vars_plot):
    ts_data = ds[key][var].groupby("time.year").mean('time')
    axs[idx_ax].plot(ts_data.year, ts_data, label = var)
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"]).groupby("time.year").mean('time')
ax_tmp.plot(ts_data.year, ts_data, 'k+', label = "FATES_FROOTC/FATES_LEAFC")
ax_tmp.legend(loc = 1)

# Subplot 6 <<<

# Subplot 7 >>>
idx_ax = 6
vars_plot = (
"""
FATES_NO3UPTAKE
FATES_NH4UPTAKE
"""
).split('\n')
vars_plot = vars_plot[1:-1]
var = vars_plot[0]
ts_data = ds[key][var].groupby("time.year").mean('time')
axs[idx_ax].plot(ts_data.year, ts_data, label = var)
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].groupby("time.year").mean('time')
ax_tmp.plot(ts_data.year,ts_data,'k+',  label = var)
ax_tmp.legend(loc = 5)
# Subplot 7 <<<

# Subplot 8 >>>
idx_ax = 7
vars_plot = (
"""
GROSS_NMIN
NET_NMIN
"""
).split('\n')
vars_plot = vars_plot[1:-1]
for i_var,var in enumerate(vars_plot):
    ts_data = ds[key][var].groupby("time.year").mean('time')
    axs[idx_ax].plot(ts_data.year, ts_data, label = var)
axs[idx_ax].set_ylabel (f"{ds[key][var].units}")
# Subplot 8 <<<

# Subplot 9 >>>
idx_ax = 8
vars_plot = (
"""
SMIN_NO3_LEACHED
DENIT
"""
).split('\n')
vars_plot = vars_plot[1:-1]
for i_var,var in enumerate(vars_plot):
    ts_data = ds[key][var].groupby("time.year").mean('time')
    axs[idx_ax].plot(ts_data.year, ts_data, label = var)
axs[idx_ax].set_ylabel (f"{ds[key][var].units}")
# Subplot 9 <<<


for i in range(len(axs)): 
    axs[i].legend()
    axs[i].axvline(x = logging_year, color = 'r',lw=5, alpha=.2)
    
fig.savefig('./Plots/' + f'{key}.pdf',bbox_inches='tight')
fig.savefig('./Plots/' + f'{key}.png',bbox_inches='tight')

All Plots >>>>>>

/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))

pwd

'/Users/ud4/repos/GitHub/FATESFACE'

import datetime

current_time = datetime.datetime.now()
formatted_time = current_time.strftime("%H:%M:%S")
print("Formatted time:", formatted_time)

Formatted time: 13:55:15

print (f"Successful run at {formatted_time} {key}")

Successful run at 13:55:15 DUK_PIDE_RD_AgBgW_testing

Ignore the plots after this cell!