solve_carbonate.py · Scientific Agent Skills

1#!/usr/bin/env python3

2"""Solve paired seawater carbonate measurements with a recorded PyCO2SYS setup."""

3

4from __future__ import annotations

5

6import argparse

7import csv

8import hashlib

9import io

10import json

11import math

12from pathlib import Path

13import platform

14

15TYPES = {"alkalinity": 1, "dic": 2, "ph": 3, "pco2": 4, "fco2": 5}

16SCALES = {"total": 1, "seawater": 2, "free": 3, "nbs": 4}

17REQUIRED = ("par1", "par2", "salinity", "temperature", "pressure",

18 "total_phosphate", "total_silicate")

19OPTIONAL = ("temperature_out", "pressure_out", "total_ammonia", "total_sulfide")

20RESULTS = ("alkalinity", "dic", "pH_total", "pCO2", "fCO2", "aqueous_CO2",

21 "bicarbonate", "carbonate", "saturation_aragonite", "saturation_calcite",

22 "revelle_factor")

23UNCERTAINTY_RESULTS = ("pH_total", "pCO2", "saturation_aragonite")

24# Published calibration ranges from the PyCO2SYS v1 argument documentation.

25RANGES = {10: (2, 35, 19, 43), 15: (0, 45, 0, 45)}

26SETTINGS = {"opt_k_bisulfate": 1, "opt_k_fluoride": 1, "opt_gas_constant": 3,

27 "opt_buffers_mode": 1, "opt_pressured_kCO2": 0,

28 "pressure_atmosphere": 1.0, "pressure_atmosphere_out": 1.0}

29

30

31def read_samples(path: Path):

32 """Reject ambiguous schemas and invalid numbers before any computation."""

33 raw = path.read_bytes()

34 reader = csv.DictReader(io.StringIO(raw.decode("utf-8-sig")))

35 fields = reader.fieldnames or []

36 if len(fields) != len(set(fields)):

37 raise ValueError("Duplicate CSV column names")

38 missing = {"sample_id", *REQUIRED} - set(fields)

39 if missing:

40 raise ValueError(f"Missing required columns: {sorted(missing)}")

41 allowed = {"sample_id", *REQUIRED, *OPTIONAL}

42 allowed.update("u_" + key for key in (*REQUIRED, *OPTIONAL))

43 unknown = set(fields) - allowed

44 if unknown:

45 raise ValueError(f"Unknown columns: {sorted(unknown)}; use the documented schema")

46 if ("temperature_out" in fields) != ("pressure_out" in fields):

47 raise ValueError("Provide both temperature_out and pressure_out, or neither")

48 for key in fields:

49 if key.startswith("u_") and key[2:] not in fields:

50 raise ValueError(f"{key} requires the corresponding {key[2:]} column")

51 rows, seen = [], set()

52 for line, row in enumerate(reader, 2):

53 if None in row or any(value is None for value in row.values()):

54 raise ValueError(f"Row {line}: wrong number of fields")

55 identifier = row["sample_id"].strip()

56 if not identifier or identifier in seen:

57 raise ValueError(f"Row {line}: empty or duplicate sample_id {identifier!r}")

58 seen.add(identifier)

59 parsed = {"sample_id": identifier}

60 for key, value in row.items():

61 if key == "sample_id":

62 continue

63 try:

64 number = float(value)

65 except ValueError as exc:

66 raise ValueError(f"Row {line}: {key} must be numeric, got {value!r}") from exc

67 if not math.isfinite(number):

68 raise ValueError(f"Row {line}: {key} must be finite")

69 nonnegative = (key.startswith(("u_", "total_", "pressure")) or key == "salinity")

70 if nonnegative and number < 0:

71 raise ValueError(f"Row {line}: {key} must be nonnegative")

72 if key.startswith("temperature") and number <= -273.15:

73 raise ValueError(f"Row {line}: {key} must exceed absolute zero")

74 parsed[key] = number

75 rows.append(parsed)

76 if not rows:

77 raise ValueError("CSV contains no samples")

78 return rows, fields, hashlib.sha256(raw).hexdigest()

79

80

81def solve(path: Path, output_dir: Path, *, par1_type: str, par2_type: str,

82 k_carbonic: int, ph_scale: str | None = None, total_borate: int = 1):

83 """Write results only after all rows produce finite, carbon-balanced solutions."""

