"""Global bedrock age / rock-type distribution of the present land surface. 4000 random land points (uniform on the sphere, Natural Earth 1:110m land mask) were queried against the Macrostrat geologic map API. The 'tiny' scale layer is the Geological Survey of Canada world map (Chorlton 2007, Open File 5529) and is the only layer with global coverage, so it is the one analysed. Output: fraction of the land surface whose mapped bedrock is sedimentary, by age, for the globe and for the USA+Canada box (to test whether North America, where the Macrostrat columns are, is representative). """ import os, json, re import numpy as np, pandas as pd from common import RAW, DATA, save_json, A_LAND SED_WORDS = ("sediment", "sandstone", "siltstone", "shale", "carbonate", "limestone", "mudstone", "conglomerate", "clay", "sand", "gravel", "glacial deposits", "turbidite", "dolomite", "evaporite", "laterite", "unconsolidated", "greywacke", "chert", "tillite", "argillite", "alluvi", "loess", "coal", "marl", "arenite") XTAL_WORDS = ("gneiss", "schist", "metamorphic", "migmatite", "intrusive", "volcanic", "basalt", "granit", "ignimbrite", "plutonic", "ophiolite", "amphibolite", "granulite", "meta", "gabbro", "ultramafic", "rhyolite", "andesite", "diorite", "tuff") ERAS = [(0, 2.58, "Quaternary"), (2.58, 66, "Tertiary"), (66, 252, "Mesozoic"), (252, 539, "Paleozoic"), (539, 2500, "Proterozoic"), (2500, 4600, "Archean")] def classify(lith): l = (lith or "").lower() if not l: return "unmapped" s = any(w in l for w in SED_WORDS) x = any(w in l for w in XTAL_WORDS) if s and not x: return "sed" if s and x: return "mixed" return "xtal" def main(): files = [f for f in os.listdir(os.path.join(RAW, "mapsample")) if f.endswith(".jsonl")] recs = {} for f in files: for line in open(os.path.join(RAW, "mapsample", f)): try: r = json.loads(line) except Exception: continue if r.get("units") is None: continue recs[r["i"]] = r rows = [] for i, r in sorted(recs.items()): tiny = [u for u in r["units"] if u.get("scale") == "tiny"] if not tiny: rows.append({"i": i, "lat": r["lat"], "lng": r["lng"], "cls": "nodata", "t_age": np.nan, "b_age": np.nan, "lith": ""}) continue u = tiny[0] rows.append({"i": i, "lat": r["lat"], "lng": r["lng"], "cls": classify(u.get("lith")), "t_age": u.get("t_age"), "b_age": u.get("b_age"), "lith": u.get("lith") or ""}) d = pd.DataFrame(rows) d["t_age"] = pd.to_numeric(d["t_age"], errors="coerce"); d["b_age"] = pd.to_numeric(d["b_age"], errors="coerce") d["dur"] = d["b_age"] - d["t_age"] d["usca"] = (d["lat"] > 25) & (d["lat"] < 84) & (d["lng"] > -170) & (d["lng"] < -52) & ~((d["lng"] > -75) & (d["lat"] > 59.5) & (d["lng"] > -73)) d["antarctica"] = d["lat"] < -60 d.to_csv(os.path.join(DATA, "mapsample_points.csv"), index=False) n = len(d) out = {"n_points": int(n), "class_fractions_all_land": (d["cls"].value_counts() / n).to_dict(), "n_usca_box": int(d["usca"].sum())} icefree = d[~d["cls"].isin(["unmapped", "nodata"])] out["n_mapped"] = int(len(icefree)) out["mapped_fraction_of_land"] = float(len(icefree) / n) out["sed_fraction_of_mapped"] = float(((icefree["cls"] == "sed") + 0.5 * (icefree["cls"] == "mixed")).sum() / len(icefree)) out["median_map_unit_duration_myr_sed"] = float(icefree.loc[icefree["cls"] == "sed", "dur"].median()) # fraction of land surface underlain by sedimentary bedrock whose mapped age range spans age t grid = np.concatenate([np.arange(0.5, 100, 1.0), np.arange(102.5, 3500, 5.0)]) def curve(sub, denom): w = (sub["cls"] == "sed") * 1.0 + (sub["cls"] == "mixed") * 0.5 t = sub["t_age"].values; b = sub["b_age"].values f_span = np.array([np.sum(w[(t <= g) & (b > g)]) for g in grid]) / denom # 'slice' estimate: share of each map unit's outcrop attributed to a 1-Myr slice # (uniform-in-time assumption within the mapped age range) dur = np.maximum(b - t, 0.01) f_slice = np.array([np.sum((w / dur)[(t <= g) & (b > g)]) for g in grid]) / denom return f_span, f_slice fs_g, fl_g = curve(d, n) us = d[d["usca"]] fs_u, fl_u = curve(us, len(us)) pd.DataFrame({"age_ma": grid, "global_sed_outcrop_frac_spanning": fs_g, "global_sed_outcrop_frac_per_myr_slice": fl_g, "usca_sed_outcrop_frac_spanning": fs_u, "usca_sed_outcrop_frac_per_myr_slice": fl_u}).to_csv( os.path.join(DATA, "mapsample_outcrop_by_age.csv"), index=False, float_format="%.5g") # era table (mid-age assignment) mid = 0.5 * (d["t_age"] + d["b_age"]) era_rows = [] for lo, hi, lab in ERAS: r = {"era": lab} for nm, sub in (("global", d), ("usca", us)): m = mid.loc[sub.index] sel = (m >= lo) & (m < hi) ns = ((sub["cls"] == "sed") & sel).sum() + 0.5 * ((sub["cls"] == "mixed") & sel).sum() nx = ((sub["cls"] == "xtal") & sel).sum() + 0.5 * ((sub["cls"] == "mixed") & sel).sum() r[f"{nm}_sed_frac_of_land"] = float(ns / len(sub)) r[f"{nm}_sed_frac_err"] = float(np.sqrt(max(ns, 1)) / len(sub)) r[f"{nm}_xtal_frac_of_land"] = float(nx / len(sub)) r["global_sed_area_km2"] = r["global_sed_frac_of_land"] * A_LAND era_rows.append(r) era = pd.DataFrame(era_rows) era.to_csv(os.path.join(DATA, "mapsample_era_table.csv"), index=False, float_format="%.4g") out["era_table"] = era_rows # ratio USA+Canada / globe for pre-Quaternary sedimentary outcrop, a representativeness check g = sum(r["global_sed_frac_of_land"] for r in era_rows[1:4]); u = sum(r["usca_sed_frac_of_land"] for r in era_rows[1:4]) out["phanerozoic_preQ_sed_outcrop_frac_global"] = float(g) out["phanerozoic_preQ_sed_outcrop_frac_usca"] = float(u) out["ratio_usca_to_global"] = float(u / g) if g > 0 else None # time-averaged outcrop fraction per Myr slice def tavg(f, t0, t1): sel = (grid >= t0) & (grid < t1) wts = np.gradient(grid)[sel] return float(np.sum(f[sel] * wts) / np.sum(wts)) out["global_sed_outcrop_frac_per_myr_slice"] = {"last_66": tavg(fl_g, 2.58, 66), "66_252": tavg(fl_g, 66, 252), "252_539": tavg(fl_g, 252, 539), "539_2400": tavg(fl_g, 539, 2400), "2.58_500": tavg(fl_g, 2.58, 500), "2.58_2400": tavg(fl_g, 2.58, 2400)} save_json(out, "mapsample_summary.json") print(json.dumps({k: v for k, v in out.items() if k != "era_table"}, indent=1)) print(era.to_string()) if __name__ == "__main__": main()