"""Download every raw input (idempotent; cached under RAW). Sources Macrostrat API v2 (Peters, Husson & Czaplewski 2018): columns, units, sections for each column project; global bedrock map sampled at random land points. Paleobiology Database data service 1.2: all fossil collections. PANGAEA 917660: CENOGRID benthic d13C/d18O splice (Westerhold et al. 2020). Zenodo 14771859: Deep-Time Marine Sedimentary Element Database (DM-SED; Lai et al. 2025). Natural Earth 1:110m land polygons (land mask for the random points). """ import os, sys, json, time, urllib.parse, http.client, ssl, threading import numpy as np from common import RAW, SEED, UA, fetch MS = "https://macrostrat.org/api/v2" # Macrostrat column projects used. 3 = eODP (every offshore drilling hole). MS_PROJECTS = {1: "north_america", 7: "caribbean", 8: "south_america", 5: "new_zealand", 9: "africa", 4: "deep_sea", 10: "na_ediacaran", 3: "eodp"} PBDB = "https://paleobiodb.org/data1.2" PBDB_INTERVALS = ["Quaternary", "Neogene", "Paleogene", "Cretaceous", "Jurassic", "Triassic", "Permian", "Carboniferous", "Devonian", "Silurian", "Ordovician", "Cambrian", "Proterozoic", "Archean"] N_MAP_POINTS = 4000 def macrostrat(): d = os.path.join(RAW, "macrostrat") fetch(f"{MS}/defs/projects?all&format=csv", os.path.join(d, "projects.csv")) fetch(f"{MS}/defs/intervals?timescale_id=1&format=csv", os.path.join(d, "intervals_ts1.csv")) fetch(f"{MS}/defs/intervals?all&format=csv", os.path.join(d, "intervals_all.csv")) for pid, name in MS_PROJECTS.items(): for ep in ("columns", "sections", "units"): if pid == 3 and ep == "units": continue # eODP units (92k rows) exceed the API's response limit (HTTP 500); not needed fetch(f"{MS}/{ep}?project_id={pid}&response=long&format=csv", os.path.join(d, f"{ep}_{name}.csv"), timeout=1200, min_bytes=200) print("macrostrat", ep, name, "ok", flush=True) def pbdb(): d = os.path.join(RAW, "pbdb") for iv in PBDB_INTERVALS: fetch(f"{PBDB}/colls/list.csv?interval={iv}&show=loc,time,strat,stratext,geo", os.path.join(d, f"colls_{iv}.csv"), timeout=1800, min_bytes=200) print("pbdb", iv, "ok", flush=True) def cenogrid(): d = os.path.join(RAW, "cenogrid") fetch("https://doi.pangaea.de/10.1594/PANGAEA.917660?format=textfile", os.path.join(d, "cenogrid_S33.tab"), min_bytes=1_000_000) print("cenogrid ok", flush=True) def dmsed(): """Deep-Time Marine Sedimentary Element Database v0.0.1 (Lai et al. 2025; built on SGP, Farrell et al. 2021), Zenodo doi:10.5281/zenodo.14771859, CC-BY 4.0.""" d = os.path.join(RAW, "dmsed") fetch("https://zenodo.org/api/records/14771859/files/DM-SED%200.0.1.csv/content", os.path.join(d, "DM-SED_0.0.1.csv"), min_bytes=10_000_000) print("dmsed ok", flush=True) def land_points(n, seed=SEED): """n random points uniform on the sphere that fall on Natural Earth 110m land.""" from matplotlib.path import Path p = fetch("https://raw.githubusercontent.com/nvkelso/natural-earth-vector/master/geojson/ne_110m_land.geojson", os.path.join(RAW, "naturalearth", "ne_110m_land.geojson"), min_bytes=10000) gj = json.load(open(p)) polys = [] for f in gj["features"]: g = f["geometry"] rings = [g["coordinates"]] if g["type"] == "Polygon" else g["coordinates"] for poly in rings: polys.append((Path(np.asarray(poly[0])), [Path(np.asarray(h)) for h in poly[1:]])) rng = np.random.default_rng(seed) out = [] while len(out) < n: m = 4 * n lng = rng.uniform(-180, 180, m) lat = np.degrees(np.arcsin(rng.uniform(-1, 1, m))) pts = np.column_stack([lng, lat]) inside = np.zeros(m, bool) for outer, holes in polys: ins = outer.contains_points(pts) for h in holes: ins &= ~h.contains_points(pts) inside |= ins for a, b in pts[inside]: out.append((round(float(a), 4), round(float(b), 4))) return out[:n] def map_sample(n=N_MAP_POINTS, nthreads=2): """Query the Macrostrat geologic map (all scales) at n random land points.""" d = os.path.join(RAW, "mapsample") os.makedirs(d, exist_ok=True) dest = os.path.join(d, f"map_points_{n}.jsonl") pts = land_points(n) done = {} if os.path.exists(dest): for line in open(dest): try: r = json.loads(line) done[r["i"]] = r except Exception: pass todo = [i for i in range(n) if i not in done] print("map sample: have", len(done), "todo", len(todo), flush=True) if not todo: return lock = threading.Lock() fout = open(dest, "a") def worker(idx): ctx = ssl.create_default_context() import urllib.request for i in idx: lng, lat = pts[i] url = f"{MS}/geologic_units/map?lat={lat}&lng={lng}&format=json" rec = None for k in range(4): try: req = urllib.request.Request(url, headers={"User-Agent": UA}) with urllib.request.urlopen(req, timeout=60) as r: js = json.loads(r.read().decode("utf-8")) units = js.get("success", {}).get("data", []) keep = [{k2: u.get(k2) for k2 in ("scale", "source_id", "name", "lith", "descrip", "t_age", "b_age", "t_int_name", "b_int_name", "best_int_name", "liths")} for u in units] rec = {"i": i, "lng": lng, "lat": lat, "units": keep} break except Exception as e: # noqa time.sleep(2 + 3 * k) if rec is None: rec = {"i": i, "lng": lng, "lat": lat, "units": None, "error": True} with lock: fout.write(json.dumps(rec) + "\n") fout.flush() time.sleep(0.05) ths = [threading.Thread(target=worker, args=(todo[k::nthreads],)) for k in range(nthreads)] [t.start() for t in ths] [t.join() for t in ths] fout.close() print("map sample ok", flush=True) def main(which=None): steps = {"macrostrat": macrostrat, "pbdb": pbdb, "cenogrid": cenogrid, "dmsed": dmsed, "map": map_sample} for k, f in steps.items(): if which and k not in which: continue f() if __name__ == "__main__": main(sys.argv[1:] or None)