84 if output_dir.exists():

85 raise FileExistsError(f"Output directory already exists: {output_dir}")

86 if par1_type not in TYPES or par2_type not in TYPES:

87 raise ValueError("Unsupported carbonate parameter type")

88 if par1_type == par2_type or {par1_type, par2_type} == {"pco2", "fco2"}:

89 raise ValueError("Provide two independent carbonate parameters")

90 if "ph" in (par1_type, par2_type) and ph_scale is None:

91 raise ValueError("pH input requires an explicit --ph-scale")

92 if ph_scale is not None and ph_scale not in SCALES:

93 raise ValueError("Unsupported pH scale")

94 if k_carbonic not in RANGES or total_borate not in (1, 2):

95 raise ValueError("Unsupported equilibrium-constant or total-borate option")

96 rows, fields, digest = read_samples(path)

97 for row in rows:

98 for key, kind in (("par1", par1_type), ("par2", par2_type)):

99 value = row[key]

100 if kind in ("dic", "pco2", "fco2") and value <= 0:

101 raise ValueError(f"{row['sample_id']}: {kind} must be positive")

102 if kind == "ph" and not 0 < value < 14:

103 raise ValueError(f"{row['sample_id']}: helper supports only 0 < pH < 14")

104

105 import numpy as np

106 import PyCO2SYS as pyco2

107

108 if pyco2.__version__ != "1.8.3.4":

109 raise ValueError("This helper targets PyCO2SYS==1.8.3.4; use the pinned environment")

110 settings = {**SETTINGS, "opt_k_carbonic": k_carbonic,

111 "opt_total_borate": total_borate, "opt_pH_scale": SCALES[ph_scale or "total"],

112 "par1_type": TYPES[par1_type], "par2_type": TYPES[par2_type]}

113 numeric_fields = [key for key in fields if key != "sample_id" and not key.startswith("u_")]

114 kwargs = {key: np.array([row[key] for row in rows]) for key in numeric_fields}

115 # Explicitly record the two optional zero assumptions, not library defaults.

116 for key in ("total_ammonia", "total_sulfide"):

117 kwargs.setdefault(key, np.zeros(len(rows)))

118 uncertainties = {key[2:]: np.array([row[key] for row in rows])

119 for key in fields if key.startswith("u_")}

120 has_output = "temperature_out" in fields

121 result_keys = list(RESULTS)

122 uncertainty_keys = list(UNCERTAINTY_RESULTS)

123 if has_output:

124 result_keys.extend(key + "_out" for key in RESULTS if key not in ("alkalinity", "dic"))

125 uncertainty_keys.extend(key + "_out" for key in UNCERTAINTY_RESULTS)

126 uncertainty_options = ({"uncertainty_from": uncertainties, "uncertainty_into": uncertainty_keys}

127 if uncertainties else {})

128 calculated = pyco2.sys(**kwargs, **settings, **uncertainty_options)

129 if uncertainties:

130 result_keys.extend("u_" + key for key in uncertainty_keys)

131 arrays = {key: np.broadcast_to(np.asarray(calculated[key], dtype=float), (len(rows),))

132 for key in result_keys}

133 for key, values in arrays.items():

134 if not np.all(np.isfinite(values)):

135 bad = [rows[i]["sample_id"] for i in np.flatnonzero(~np.isfinite(values))]

136 raise ValueError(f"Non-finite result {key} for {bad}; inspect the input pair and conditions")

137 for suffix in ("", "_out") if has_output else ("",):

138 components = sum(arrays[key + suffix] for key in ("aqueous_CO2", "bicarbonate", "carbonate"))

139 if not np.allclose(components, arrays["dic"], rtol=1e-8, atol=1e-6):

140 raise ValueError("Dissolved inorganic carbon mass balance failed")

141 if np.any(arrays["dic"] <= 0) or any(

142 np.any(arrays[key + suffix] < 0)

143 for key in ("aqueous_CO2", "bicarbonate", "carbonate")

144 ):

145 raise ValueError("Solution has nonphysical carbon concentrations")

146

147 warnings = []

148 tmin, tmax, smin, smax = RANGES[k_carbonic]

149 result_rows = []

150 for index, row in enumerate(rows):

151 flags = []

152 for suffix in ("", "_out") if has_output else ("",):

153 label = "output" if suffix else "input"

