"""How much of ancient land-surface space-time have palaeontologists actually looked at? Every Paleobiology Database collection is treated as one close visual inspection of a small patch of bedding of known place and age. An installation of area A that stood for tau years at age t counts as 'sampled' when at least one collection lies inside its footprint and the collection's (uncertain) age window overlaps the occupation. For a collection whose age is only known to lie in an interval of width W, the overlap probability is min(1, (tau + h) / W), h being the time one collection represents (time-averaging plus the thickness collected). P_sampled(t; A, tau) = sum over equal-area cells of side sqrt(A) of [1 - prod_i (1 - p_i)] * A / A_land Because fossil collections exist only where sediment was deposited, survived and is exposed, this single empirical number folds together preservation, exposure and examination. Outputs: pbdb_density_1myr.csv, pbdb_psampled_slices.csv, pbdb_period_table.csv, pbdb_summary.json """ import os, glob, json import numpy as np, pandas as pd from common import RAW, DATA, save_json, A_LAND, A_EARTH TERR = {"terrestrial indet.", "fluvial indet.", "fluvial-lacustrine indet.", "lacustrine indet.", "lacustrine - large", "lacustrine - small", "cave", '"floodplain"', '"channel"', "wet floodplain", "dry floodplain", "mire/swamp", "pond", "fissure fill", "crevasse splay", "fine channel fill", "coarse channel fill", "channel lag", "crater lake", "glacial", "eolian indet.", "alluvial fan", "karst indet.", "lacustrine deltaic indet.", "dune", "sinkhole", "interdune", "tar", "loess", "lacustrine delta plain", "spring", "levee", "lacustrine delta front", "lacustrine prodelta", "lacustrine interdistributary bay", "fluvial-deltaic indet."} R_E = 6371.0 SCALES_KM = [1.0, 10.0, 100.0] # footprint side; areas 1, 100, 10,000 km2 TAUS_MYR = [3e-4, 1e-2, 1.0] # occupation 300 yr, 10 kyr, 1 Myr H_MYR = 1e-2 # time represented by one collection (10 kyr); varied in s07 PERIODS = [(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, 2500, "Proterozoic")] def load(): fs = sorted(glob.glob(os.path.join(RAW, "pbdb", "colls_*.csv"))) use = ["collection_no", "lng", "lat", "n_occs", "max_ma", "min_ma", "cc", "formation", "stratscale", "geogscale", "environment"] d = pd.concat([pd.read_csv(f, usecols=use, low_memory=False) for f in fs]) n_raw = len(d) d = d.drop_duplicates("collection_no") n_unique = len(d) d = d.dropna(subset=["lng", "lat", "max_ma", "min_ma"]) d = d[(d["max_ma"] > d["min_ma"])].copy() d["W"] = d["max_ma"] - d["min_ma"] d["terr"] = d["environment"].isin(TERR) d["env_known"] = d["environment"].notna() return d.reset_index(drop=True), n_raw, n_unique def cell_id(lat, lng, s_km): """Equal-area (Lambert cylindrical) cells of side s_km.""" x = np.floor(R_E * np.radians(lng) / s_km).astype(np.int64) y = np.floor(R_E * np.sin(np.radians(lat)) / s_km).astype(np.int64) return x * 100003 + y def p_sampled(d, t, s_km, tau, h=H_MYR, by_formation=False, area_norm=A_LAND): sel = (d["min_ma"].values <= t) & (d["max_ma"].values > t) if sel.sum() == 0: return 0.0, 0 sub = d[sel] p = np.minimum(1.0, (tau + h) / sub["W"].values) cid = cell_id(sub["lat"].values, sub["lng"].values, s_km) if by_formation: # collections from the same formation in the same cell count once (they revisit the same beds) fm = sub["formation"].fillna("?").astype(str).values g = pd.DataFrame({"c": cid, "f": fm, "p": p}).groupby(["c", "f"])["p"].max().reset_index() cid = g["c"].values; p = g["p"].values lg = np.log1p(-np.minimum(p, 1 - 1e-12)) s = pd.Series(lg).groupby(cid).sum() pc = 1.0 - np.exp(s.values) return float(min(1.0, pc.sum() * s_km * s_km / area_norm)), int(len(pc)) def main(): d, n_raw, n_unique = load() out = {"n_rows_downloaded": int(n_raw), "n_unique_collections": int(n_unique), "n_used": int(len(d)), "n_terrestrial": int(d["terr"].sum()), "n_env_unknown": int((~d["env_known"]).sum()), "median_age_window_myr": float(d["W"].median()), "h_myr": H_MYR} # --- density of collections per Myr of geological time ------------------------------ grid = np.arange(0, 600) + 0.5 def dens(sub): n = np.zeros(len(grid)) lo = np.searchsorted(grid, sub["min_ma"].values, "left"); hi = np.searchsorted(grid, sub["max_ma"].values, "left") w = 1.0 / np.maximum(sub["W"].values, 1.0) # spread over >= 1 Myr on this 1-Myr grid for a, b, ww, mn in zip(lo, hi, w, sub["min_ma"].values): if b > a: n[a:b] += ww * 1.0 elif a < len(n): n[min(a, len(n) - 1)] += 1.0 # window narrower than the grid step return n dd = pd.DataFrame({"age_ma": grid, "colls_per_myr_all": dens(d), "colls_per_myr_terrestrial": dens(d[d["terr"]])}) dd.to_csv(os.path.join(DATA, "pbdb_density_1myr.csv"), index=False, float_format="%.5g") # --- P_sampled on time slices -------------------------------------------------------- slices = np.concatenate([[0.005, 0.05, 0.3, 1.0, 2.0], np.arange(3.0, 540.0, 2.0)]) wts = np.gradient(slices); wts[0] = 0.0117; wts[1] = 0.1; wts[2] = 0.5; wts[3] = 0.9; wts[4] = 1.0 rows = [] terr = d[d["terr"]] for t in slices: r = {"age_ma": float(t), "n_terr_spanning": int(((terr["min_ma"] <= t) & (terr["max_ma"] > t)).sum()), "n_all_spanning": int(((d["min_ma"] <= t) & (d["max_ma"] > t)).sum())} st = terr[(terr["min_ma"] <= t) & (terr["max_ma"] > t)] sa = d[(d["min_ma"] <= t) & (d["max_ma"] > t)] for tau in TAUS_MYR: # expected number of collections worldwide whose age window overlaps the occupation r[f"Nov_terr_tau{tau:g}"] = float(np.minimum(1.0, (tau + H_MYR) / st["W"].values).sum()) r[f"Nov_all_tau{tau:g}"] = float(np.minimum(1.0, (tau + H_MYR) / sa["W"].values).sum()) for s in SCALES_KM: for tau in TAUS_MYR: key = f"A{int(s*s)}_tau{tau:g}" r[f"Pterr_{key}"], r[f"ncell_terr_{int(s)}km"] = p_sampled(terr, t, s, tau) r[f"Pterr_fm_{key}"], _ = p_sampled(terr, t, s, tau, by_formation=True) r[f"Pall_{key}"], _ = p_sampled(d, t, s, tau, area_norm=A_EARTH) rows.append(r) sl = pd.DataFrame(rows) sl["weight_myr"] = wts sl.to_csv(os.path.join(DATA, "pbdb_psampled_slices.csv"), index=False, float_format="%.5g") # --- table by period --------------------------------------------------------------------- prow = [] for lo, hi, lab in PERIODS: mid = 0.5 * (d["max_ma"] + d["min_ma"]) selp = (mid >= lo) & (mid < hi) r = {"period": lab, "age_top_ma": lo, "age_base_ma": hi, "n_collections": int(selp.sum()), "n_terrestrial": int((selp & d["terr"]).sum()), "terr_colls_per_myr": float((selp & d["terr"]).sum() / (hi - lo)), "terr_cells_100km": int(len(np.unique(cell_id(d.loc[selp & d["terr"], "lat"].values, d.loc[selp & d["terr"], "lng"].values, 100.0)))), "terr_cells_10km": int(len(np.unique(cell_id(d.loc[selp & d["terr"], "lat"].values, d.loc[selp & d["terr"], "lng"].values, 10.0))))} r["frac_land_cells_100km_with_terr_coll"] = r["terr_cells_100km"] * 1e4 / A_LAND ss = sl[(sl["age_ma"] >= lo) & (sl["age_ma"] < hi)] if len(ss): for c in sl.columns: if c.startswith("Pterr_"): r[c] = float(np.average(ss[c], weights=ss["weight_myr"])) prow.append(r) pt = pd.DataFrame(prow) pt.to_csv(os.path.join(DATA, "pbdb_period_table.csv"), index=False, float_format="%.4g") # --- time averages ----------------------------------------------------------------------- def tavg(col, t0, t1): ss = sl[(sl["age_ma"] >= t0) & (sl["age_ma"] < t1)] return float(np.average(ss[col], weights=ss["weight_myr"])) av = {} for c in [c for c in sl.columns if c.startswith("P") or c.startswith("Nov_")]: av[c] = {"3_500Myr": tavg(c, 3, 500), "0_3Myr": tavg(c, 0, 3), "3_66Myr": tavg(c, 3, 66), "66_252Myr": tavg(c, 66, 252), "252_539Myr": tavg(c, 252, 539)} out["time_averages"] = av out["terr_colls_per_myr"] = {"3_66Myr": float(dd.loc[(grid > 3) & (grid < 66), "colls_per_myr_terrestrial"].mean()), "66_252Myr": float(dd.loc[(grid > 66) & (grid < 252), "colls_per_myr_terrestrial"].mean()), "252_539Myr": float(dd.loc[(grid > 252) & (grid < 539), "colls_per_myr_terrestrial"].mean())} out["n_precambrian_collections"] = int((d["min_ma"] >= 538.8).sum()) out["n_precambrian_terrestrial"] = int(((d["min_ma"] >= 538.8) & d["terr"]).sum()) save_json(out, "pbdb_summary.json") pd.set_option("display.width", 250) print(json.dumps({k: v for k, v in out.items() if k != "time_averages"}, indent=1)) print(pt[["period", "n_collections", "n_terrestrial", "terr_colls_per_myr", "terr_cells_100km", "frac_land_cells_100km_with_terr_coll", "Pterr_A10000_tau0.0003", "Pterr_A10000_tau0.01", "Pterr_A10000_tau1", "Pterr_fm_A10000_tau1", "Pterr_A100_tau0.01", "Pterr_A1_tau0.01"]].to_string()) for c in ("Pterr_A10000_tau0.0003", "Pterr_A10000_tau0.01", "Pterr_A10000_tau1", "Pterr_fm_A10000_tau1", "Pterr_A100_tau0.01", "Pterr_A100_tau1", "Pterr_A1_tau0.01", "Pterr_A1_tau1"): print(c, av[c]) if __name__ == "__main__": main()