"""How densely has deep-time marine sediment actually been sampled geochemically? Deep-Time Marine Sedimentary Element Database v0.0.1 (Lai et al. 2025, Zenodo 14771859), which extends the Sedimentary Geochemistry and Paleoenvironments Project compilation (Farrell et al. 2021): 63,627 samples with ages, sites and element/isotope values. Measures, per geological interval: * samples per Myr and distinct sites; * expected number of existing samples inside a random time window of length D; * fraction of time covered by >= 1 (and >= 2) sites whose mean sample spacing is <= X, for all samples, for carbon (TOC or d13C) and for lead (exact interval union). These are statements about this one compilation and are lower bounds on the true state of sampling; s04 supplies the CENOGRID spacing coverage and s07 combines the two floors. Outputs: dmsed_bins.csv, dmsed_coverage.csv, dmsed_summary.json """ import os, json import numpy as np, pandas as pd from common import RAW, DATA, save_json BINS = [(0, 2.58, "Quaternary"), (2.58, 23.03, "Neogene"), (23.03, 66, "Paleogene"), (66, 145, "Cretaceous"), (145, 201.4, "Jurassic"), (201.4, 251.9, "Triassic"), (251.9, 298.9, "Permian"), (298.9, 358.9, "Carboniferous"), (358.9, 419.2, "Devonian"), (419.2, 443.8, "Silurian"), (443.8, 485.4, "Ordovician"), (485.4, 538.8, "Cambrian"), (538.8, 1000, "Neoproterozoic"), (1000, 1600, "Mesoproterozoic"), (1600, 2500, "Paleoproterozoic"), (2500, 3500, "Archean (to 3.5 Ga)")] SPACINGS_KYR = [0.1, 1, 3.33, 10, 50, 200, 333] H_SAMPLE_MYR = 1e-3 # time represented by one geochemical sample (1-2 cm at ~1-2 cm/kyr) def load(): d = pd.read_csv(os.path.join(RAW, "dmsed", "DM-SED_0.0.1.csv"), low_memory=False, encoding_errors="replace") cols = {c: c for c in d.columns} def find(sub): return [c for c in d.columns if sub in c][0] d["age"] = pd.to_numeric(d["Age (Ma)"], errors="coerce") d["has_c"] = (pd.to_numeric(d[find("TOC (wt%)")], errors="coerce").notna() | pd.to_numeric(d[find("13CTOC")], errors="coerce").notna() | pd.to_numeric(d[find("13Ccarb")], errors="coerce").notna()) d["has_d13c"] = (pd.to_numeric(d[find("13CTOC")], errors="coerce").notna() | pd.to_numeric(d[find("13Ccarb")], errors="coerce").notna()) d["has_pb"] = pd.to_numeric(d[find("Pb(ppm)")], errors="coerce").notna() d["site"] = d["SiteName"].astype(str) + "|" + d["ModLat"].round(2).astype(str) + "|" + d["ModLon"].round(2).astype(str) return d.dropna(subset=["age"]) def site_table(d, flag=None): s = d if flag is None else d[d[flag]] g = s.groupby("site")["age"].agg(["count", "min", "max"]) g["span"] = g["max"] - g["min"] g = g[(g["count"] >= 5) & (g["span"] > 0)].copy() g["spacing_kyr"] = g["span"] / (g["count"] - 1) * 1000.0 return g def coverage(g, max_spacing_kyr, lo, hi, nmin=1): """Fraction of [lo, hi] Ma covered by >= nmin sites with mean spacing <= max_spacing_kyr. Exact interval arithmetic (sweep over the clipped site age ranges). The first-pass version evaluated coverage on a 4,000-point grid, which overstated very small coverages (0.075% instead of 0.021% for <= 1 kyr in 2.58-500 Ma; referee point 1).""" s = g[g["spacing_kyr"] <= max_spacing_kyr] if len(s) == 0: return 0.0 a = np.clip(s["min"].values.astype(float), lo, hi); b = np.clip(s["max"].values.astype(float), lo, hi) ok = b > a if not ok.any(): return 0.0 ev = sorted([(x, 1) for x in a[ok]] + [(x, -1) for x in b[ok]]) depth = 0; covered = 0.0; prev = lo for x, dlt in ev: if depth >= nmin: covered += x - prev depth += dlt; prev = x return float(covered / (hi - lo)) def main(): d = load() out = {"n_samples": int(len(d)), "n_sites": int(d["site"].nunique()), "n_refs": int(d["Reference"].nunique()), "n_with_carbon": int(d["has_c"].sum()), "n_with_d13c": int(d["has_d13c"].sum()), "n_with_pb": int(d["has_pb"].sum())} G = {"all": site_table(d), "carbon": site_table(d, "has_c"), "d13c": site_table(d, "has_d13c"), "pb": site_table(d, "has_pb")} out["n_sites_with_age_span_ge5samples"] = {k: int(len(v)) for k, v in G.items()} out["median_site_spacing_kyr"] = {k: float(v["spacing_kyr"].median()) for k, v in G.items()} st_all = d.groupby("site")["age"].agg(["count", "min", "max"]) out["frac_sites_ge5_samples_with_single_age"] = float(((st_all["count"] >= 5) & (st_all["max"] == st_all["min"])).sum() / (st_all["count"] >= 5).sum()) single = (st_all["count"] >= 5) & (st_all["max"] == st_all["min"]) out["n_sites_ge5_samples"] = int((st_all["count"] >= 5).sum()) out["n_sites_ge5_samples_with_single_age"] = int(single.sum()) out["n_samples_in_single_age_sites"] = int(st_all.loc[single, "count"].sum()) out["scope_note"] = ("DM-SED v0.0.1 is one compilation of (mostly shale and carbonate) element geochemistry. It is NOT a census " "of geochemical sampling: it holds no ocean-drilling stable-isotope or core-scanning data, and sites whose " "samples all carry one age (dense sections without an age model) cannot be used for spacing. Every coverage " "number derived from it is a lower bound on what exists.") for flag, nm in (("has_d13c", "d13c"), ("has_c", "carbon")): for lo_, hi_, lab_ in ((2.58, 66.0, "2.58-66Ma"), (66.0, 252.0, "66-252Ma"), (252.0, 539.0, "252-539Ma")): out[f"n_{nm}_samples_{lab_}"] = int((d[flag] & (d["age"] >= lo_) & (d["age"] < hi_)).sum()) d13carb = pd.to_numeric(d[[c for c in d.columns if "13Ccarb" in c][0]], errors="coerce").notna() out["n_d13Ccarb_samples_2.58-66Ma"] = int((d13carb & (d["age"] >= 2.58) & (d["age"] < 66)).sum()) rows = []; crow = [] for lo, hi, lab in BINS + [(2.58, 500, "2.58-500 Ma"), (2.58, 2400, "2.58-2400 Ma"), (66, 252, "Mesozoic"), (252, 539, "Paleozoic")]: s = d[(d["age"] >= lo) & (d["age"] < hi)] r = {"bin": lab, "age_top_ma": lo, "age_base_ma": hi, "n_samples": int(len(s)), "n_sites": int(s["site"].nunique()), "samples_per_myr": len(s) / (hi - lo), "pb_samples_per_myr": s["has_pb"].sum() / (hi - lo), "carbon_samples_per_myr": s["has_c"].sum() / (hi - lo)} for D, nm in ((3e-4, "300yr"), (1e-2, "10kyr"), (1.0, "1Myr")): r[f"expected_samples_in_window_{nm}"] = r["samples_per_myr"] * (D + H_SAMPLE_MYR) r[f"expected_pb_samples_in_window_{nm}"] = r["pb_samples_per_myr"] * (D + H_SAMPLE_MYR) rows.append(r) for key, g in G.items(): for X in SPACINGS_KYR: crow.append({"bin": lab, "age_top_ma": lo, "age_base_ma": hi, "set": key, "max_spacing_kyr": X, "frac_time_ge1_site": coverage(g, X, lo, hi, 1), "frac_time_ge2_sites": coverage(g, X, lo, hi, 2)}) pd.DataFrame(rows).to_csv(os.path.join(DATA, "dmsed_bins.csv"), index=False, float_format="%.4g") cov = pd.DataFrame(crow) cov.to_csv(os.path.join(DATA, "dmsed_coverage.csv"), index=False, float_format="%.4g") out["bins"] = rows save_json(out, "dmsed_summary.json") pd.set_option("display.width", 250) print(json.dumps({k: v for k, v in out.items() if k != "bins"}, indent=1)) print(pd.DataFrame(rows)[["bin", "n_samples", "n_sites", "samples_per_myr", "pb_samples_per_myr", "expected_samples_in_window_300yr", "expected_samples_in_window_10kyr"]].to_string()) pv = cov[cov["set"].isin(["all", "carbon", "pb"])].pivot_table(index=["bin"], columns=["set", "max_spacing_kyr"], values="frac_time_ge1_site", sort=False) print(pv.round(3).to_string()) if __name__ == "__main__": main()