"""Rock-record preservation statistics from Macrostrat (Peters, Husson & Czaplewski 2018). Outputs (data/): macrostrat_coverage_1myr.csv per-Myr coverage of sedimentary rock, by region and environment macrostrat_bins.csv the same in geological bins plus resolution (yr/cm), unit durations macrostrat_hiatus.csv hiatus (gap) duration distribution per region macrostrat_sadler.csv accumulation rate vs. time span (Sadler scaling) fits macrostrat_summary.json """ import os, json import numpy as np, pandas as pd from common import RAW, DATA, save_json, A_LAND MS = os.path.join(RAW, "macrostrat") CONT = ["north_america", "caribbean", "south_america", "new_zealand", "africa"] TMAX = 3500 GRID = np.arange(0, TMAX) + 0.5 # mid-points of 1-Myr bins BINS = [(0, 0.0117, "Holocene"), (0.0117, 2.58, "Pleistocene"), (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)")] def parse_lith(s): """Return (sedimentary fraction, metasedimentary fraction) from a Macrostrat lith string.""" if not isinstance(s, str) or not s.strip(): return np.nan, np.nan sed = meta = tot = 0.0 for part in s.split("|"): if "~" not in part: continue desc, prop = part.rsplit("~", 1) try: p = float(prop) except ValueError: continue toks = desc.split() cls = toks[-1] if toks else "" typ = toks[-2] if len(toks) > 1 else "" tot += p if cls == "sedimentary": sed += p if typ == "metasedimentary": meta += p if tot <= 0: return np.nan, np.nan return sed / tot, meta / tot def parse_env(s): """Return (has_marine, has_nonmarine) from a Macrostrat environ string.""" if not isinstance(s, str) or not s.strip(): return False, False mar = non = False for part in s.split("|"): toks = part.split() if not toks: continue if toks[-1] == "non-marine": non = True elif toks[-1] == "marine": mar = True return mar, non def load_units(name): u = pd.read_csv(os.path.join(MS, f"units_{name}.csv"), low_memory=False) u["region"] = name sf = u["lith"].map(parse_lith) u["sed_frac"] = [a for a, b in sf] u["meta_frac"] = [b for a, b in sf] ev = u["environ"].map(parse_env) u["marine"] = [a for a, b in ev] u["nonmarine"] = [b for a, b in ev] u["dur"] = u["b_age"] - u["t_age"] return u def load_cols(name): c = pd.read_csv(os.path.join(MS, f"columns_{name}.csv"), low_memory=False) c["region"] = name return c def coverage(units, cols): """Per-Myr fraction of columns (by count and by area) holding a unit that spans each age.""" ids = cols["col_id"].values idx = {c: i for i, c in enumerate(ids)} area = cols["col_area"].values.astype(float) M = np.zeros((len(ids), TMAX), bool) for cid, t, b in zip(units["col_id"].values, units["t_age"].values, units["b_age"].values): if cid not in idx or not (b > t): continue i0 = int(np.floor(t)); i1 = int(np.ceil(b)) # a unit covers the bin mid-point t_mid if t <= t_mid < b lo = int(np.searchsorted(GRID, t, "left")); hi = int(np.searchsorted(GRID, b, "left")) if hi > lo: M[idx[cid], lo:hi] = True elif i0 < TMAX: # unit shorter than the grid step and not containing a mid-point: assign its bin M[idx[cid], min(i0, TMAX - 1)] = True n = M.sum(0) return M, n / len(ids), (M * area[:, None]).sum(0) / area.sum(), n def wquant(x, w, q): o = np.argsort(x); x = np.asarray(x)[o]; w = np.asarray(w)[o] cw = np.cumsum(w) - 0.5 * w cw /= w.sum() return np.interp(q, cw, x) def main(): out = {} U = {n: load_units(n) for n in CONT + ["deep_sea"]} C = {n: load_cols(n) for n in CONT + ["deep_sea"]} out["n_columns"] = {n: int(len(C[n])) for n in C} out["column_area_km2"] = {n: float(C[n]["col_area"].sum()) for n in C} out["n_units"] = {n: int(len(U[n])) for n in U} out["median_column_area_km2"] = {n: float(C[n]["col_area"].median()) for n in C} cov = pd.DataFrame({"age_ma": GRID}) Msed = {} for n in CONT + ["deep_sea"]: u = U[n]; c = C[n] sed = u[(u["sed_frac"] >= 0.5)] out.setdefault("n_sed_units", {})[n] = int(len(sed)) M, fc, fa, cnt = coverage(sed, c) Msed[n] = M cov[f"{n}_sed_fcols"] = fc; cov[f"{n}_sed_farea"] = fa; cov[f"{n}_sed_ncols"] = cnt for lab, sel in (("marine", sed[sed["marine"]]), ("nonmarine", sed[sed["nonmarine"]]), ("envunknown", sed[~sed["marine"] & ~sed["nonmarine"]])): _, fc2, fa2, cnt2 = coverage(sel, c) cov[f"{n}_{lab}_fcols"] = fc2; cov[f"{n}_{lab}_farea"] = fa2 # any rock at all (igneous, metamorphic included) _, fc3, fa3, _ = coverage(u, c) cov[f"{n}_anyrock_farea"] = fa3 # sedimentary rock that crops out at the surface (Macrostrat 'outcrop' = surface or both) oc = sed[sed["outcrop"].isin(["surface", "both"])] _, _, fa4, _ = coverage(oc, c) cov[f"{n}_sed_outcrop_farea"] = fa4 # pooled continental (area-weighted over all continental column projects) atot = sum(out["column_area_km2"][n] for n in CONT) for lab in ("sed", "marine", "nonmarine"): cov[f"continental_{lab}_farea"] = sum(cov[f"{n}_{lab}_farea"] * out["column_area_km2"][n] for n in CONT) / atot cov.to_csv(os.path.join(DATA, "macrostrat_coverage_1myr.csv"), index=False, float_format="%.5g") # ---- binned table with resolution ------------------------------------------------- rows = [] na = U["north_america"]; ds = U["deep_sea"] for lo, hi, lab in BINS: r = {"bin": lab, "age_top_ma": lo, "age_base_ma": hi} sel = (GRID >= lo) & (GRID < hi) if sel.sum() == 0: # bins shorter than 1 Myr: use the first grid point sel = np.zeros(TMAX, bool); sel[0] = True for n in ("north_america", "south_america", "caribbean", "deep_sea"): r[f"{n}_sed_farea"] = float(cov.loc[sel, f"{n}_sed_farea"].mean()) r["na_marine_farea"] = float(cov.loc[sel, "north_america_marine_farea"].mean()) r["na_nonmarine_farea"] = float(cov.loc[sel, "north_america_nonmarine_farea"].mean()) r["na_envunknown_farea"] = float(cov.loc[sel, "north_america_envunknown_farea"].mean()) r["na_sed_outcrop_farea"] = float(cov.loc[sel, "north_america_sed_outcrop_farea"].mean()) r["na_sed_ncols_mean"] = float(cov.loc[sel, "north_america_sed_ncols"].mean()) r["na_sed_area_km2"] = r["north_america_sed_farea"] * out["column_area_km2"]["north_america"] for nm, u in (("na", na), ("ds", ds)): s = u[(u["sed_frac"] >= 0.5) & (u["dur"] > 0)] ov = np.minimum(s["b_age"], hi) - np.maximum(s["t_age"], lo) s = s[ov > 0]; ov = ov[ov > 0] r[f"{nm}_n_units"] = int(len(s)) if len(s) == 0: continue r[f"{nm}_unit_dur_med_myr"] = float(wquant(s["dur"].values, ov.values, 0.5)) # overlap-time-weighted mids = 0.5 * (s["t_age"].values + s["b_age"].values) if ((mids >= lo) & (mids < hi)).sum() >= 5: r[f"{nm}_unit_dur_med_myr_unweighted"] = float(np.median(s["dur"].values[(mids >= lo) & (mids < hi)])) for envlab, ss, oo in (("all", s, ov), ("marine", s[s["marine"] & ~s["nonmarine"]], ov[s["marine"] & ~s["nonmarine"]]), ("nonmarine", s[s["nonmarine"] & ~s["marine"]], ov[s["nonmarine"] & ~s["marine"]])): th = ss["max_thick"].values.astype(float) ok = th > 0 if ok.sum() < 5: continue ypc = ss["dur"].values[ok] * 1e6 / (th[ok] * 100.0) # years per cm w = oo.values[ok] q = wquant(ypc, w, [0.16, 0.5, 0.84]) # unweighted (per-unit) median over units whose mid-age falls in the bin midu = 0.5 * (ss["t_age"].values + ss["b_age"].values) inb = ok & (midu >= lo) & (midu < hi) if inb.sum() >= 5: r[f"{nm}_{envlab}_yr_per_cm_med_unweighted"] = float(np.median(ss["dur"].values[inb] * 1e6 / (th[inb] * 100.0))) r[f"{nm}_{envlab}_yr_per_cm_p16"] = float(q[0]) r[f"{nm}_{envlab}_yr_per_cm_med"] = float(q[1]) r[f"{nm}_{envlab}_yr_per_cm_p84"] = float(q[2]) r[f"{nm}_{envlab}_thick_med_m"] = float(wquant(th[ok], w, 0.5)) # same with min_thick where given (robustness) mt = ss["min_thick"].values.astype(float) ok2 = mt > 0 if ok2.sum() >= 5: r[f"{nm}_{envlab}_yr_per_cm_med_minthick"] = float( wquant(ss["dur"].values[ok2] * 1e6 / (mt[ok2] * 100.0), oo.values[ok2], 0.5)) rows.append(r) bins = pd.DataFrame(rows) bins.to_csv(os.path.join(DATA, "macrostrat_bins.csv"), index=False, float_format="%.5g") # ---- hiatuses: gaps between sedimentary packages within a column -------------------- hrows = [] hi_summary = {} for n in ("north_america", "south_america", "deep_sea"): M = Msed[n] gaps = []; frac_phan = []; frac_500 = []; frac_2400 = [] for i in range(M.shape[0]): row = M[i] if not row.any(): frac_phan.append(0.0); frac_500.append(0.0); frac_2400.append(0.0) continue frac_phan.append(row[:539].mean()); frac_500.append(row[:500].mean()); frac_2400.append(row[:2400].mean()) first, last = np.argmax(row), len(row) - 1 - np.argmax(row[::-1]) seg = row[first:last + 1] # run lengths of False inside [first, last] d = np.diff(np.concatenate([[1], seg.astype(int), [1]])) starts = np.where(d == -1)[0]; ends = np.where(d == 1)[0] for s0, e0 in zip(starts, ends): gaps.append((first + s0 + 0.5, e0 - s0)) gaps = np.array(gaps) if gaps else np.zeros((0, 2)) gp = gaps[gaps[:, 0] < 539] if len(gaps) else gaps hi_summary[n] = { "n_columns": int(M.shape[0]), "frac_time_with_sediment_last500Myr_median_column": float(np.median(frac_500)), "frac_time_with_sediment_last500Myr_mean_column": float(np.mean(frac_500)), "frac_time_with_sediment_last500Myr_p16": float(np.percentile(frac_500, 16)), "frac_time_with_sediment_last500Myr_p84": float(np.percentile(frac_500, 84)), "frac_time_with_sediment_last2400Myr_mean_column": float(np.mean(frac_2400)), "frac_time_with_sediment_last2400Myr_median_column": float(np.median(frac_2400)), "n_hiatus_phanerozoic": int(len(gp)), "hiatus_median_myr_phanerozoic": float(np.median(gp[:, 1])) if len(gp) else None, "hiatus_p16_myr": float(np.percentile(gp[:, 1], 16)) if len(gp) else None, "hiatus_p84_myr": float(np.percentile(gp[:, 1], 84)) if len(gp) else None, "hiatus_mean_myr": float(np.mean(gp[:, 1])) if len(gp) else None, "hiatuses_per_column_phanerozoic": float(len(gp) / M.shape[0]), } for a, dlen in gaps: hrows.append((n, a, dlen)) pd.DataFrame(hrows, columns=["region", "gap_top_age_ma", "gap_duration_myr"]).to_csv( os.path.join(DATA, "macrostrat_hiatus.csv"), index=False) out["hiatus"] = hi_summary # ---- Sadler scaling: accumulation rate vs. time span -------------------------------- srows = [] for nm, u, sel in (("na_marine", na, na["marine"] & ~na["nonmarine"]), ("na_nonmarine", na, na["nonmarine"] & ~na["marine"]), ("deep_sea", ds, np.ones(len(ds), bool))): s = u[sel & (u["sed_frac"] >= 0.5) & (u["dur"] > 0) & (u["max_thick"] > 0)] x = np.log10(s["dur"].values * 1e6) # yr y = np.log10(s["max_thick"].values * 100.0 / (s["dur"].values * 1e6)) # cm / yr A = np.vstack([x, np.ones_like(x)]).T coef, res, *_ = np.linalg.lstsq(A, y, rcond=None) rng = np.random.default_rng(1) bs = [] for _ in range(300): k = rng.integers(0, len(x), len(x)) bs.append(np.linalg.lstsq(A[k], y[k], rcond=None)[0][0]) srows.append({"set": nm, "n_units": int(len(s)), "slope_lograte_vs_logspan": float(coef[0]), "slope_err_boot": float(np.std(bs)), "intercept": float(coef[1]), "span_p05_yr": float(10 ** np.percentile(x, 5)), "span_p95_yr": float(10 ** np.percentile(x, 95)), "median_rate_cm_per_kyr": float(10 ** np.median(y) * 1e3)}) pd.DataFrame(srows).to_csv(os.path.join(DATA, "macrostrat_sadler.csv"), index=False, float_format="%.5g") out["sadler"] = srows # ---- headline scalars -------------------------------------------------------------- def tavg(col, t0, t1): sel = (GRID >= t0) & (GRID < t1) return float(cov.loc[sel, col].mean()) out["mean_coverage"] = { k: {"last_500Myr": tavg(k, 0, 500), "last_2400Myr": tavg(k, 0, 2400), "500_2400Myr": tavg(k, 500, 2400), "last_66Myr": tavg(k, 0, 66), "66_252Myr": tavg(k, 66, 252), "252_539Myr": tavg(k, 252, 539)} for k in ("north_america_sed_farea", "north_america_marine_farea", "north_america_nonmarine_farea", "north_america_envunknown_farea", "north_america_sed_outcrop_farea", "south_america_sed_farea", "south_america_nonmarine_farea", "caribbean_sed_farea", "continental_sed_farea", "continental_nonmarine_farea", "deep_sea_sed_fcols")} out["continental_column_area_total_km2"] = atot out["continental_column_area_fraction_of_land"] = atot / A_LAND na_sed = na[na["sed_frac"] >= 0.5] out["na_sed_unit_duration_median_myr"] = float(na_sed["dur"].median()) # ---- unweighted (per-unit) medians, to set beside the time-weighted ones in macrostrat_bins.csv. # The time-weighted median answers "how long is the unit that holds a random moment of preserved # time?"; the unweighted one answers "how long is a typical unit?". Both are reported (referee point 6). unw = {} mid_na = 0.5 * (na_sed["t_age"] + na_sed["b_age"]) ds_sed = ds[ds["sed_frac"] >= 0.5] mid_ds = 0.5 * (ds_sed["t_age"] + ds_sed["b_age"]) for lo, hi, lab in ((0, 66, "0-66 Ma"), (66, 252, "66-252 Ma"), (252, 539, "252-539 Ma"), (539, 2400, "539-2400 Ma"), (0, 539, "0-539 Ma")): r = {} for nm, u, mid in (("na", na_sed, mid_na), ("ds", ds_sed, mid_ds)): s = u[(mid >= lo) & (mid < hi) & (u["dur"] > 0)] r[f"{nm}_n_units"] = int(len(s)) if len(s) < 5: continue r[f"{nm}_unit_dur_median_myr_unweighted"] = float(s["dur"].median()) for envlab, ss in (("all", s), ("marine", s[s["marine"] & ~s["nonmarine"]]), ("nonmarine", s[s["nonmarine"] & ~s["marine"]])): ok = ss["max_thick"].values.astype(float) > 0 if ok.sum() < 5: continue ypc = ss["dur"].values[ok] * 1e6 / (ss["max_thick"].values.astype(float)[ok] * 100.0) r[f"{nm}_{envlab}_yr_per_cm_median_unweighted"] = float(np.median(ypc)) r[f"{nm}_{envlab}_n_with_thickness"] = int(ok.sum()) # time-weighted counterparts over the same era (weights = overlap with the era) s2 = u[(u["dur"] > 0)] ov = np.minimum(s2["b_age"], hi) - np.maximum(s2["t_age"], lo) s2 = s2[ov > 0]; ov = ov[ov > 0] if len(s2) >= 5: r[f"{nm}_unit_dur_median_myr_timeweighted"] = float(wquant(s2["dur"].values, ov.values, 0.5)) ok = s2["max_thick"].values.astype(float) > 0 if ok.sum() >= 5: ypc = s2["dur"].values[ok] * 1e6 / (s2["max_thick"].values.astype(float)[ok] * 100.0) r[f"{nm}_all_yr_per_cm_median_timeweighted"] = float(wquant(ypc, ov.values[ok], 0.5)) unw[lab] = r out["unit_statistics_weighted_vs_unweighted"] = unw out["na_fraction_sed_units_env_unknown"] = float((~na_sed["marine"] & ~na_sed["nonmarine"]).mean()) save_json(out, "macrostrat_summary.json") print(json.dumps(out["mean_coverage"], indent=1)) print(json.dumps(out["hiatus"], indent=1)) print(json.dumps(out["unit_statistics_weighted_vs_unweighted"], indent=1)) print(pd.DataFrame(srows).to_string()) pd.set_option("display.width", 250) print(bins[["bin", "north_america_sed_farea", "na_marine_farea", "na_nonmarine_farea", "na_sed_outcrop_farea", "south_america_sed_farea", "deep_sea_sed_farea", "na_unit_dur_med_myr", "na_all_yr_per_cm_med", "na_nonmarine_yr_per_cm_med", "ds_all_yr_per_cm_med", "ds_unit_dur_med_myr"]].to_string()) if __name__ == "__main__": main()