before coffee

OH2 1 scale reference.py

@@ -1,8 +1,8 @@
-# This code is used to scale the reference gases
+# The following code is used to scale the reference gases
 # based on the "accepted" values for IAEA-603 and NBS-18
 # from Wostbrock et al. (2020)
 
-# INPUT: OH2 Table S1.csv
+# INPUT: OH2 Table S2.csv
 # OUTPUT: OH2 Figure S1.png
 
 # >>>>>>>>>
@@ -14,16 +14,16 @@
 import matplotlib.pyplot as plt
 
 # Plot parameters
-plt.rcParams["legend.loc"] = "best"
 plt.rcParams.update({'font.size': 7})
 plt.rcParams['scatter.edgecolors'] = "k"
 plt.rcParams['scatter.marker'] = "o"
 plt.rcParams["lines.linewidth"] = 0.5
 plt.rcParams["patch.linewidth"] = 0.5
 plt.rcParams["figure.figsize"] = (9, 4)
-plt.rcParams["savefig.dpi"] = 600
+plt.rcParams["savefig.dpi"] = 800
 plt.rcParams["savefig.bbox"] = "tight"
 
+
 # Functions that make life easier
 def prime(delta):
     dprime = 1000 * np.log(delta/1000 + 1)
@@ -36,40 +36,37 @@ def unprime(dprime):
 
 
 def Dp17O(d17O, d18O):
-    # value in ppm
     return (prime(d17O) - 0.528 * prime(d18O)) * 1000
 
 
 def to_VSMOW(VPDB):
     return (VPDB*1.03092)+30.92
 
 
-# Read data from Table S1
-df = pd.read_csv(sys.path[0] + "/" + "OH2 Table S1.csv")
-df = df[df["MeasurementPeriod"].str.contains("2023-09-20 to 2023-10-08")]
-df = df[df["SampleName"].str.contains("light|heavy|NBS18|IAEA")]
+def print_info(df, d18O_col, Dp17O_col, sample_name):
+    gas_subset = df[df["SampleName"].str.contains(sample_name)].copy()
+    d18O_mean = gas_subset[d18O_col].mean()
+    d18O_std = gas_subset[d18O_col].std()
+    Dp17O_mean = gas_subset[Dp17O_col].mean()
+    Dp17O_std = gas_subset[Dp17O_col].std()
+    N_gas = len(gas_subset)
+    print(f"{sample_name}, N = {N_gas}, d18O = {d18O_mean:.3f}(±{d18O_std:.3f})‰, ∆'17O = {Dp17O_mean:.0f}(±{Dp17O_std:.0f}) ppm", end="")
 
-heavy_d18O_SD = df[df["SampleName"].str.contains("heavy")]["d18O"].std()
-heavy_Dp17O_SD = df[df["SampleName"].str.contains("heavy")]["Dp17O"].std()
-N_heavy = len(df[df["SampleName"].str.contains("heavy")])
-print(f"\nHeavy reference gas, N = {N_heavy}, d18O SD = ±{heavy_d18O_SD:.3f}‰, ∆'17O = ±{heavy_Dp17O_SD:.0f} ppm")
 
-light_d18O_SD = df[df["SampleName"].str.contains("light")]["d18O"].std()
-light_Dp17O_SD = df[df["SampleName"].str.contains("light")]["Dp17O"].std()
-N_light = len(df[df["SampleName"].str.contains("light")])
-print(f"Light reference gas, N = {N_light}, d18O: ±{light_d18O_SD:.3f}‰, ∆'17O: ±{light_Dp17O_SD:.0f} ppm")
+# Read data from Table S1
+df = pd.read_csv(sys.path[0] + "/" + "OH2 Table S2.csv")
+df = df[df["measurementPeriod"].str.contains("2023-09-20 to 2023-10-08")]
+df = df[df["SampleName"].str.contains("light|heavy|NBS18|IAEA")]
 
