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_PIDG_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_PIDG_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"

PIDG

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

#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.00018634275 0.0002601458746008575

# 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.005
 pid_kd : 0.5
 pid_ki : 0.0 
pid_kp : 0.005
 pid_kd : 0.5
 pid_ki : 0.0 

Printing Variables that go in L2FR calculations

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

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

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

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

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

filenames_sim_out

['/Users/ud4/FATESMDS_analysis/outputs/runs/tests_alp/230309/PIDG_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 53290

	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.006661	0.567606	1.006684	0.006661	0.006661	0.000333	0.006661	0.000200	0.687673	...	0.076460	0.000036	0.000065	0.01	0.000570	0.000570	0.000036	0.000065	0.390608	0
1	1850-01-01	0.055128	0.366171	1.056676	0.570345	0.000751	0.002389	0.055128	0.001470	0.693239	...	0.082994	0.000040	0.000109	0.01	0.000619	0.000589	0.000040	0.000109	0.404701	0
2	1850-01-01	0.054447	0.366862	1.055957	2.994859	0.000448	0.002314	0.054447	0.001429	0.696358	...	0.097365	0.000046	0.000127	0.01	0.000726	0.000692	0.000046	0.000127	0.454223	0
3	1850-01-01	0.052606	0.369313	1.054014	6.943459	0.000430	0.002236	0.052606	0.001381	0.701776	...	0.121405	0.000058	0.000158	0.01	0.000905	0.000864	0.000058	0.000158	0.525736	0
4	1850-01-01	0.050860	0.371754	1.052175	10.818777	0.000413	0.002163	0.050860	0.001336	0.706703	...	0.149802	0.000072	0.000193	0.01	0.001117	0.001068	0.000072	0.000193	0.603827	0
5	1850-01-01	0.049195	0.374293	1.050426	14.463816	0.000397	0.002092	0.049195	0.001292	0.710440	...	0.182888	0.000087	0.000235	0.01	0.001364	0.001306	0.000087	0.000235	0.688006	0
6	1850-01-01	0.047590	0.376822	1.048740	18.354883	0.000383	0.002025	0.047590	0.001250	0.714003	...	0.221559	0.000106	0.000282	0.01	0.001653	0.001584	0.000106	0.000282	0.779220	0
7	1850-01-01	0.046025	0.379426	1.047100	22.281203	0.000370	0.001959	0.046025	0.001209	0.716954	...	0.266642	0.000127	0.000337	0.01	0.001990	0.001910	0.000127	0.000337	0.877790	0
8	1850-01-01	0.044631	0.382143	1.045642	26.360688	0.000342	0.001895	0.044631	0.001170	0.719489	...	0.319381	0.000153	0.000401	0.01	0.002383	0.002291	0.000153	0.000401	0.984702	0
9	1850-01-01	0.005353	0.382341	1.005368	30.711033	-0.000403	0.000035	0.005353	0.000045	0.718935	...	0.371038	0.000177	0.000466	0.01	0.002769	0.002768	0.000177	0.000466	1.102583	0
10	1850-01-01	0.005448	0.381409	1.005463	33.842410	-0.000403	0.000037	0.005448	0.000046	0.711632	...	0.444401	0.000212	0.000559	0.01	0.003316	0.003315	0.000212	0.000559	1.289952	0
11	1850-01-01	0.005361	0.383091	1.005375	37.943491	-0.000383	0.000037	0.005361	0.000045	0.708947	...	0.521534	0.000249	0.000653	0.01	0.003892	0.003891	0.000249	0.000653	1.444714	0
12	1850-01-01	0.005421	0.383341	1.005436	42.072525	-0.000374	0.000038	0.005421	0.000046	0.696998	...	0.606267	0.000290	0.000759	0.01	0.004524	0.004522	0.000290	0.000759	1.643261	0
13	1850-01-01	0.005402	0.384680	1.005416	51.483147	-0.000352	0.000039	0.005402	0.000047	0.686176	...	0.762162	0.000364	0.000951	0.01	0.005688	0.005685	0.000364	0.000951	1.965826	0
14	1850-01-01	0.005288	0.396773	1.005302	42.534206	-0.000280	0.000042	0.005288	0.000047	0.609192	...	0.868728	0.000415	0.001051	0.01	0.006483	0.006481	0.000415	0.001051	2.142658	0
15	1850-01-01	0.004446	0.413444	1.004455	0.055925	-0.000092	0.000050	0.004446	0.000047	0.488974	...	8.505434	0.004066	0.009875	0.01	0.063529	0.063508	0.004066	0.009875	12.746414	0
16	1850-01-02	0.006985	0.098021	1.007010	0.006985	0.006985	0.000349	0.006985	0.000210	0.687718	...	0.076514	0.000036	0.000375	0.01	0.000570	0.000570	0.000036	0.000375	0.390815	1
17	1850-01-02	0.054063	0.091456	1.055551	0.620221	-0.000706	0.002220	0.054063	0.001380	0.694578	...	0.082947	0.000040	0.000435	0.01	0.000618	0.000589	0.000040	0.000435	0.404597	1
18	1850-01-02	0.053735	0.092001	1.055205	3.048594	-0.000712	0.002163	0.053735	0.001350	0.697708	...	0.097368	0.000046	0.000508	0.01	0.000726	0.000692	0.000046	0.000508	0.454223	1
19	1850-01-02	0.051906	0.093929	1.053276	6.995365	-0.000700	0.002090	0.051906	0.001304	0.703080	...	0.121408	0.000058	0.000620	0.01	0.000905	0.000865	0.000058	0.000620	0.525736	1
20	1850-01-02	0.050173	0.095901	1.051453	10.868950	-0.000687	0.002020	0.050173	0.001261	0.707964	...	0.149805	0.000072	0.000750	0.01	0.001117	0.001069	0.000072	0.000750	0.603827	1
21	1850-01-02	0.048524	0.098010	1.049721	14.512340	-0.000671	0.001954	0.048524	0.001220	0.711660	...	0.182893	0.000087	0.000896	0.01	0.001364	0.001307	0.000087	0.000896	0.688006	1
22	1850-01-02	0.046934	0.100173	1.048053	18.401817	-0.000655	0.001891	0.046934	0.001180	0.715183	...	0.221564	0.000106	0.001062	0.01	0.001653	0.001586	0.000106	0.001062	0.779220	1
23	1850-01-02	0.045387	0.102467	1.046433	22.326590	-0.000638	0.001829	0.045387	0.001141	0.718095	...	0.266648	0.000127	0.001249	0.01	0.001990	0.001911	0.000127	0.001249	0.877790	1
24	1850-01-02	0.044018	0.104905	1.045001	26.404706	-0.000614	0.001769	0.044018	0.001105	0.720594	...	0.319390	0.000153	0.001461	0.01	0.002383	0.002292	0.000153	0.001461	0.984702	1
25	1850-01-02	0.005288	0.105315	1.005302	30.716322	-0.000065	0.000030	0.005288	0.000042	0.718977	...	0.371150	0.000177	0.001692	0.01	0.002770	0.002770	0.000177	0.001692	1.102820	1
26	1850-01-02	0.005382	0.104448	1.005396	33.847792	-0.000066	0.000031	0.005382	0.000043	0.711675	...	0.444536	0.000212	0.002043	0.01	0.003317	0.003317	0.000212	0.002043	1.290229	1
27	1850-01-02	0.005282	0.105995	1.005296	37.948773	-0.000079	0.000031	0.005282	0.000042	0.708989	...	0.521684	0.000249	0.002363	0.01	0.003893	0.003893	0.000249	0.002363	1.445006	1
28	1850-01-02	0.005335	0.106208	1.005350	42.077860	-0.000086	0.000032	0.005335	0.000043	0.697041	...	0.606437	0.000290	0.002741	0.01	0.004525	0.004525	0.000290	0.002741	1.643583	1
29	1850-01-02	0.005300	0.107452	1.005314	51.488448	-0.000101	0.000032	0.005300	0.000043	0.686219	...	0.762364	0.000364	0.003406	0.01	0.005689	0.005689	0.000364	0.003406	1.966184	1

30 rows × 21 columns

Number of Cohorts over time

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

Maximum Cohorts: 32

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

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

FATES_NPP just for reference

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

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

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

Plot of L2FR from print file

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

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

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

Maximum Cohorts: 25

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

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

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

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

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

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

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

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_PIDG_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_PIDG_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 0x17b253340>

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

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

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

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

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

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

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

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

FATES_NPP (for C-only at that site and RD/ECA) [done] FATES_NPP [done] FATES_STOREC FATES_STOREC_TF FATES_STOREN FATES_STOREN_TF FATES_STOREC_TF / FATES_STOREN_TF FATES_L2FR FATES_FROOTC_ALLOC / FATES_LEAFC_ALLOC FATES_FROOTC / FATES_LEAFC FATES_FROOTC FATES_LEAFC FATES_NO3UPTAKE FATES_NH4UPTAKE NET_NMIN GROSS_NMIN It's also be good to look at N leaching (NLEACH?) and denitrification (DENIT?). Not sure those are the variable names but the variables (however named) should be in the output files I hope. And let's take a look at PID-D for both Duke and ORNL, and PID-B for Duke.

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

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

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

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

# Subplot 1 <<<

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

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

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

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

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

# Subplot 6 <<<

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

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

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


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

All Plots >>>>>>

/Users/ud4/opt/anaconda3/envs/pyces/lib/python3.9/site-packages/dask/core.py:121: RuntimeWarning: invalid value encountered in true_divide
  return func(*(_execute_task(a, cache) for a in args))

pwd

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

import datetime

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

Formatted time: 13:30:11

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

Successful run at 13:30:11 DUK_PIDG_RD_AgBgW

Ignore the plots after this cell!