"""Quality checks to run *before* any mixing analysis.
A mixing calculation cannot tell you that an analysis is wrong -- it will happily
apportion a charge-imbalanced sample among your end-members and return confident
numbers. These checks are the cheapest protection available, and they belong at the
top of every notebook.
"""
from __future__ import annotations
from typing import TYPE_CHECKING, Any
import numpy as np
from mescla._util import as_chem
from mescla.prep.units import lookup_ion, to_meq_per_l
if TYPE_CHECKING: # pragma: no cover
import pandas as pd
__all__ = [
"charge_balance_error",
"charge_balance_report",
"detection_limit_report",
"ec_check",
"ec_from_ions",
"missing_data_report",
]
#: Empirical factor relating the cation sum (meq/L) to electrical conductivity
#: (uS/cm) in ordinary fresh waters. Good to roughly +/-10%; degrades badly in
#: brines and in very dilute waters.
EC_PER_MEQ = 100.0
def _meq_split(table: Any) -> tuple[np.ndarray, np.ndarray, list[str], list[str]]:
"""Return cation and anion sums in meq/L, plus the species used in each."""
# Non-ionic species (EC, SiO2, isotopes) carry no charge and are excluded from
# the sums by design, so skipping them is the intended policy rather than an
# oversight -- hence on_unknown="skip" instead of the stricter default.
chem = to_meq_per_l(as_chem(table), on_unknown="skip")
cations = np.zeros(chem.n_rows)
anions = np.zeros(chem.n_rows)
cat_names: list[str] = []
an_names: list[str] = []
for j, sp in enumerate(chem.species):
ion = lookup_ion(sp)
if ion is None or ion.charge == 0:
continue
col = np.nan_to_num(chem.data[:, j], nan=0.0)
if ion.charge > 0:
cations += col
cat_names.append(sp)
else:
anions += col
an_names.append(sp)
return cations, anions, cat_names, an_names
[docs]
def charge_balance_error(table: Any, convention: str = "sum") -> np.ndarray:
"""Charge balance error in percent, per sample.
Two conventions are in circulation and they differ by exactly a factor of
two, so a CBE is meaningless unless the convention travels with it.
``convention="sum"`` (the default) normalises by the *sum* of the two sums,
.. math::
\\mathrm{CBE} = 100 \\,
\\frac{\\sum \\mathrm{cations} - \\sum \\mathrm{anions}}
{\\sum \\mathrm{cations} + \\sum \\mathrm{anions}}
which is Freeze & Cherry (1979) eq. 3.21 and what PHREEQC prints as
``Percent error, 100*(Cat-|An|)/(Cat+|An|)``.
``convention="mean"`` normalises by their *mean*, i.e. the same expression
with 200 in place of 100. This is the analytical-laboratory "ion balance
percent difference" of APHA Standard Methods 1030E, and the form used by
much of the water-quality QA literature.
Neither is wrong; ``"mean"`` simply reports twice the number for the same
water. The ±5% acceptance rule is quoted for *both*, which is the real trap:
screening at 5% under ``"mean"`` is twice as strict as under ``"sum"``.
Relax to 10% for very dilute waters, where a small absolute error is a large
relative one.
A large positive CBE usually means a missing anion (organic acids, nitrate,
unmeasured alkalinity); a large negative one, a missing cation.
Parameters
----------
convention : {"sum", "mean"}, default "sum"
Denominator convention, as above.
"""
if convention not in ("sum", "mean"):
raise ValueError(f"convention must be 'sum' or 'mean'; got {convention!r}")
cations, anions, cat, an = _meq_split(table)
if not cat or not an:
raise ValueError(
f"charge balance needs both cations and anions; found cations={cat}, anions={an}"
)
factor = 100.0 if convention == "sum" else 200.0
total = cations + anions
with np.errstate(divide="ignore", invalid="ignore"):
return np.where(total > 0, factor * (cations - anions) / total, np.nan)
[docs]
def charge_balance_report(
table: Any, tolerance: float = 5.0, convention: str = "sum"
) -> pd.DataFrame:
"""Per-sample charge balance with a pass/fail flag.
Parameters
----------
tolerance : float, default 5.0
Acceptance threshold on ``|CBE|`` in percent. It is only meaningful together
with ``convention`` -- see :func:`charge_balance_error`.
convention : {"sum", "mean"}, default "sum"
Denominator convention for the CBE.
Returns
-------
DataFrame
Columns ``cations_meq``, ``anions_meq``, ``cbe_percent``, ``pass``. The
species that contributed to each sum are recorded in ``.attrs["species_used"]``,
and the convention used in ``.attrs["convention"]``.
"""
import pandas as pd
chem = as_chem(table)
cations, anions, cat, an = _meq_split(chem)
cbe = charge_balance_error(chem, convention=convention)
report = pd.DataFrame(
{
"cations_meq": cations,
"anions_meq": anions,
"cbe_percent": cbe,
"pass": np.abs(cbe) <= tolerance,
},
index=list(chem.index),
)
report.attrs["species_used"] = {"cations": cat, "anions": an}
report.attrs["tolerance"] = tolerance
report.attrs["convention"] = convention
return report
[docs]
def ec_from_ions(table: Any) -> np.ndarray:
"""Estimate electrical conductivity (uS/cm) from the cation sum.
Approximate by construction -- see :data:`EC_PER_MEQ`. Its value is as a
*consistency* check against a measured EC, not as a substitute for one.
"""
cations, _, _, _ = _meq_split(table)
return cations * EC_PER_MEQ
[docs]
def ec_check(table: Any, ec_species: str = "EC", tolerance: float = 15.0) -> pd.DataFrame:
"""Compare measured EC with EC estimated from the ion sum.
A mismatch beyond ``tolerance`` percent points to a major ion that was not
analysed, or to a unit error.
"""
import pandas as pd
chem = as_chem(table)
if ec_species not in chem.species:
raise KeyError(f"{ec_species!r} not among species: {chem.species}")
measured = chem.data[:, chem.species.index(ec_species)]
estimated = ec_from_ions(chem.drop_species(ec_species))
with np.errstate(divide="ignore", invalid="ignore"):
diff = np.where(measured > 0, 100.0 * (estimated - measured) / measured, np.nan)
return pd.DataFrame(
{
"ec_measured": measured,
"ec_from_ions": estimated,
"difference_percent": diff,
"pass": np.abs(diff) <= tolerance,
},
index=list(chem.index),
)
[docs]
def detection_limit_report(table: Any, limits: dict[str, float]) -> pd.DataFrame:
"""Count values at or below the detection limit, per species.
Species with many non-detects are poor tracers: near the detection limit the
error distribution is neither Gaussian nor homoscedastic, which breaks the
weighting assumed by every least-squares and maximum-likelihood estimator here.
Carry such species with a large sigma instead of dropping them outright.
"""
import pandas as pd
chem = as_chem(table)
rows = []
for j, sp in enumerate(chem.species):
lod = limits.get(sp)
col = chem.data[:, j]
n_below = int(np.sum(col <= lod)) if lod is not None else 0
rows.append(
{
"species": sp,
"detection_limit": lod,
"n_at_or_below": n_below,
"fraction": n_below / chem.n_rows if chem.n_rows else np.nan,
"min": np.nanmin(col) if col.size else np.nan,
}
)
return pd.DataFrame(rows).set_index("species")
[docs]
def missing_data_report(table: Any) -> pd.DataFrame:
"""Count missing values per species, with the count of complete samples."""
import pandas as pd
chem = as_chem(table)
nan = np.isnan(chem.data)
out = pd.DataFrame(
{
"n_missing": nan.sum(axis=0),
"fraction_missing": nan.mean(axis=0),
},
index=list(chem.species),
)
out.attrs["n_complete_samples"] = int(np.sum(~nan.any(axis=1)))
return out