-IAEA603_d18O_SD = df[df["SampleName"].str.contains("IAEA603")]["d18O"].std()
-IAEA603_Dp17O_SD = df[df["SampleName"].str.contains("IAEA603")]["Dp17O"].std()
-N_IAEA603 = len(df[df["SampleName"].str.contains("IAEA603")])
-print(f"IAEA-603, N = {N_IAEA603}, d18O: ±{IAEA603_d18O_SD:.3f}‰, ∆'17O: ±{IAEA603_Dp17O_SD:.0f} ppm")
+print("\nNumber of replicates and standard deviations for the references:")
+print_info(df, "d18O", "Dp17O", "light"); print("\t<--- unscaled")
+print_info(df, "d18O", "Dp17O", "heavy"); print("\t<--- unscaled")
+print_info(df, "d18O", "Dp17O", "IAEA603"); print("\t<--- unscaled")
+print_info(df, "d18O", "Dp17O", "NBS18"); print("\t<--- unscaled")
 
-NBS18_d18O_SD = df[df["SampleName"].str.contains("NBS18")]["d18O"].std()
-NBS18_Dp17O_SD = df[df["SampleName"].str.contains("NBS18")]["Dp17O"].std()
-N_NBS18 = len(df[df["SampleName"].str.contains("NBS18")])
-print(f"NBS-18, N = {N_NBS18}, d18O: ±{NBS18_d18O_SD:.3f}‰, ∆'17O: ±{NBS18_Dp17O_SD:.0f} ppm")
+# Scale values
 
-# Do the scaling here
+# Measured CO2 values
 IAEA603_d18O_measured = df[df["SampleName"].str.contains("IAEA603")]["d18O"].mean()
 IAEA603_d17O_measured = df[df["SampleName"].str.contains("IAEA603")]["d17O"].mean()
 
@@ -80,45 +77,46 @@ def to_VSMOW(VPDB):
 IAEA603_d18O_accepted = (to_VSMOW(-2.37) + 1000) * 1.01025 - 1000
 IAEA603_Dp17O_accepted = -147 # from Wostbrock et al. (2020)
 IAEA603_d17O_accepted = unprime(IAEA603_Dp17O_accepted/1000 + 0.528 * prime(IAEA603_d18O_accepted))
-print(f"\nAccepted for IAEA603 d18O = {IAEA603_d18O_accepted:.3f}‰, ∆'17O = {IAEA603_Dp17O_accepted:.0f} ppm")
+print("\nAccepted values:")
+print(f"IAEA-603: d18O = {IAEA603_d18O_accepted:.3f}‰, ∆'17O = {IAEA603_Dp17O_accepted:.0f} ppm")
 
 NBS18_d18O_accepted = (to_VSMOW(-23.2) + 1000) * 1.01025 - 1000
 NBS18_Dp17O_accepted = -100 # from Wostbrock et al. (2020)
 NBS18_d17O_accepted = unprime(NBS18_Dp17O_accepted/1000 + 0.528 * prime(NBS18_d18O_accepted))
-print(f"Accepted for NBS18 d18O = {NBS18_d18O_accepted:.3f}‰, ∆'17O = {NBS18_Dp17O_accepted:.0f} ppm")
+print(f"NBS-18: d18O = {NBS18_d18O_accepted:.3f}‰, ∆'17O = {NBS18_Dp17O_accepted:.0f} ppm")
 
+# Calculate the scaling factors
 slope_d18O = (NBS18_d18O_accepted - IAEA603_d18O_accepted) / (NBS18_d18O_measured - IAEA603_d18O_measured)
 intercept_d18O = IAEA603_d18O_accepted - slope_d18O * IAEA603_d18O_measured
 
 slope_d17O = (NBS18_d17O_accepted - IAEA603_d17O_accepted) / (NBS18_d17O_measured - IAEA603_d17O_measured)
 intercept_d17O = IAEA603_d17O_accepted - slope_d17O * IAEA603_d17O_measured
 
+# Scale the measured values
 df["d18O_scaled"] = slope_d18O*df['d18O']+intercept_d18O
 df["d17O_scaled"] = slope_d17O*df['d17O']+intercept_d17O
