"""Stage `environment`: Holdridge life zone + seasonality tag + lithology + landform + ground.""" from __future__ import annotations import numpy as np from .config import params from .crust import CRATON, LIP, POST_OROGEN, PRE_OROGEN, RIFT, SCAR, VOLCANO from .graph import nbr_max, nbr_min, steepest_receivers from . import minerals as MN from .noise import fbm REGIONS = ["polar", "subpolar", "boreal", "cool temperate", "warm temperate", "subtropical", "tropical"] ZONES = { "polar": ["desert"], "subpolar": ["dry tundra", "moist tundra", "wet tundra", "rain tundra"], "boreal": ["desert", "dry scrub", "moist forest", "wet forest", "rain forest"], "cool temperate": ["desert", "desert scrub", "steppe", "moist forest", "wet forest", "rain forest"], "warm temperate": ["desert", "desert scrub", "thorn steppe", "dry forest", "moist forest", "wet forest", "rain forest"], "subtropical": ["desert", "desert scrub", "thorn woodland", "dry forest", "moist forest", "wet forest", "rain forest"], "tropical": ["desert", "desert scrub", "thorn woodland", "very dry forest", "dry forest", "moist forest", "wet forest", "rain forest"], } THRESH = { # annual precipitation (mm) upper bounds between successive zones "polar": [], "subpolar": [125, 250, 500], "boreal": [125, 250, 500, 1000], "cool temperate": [125, 250, 500, 1000, 2000], "warm temperate": [125, 250, 500, 1000, 2000, 4000], "subtropical": [125, 250, 500, 1000, 2000, 4000], "tropical": [125, 250, 500, 1000, 2000, 4000, 8000], } HOLDRIDGE_NAMES = [f"{r} {z}" for r in REGIONS for z in ZONES[r]] SEAS_F, SEAS_S, SEAS_W, SEAS_M = 0, 1, 2, 3 SEASONALITY_NAMES = ["even rainfall", "dry summer", "dry winter / summer monsoon", "tropical monsoon"] LI_OCEAN, LI_GRANITE, LI_LIMESTONE, LI_BASALT, LI_ANDESITE, LI_SANDSTONE, LI_METAMORPHIC = range(7) LITHOLOGY_NAMES = ["oceanic basalt", "granite/gneiss", "limestone", "basalt", "andesite", "sandstone/shale", "metamorphic"] (LF_OCEAN, LF_PLAIN, LF_HILLS, LF_MOUNTAINS, LF_PLATEAU, LF_RIFT, LF_ESCARPMENT, LF_ARC, LF_MASSIF, LF_BASALT_PLATEAU, LF_DUNES, LF_BADLANDS) = range(12) LANDFORM_NAMES = ["ocean", "plain", "hills", "mountains", "plateau", "rift valley", "escarpment", "volcanic arc", "volcanic massif", "basalt plateau", "dunes", "badlands"] (GR_NONE, GR_WETLAND, GR_BOG, GR_FLOODPLAIN, GR_DELTA, GR_MANGROVE, GR_SALT_FLAT, GR_PERMAFROST, GR_KARST) = range(9) GROUND_NAMES = ["—", "wetland", "bog", "floodplain", "delta", "mangrove", "salt flat", "permafrost", "karst"] LEGENDS = {"deposits": MN.LEGEND, "holdridge": HOLDRIDGE_NAMES, "seasonality": SEASONALITY_NAMES, "landform": LANDFORM_NAMES, "lithology": LITHOLOGY_NAMES, "ground": GROUND_NAMES} DEFAULTS = {"relief_km": 100.0, "coastal_km": 60.0, "dry_ratio_s": 3.0, "dry_ratio_w": 4.0, "monsoon_min_mm": 1500.0, "flat_m_per_km": 1.0, "wet_min_mm": 700.0, "karst_min_mm": 800.0, "small_river_km3_yr": 5.0, "life": True} def holdridge(bio, p, tmin): region = np.select([bio < 1.5, bio < 3, bio < 6, bio < 12, (bio < 24) & (tmin < 0), bio < 24], [0, 1, 2, 3, 4, 5], default=6).astype(np.int8) zone = np.zeros(len(bio), np.int16) offset = 0 for ri, r in enumerate(REGIONS): sel = region == ri zone[sel] = offset + np.searchsorted(THRESH[r], p[sel], side="right") offset += len(ZONES[r]) return zone, region def seasonality(p_jun, p_dec, lat, tmin, p_ann, P): summer = np.where(lat >= 0, p_jun, p_dec) winter = np.where(lat >= 0, p_dec, p_jun) tag = np.full(len(lat), SEAS_F, np.int8) tag[summer < winter / P["dry_ratio_s"]] = SEAS_S w = winter < summer / P["dry_ratio_w"] tag[w] = SEAS_W tag[w & (tmin >= 18.0) & (p_ann >= P["monsoon_min_mm"])] = SEAS_M return tag def lithology(g, z, continental, age, d_over, seed): nz = fbm(g.xyz, seed + 51, 4, 5.0) lit = np.full(g.n, LI_GRANITE, np.int8) lit[continental & (z < 600) & (nz > 0.1)] = LI_SANDSTONE lit[continental & (z < 300) & (nz < -0.1)] = LI_LIMESTONE lit[np.isin(age, [POST_OROGEN, PRE_OROGEN]) & (z > 2000)] = LI_METAMORPHIC lit[~continental] = LI_OCEAN lit[(d_over > 100) & (d_over < 320)] = LI_ANDESITE lit[np.isin(age, [LIP, VOLCANO, SCAR])] = LI_BASALT return lit def landform(g, z, age, lit, d_over, p_ann, ocean=None, relief_km=100.0): hi, lo = np.asarray(z, dtype=np.float64), np.asarray(z, dtype=np.float64) for _ in range(max(1, int(round(relief_km / g.spacing_km)))): # window ≈ relief_km at any resolution hi, lo = nbr_max(g, hi), nbr_min(g, lo) relief = hi - lo lf = np.full(g.n, LF_PLAIN, np.int8) lf[relief >= 300] = LF_HILLS lf[(z > 1000) & (relief < 500)] = LF_PLATEAU lf[(relief >= 1500) | (z > 2500)] = LF_MOUNTAINS lf[(p_ann < 250) & (relief < 300) & (lit == LI_SANDSTONE)] = LF_DUNES lf[(p_ann >= 250) & (p_ann < 500) & (relief >= 200) & (relief < 800) & (lit == LI_SANDSTONE)] = LF_BADLANDS lf[age == RIFT] = LF_RIFT lf[age == LIP] = LF_BASALT_PLATEAU ocean = z <= 0 if ocean is None else ocean lf[(d_over > 100) & (d_over < 320) & ~ocean] = LF_ARC lf[age == VOLCANO] = LF_MASSIF lf[age == SCAR] = LF_ESCARPMENT lf[ocean] = LF_OCEAN return lf, relief def ground(g, z, lit, t_mean, t_min, p_ann, dist_ocean, river, strahler, discharge, lake, salt_flat, P, ocean=None): land = ~(z <= 0 if ocean is None else ocean) _, slope, _ = steepest_receivers(g, z) flat = slope < P["flat_m_per_km"] coastal = land & (dist_ocean < P["coastal_km"]) wet = land & flat & ~lake & (p_ann > P["wet_min_mm"]) & (discharge > P["small_river_km3_yr"]) gr = np.zeros(g.n, np.int8) gr[land & (t_mean < -2.0)] = GR_PERMAFROST gr[land & (lit == LI_LIMESTONE) & (p_ann > P["karst_min_mm"])] = GR_KARST gr[wet & (t_mean < 5.0)] = GR_BOG gr[wet & (t_mean >= 5.0)] = GR_WETLAND gr[river & (strahler >= 3) & flat] = GR_FLOODPLAIN gr[coastal & flat & (t_min >= 20.0) & (p_ann > 1000.0)] = GR_MANGROVE gr[river & (strahler >= 4) & coastal] = GR_DELTA gr[salt_flat] = GR_SALT_FLAT if not P["life"]: gr[np.isin(gr, [GR_WETLAND, GR_BOG, GR_MANGROVE])] = GR_NONE return gr def run(ctx) -> dict: g = ctx.grid P = params(ctx.cfg, "environment", DEFAULTS) (z, cont, age, d_over, t_mean, t_min, p_ann, p_jun, p_dec, bio, dist_ocean, river, strahler, q, lake, salt) = ctx.need("elevation_eroded_m", "continental", "age_class", "d_over_km", "T_mean", "T_min", "P_ann", "P_jun", "P_dec", "biotemp", "dist_ocean_km", "river", "strahler", "discharge_km3_yr", "lake", "salt_flat") z = z.astype(np.float64) bed = z if "lake_level_m" in ctx.data: # the ground's shape: lakes at their water, not their beds lev = np.asarray(ctx.data["lake_level_m"], dtype=np.float64) z = np.where(np.isfinite(lev), np.maximum(z, lev), z) ocean = np.asarray(ctx.data["ocean"]) if "ocean" in ctx.data else z <= 0 zone, region = holdridge(bio, p_ann, t_min) lit = lithology(g, z, cont, age, d_over, ctx.seed) lf, relief = landform(g, z, age, lit, d_over, p_ann, ocean, P["relief_km"]) gr = ground(g, z, lit, t_mean, t_min, p_ann, dist_ocean, river, strahler, q, lake, salt, P, ocean) coal = (lit == LI_SANDSTONE) & (p_ann > 600) & ~ocean & bool(P["life"]) get = lambda k, v: np.asarray(ctx.data[k]) if k in ctx.data else np.full(g.n, v) deposits, main = MN.place(g, {"land": ~ocean & ~lake, "z": bed, "lit": lit, "age": age, "lf": lf, "relief": relief, "t_mean": t_mean, "p_ann": p_ann, "d_coll": get("d_coll_km", np.inf), "salt": salt, "endo": get("endorheic", False), "dist_ocean": dist_ocean, "coal": coal}, ctx.seed) iron = np.uint32(MN.BIT["bog_iron"] | MN.BIT["ironstone"] | MN.BIT["iron_high"]) return {"deposits": deposits, "deposit_main": main,"holdridge": zone, "hold_region": region, "seasonality": seasonality(p_jun, p_dec, g.lat, t_min, p_ann, P), "landform": lf, "relief_m": relief, "lithology": lit, "ground": gr, "coal_potential": coal, "iron_potential": (deposits & iron) > 0}