1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
|
"""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}
|
