worldgen

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

master

raw · 1812 bytes

"""Stage `ice`: ice sheets, mountain glaciers, seasonal/perennial sea ice."""
from __future__ import annotations

import numpy as np

from .config import params
from .graph import components, distance_to

ICE_NONE, ICE_SHEET, ICE_GLACIER, ICE_SEA_SEASONAL, ICE_SEA_PERENNIAL = range(5)
ICE_NAMES = ["none", "ice sheet", "glacier", "seasonal sea ice", "perennial sea ice"]
DEFAULTS = {"melt_summer_c": 0.0, "sheet_min_km2": 250000.0, "sheet_min_p_mm": 100.0, "sea_ice_t_c": -1.8, "sheet_max_m": 3000.0,
            "sheet_edge_m": 200.0, "sheet_growth_m_per_km": 3.0}


def run(ctx) -> dict:
    g = ctx.grid
    P = params(ctx.cfg, "ice", DEFAULTS)
    z, t_mean, t_jun, t_dec, p_ann = ctx.need("elevation_eroded_m", "T_mean", "T_jun", "T_dec", "P_ann")
    z = z.astype(np.float64)
    land = ~np.asarray(ctx.data["ocean"]) if "ocean" in ctx.data else z > 0
    summer = np.where(g.lat >= 0, t_jun, t_dec)
    winter = np.where(g.lat >= 0, t_dec, t_jun)
    cold = land & (summer < P["melt_summer_c"])                  # snow survives the summer → permanent ice
    lab = components(g, cold)
    area = np.bincount(lab[cold], weights=g.area_km2[cold]) if cold.any() else np.zeros(1)
    big = cold & (area[np.maximum(lab, 0)] >= P["sheet_min_km2"])
    sheet = big & (p_ann > P["sheet_min_p_mm"])
    thick = np.zeros(g.n)
    if sheet.any():
        d_edge = distance_to(g, ~sheet)
        thick = np.where(sheet, np.minimum(P["sheet_max_m"], P["sheet_edge_m"] + P["sheet_growth_m_per_km"] * d_edge), 0.0)
    ice = np.zeros(g.n, np.int8)
    ice[cold & ~sheet] = ICE_GLACIER
    ice[sheet] = ICE_SHEET
    ice[~land & (winter < P["sea_ice_t_c"])] = ICE_SEA_SEASONAL
    ice[~land & (summer < P["sea_ice_t_c"])] = ICE_SEA_PERENNIAL
    return {"ice": ice, "ice_thickness_m": thick, "z_surface_m": z + thick}