154 if not (tmin < row["temperature" + suffix] < tmax and smin < row["salinity"] < smax):

155 flags.append(f"outside_k_carbonic_{label}_range")

156 if row["pressure" + suffix] > 0:

157 flags.append(f"gas_pressure_correction_disabled_{label}")

158 if flags:

159 warnings.append({"sample_id": row["sample_id"], "flags": flags})

160 result_rows.append({"sample_id": row["sample_id"],

161 **{"input_" + key: row[key] for key in fields if key != "sample_id"},

162 **{key: float(values[index]) for key, values in arrays.items()},

163 "qc_flags": ";".join(flags)})

164 assumptions = [f"{key}=0 micromol/kg seawater" for key in ("total_ammonia", "total_sulfide")

165 if key not in fields]

166 manifest = {

167 "schema_version": "1.0", "skill_version": "1.1",

168 "software": {"PyCO2SYS": pyco2.__version__, "numpy": np.__version__,

169 "python": platform.python_version()},

170 "input_file": path.name, "input_sha256": digest, "sample_count": len(rows),

171 "parameter_types": {"par1": par1_type, "par2": par2_type}, "settings": settings,

172 "output_conditions": "provided per row" if has_output else "input conditions only",

173 "gas_pressure_convention": (

174 "CO2 solubility and fugacity use 1 atm atmospheric pressure without hydrostatic "

175 "gas corrections (opt_pressured_kCO2=0). Other pressure-dependent equilibria "

176 "still use the supplied sea pressure; nonzero-pressure pCO2/fCO2 must not be "

177 "interpreted as fully pressure-corrected in situ gas values."),

178 "units": {"carbon_and_nutrients": "micromol/kg seawater", "pCO2_and_fCO2": "microatm",

179 "temperature": "degrees Celsius", "pressure": "dbar, sea pressure",

180 "salinity": "Practical Salinity, dimensionless", "pH_total": "total scale",

181 "saturation_and_revelle": "dimensionless",

182 "uncertainties": "absolute 1-sigma, in the corresponding parameter units"},

183 "uncertainty_sources": sorted(uncertainties),

184 "uncertainty_scope": ("Independent input uncertainties only; unspecified inputs and "

185 "equilibrium constants treated as exact; no covariance or model error."),

186 "assumptions": assumptions, "warnings": warnings,

187 }

188 output_dir.mkdir(parents=True, exist_ok=False)

189 with (output_dir / "carbonate.csv").open("w", newline="", encoding="utf-8") as handle:

190 writer = csv.DictWriter(handle, fieldnames=list(result_rows[0]))

191 writer.writeheader()

192 writer.writerows(result_rows)

193 (output_dir / "provenance.json").write_text(json.dumps(manifest, indent=2, allow_nan=False) + "\n")

194 return manifest

195

196

197def main():

198 parser = argparse.ArgumentParser(description=__doc__)

199 parser.add_argument("input_csv", type=Path)

200 parser.add_argument("--output-dir", type=Path, required=True, help="New directory for CSV and provenance")

201 parser.add_argument("--par1-type", choices=TYPES, required=True)

202 parser.add_argument("--par2-type", choices=TYPES, required=True)

203 parser.add_argument("--ph-scale", choices=SCALES, help="Required when either input is pH")

204 parser.add_argument("--k-carbonic", type=int, choices=RANGES, required=True,

205 help="10: Lueker et al. 2000; 15: Waters et al. 2014")

206 parser.add_argument("--total-borate", type=int, choices=(1, 2), default=1,

207 help="1: Uppstrom 1974 (default); 2: Lee et al. 2010")

208 args = parser.parse_args()

209 try:

210 manifest = solve(args.input_csv, args.output_dir, par1_type=args.par1_type,

211 par2_type=args.par2_type, ph_scale=args.ph_scale,

212 k_carbonic=args.k_carbonic, total_borate=args.total_borate)

213 except (ValueError, OSError, UnicodeError) as exc:

214 parser.exit(2, f"error: {exc}\n")

215 print(f"Solved {manifest['sample_count']} samples; {len(manifest['warnings'])} with QC flags.")

216 for entry in manifest["warnings"]:

217 print(f"WARNING {entry['sample_id']}: {', '.join(entry['flags'])}")

218 print(f"Results: {args.output_dir / 'carbonate.csv'}")

219

220

221if __name__ == "__main__":

222 main()