Detectability limits for OAE experiment

PROJECT_ROOT = "/anvil/projects/x-ees250129/x-uheede/C-Star-in-Hvalfjordur"
obs_ctd_mooring_path_july = f"{PROJECT_ROOT}/data/staged/seanoe_hvalfjordur/2026-04-27/seanoe_113246/ctd_mooring_qc/Mooring_data_final_.csv/HVIN1_CTDphox_60m_20240404_20240814.csv"
obs_ctd_mooring_path_nov = f"{PROJECT_ROOT}/data/staged/seanoe_hvalfjordur/2026-04-27/seanoe_113246/ctd_mooring_qc/Mooring_data_final_.csv/HVIN2_CTDphox_57m_20240814_20250217.csv"
obs_bgc_path = f"{PROJECT_ROOT}/data/staged/seanoe_hvalfjordur/2026-04-27/seanoe_110401/discrete_samples_qc/discrete_samples_qc.xlsx"

import pandas as pd


# If wildcard still present, resolve it once
import glob
file1 = obs_ctd_mooring_path_july
file2 = obs_ctd_mooring_path_nov

# Read data
df = pd.read_csv(file)

# -----------------------------
# Convert time to datetime
# -----------------------------
df["time"] = pd.to_datetime("1950-01-01") + pd.to_timedelta(
    df["Time[days_since_1950-01-01T00:00:00Z]"], unit="D"
)

df = df.set_index("time")

# -----------------------------
# Select July 2024
# -----------------------------
df_july = df.loc["2024-07-01":"2024-07-31"]
df_nov = df.loc["2024-07-01":"2024-07-31"]
# -----------------------------
# Standard deviation of pH
# -----------------------------
ph_std_july = df_july["pH_qc[mL/L]"].std()
ph_std_nov = df_nov["pH_qc[mL/L]"].std()


model_mean_ph_july=8.088
model_mean_ph_nov=8.123

print("July detectability limit:", ph_std_july + model_mean_ph_july)
print("July detectability limit double:", ph_std_july*2 + model_mean_ph_july)
print("November detectability limit:", ph_std_nov + model_mean_ph_nov)
print("November detectability limit double:", ph_std_nov*2 + model_mean_ph_nov)

ph_std_july

import numpy as np
import xarray as xr
import gsw as sw

xls = pd.ExcelFile(obs_bgc_path)
combo = pd.read_excel(xls, 'Sheet1', decimal='.', header=27, skiprows=[28])

# We use .str.replace with regex=False for a simple string swap
combo['HV'] = combo['Station_name'].str.replace('HV', '', regex=False).astype(int)

# 2. Create the 'time' column
combo['time'] = pd.to_datetime(combo[['Year_UTC', 'Month_UTC', 'Day_UTC']].rename(
    columns={'Year_UTC': 'year', 'Month_UTC': 'month', 'Day_UTC': 'day'}
))

# -----------------------------
# CREATE DEPTH BINS (0–9.9, 10–19.9, ...)
# -----------------------------
combo["depth_bin"] = (np.floor(combo["Depth"] / 10) * 10).astype(int)

obs=xr.Dataset.from_dataframe(combo)
obs=obs.set_index(index=['HV','depth_bin','time'])
obs=obs.drop_duplicates('index')
obs=obs.unstack('index')
obs=obs.rename(name_dict={'Latitude':'lat','Longitude':'lon'})



density=sw.density.sigma0(obs['Salinity_CTD'],obs['Temperature_CTD'])

obs['Alk']=obs['TA'] * (1 / 1000) * (density+1000)
obs['DIC']=obs['DIC'] * (1 / 1000) * (density+1000)

# -----------------------------
# QC MASK (all variables must pass)
# -----------------------------
qc_mask = (
    (obs["TEMP_flag"] == 2) &
    (obs["Salinity_CTD_flag"] == 2) &
    (obs["TA_flag"] == 2) &
    (obs["DIC_flag"] == 2)
)

# -----------------------------
# APPLY MASK INTO obs (IMPORTANT)
# -----------------------------
obs["Alk_qc"]      = obs["Alk"].where(qc_mask)
obs["DIC_qc"]      = obs["DIC"].where(qc_mask)
obs["Salinity_qc"] = obs["Salinity_CTD"].where(qc_mask)

# -----------------------------
# Constants
# -----------------------------
TA0 = 282
Sref = 34
DIC0 = 313

# -----------------------------
# NORMALIZATION (QC-CONSISTENT)
# -----------------------------
obs["nTA"] = ((obs["Alk_qc"] - TA0) / obs["Salinity_qc"]) * Sref + TA0
obs["nDIC"] = ((obs["DIC_qc"] - DIC0) / obs["Salinity_qc"]) * Sref + DIC0

obs_hv10=obs.sel(HV=10)
TA_mmol_m3=obs_hv10['nTA']
DIC_mmol_m3=obs_hv10['nDIC']
# -----------------------------
# Select months
# -----------------------------
july_mask = obs_hv10['time'].dt.month == 7
nov_mask  = obs_hv10['time'].dt.month == 11

# -----------------------------
# JULY
# -----------------------------
TA_july_std  = TA_mmol_m3.where(july_mask).std(dim=['time', 'depth_bin'], skipna=True)
DIC_july_std = DIC_mmol_m3.where(july_mask).std(dim=['time', 'depth_bin'], skipna=True)

# -----------------------------
# NOVEMBER
# -----------------------------
TA_nov_std  = TA_mmol_m3.where(nov_mask).std(dim=['time', 'depth_bin'], skipna=True)
DIC_nov_std = DIC_mmol_m3.where(nov_mask).std(dim=['time', 'depth_bin'], skipna=True)

# -----------------------------
# Print results
# -----------------------------
print("July TA std (mmol/m³):", TA_july_std.values)
print("July DIC std (mmol/m³):", DIC_july_std.values)

print("November TA std (mmol/m³):", TA_nov_std.values)
print("November DIC std (mmol/m³):", DIC_nov_std.values)

import pandas as pd

# File path
file_path = f"{PROJECT_ROOT}/data/2024-11-05_2024-12-10+pCO2_moored.txt"

# Column names
cols = [
    "Mode", "Year", "Month", "Day", "Hour", "Minute", "Second",
    "SampleNum", "RecordNum", "CO2_ppm", "TempWater_C",
    "DetV_V", "LampV_V", "Status", "TempInternal_C",
    "TempBench_C", "SupplyV_V"
]

# Read file (adjust delimiter if needed, often whitespace or comma)
df = pd.read_csv(file_path, sep=",", names=cols)

# Create datetime column
df["time"] = pd.to_datetime(df[["Year", "Month", "Day", "Hour", "Minute", "Second"]])

# Set as index (optional but cleaner)
df = df.set_index("time")

# Filter for November
df_nov = df[df.index.month == 11]

# Compute standard deviation
pco2_std_nov = df_nov["CO2_ppm"].std()

print(f"November pCO2 standard deviation: {pco2_std_nov:.3f} ppm")