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_PIDC_RD_AgBgW

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

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"

PIDC

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.0009045469760894776

#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])

-0.00016459823 -0.0008314407736063004

# 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.0001
 pid_kd : 0.01
 pid_ki : 0.0 
pid_kp : 0.0001
 pid_kd : 0.01
 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/*PIDC*l2fr.txt'

filenames_sim_out

['/Users/ud4/FATESMDS_analysis/outputs/runs/tests_alp/230309/PIDC_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.212186	0.414966	0.808814	-352.915922	-0.000145	-0.000036	-0.212186	-2.158065e-05	0.473836	...	12.570138	0.006023	0.014541	0.01	0.094106	0.116567	0.006023	0.014541	17.296259	0
1	1850-01-01	0.004225	0.409707	1.004234	1.716733	-0.000058	0.000047	0.004225	8.876873e-07	0.515037	...	13.315416	0.006367	0.015600	0.01	0.099478	0.099445	0.006367	0.015600	18.046229	0
2	1850-01-02	-0.211983	0.143543	0.808979	-353.127904	0.000203	-0.000024	-0.211983	-2.144075e-05	0.473815	...	12.573126	0.006023	0.042046	0.01	0.094111	0.116572	0.006023	0.042046	17.297046	1
3	1850-01-02	0.003892	0.135507	1.003900	1.720625	-0.000333	0.000028	0.003892	6.645410e-07	0.515037	...	13.315448	0.006367	0.047167	0.01	0.099479	0.099479	0.006367	0.047167	18.046353	1
4	1850-01-03	0.000616	0.126978	1.000616	0.000616	0.000616	0.000031	0.000616	3.696501e-07	0.510021	...	0.076056	0.000036	0.000288	0.01	0.000570	0.000570	0.000036	0.000288	0.390669	2
5	1850-01-03	-0.215707	0.143058	0.805972	-353.343611	-0.003724	-0.000209	-0.215707	-2.366297e-05	0.473791	...	12.532482	0.006023	0.042050	0.01	0.094111	0.116621	0.006023	0.042050	17.297046	2
6	1850-01-03	0.000805	0.135103	1.000805	1.721429	-0.003087	-0.000128	0.000805	-1.201721e-06	0.515036	...	13.274516	0.006367	0.047165	0.01	0.099485	0.099485	0.006367	0.047165	18.047260	2
7	1850-01-04	0.003652	0.127349	1.003658	0.003652	0.003652	0.000183	0.003652	2.191040e-06	0.510022	...	0.076259	0.000036	0.000288	0.01	0.000570	0.000570	0.000036	0.000288	0.390739	3
8	1850-01-04	0.003883	0.113069	1.003890	0.004499	0.003267	0.000193	0.003883	2.314359e-06	0.510023	...	0.076292	0.000037	0.000324	0.01	0.000570	0.000570	0.000037	0.000324	0.390870	3
9	1850-01-04	-0.212930	0.143515	0.808213	-353.550876	0.002774	-0.000060	-0.212930	-2.189389e-05	0.473770	...	12.569492	0.006023	0.042040	0.01	0.094111	0.116642	0.006023	0.042040	17.297075	3
10	1850-01-04	0.003675	0.135491	1.003681	1.725104	0.002870	0.000022	0.003675	5.845066e-07	0.515037	...	13.313529	0.006367	0.047168	0.01	0.099492	0.099492	0.006367	0.047168	18.048165	3
11	1850-01-05	0.003485	0.127328	1.003491	0.003485	0.003485	0.000174	0.003485	2.091154e-06	0.510022	...	0.076246	0.000036	0.000288	0.01	0.000570	0.000570	0.000036	0.000288	0.390741	4
12	1850-01-05	0.003485	0.112997	1.003491	0.007309	-0.000214	0.000165	0.003485	1.995468e-06	0.510024	...	0.076286	0.000037	0.000324	0.01	0.000570	0.000570	0.000037	0.000324	0.390898	4
13	1850-01-05	-0.213074	0.143493	0.808096	-353.758336	-0.000148	-0.000064	-0.213074	-2.195215e-05	0.473749	...	12.568642	0.006023	0.042045	0.01	0.094115	0.116655	0.006023	0.042045	17.297652	4
14	1850-01-05	0.003502	0.135468	1.003508	1.728606	-0.000173	0.000012	0.003502	4.698207e-07	0.515037	...	13.312080	0.006368	0.047171	0.01	0.099498	0.099498	0.006368	0.047171	18.049082	4
15	1850-01-06	0.004162	0.127415	1.004171	0.004162	0.004162	0.000208	0.004162	2.497379e-06	0.510023	...	0.076298	0.000036	0.000288	0.01	0.000570	0.000570	0.000036	0.000288	0.390741	5
16	1850-01-06	0.004162	0.113073	1.004171	0.009692	0.000677	0.000195	0.004162	2.361528e-06	0.510026	...	0.076352	0.000037	0.000324	0.01	0.000571	0.000571	0.000037	0.000324	0.390958	5
17	1850-01-06	-0.212482	0.143582	0.808575	-353.964948	0.000588	-0.000032	-0.212482	-2.156673e-05	0.473728	...	12.576053	0.006024	0.042045	0.01	0.094121	0.116658	0.006024	0.042045	17.298532	5
18	1850-01-06	0.004255	0.135570	1.004265	1.732861	0.000754	0.000049	0.004255	9.161100e-07	0.515038	...	13.322986	0.006368	0.047175	0.01	0.099505	0.099505	0.006368	0.047175	18.050007	5
19	1850-01-07	0.003223	0.127296	1.003228	0.003223	0.003223	0.000161	0.003223	1.933707e-06	0.510022	...	0.076226	0.000036	0.000288	0.01	0.000570	0.000570	0.000036	0.000288	0.390743	6
20	1850-01-07	0.003222	0.112968	1.003228	0.010693	-0.000940	0.000143	0.003222	1.752321e-06	0.510026	...	0.076292	0.000037	0.000324	0.01	0.000571	0.000571	0.000037	0.000324	0.391005	6
21	1850-01-07	-0.213245	0.143486	0.807958	-354.172182	-0.000766	-0.000069	-0.213245	-2.201013e-05	0.473706	...	12.568605	0.006024	0.042047	0.01	0.094125	0.116674	0.006024	0.042047	17.299050	6
22	1850-01-07	0.002996	0.135399	1.003000	1.735857	-0.001260	-0.000016	0.002996	1.356942e-07	0.515038	...	13.307086	0.006369	0.047178	0.01	0.099511	0.099511	0.006369	0.047178	18.050941	6
23	1850-01-08	-0.212742	0.143567	0.808364	-354.378717	0.000499	-0.000040	-0.212742	-2.167617e-05	0.473685	...	12.575341	0.006024	0.042047	0.01	0.094131	0.116681	0.006024	0.042047	17.299843	7
24	1850-01-08	0.002983	0.135398	1.002988	1.738840	-0.000013	-0.000016	0.002983	1.363489e-07	0.515038	...	13.307801	0.006369	0.047181	0.01	0.099518	0.099518	0.006369	0.047181	18.051879	7
25	1850-01-09	-0.212307	0.143626	0.808717	-354.584752	0.000432	-0.000017	-0.212307	-2.139657e-05	0.473665	...	12.580515	0.006025	0.042047	0.01	0.094137	0.116680	0.006025	0.042047	17.300831	8
26	1850-01-09	0.003915	0.135524	1.003922	1.742754	0.000931	0.000031	0.003915	7.033159e-07	0.515039	...	13.321089	0.006370	0.047184	0.01	0.099524	0.099524	0.006370	0.047184	18.052818	8
27	1850-01-10	-0.212255	0.143626	0.808759	-354.790693	0.000048	-0.000013	-0.212255	-2.135902e-05	0.473644	...	12.581122	0.006025	0.042049	0.01	0.094144	0.116679	0.006025	0.042049	17.301804	9
28	1850-01-10	0.004839	0.135649	1.004851	1.747593	0.000924	0.000076	0.004839	1.242277e-06	0.515040	...	13.334296	0.006370	0.047187	0.01	0.099531	0.099531	0.006370	0.047187	18.053760	9
29	1850-01-11	-0.213210	0.143501	0.807987	-354.997547	-0.000959	-0.000061	-0.213210	-2.192723e-05	0.473623	...	12.571440	0.006025	0.042052	0.01	0.094148	0.116697	0.006025	0.042052	17.302314	10

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: 30

plt.figure(figsize=(15,6))
plt.plot(np.unique(df_vals['timex'], return_counts=1)[0], np.unique(df_vals['timex_code'], return_counts=1)[1])
plt.ylabel("Number of cohorts", fontsize =14)
plt.axvline(x = dt.date(1855,1,1), color = 'r',lw=5, alpha=.2)
plt.title( "Cohorts Vs Time (daily)\n", fontsize  = 16)

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: 10

# 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: 27

# 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 0x17d643280>

#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 0x182b7b2e0>

#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 0x182e34700>

#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 0x17fa942e0>

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 0x1822cc310>

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 0x17f6592b0>

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_PIDC_RD_AgBgW

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

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 0x181cf7ee0>

ts_data.year

<xarray.DataArray 'year' (year: 146)>
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])
Coordinates:
  * year     (year) int64 1850 1851 1852 1853 1854 ... 1991 1992 1993 1994 1995

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 0x181cf7310>

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 0x1829c97c0>

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 0x181302c40>

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 0x1806620d0>

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 0x18027bb50>

### *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 0x18061ca00>

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 0x1811206a0>

# 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 0x17e7edaf0>

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: 14:01:12

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

Successful run at 14:01:12 DUK_PIDC_RD_AgBgW

Ignore the plots after this cell!