worldgen

git clone https://git.godosa.eu/worldgen

master

history zip tar.gz

raw · 8431 bytes

  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}