-df["Dp17O_scaled"] = (prime(df["d17O_scaled"]) - 0.528 * prime(df["d18O_scaled"]))*1000
+df["Dp17O_scaled"] = Dp17O(df["d17O_scaled"], df["d18O_scaled"])
+
 
-# Calculate the mean values from the replicate measurements
-df_scaled = df.loc[:, ["SampleName", "d18O_scaled", "Dp17O_scaled"]]
-df_mean = df_scaled.groupby('SampleName').mean().reset_index()
-df_mean = df_mean.rename(columns={'d18O_scaled': 'd18O_CO2', 'Dp17O_scaled': 'Dp17O_CO2'})
+# Calculate the mean and SD of the replicate measurements
+df_gases = df.loc[:, ["SampleName", "d18O_scaled", "Dp17O_scaled"]]
+df_gases = df_gases[df_gases["SampleName"].str.contains("light|heavy")]
+df_gases_A = df_gases.groupby('SampleName').mean().reset_index()
+df_gases_A = df_gases_A.rename(columns={'d18O_scaled': 'd18O_CO2', 'Dp17O_scaled': 'Dp17O_CO2'})
+df_gases_SD = df_gases.groupby('SampleName').std().reset_index()
+df_gases_SD = df_gases_SD.rename(columns={'d18O_scaled': 'd18O_SD', 'Dp17O_scaled': 'Dp17O_SD'})
 
-# Calculate the standard deviation from the replicate measurements
-df_std = df_scaled.groupby('SampleName').std().reset_index()
-df_std = df_std.rename(columns={'d18O_scaled': 'd18O_error',  'Dp17O_scaled': 'Dp17O_error'})
+print("\nScaled values:")
+print(df_gases_A.merge(df_gases_SD, on='SampleName').round({"Dp17O_CO2": 0, "d18O_CO2": 3, "Dp17O_SD": 0, "d18O_SD": 3}))
 
-dfMerged = df_mean.merge(df_std, on='SampleName')
-dfMerged['Replicates'] = df_scaled.groupby('SampleName').size().reset_index(name='counts')['counts']
 
-print("\nScaled values:\n", dfMerged.round({"Dp17O_CO2": 0, "Dp17O_error": 0}).round(3))
+# Create Figure S1
+fig, (ax1, ax2) = plt.subplots(1, 2)
 
 # Create new dataframes for IAEA-603 and NBS-18
 df_IAEA603 = df[df["SampleName"].str.contains("IAEA603")]
 df_NBS18 = df[df["SampleName"].str.contains("NBS18")]
 
-
-# Create Figure S1
-fig, (ax1, ax2) = plt.subplots(1, 2)
-
 # Plot IAEA-603 data
 ax1.scatter(df_IAEA603["d18O"], df_IAEA603["Dp17O"],
             label="IAEA-603", c="#1455C0", marker="s")
@@ -129,7 +127,7 @@ def to_VSMOW(VPDB):
 ax1.text(0.98, 0.98, "A", size=14, ha="right", va="top",
          transform=ax1.transAxes, fontweight="bold")
 ax1.set_title("IAEA-603")
-ax1.set_xlabel("$\delta^{18}$O (‰, VSMOW)")
+ax1.set_xlabel("$\delta^{18}$O (‰, unscaled)")
 ax1.set_ylabel("$\Delta^{\prime 17}$O (ppm, unscaled)")
 
 # Plot NBS-18 data
@@ -142,7 +140,7 @@ def to_VSMOW(VPDB):
 ax2.text(0.98, 0.98, "B", size=14, ha="right", va="top",
          transform=ax2.transAxes, fontweight="bold")
 ax2.set_title("NBS-18")
-ax2.set_xlabel("$\delta^{18}$O (‰, VSMOW)")
+ax2.set_xlabel("$\delta^{18}$O (‰, unscaled)")
 ax2.set_ylabel("$\Delta^{\prime 17}$O (ppm, unscaled)")
 
 plt.savefig(sys.path[0] + "/" + "OH2 Figure S1.png")