worldgen: initial public history

31 files changed, 3877 insertions, 0 deletions

diff --git a/tests/__init__.py b/tests/__init__.py
new file mode 100644
index 0000000..ccdaff7
--- /dev/null
+++ b/tests/__init__.py
@@ -0,0 +1,10 @@
+import atexit
+import os
+import shutil
+import tempfile
+
+if "WORLDGEN_TEST_TMP" not in os.environ:
+    _root = tempfile.mkdtemp(prefix="worldgen-tests-")
+    os.environ["WORLDGEN_TEST_TMP"] = os.environ["TMPDIR"] = tempfile.tempdir = _root
+    _pid = os.getpid()
+    atexit.register(lambda: os.getpid() == _pid and shutil.rmtree(_root, True))
diff --git a/tests/data/climate_golden_r2.npz b/tests/data/climate_golden_r2.npz
new file mode 100644
index 0000000..1346781
--- /dev/null
+++ b/tests/data/climate_golden_r2.npz
diff --git a/tests/helpers.py b/tests/helpers.py
new file mode 100644
index 0000000..1bb865b
--- /dev/null
+++ b/tests/helpers.py
@@ -0,0 +1,39 @@
+import copy
+import functools
+from pathlib import Path
+
+ROOT = Path(__file__).resolve().parents[1]
+
+BASE_CFG = {
+    "planet": {"radius_km": 12742.0, "gravity_g": 1.05, "day_hours": 31.149, "year_days": 216,
+               "tilt_deg": 20.0, "sea_level_pressure_bar": 1.0, "scale_height_km": 8.0, "o2_fraction": 0.21},
+    "build": {"seed": 7, "res_dev": 2, "res_final": 2, "raster_width": 256, "preview_width": 128,
+              "land_fraction": 0.22},
+}
+
+
+@functools.lru_cache(maxsize=4)
+def small_grid(res: int = 2):
+    from mapgen.grid import build_grid
+    return build_grid(res, 12742.0)
+
+
+def make_ctx(res=2, data=None, cfg=None, tect=None, root=None):
+    from mapgen.pipeline import Ctx
+    c = copy.deepcopy(BASE_CFG)
+    for sec, vals in (cfg or {}).items():
+        c.setdefault(sec, {}).update(vals)
+    ctx = Ctx(root or ROOT, c, tect or {"plate": []}, res)
+    ctx.grid = small_grid(res)
+    ctx.data = dict(data or {})
+    return ctx
+
+
+def big_river(a) -> int:
+    """The cell of the biggest river that flows on over land (not into a lake at its own level, which cuts no
+    valley, nor off a coast, like a basin lake's spillway): the fixtures' river."""
+    import numpy as np
+    river, lake = np.asarray(a["river"]).astype(bool), np.asarray(a["lake"]).astype(bool)
+    recv = np.asarray(a["recv"])
+    wet = lake | np.asarray(a["ocean"]).astype(bool)
+    return int(np.argmax(np.asarray(a["discharge_km3_yr"]) * (river & ~wet & ~wet[recv])))
diff --git a/tests/test_check.py b/tests/test_check.py
new file mode 100644
index 0000000..a4ba179
--- /dev/null
+++ b/tests/test_check.py
@@ -0,0 +1,165 @@
+import unittest
+
+import numpy as np
+
+from mapgen import check as CK
+from tests.helpers import small_grid
+
+
+class CheckTest(unittest.TestCase):
+    def setUp(self):
+        self.g = small_grid(2)
+
+    def test_land_fraction(self):
+        a = np.ones(100)
+        z = np.where(np.arange(100) < 22, 100.0, -100.0)
+        self.assertIsNone(CK.check_land_fraction(z, a, 0.22))
+        self.assertIn("land fraction", CK.check_land_fraction(z, a, 0.40))
+
+    def test_hypsometry(self):
+        rng = np.random.default_rng(0)
+        z = np.concatenate([rng.normal(-4000, 600, 7000), rng.normal(500, 400, 3000)])
+        self.assertIsNone(CK.check_hypsometry(z, np.ones(len(z))))
+        self.assertIn("bimodal", CK.check_hypsometry(rng.normal(0, 3000, 10000), np.ones(10000)))
+
+    def test_max_elevation(self):
+        self.assertIsNone(CK.check_max_elevation(np.array([11000.0]), np.array([1.0])))
+        self.assertIsNone(CK.check_max_elevation(np.array([20000.0]), np.array([0.5])))
+        self.assertIn("12000", CK.check_max_elevation(np.array([13000.0]), np.array([1.0])))
+
+    def test_drainage_and_uphill(self):
+        recv = np.array([0, 0, 1, 3])
+        z = np.array([-10.0, 5.0, 8.0, 20.0])
+        endo = np.array([False, False, False, False])
+        self.assertIn("sink", CK.check_drainage(recv, z, endo))
+        endo[3] = True
+        self.assertIsNone(CK.check_drainage(recv, z, endo))
+        self.assertIsNone(CK.check_no_uphill(z, recv, z, endo))
+        self.assertIn("uphill", CK.check_no_uphill(np.array([-10.0, 9.0, 8.0, 20.0]), recv, z, endo))
+
+    def test_temperature_and_rain_bands(self):
+        g = self.g
+        self.assertIsNone(CK.check_temperature(g, 30 - 0.6 * np.abs(g.lat)))
+        self.assertIn("poleward", CK.check_temperature(g, 0.3 * np.abs(g.lat)))
+        wet = 2000 * np.exp(-(g.lat / 8) ** 2) + 900 * np.exp(-((np.abs(g.lat) - 40) / 8) ** 2) + 100
+        self.assertIsNone(CK.check_rain_bands(g, wet, 20.0))
+        self.assertIn("subtropic", CK.check_rain_bands(g, np.full(g.n, 800.0), 20.0))
+
+    def test_drainage_terminals_must_hold_water_or_salt(self):
+        recv = np.array([0, 0, 1, 3])
+        z = np.array([-10.0, 5.0, 8.0, 20.0])
+        endo = np.array([False, False, False, True])
+        dry = np.zeros(4, bool)
+        self.assertIn("lake", CK.check_drainage(recv, z, endo, terminal_ok=dry))
+        wet = dry.copy()
+        wet[3] = True
+        self.assertIsNone(CK.check_drainage(recv, z, endo, terminal_ok=wet))
+
+    def test_lake_surface_routing_is_exempt(self):
+        recv = np.array([0, 0, 1, 2])
+        zf = np.array([-10.0, 5.0, 5.02, 5.01])      # cell 3 → 2 is ε-uphill across a flat lake
+        z = np.array([-10.0, 5.0, 4.0, 4.5])
+        endo = np.zeros(4, bool)
+        lake = np.array([False, False, True, True])
+        self.assertIn("uphill", CK.check_no_uphill(zf, recv, z, endo))
+        self.assertIsNone(CK.check_no_uphill(zf, recv, z, endo, lake=lake))
+
+
+class RevisionCheckTest(unittest.TestCase):
+    P1 = {"name": "east-flank", "center": [-39.9, -18.0], "area_km2": 3.0e6, "elongation": 1.9, "azimuth_deg": 30.0,
+          "top_m": [1500.0, 3000.0]}
+    P2 = {"name": "plateau-02", "center": [-0.4, 176.6], "area_km2": 1.0e6, "top_m": [1000.0, 1500.0],
+          "islands": True}
+
+    def world(self):
+        from mapgen import plateaus as PL
+        g = small_grid(4)
+        ids = PL.cell_ids(g.xyz, [self.P1, self.P2], 7, g.radius_km)
+        z = PL.apply(g, np.full(g.n, -6000.0), [self.P1, self.P2], ids, 7)
+        return g, {"elevation_eroded_m": z, "ocean": z <= 0, "plateau_id": ids}
+
+    def test_plateaus_pass_and_fail(self):
+        g, d = self.world()
+        self.assertEqual(CK.check_plateaus(g, d, [self.P1, self.P2], 7), [])
+        z = d["elevation_eroded_m"].copy()
+        z[np.flatnonzero(d["plateau_id"] == 0)[0]] = -100.0
+        z[z > 0] = -700.0
+        bad = CK.check_plateaus(g, {**d, "elevation_eroded_m": z, "ocean": z <= 0}, [self.P1, self.P2], 7)
+        self.assertTrue(any("east-flank" in m and "−500" in m for m in bad), bad)
+        self.assertTrue(any("plateau-02" in m and "island" in m for m in bad), bad)
+
+    def test_zones_gradual_rule(self):
+        from mapgen import zones as ZN
+        g = small_grid(4)
+        m = ZN.contribution(g.xyz, {"name": "R1", "field": "o2", "center": [0.0, 0.0], "radius_km": 3200.0,
+                                    "v": 1.0}, 7, g.radius_km)
+        smooth = {"o2_fraction": 0.21 * (1 + 0.5 * m), "gravity_g": np.full(g.n, 1.05)}
+        self.assertEqual(CK.check_zones(g, smooth), [])
+        step = {"o2_fraction": np.where(g.lon > 0, 0.35, 0.21), "gravity_g": np.full(g.n, 1.05)}
+        self.assertTrue(any("o2_fraction" in m for m in CK.check_zones(g, step)))
+
+    def test_compare_skips_o2_fields_inside_o2_zones(self):
+        g = small_grid(3)
+        tect = {"plateau": [], "land_patch": [],
+                "zone": [{"name": "r", "field": "o2", "center": [0.0, 0.0], "radius_km": 2000.0, "v": 1.0}]}
+        o2 = CK.o2_areas(g.xyz, tect, g.radius_km)
+        self.assertTrue(o2[g.cell_index(0.0, 0.0)])
+        self.assertFalse(o2[g.cell_index(0.0, 90.0)])
+        ex = CK.edited_areas(g.xyz, tect, g.radius_km)
+        self.assertFalse(ex.any(), "O₂ zones don't change relief")
+        old = {"po2_bar": np.zeros(g.n), "T_mean": np.zeros(g.n)}
+        new = {"po2_bar": np.where(o2, 0.1, 0.0), "T_mean": np.where(o2, 1.0, 0.0)}
+        lines = CK.compare(old, new, ex, {k: o2 for k in CK.O2_FIELDS})
+        self.assertIn("po2_bar: 0.0% changed, max 0, p99 0", lines)
+        self.assertTrue(any(s.startswith("T_mean:") and not s.startswith("T_mean: 0.0%") for s in lines), lines)
+
+    def test_compare_command_needs_a_world_cells_file(self):
+        import importlib.util
+        import tempfile
+        from pathlib import Path
+        from mapgen.testing import built_world
+        spec = importlib.util.spec_from_file_location("mapgen_cli", Path(CK.__file__).parents[1] / "mapgen.py")
+        cli = importlib.util.module_from_spec(spec)
+        spec.loader.exec_module(cli)
+        old = Path(tempfile.mkdtemp()) / "other.npz"
+        np.savez(old, a=np.zeros(3))
+        self.assertEqual(cli.main(["compare", "--res", "2", "--old", str(old)], root=built_world()), 2)
+
+    def test_compare_ignores_edited_areas(self):
+        g = small_grid(3)
+        tect = {"plateau": [self.P1], "land_patch": [], "zone": []}
+        ex = CK.edited_areas(g.xyz, tect, g.radius_km)
+        self.assertTrue(ex[g.cell_index(-39.9, -18.0)])
+        self.assertFalse(ex[g.cell_index(40.0, 100.0)])
+        old = {"a": np.zeros(g.n), "b": np.zeros(g.n), "g_ids": g.ids}
+        new = {"a": np.where(ex, 5.0, 0.0), "b": np.where(ex, 0.0, 1.0), "g_ids": g.ids}
+        lines = CK.compare(old, new, ex)
+        self.assertIn("a: 0.0% changed, max 0, p99 0", lines)
+        self.assertTrue(any(line.startswith("b: 100.0% changed, max 1") for line in lines), lines)
+
+    def test_site_report_warns_about_land_nearby(self):
+        from mapgen import plateaus as PL
+        g, d = self.world()
+        land = np.abs(g.lat - (-39.9)) < 2.0
+        near = land & (np.abs(g.lon - (-18.0)) < 30.0) & (d["plateau_id"] < 0)
+        data = {**d, "ocean": ~near, "bnd_type": np.zeros(g.n, np.int8)}
+        msgs = PL.site_report(g, data, [self.P1, self.P2])
+        self.assertTrue(any("east-flank" in m and "from land" in m for m in msgs), msgs)
+
+    def test_compare_command(self):
+        import importlib.util
+        import shutil
+        import tempfile
+        from pathlib import Path
+        from mapgen.testing import built_world
+        spec = importlib.util.spec_from_file_location("mapgen_cli", Path(CK.__file__).parents[1] / "mapgen.py")
+        cli = importlib.util.module_from_spec(spec)
+        spec.loader.exec_module(cli)
+        old = Path(tempfile.mkdtemp()) / "cells.npz"
+        shutil.copy(built_world() / "out" / "r2" / "cells.npz", old)
+        self.assertEqual(cli.main(["compare", "--res", "2", "--old", str(old)], root=built_world()), 0)
+
+    def test_compare_equal_infinities_are_unchanged(self):
+        ex = np.zeros(3, bool)
+        lines = CK.compare({"d": np.array([np.inf, 1.0, 2.0])}, {"d": np.array([np.inf, 1.0, 3.0])}, ex)
+        self.assertIn("d: 33.3% changed, max 1, p99 0.98", lines)
diff --git a/tests/test_climate_ocean.py b/tests/test_climate_ocean.py
new file mode 100644
index 0000000..8132c6e
--- /dev/null
+++ b/tests/test_climate_ocean.py
@@ -0,0 +1,81 @@
+import unittest
+from pathlib import Path
+
+import numpy as np
+
+from mapgen import climate as CL
+from tests.helpers import make_ctx
+
+GOLDEN = Path(__file__).parent / "data" / "climate_golden_r2.npz"
+NEW = ("current", "current_speed", "sst", "upwelling", "productivity")
+
+
+def golden_ctx(cfg=None):
+    ctx = make_ctx(2, cfg=cfg)
+    g = ctx.grid
+    land = (np.abs(g.lat) < 40) & (np.abs(g.lon) < 50)
+    ctx.data["elevation_eroded_m"] = np.where(land, np.where(np.abs(g.lon) < 10, 2500.0, 300.0),
+                                              -4000.0).astype(np.float32)
+    return ctx, land
+
+
+class OceanDisabledTest(unittest.TestCase):
+    def test_disabled_reproduces_pre_ocean_climate(self):
+        ctx, land = golden_ctx({"ocean": {"enabled": False}})
+        out = CL.run(ctx)
+        ref = np.load(GOLDEN)
+        for k in ref.files:
+            np.testing.assert_array_equal(out[k], ref[k], err_msg=k)
+        self.assertTrue(np.all(out["current"] == 0) and np.all(out["productivity"] == 0))
+        np.testing.assert_array_equal(out["sst"], out["T_mean"].astype(np.float32))
+
+
+class OceanCoupledTest(unittest.TestCase):
+    @classmethod
+    def setUpClass(cls):
+        cls.ctx, cls.land = golden_ctx()
+        cls.out = CL.run(cls.ctx)
+
+    def test_new_fields(self):
+        g, out = self.ctx.grid, self.out
+        for k in NEW:
+            self.assertIn(k, out)
+            self.assertEqual(out[k].dtype, np.float32, k)
+            self.assertTrue(np.all(np.isfinite(out[k])), k)
+        self.assertEqual(out["current"].shape, (g.n, 3))
+        self.assertTrue(np.all(out["current"][self.land] == 0))
+        np.testing.assert_allclose(out["current_speed"], np.linalg.norm(out["current"], axis=1), rtol=1e-5)
+        sea = ~self.land
+        self.assertTrue(np.all((out["sst"][sea] >= -1.8) & (out["sst"][sea] < 40)))            # sea water freezes at −1.8 °C
+
+    def test_currents_change_the_climate(self):
+        ref = np.load(GOLDEN)
+        self.assertGreater(np.abs(self.out["T_mean"] - ref["T_mean"]).max(), 0.5)   # the old gyre rule is gone
+
+    def test_winds_come_from_current_free_temperatures(self):
+        ctx, _ = golden_ctx({"ocean": {"relax_days": 30.0}})                      # other SST, same winds
+        np.testing.assert_array_equal(CL.run(ctx)["wind_jun"], self.out["wind_jun"])
+
+    def test_unknown_ocean_key_is_an_error(self):
+        from mapgen.config import ConfigError
+        ctx, _ = golden_ctx({"ocean": {"frictoin_days": 3.0}})
+        with self.assertRaises(ConfigError):
+            CL.run(ctx)
+
+    def test_all_land_world(self):
+        ctx = make_ctx(2)
+        ctx.data["elevation_eroded_m"] = np.full(ctx.grid.n, 500.0, np.float32)
+        out = CL.run(ctx)
+        self.assertTrue(np.all(out["current"] == 0))
+
+    def test_locked_world_has_still_ocean(self):
+        ctx, land = golden_ctx({"climate": {"lock": True}})
+        out = CL.run(ctx)
+        self.assertTrue(np.all(out["current"] == 0) and np.all(out["upwelling"] == 0))
+        np.testing.assert_array_equal(out["sst"], out["T_mean"].astype(np.float32))
+
+    def test_upwelling_and_productivity_filled(self):
+        sea = ~self.land
+        self.assertGreater(np.abs(self.out["upwelling"][sea]).max(), 1.0)
+        self.assertGreater(self.out["productivity"][sea].max(), 0.05)
+        self.assertTrue(np.all(self.out["productivity"][self.land] == 0))
diff --git a/tests/test_climate_rain.py b/tests/test_climate_rain.py
new file mode 100644
index 0000000..38afed5
--- /dev/null
+++ b/tests/test_climate_rain.py
@@ -0,0 +1,117 @@
+import unittest
+
+import numpy as np
+
+from mapgen import climate as CL
+from mapgen.sphere import east_north
+from tests.helpers import make_ctx, small_grid
+
+
+class WindTest(unittest.TestCase):
+    def test_bands(self):
+        g = small_grid(3)
+        w = CL.band_winds(g, 0.0, CL.DEFAULTS)
+        e, n = east_north(g.xyz)
+        u = np.sum(w * e, axis=1)
+        v = np.sum(w * n, axis=1)
+        self.assertLess(u[np.abs(g.lat - 10) < 2].mean(), -3)     # trades from the east
+        self.assertGreater(u[np.abs(g.lat - 40) < 2].mean(), 3)   # westerlies
+        self.assertLess(v[np.abs(g.lat - 10) < 2].mean(), 0)      # NH trades flow equatorward
+        self.assertGreater(v[np.abs(g.lat + 10) < 2].mean(), 0)   # SH trades flow equatorward
+
+    def test_monsoon_reverses_onshore_flow(self):
+        g = small_grid(3)
+        land = (g.lat > 10) & (g.lat < 40) & (np.abs(g.lon) < 50)
+        z = np.where(land, 200.0, -4000.0)
+        T, _ = CL.temperatures(g, z, land, CL.DEFAULTS, 20.0)
+        _, n = east_north(g.xyz)
+        south = (g.lat > 0) & (g.lat < 8) & (np.abs(g.lon) < 40)
+        vj = np.sum(CL.monsoon_winds(g, T["jun"], land, CL.DEFAULTS) * n, axis=1)[south].mean()
+        vd = np.sum(CL.monsoon_winds(g, T["dec"], land, CL.DEFAULTS) * n, axis=1)[south].mean()
+        self.assertGreater(vj, 0.5)
+        self.assertLess(vd, -0.5)
+
+
+class RainTest(unittest.TestCase):
+    def test_itcz_wetter_than_subtropics_on_aquaplanet(self):
+        g = small_grid(3)
+        z = np.full(g.n, -4000.0)
+        T, _ = CL.temperatures(g, z, z > 0, CL.DEFAULTS, 20.0)
+        w = CL.band_winds(g, 0.0, CL.DEFAULTS)
+        p = CL.precipitation(g, w, T["eq"], z, z > 0, 0.0, CL.DEFAULTS)
+        self.assertGreater(p[np.abs(g.lat) < 5].mean(), 1.5 * p[np.abs(np.abs(g.lat) - 21) < 3].mean())
+
+    def test_windward_wetter_than_lee(self):
+        g = small_grid(3)
+        land = (g.lat > 30) & (g.lat < 50) & (np.abs(g.lon) < 60)
+        z = np.where(land, 200.0, -4000.0)
+        z = np.where(land & (np.abs(g.lon) < 5), 3000.0, z)
+        T, _ = CL.temperatures(g, z, land, CL.DEFAULTS, 20.0)
+        w = CL.band_winds(g, 0.0, CL.DEFAULTS)
+        p = CL.precipitation(g, w, T["eq"], z, land, 0.0, CL.DEFAULTS)
+        band = (g.lat > 35) & (g.lat < 45)
+        west = p[band & (g.lon > -12) & (g.lon < -3)].mean()
+        east = p[band & (g.lon > 3) & (g.lon < 12)].mean()
+        self.assertGreater(west, 1.3 * east)
+
+    def test_stage_outputs(self):
+        ctx = make_ctx(2)
+        g = ctx.grid
+        land = (np.abs(g.lat) < 40) & (np.abs(g.lon) < 50)
+        ctx.data["elevation_eroded_m"] = np.where(land, 300.0, -4000.0).astype(np.float32)
+        out = CL.run(ctx)
+        for k in ("T_jun", "T_dec", "T_eq", "T_mean", "T_range", "T_min", "P_jun", "P_dec", "P_eq", "P_ann",
+                  "wind_jun", "biotemp", "PET", "dist_ocean_km"):
+            self.assertIn(k, out)
+            self.assertTrue(np.all(np.isfinite(out[k])), k)
+        mean = np.sum(out["P_ann"] * g.area_km2) / g.area_km2.sum()
+        self.assertAlmostEqual(mean, 1000.0, delta=1.0)
+        self.assertTrue(np.all(out["P_ann"] >= 0))
+
+
+class StormTrackTest(unittest.TestCase):
+    def test_midlatitudes_wetter_than_subtropics(self):
+        g = small_grid(3)
+        z = np.full(g.n, -4000.0)
+        T, _ = CL.temperatures(g, z, z > 0, CL.DEFAULTS, 20.0)
+        p = CL.precipitation(g, CL.band_winds(g, 0.0, CL.DEFAULTS), T["eq"], z, z > 0, 0.0, CL.DEFAULTS)
+        a = np.abs(g.lat)
+        h = CL.DEFAULTS["hadley_edge_deg"]        # storm track ≈ h+12..h+25, dry belt ≈ h−2..h+8
+        self.assertGreater(p[(a > h + 12) & (a < h + 25)].mean(), 1.2 * p[(a > h - 2) & (a < h + 8)].mean())
+
+
+class LandMoistureTest(unittest.TestCase):
+    def test_equatorial_land_stays_wet(self):
+        g = small_grid(3)
+        land = (np.abs(g.lat) < 12) & (np.abs(g.lon) < 40)
+        z = np.where(land, 300.0, -4000.0)
+        T, _ = CL.temperatures(g, z, land, CL.DEFAULTS, 20.0)
+        p = CL.precipitation(g, CL.band_winds(g, 0.0, CL.DEFAULTS), T["eq"], z, land, 0.0, CL.DEFAULTS)
+        eq = np.abs(g.lat) < 8
+        self.assertGreater(p[land & eq].mean(), 0.5 * p[~land & eq].mean())
+
+
+class OceanMaskUseTest(unittest.TestCase):
+    def test_inland_basin_is_land_for_climate(self):
+        ctx = make_ctx(3)
+        g = ctx.grid
+        land = (np.abs(g.lat) < 40) & (np.abs(g.lon) < 60)
+        z = np.where(land, 300.0, -4000.0)
+        basin = (np.abs(g.lat) < 8) & (np.abs(g.lon) < 8)
+        z[basin] = -30.0
+        ctx.data.update({"elevation_eroded_m": z.astype(np.float32), "ocean": ~land})
+        out = CL.run(ctx)
+        self.assertTrue(np.all(out["dist_ocean_km"][basin] > 1000))
+
+
+class InteriorRainTest(unittest.TestCase):
+    def test_continental_interior_keeps_a_third_of_coastal_rain(self):
+        from mapgen.graph import distance_to
+        g = small_grid(3)
+        land = (g.lat > 10) & (g.lat < 50) & (np.abs(g.lon) < 70)
+        z = np.where(land, 300.0, -4000.0)
+        T, _ = CL.temperatures(g, z, land, CL.DEFAULTS, 20.0)
+        p = CL.precipitation(g, CL.band_winds(g, 0.0, CL.DEFAULTS), T["eq"], z, land, 0.0, CL.DEFAULTS)
+        d = distance_to(g, ~land)
+        coast, interior = p[land & (d < 500)].mean(), p[land & (d > 1500)].mean()
+        self.assertGreater(interior, 0.33 * coast)          # was ≈0.25 before the desert retune
diff --git a/tests/test_climate_temp.py b/tests/test_climate_temp.py
new file mode 100644
index 0000000..2a61fc3
--- /dev/null
+++ b/tests/test_climate_temp.py
@@ -0,0 +1,93 @@
+import unittest
+
+import numpy as np
+
+from mapgen import climate as CL
+from tests.helpers import small_grid
+
+
+class InsolationTest(unittest.TestCase):
+    def test_insolation_finite_at_poles(self):
+        q = CL.insolation(np.array([90.0, -90.0, 89.999]), 20.0)
+        self.assertTrue(np.all(np.isfinite(q)))
+        self.assertAlmostEqual(q[1], 0.0, places=6)          # polar night
+        self.assertGreater(q[0], 400)                          # polar day
+
+    def test_symmetry(self):
+        lat = np.linspace(-80, 80, 17)
+        np.testing.assert_allclose(CL.insolation(lat, 20.0), CL.insolation(-lat, -20.0), rtol=1e-12)
+
+
+class TemperatureTest(unittest.TestCase):
+    def setUp(self):
+        self.g = small_grid(3)
+        self.P = dict(CL.DEFAULTS)
+
+    def test_aquaplanet_zonal_means_fall_poleward(self):
+        g = self.g
+        z = np.full(g.n, -4000.0)
+        T, _ = CL.temperatures(g, z, z > 0, self.P, 20.0)
+        tm = (T["jun"] + T["dec"] + 2 * T["eq"]) / 4
+        bands = [tm[(np.abs(g.lat) >= a) & (np.abs(g.lat) < a + 10)].mean() for a in range(0, 90, 10)]
+        self.assertTrue(all(b1 > b2 for b1, b2 in zip(bands, bands[1:])), bands)
+        self.assertTrue(20 < bands[0] < 32 and bands[-1] < -5, bands)
+
+    def test_lapse_rate(self):
+        g = self.g
+        land = (np.abs(g.lat) < 30) & (np.abs(g.lon) < 60)
+        flat = np.where(land, 10.0, -4000.0)
+        high = np.where(land & (np.abs(g.lon) < 20), 3000.0, flat)
+        T0, _ = CL.temperatures(g, flat, land, self.P, 20.0)
+        T1, _ = CL.temperatures(g, high, land, self.P, 20.0)
+        i = g.cell_index(0.0, 0.0)
+        self.assertAlmostEqual(T0["eq"][i] - T1["eq"][i], 6.5 * 2.99, delta=0.5)
+
+    def test_cold_west_coast_warm_east_coast(self):
+        g = self.g
+        land = (g.lat > 15) & (g.lat < 45) & (np.abs(g.lon) < 30)
+        a = CL.current_anomaly(g, land, self.P)
+        self.assertLess(a[g.cell_index(30.0, -33.0)], -0.5)
+        self.assertGreater(a[g.cell_index(30.0, 33.0)], 0.5)
+
+    def test_continental_interior_has_bigger_seasons(self):
+        g = self.g
+        land = (g.lat > 20) & (g.lat < 70) & (np.abs(g.lon) < 90)
+        z = np.where(land, 200.0, -4000.0)
+        T, _ = CL.temperatures(g, z, land, self.P, 20.0)
+        rng = np.abs(T["jun"] - T["dec"])
+        self.assertGreater(rng[g.cell_index(50.0, 0.0)], rng[g.cell_index(50.0, 150.0)] + 10)
+
+    def test_biotemperature(self):
+        bio = CL.biotemperature(np.array([10.0, -5.0, 35.0, 0.0]), np.array([0.0, 0.0, 0.0, 20.0]))
+        np.testing.assert_allclose(bio[:3], [10.0, 0.0, 30.0])
+        self.assertTrue(2.0 < bio[3] < 4.0)
+
+
+class CalibrationTest(unittest.TestCase):
+    def test_earthlike_midlatitudes_and_sea_ice_edge(self):
+        g = small_grid(3)
+        z = np.full(g.n, -4000.0)
+        T, _ = CL.temperatures(g, z, z > 0, CL.DEFAULTS, 20.0)
+        tm = (T["jun"] + T["dec"] + 2 * T["eq"]) / 4
+        band = lambda a: np.abs(np.abs(g.lat) - a) < 2
+        self.assertTrue(10 < tm[band(45)].mean() < 18, tm[band(45)].mean())
+        self.assertTrue(-1 < tm[band(60)].mean() < 8, tm[band(60)].mean())
+        winter = np.where(g.lat >= 0, T["dec"], T["jun"])
+        self.assertGreater(winter[band(48)].mean(), -1.8)     # no sea ice at 48°
+        self.assertLess(winter[band(66)].mean(), -1.8)        # seasonal sea ice by 66°
+
+    def test_subtropical_land_heats_in_summer(self):
+        g = small_grid(3)
+        land = (g.lat > 10) & (g.lat < 40) & (np.abs(g.lon) < 50)
+        z = np.where(land, 200.0, -4000.0)
+        T, _ = CL.temperatures(g, z, land, CL.DEFAULTS, 20.0)
+        anom = T["jun"] - CL.zonal_mean(g, T["jun"])
+        self.assertGreater(anom[land].mean(), 2.5)
+        self.assertLess(T["jun"][land].max(), 42.0)
+
+    def test_ocean_seasonal_range_earthlike(self):
+        g = small_grid(3)
+        z = np.full(g.n, -4000.0)
+        T, _ = CL.temperatures(g, z, z > 0, CL.DEFAULTS, 20.0)
+        rng = np.abs(T["jun"] - T["dec"])[np.abs(np.abs(g.lat) - 50) < 2].mean()
+        self.assertTrue(4 < rng < 10, rng)
diff --git a/tests/test_config.py b/tests/test_config.py
new file mode 100644
index 0000000..5aafbef
--- /dev/null
+++ b/tests/test_config.py
@@ -0,0 +1,182 @@
+import shutil
+import tempfile
+import unittest
+from pathlib import Path
+
+from mapgen import config as C
+
+from mapgen.testing import FIXTURE_TOML
+
+TECT = """
+[[plate]]
+id = "a"
+seed = [0.0, 0.0]
+kind = "continental"
+motion = [90.0, 3.0]
+[[plate]]
+id = "b"
+seed = [0.0, 90.0]
+kind = "oceanic"
+motion = [270.0, 3.0]
+"""
+
+
+class ConfigTest(unittest.TestCase):
+    def setUp(self):
+        self.tmp = Path(tempfile.mkdtemp())
+        (self.tmp / "config").mkdir()
+        shutil.copy(FIXTURE_TOML, self.tmp / "config" / "world.toml")
+        (self.tmp / "config" / "tectonics.toml").write_text(TECT)
+
+    def tearDown(self):
+        shutil.rmtree(self.tmp)
+
+    def test_loads_repo_world(self):
+        cfg, tect = C.load(self.tmp)
+        self.assertEqual(cfg["planet"]["radius_km"], 12742.0)
+        self.assertEqual(len(tect["plate"]), 2)
+
+    def test_out_of_range_is_error(self):
+        p = self.tmp / "config" / "world.toml"
+        p.write_text(p.read_text().replace("tilt_deg = 20.0", "tilt_deg = 120.0"))
+        with self.assertRaisesRegex(C.ConfigError, "tilt_deg"):
+            C.load(self.tmp)
+
+    def test_missing_key_is_error(self):
+        p = self.tmp / "config" / "world.toml"
+        p.write_text(p.read_text().replace("seed = 1296\n", ""))
+        with self.assertRaisesRegex(C.ConfigError, "seed"):
+            C.load(self.tmp)
+
+    def test_params_merge_and_typo_guard(self):
+        cfg = {"erosion": {"k": 0.5}}
+        self.assertEqual(C.params(cfg, "erosion", {"k": 0.1, "m": 0.5}), {"k": 0.5, "m": 0.5})
+        with self.assertRaisesRegex(C.ConfigError, "kk"):
+            C.params({"erosion": {"kk": 1}}, "erosion", {"k": 0.1})
+
+    def test_tectonics_validation(self):
+        (self.tmp / "config" / "tectonics.toml").write_text(TECT.replace('kind = "oceanic"', 'kind = "lava"'))
+        with self.assertRaisesRegex(C.ConfigError, "kind"):
+            C.load(self.tmp)
+        (self.tmp / "config" / "tectonics.toml").write_text(TECT.replace('id = "b"', 'id = "a"'))
+        with self.assertRaisesRegex(C.ConfigError, "duplicate"):
+            C.load(self.tmp)
+
+
+REV = """
+[[plateau]]
+name = "p1"
+center = [-40.0, -18.0]
+area_km2 = 3.0e6
+elongation = 1.9
+azimuth_deg = 30.0
+top_m = [1500.0, 3000.0]
+[[plateau]]
+name = "p2"
+center = [0.0, 176.0]
+area_km2 = 1.0e6
+top_m = [1000.0, 1500.0]
+islands = true
+[[land_patch]]
+name = "fill"
+center = [-10.0, 102.0]
+radius_km = 2850.0
+edge_noise = 0.4
+[[zone]]
+name = "R1"
+field = "o2"
+center = [-5.6, -118.8]
+radius_km = 3200.0
+v = 1.0
+"""
+ERAS = """
+[[event]]
+name = "cut"
+kind = "disintegrate"
+center = [-10.0, 102.0]
+radius_km = 3000.0
+depth_m = 3000.0
+[[event]]
+name = "aura"
+kind = "zone"
+shape = "landmass"
+seed = [-12.7, 121.5]
+reach_km = [500.0, 1000.0]
+edge_km = 5.0
+fields = { gravity_g = 0.35, pressure_bar = 2.0, o2_fraction = 0.35, fire_reactivity = 0.5 }
+[eras]
+order = ["before", "after"]
+default = "after"
+[eras.before]
+label = "Before"
+events = []
+[eras.after]
+label = "After"
+events = ["cut", "aura"]
+"""
+
+
+class RevisionConfigTest(unittest.TestCase):
+    def setUp(self):
+        self.tmp = Path(tempfile.mkdtemp())
+        (self.tmp / "config").mkdir()
+        shutil.copy(FIXTURE_TOML, self.tmp / "config" / "world.toml")
+        self._set()
+
+    def tearDown(self):
+        shutil.rmtree(self.tmp)
+
+    def _set(self, rev=REV, eras=ERAS):
+        (self.tmp / "config" / "tectonics.toml").write_text(TECT + rev)
+        (self.tmp / "config" / "eras.toml").write_text(eras)
+
+    def test_valid_revision_config_loads(self):
+        _, t = C.load(self.tmp)
+        self.assertEqual([p["name"] for p in t["plateau"]], ["p1", "p2"])
+        self.assertEqual([e["name"] for e in C.era_events(t, "after")], ["cut", "aura"])
+        self.assertEqual(C.era_events(t, "before"), [])
+
+    def test_era_events_are_cumulative(self):
+        eras = ERAS.replace('order = ["before", "after"]', 'order = ["before", "mid", "after"]').replace(
+            '[eras.after]\nlabel = "After"\nevents = ["cut", "aura"]',
+            '[eras.mid]\nlabel = "Mid"\nevents = ["aura"]\n[eras.after]\nlabel = "After"\nevents = ["cut"]')
+        self._set(REV, eras)
+        _, t = C.load(self.tmp)
+        self.assertEqual([e["name"] for e in C.era_events(t, "mid")], ["aura"])
+        self.assertEqual([e["name"] for e in C.era_events(t, "after")], ["aura", "cut"])
+        with self.assertRaises(C.ConfigError):
+            C.era_events(t, "never")
+
+    def test_errors_name_the_problem(self):
+        dup = '[[plateau]]\nname = "p1"\ncenter = [60.0, 60.0]\narea_km2 = 1e6\ntop_m = [1000.0, 2000.0]\n'
+        cases = [
+            (REV.replace("top_m = [1500.0, 3000.0]", "top_m = [3000.0, 1500.0]"), ERAS, "top_m"),
+            (REV.replace("center = [0.0, 176.0]", "center = [-41.0, -10.0]"), ERAS, "apart"),
+            (REV + dup, ERAS, "duplicate"),
+            (REV.replace('field = "o2"', 'field = "heat"'), ERAS, "field"),
+            (REV, ERAS.replace("fire_reactivity = 0.5", "sparkle = 1.0"), "sparkle"),
+            (REV, ERAS.replace('default = "after"', 'default = "later"'), "default"),
+            (REV, ERAS.replace('events = ["cut", "aura"]', 'events = ["cut", "nope"]'), "nope"),
+            (REV, ERAS.replace('kind = "disintegrate"', 'kind = "flood"'), "kind"),
+            (REV, ERAS.replace("edge_km = 5.0\n", ""), "edge_km"),
+        ]
+        for rev, eras, needle in cases:
+            with self.subTest(needle=needle):
+                self._set(rev, eras)
+                with self.assertRaises(C.ConfigError) as cm:
+                    C.load(self.tmp)
+                self.assertIn(needle, str(cm.exception))
+
+    def test_events_belong_in_eras_toml(self):
+        (self.tmp / "config" / "tectonics.toml").write_text(TECT + REV + ERAS)
+        with self.assertRaises(C.ConfigError) as cm:
+            C.load(self.tmp)
+        self.assertIn("eras.toml", str(cm.exception))
+
+    def test_eras_toml_is_not_part_of_the_base_inputs_key(self):
+        from mapgen import pipeline as P
+        k = P.inputs_key(self.tmp, 1)
+        self._set(REV, ERAS.replace("fire_reactivity = 0.5", "fire_reactivity = 0.6"))
+        self.assertEqual(P.inputs_key(self.tmp, 1), k, "tuning an event must not rebuild the base world")
+        self._set(REV.replace("v = 1.0", "v = 0.9"), ERAS)
+        self.assertNotEqual(P.inputs_key(self.tmp, 1), k)
diff --git a/tests/test_crust.py b/tests/test_crust.py
new file mode 100644
index 0000000..1bfca60
--- /dev/null
+++ b/tests/test_crust.py
@@ -0,0 +1,114 @@
+import unittest
+
+import numpy as np
+
+from mapgen import crust as CR
+from mapgen.plates import CONV, DIV
+from tests.helpers import make_ctx
+
+TECT = {"plate": [],
+        "lip": [{"name": "l", "center": [10.0, -30.0], "radius_km": 1500.0}],
+        "volcano": [{"name": "v", "center": [-10.0, 20.0], "radius_km": 700.0, "height_m": 7000.0,
+                     "scar_azimuths": [90.0]}]}
+
+
+def ctx3():
+    ctx = make_ctx(3, tect=TECT, cfg={"crust": {"edge_noise": 0.0}})
+    g = ctx.grid
+    land = ((np.abs(g.lat) < 35) & (np.abs(g.lon) < 60)).astype(float)
+    bt = np.zeros(g.n, np.int8)
+    bt[np.abs(g.lon - 120) < 1.0] = DIV   # a mid-ocean ridge along lon 120
+    bt[np.abs(g.lon - 0) < 1.0] = CONV
+    ctx.data.update({"sk_land": land, "m_land_hint": np.zeros(g.n), "bnd_type": bt})
+    return ctx
+
+
+class CrustTest(unittest.TestCase):
+    def test_classes(self):
+        ctx = ctx3()
+        out = CR.run(ctx)
+        g = ctx.grid
+        age = out["age_class"]
+        self.assertTrue(out["continental"][g.cell_index(20.0, -50.0)])
+        self.assertFalse(out["continental"][g.cell_index(0.0, 150.0)])
+        self.assertEqual(age[g.cell_index(0.0, 150.0)], CR.OCEANIC)
+        self.assertEqual(age[g.cell_index(10.0, -30.0)], CR.LIP)
+        self.assertEqual(age[g.cell_index(-10.0, 20.0)], CR.VOLCANO)
+        self.assertEqual(age[g.cell_index(-10.0, 24.0)], CR.SCAR)      # east sector, ~900 km out
+        self.assertEqual(age[g.cell_index(20.0, 3.0)], CR.POST_OROGEN)
+
+    def test_ocean_age_grows_from_ridge(self):
+        ctx = ctx3()
+        out = CR.run(ctx)
+        g = ctx.grid
+        ocean = ~out["continental"] & (np.abs(g.lat) < 20) & (np.abs(g.lon - 120) < 20)
+        r = np.corrcoef(np.abs(g.lon[ocean] - 120), out["ocean_age_myr"][ocean])[0, 1]
+        self.assertGreater(r, 0.95)
+        self.assertTrue(np.all(out["ocean_age_myr"][out["continental"]] == 0))
+
+
+class ShelfTest(unittest.TestCase):
+    def test_continental_shelf_is_resolution_independent(self):
+        from mapgen.graph import distance_to
+        for res in (3, 4):
+            ctx = make_ctx(res, tect={"plate": []}, cfg={"crust": {"edge_noise": 0.0}})
+            g = ctx.grid
+            land = ((np.abs(g.lat) < 30) & (np.abs(g.lon) < 40)).astype(float)
+            ctx.data.update({"sk_land": land, "m_land_hint": np.zeros(g.n), "bnd_type": np.zeros(g.n, np.int8)})
+            cont = CR.run(ctx)["continental"]
+            d = distance_to(g, land > 0.5)
+            self.assertTrue(np.all(cont[land > 0.5]), res)
+            self.assertGreater(cont[(d > 0) & (d < 300)].mean(), 0.8, res)     # shelf
+            self.assertEqual(cont[d > 700].sum(), 0, res)                      # but not far out
+
+
+class FractalMarginTest(unittest.TestCase):
+    @staticmethod
+    def _edge(noise):
+        from mapgen.graph import components
+        ctx = make_ctx(4, tect={"plate": []}, cfg={"crust": {"edge_noise": noise}})
+        g = ctx.grid
+        land = ((np.abs(g.lat) < 30) & (np.abs(g.lon) < 60)).astype(float)
+        ctx.data.update({"sk_land": land, "m_land_hint": np.zeros(g.n), "bnd_type": np.zeros(g.n, np.int8)})
+        cont = CR.run(ctx)["continental"]
+        lab = components(g, cont)
+        sizes = np.bincount(lab[cont])
+        return int(np.sum(cont[g.src] != cont[g.dst])), int(np.sum((sizes > 0) & (sizes < 0.05 * sizes.max())))
+
+    def test_margins_are_fractal_with_fragments(self):
+        smooth, _ = self._edge(0.0)
+        edges, fragments = self._edge(CR.DEFAULTS["edge_noise"])
+        self.assertGreater(edges, 1.4 * smooth)
+        self.assertGreaterEqual(fragments, 3)
+
+
+class RevisionCrustTest(unittest.TestCase):
+    def test_no_patch_no_plateau_keeps_crust(self):
+        out = CR.run(ctx3())
+        np.testing.assert_array_equal(out["continental"], out["continental_base"])
+        self.assertTrue(np.all(out["plateau_id"] == -1))
+
+    def test_land_patch_adds_continental_crust_only_around_its_disc(self):
+        ctx = ctx3()
+        ctx.tect = {**TECT, "land_patch": [{"name": "fill", "center": [0.0, 150.0], "radius_km": 2000.0,
+                                            "edge_noise": 0.4}]}
+        out = CR.run(ctx)
+        g = ctx.grid
+        i = g.cell_index(0.0, 150.0)
+        self.assertTrue(out["continental"][i])
+        self.assertFalse(out["continental_base"][i])
+        far = CR.center_dist(g, [0.0, 150.0]) > 4500.0
+        np.testing.assert_array_equal(out["continental"][far], out["continental_base"][far])
+        self.assertEqual(out["age_class"][i] == CR.OCEANIC, False)
+
+    def test_plateau_cells_become_continental_crust_with_ids(self):
+        ctx = ctx3()
+        ctx.tect = {**TECT, "plateau": [{"name": "p", "center": [0.0, 150.0], "area_km2": 3.0e6,
+                                         "top_m": [1500.0, 3000.0]}]}
+        out = CR.run(ctx)
+        i = ctx.grid.cell_index(0.0, 150.0)
+        self.assertEqual(out["plateau_id"][i], 0)
+        self.assertTrue(out["continental"][i])
+        self.assertFalse(out["continental_base"][i])
+        self.assertNotEqual(out["age_class"][i], CR.OCEANIC)
+        self.assertGreater(out["ocean_age_myr"][i], 0.0, "the plateau keeps the age of the sea floor around it")
diff --git a/tests/test_elevation.py b/tests/test_elevation.py
new file mode 100644
index 0000000..3a96460
--- /dev/null
+++ b/tests/test_elevation.py
@@ -0,0 +1,207 @@
+import unittest
+
+import numpy as np
+
+from mapgen import crust, elevation as EL, plates
+from tests.helpers import make_ctx
+
+
+def scenario(tect, land_fn, land_fraction, res=3, gravity=None, lock=None):
+    ctx = make_ctx(res, tect=tect, cfg={"build": {"land_fraction": land_fraction}})
+    g = ctx.grid
+    z0 = np.zeros(g.n)
+    ctx.data.update({"sk_land": land_fn(g).astype(float), "m_land_hint": z0, "sk_mountains": z0,
+                     "m_mountain_hint": z0, "m_gravity_zones": z0 if gravity is None else gravity(g),
+                     "m_lock": z0 if lock is None else lock(g)})
+    ctx.data.update(plates.run(ctx))
+    ctx.data.update(crust.run(ctx))
+    return ctx
+
+
+SUBDUCT = {"plate": [
+    {"id": "c", "seed": [0.0, -30.0], "kind": "continental", "motion": [90.0, 3.0]},
+    {"id": "o", "seed": [0.0, 10.0], "kind": "oceanic", "motion": [270.0, 5.0]}]}
+COLLIDE = {"plate": [
+    {"id": "a", "seed": [0.0, -30.0], "kind": "continental", "motion": [90.0, 4.0]},
+    {"id": "b", "seed": [0.0, 30.0], "kind": "continental", "motion": [270.0, 4.0]}]}
+
+
+class ElevationTest(unittest.TestCase):
+    def test_sea_level_hits_target(self):
+        rng = np.random.default_rng(1)
+        z = rng.normal(size=5000) * 1000
+        area = rng.uniform(1, 2, size=5000)
+        zs = EL.solve_sea_level(z, area, 0.3)
+        self.assertAlmostEqual(area[zs > 0].sum() / area.sum(), 0.3, delta=0.005)
+
+    def test_subduction_trench_and_coastal_range(self):
+        ctx = scenario(SUBDUCT, lambda g: (np.abs(g.lat) < 40) & (g.lon > -70) & (g.lon < 0), 0.12)
+        out = EL.run(ctx)
+        z, g = out["elevation_m"], ctx.grid
+        band = np.abs(g.lat) < 30
+        ocean_side = band & (ctx.data["plate"] == 1) & (out["d_sub_km"] < 300)
+        cont_side = band & (ctx.data["plate"] == 0) & (out["d_over_km"] < 600)
+        self.assertLess(z[ocean_side].min(), -7000)
+        self.assertGreater(z[cont_side].max(), 2500)
+
+    def test_collision_builds_high_range(self):
+        ctx = scenario(COLLIDE, lambda g: (np.abs(g.lat) < 30) & (np.abs(g.lon) < 60), 0.15)
+        out = EL.run(ctx)
+        near = out["d_coll_km"] < 300
+        self.assertGreater(out["elevation_m"][near].max(), 6000)
+        self.assertLessEqual(out["elevation_m"].max(), 12000)
+
+    def test_low_gravity_zone_raises_relief(self):
+        land = lambda g: (np.abs(g.lat) < 30) & (np.abs(g.lon) < 60)
+        zone = lambda g: np.where((np.abs(g.lat) < 20) & (g.lon > 10) & (g.lon < 50), -1.0, 0.0)
+        base = EL.run(scenario(COLLIDE, land, 0.15))["elevation_m"]
+        ctx = scenario(COLLIDE, land, 0.15, gravity=zone)
+        low = EL.run(ctx)["elevation_m"]
+        inside = zone(ctx.grid) < 0
+        self.assertGreater(low[inside & (base > 0)].mean(), base[inside & (base > 0)].mean() * 1.5)
+
+    def test_lock_forces_coast(self):
+        land = lambda g: (np.abs(g.lat) < 30) & (np.abs(g.lon) < 60)
+        lock = lambda g: np.ones(g.n)
+        ctx = scenario(COLLIDE, land, 0.15, lock=lock)
+        z = EL.run(ctx)["elevation_m"]
+        self.assertTrue(np.all((z > 0) == (ctx.data["sk_land"] > 0.5)))
+
+
+class SeaLevelRegressionTest(unittest.TestCase):
+    LAND = staticmethod(lambda g: (g.lat > -25) & (g.lat < 45) & (np.abs(g.lon) < 62))   # ≈ 0.195 of the sphere
+
+    def test_continents_not_lifted(self):
+        ctx = scenario(COLLIDE, self.LAND, 0.19)
+        out = EL.run(ctx)
+        z, g = out["elevation_m"], ctx.grid
+        self.assertLess(np.median(z[self.LAND(g) & (out["d_coll_km"] > 1500)]), 1500)
+        self.assertLess(np.median(z[~ctx.data["continental"]]), -3000)
+
+    def test_target_above_continental_area_errors(self):
+        from mapgen.pipeline import StageError
+        with self.assertRaisesRegex(StageError, "land_fraction"):
+            EL.run(scenario(COLLIDE, self.LAND, 0.30))
+
+
+class LockThresholdTest(unittest.TestCase):
+    def test_near_neutral_lock_does_nothing(self):
+        land = lambda g: (np.abs(g.lat) < 30) & (np.abs(g.lon) < 60)
+        base = EL.run(scenario(COLLIDE, land, 0.15))["elevation_m"]
+        tiny = EL.run(scenario(COLLIDE, land, 0.15, lock=lambda g: np.full(g.n, 1 / 255)))["elevation_m"]
+        np.testing.assert_array_equal(base, tiny)
+
+
+class MarginProfileTest(unittest.TestCase):
+    @staticmethod
+    def _profile(extra):
+        from mapgen.graph import distance_to, ocean_mask
+        land = lambda g: (np.abs(g.lat) < 30) & (np.abs(g.lon) < 60)
+        g = make_ctx(4).grid
+        frac = g.area_km2[land(g)].sum() / g.area_km2.sum()
+        ctx = make_ctx(4, tect=COLLIDE, cfg={"build": {"land_fraction": round(frac - 0.003, 4)},
+                                             "crust": {"edge_noise": 0.0}, "elevation": extra})
+        z0 = np.zeros(g.n)
+        ctx.data.update({"sk_land": land(g).astype(float), "m_land_hint": z0, "sk_mountains": z0,
+                         "m_mountain_hint": z0, "m_gravity_zones": z0, "m_lock": z0})
+        ctx.data.update(plates.run(ctx))
+        ctx.data.update(crust.run(ctx))
+        z = EL.run(ctx)["elevation_m"].astype(float)
+        sea = ocean_mask(g, z, 5.0e6)
+        d_coast_land = distance_to(g, sea)
+        d_coast_sea = distance_to(g, ~sea)
+        far = (~sea) & (d_coast_land > 1200) & (np.abs(g.lon) > 20)       # interior, away from the collision belt
+        coastal = (~sea) & (d_coast_land < 250) & (np.abs(g.lon) > 20)
+        shelf = sea & (d_coast_sea < 200)
+        return np.median(z[coastal]), np.median(z[far]), np.median(z[shelf])
+
+    def test_coastal_lowlands_and_shelves(self):
+        step = self._profile({"coast_noise_m": 0.0, "margin_km": 1.0, "slope_km": 1.0})
+        coastal, interior, shelf = self._profile({})
+        self.assertLess(coastal, interior - 300)              # coastal plains sit well below the interior
+        self.assertGreater(shelf, -800)                        # a shallow shelf fringes the coast
+        self.assertGreaterEqual(step[0], step[1] - 300)        # the old step margin had no coastal lowlands
+
+
+class ConnectedSeaLevelTest(unittest.TestCase):
+    def test_interior_pits_do_not_count_as_sea(self):
+        from mapgen.graph import ocean_mask
+        from tests.helpers import small_grid
+        g = small_grid(3)
+        land = (np.abs(g.lat) < 40) & (np.abs(g.lon) < 70)
+        z = np.where(land, 300.0 + 40.0 * (70.0 - np.abs(g.lon)), -4000.0)     # rises inland (continuous)
+        pits = land & ((g.lat % 10) < 3) & ((g.lon % 10) < 3) & (np.abs(g.lon) < 55)
+        z[pits] = -500.0                                                    # many deep interior pits
+        target = 0.8 * g.area_km2[land].sum() / g.area_km2.sum()
+        zs = EL.solve_sea_level_connected(g, z, target, 5.0e6)
+        sea = ocean_mask(g, zs, 5.0e6)
+        self.assertAlmostEqual(g.area_km2[~sea].sum() / g.area_km2.sum(), target, delta=0.01)
+
+
+class RevisionElevationTest(unittest.TestCase):
+    LAND = staticmethod(lambda g: np.abs(g.lon + 30) < 25)
+
+    def test_without_patches_the_relief_is_solved_once(self):
+        from unittest import mock
+        ctx = scenario(SUBDUCT, self.LAND, 0.1)
+        with mock.patch.object(EL, "_relief", wraps=EL._relief) as rel:
+            EL.run(ctx)
+        self.assertEqual(rel.call_count, 1)
+
+    def test_a_land_patch_keeps_the_sea_level_of_the_world_without_it(self):
+        from mapgen.crust import center_dist
+        za = EL.run(scenario(SUBDUCT, self.LAND, 0.1))["elevation_m"]
+        tect = {**SUBDUCT, "land_patch": [{"name": "fill", "center": [0.0, 120.0], "radius_km": 1500.0}]}
+        b = scenario(tect, self.LAND, 0.1)
+        zb = EL.run(b)["elevation_m"]
+        g = b.grid
+        far = (za > 0) & (center_dist(g, [0.0, 120.0]) > 5000.0)
+        self.assertGreater(far.sum(), 20)
+        np.testing.assert_allclose(zb[far], za[far], atol=1e-3)     # the other coasts don't move
+        i = g.cell_index(0.0, 120.0)
+        self.assertGreater(zb[i] - za[i], 1000.0, "the patch is continental ground now")
+
+    def test_a_land_patch_is_land_though_the_sea_level_drowns_bare_crust(self):
+        from mapgen.crust import center_dist
+        tect = {**SUBDUCT, "land_patch": [{"name": "fill", "center": [0.0, 120.0], "radius_km": 1500.0}]}
+        b = scenario(tect, self.LAND, 0.05)
+        z = EL.run(b)["elevation_m"]
+        inner = center_dist(b.grid, [0.0, 120.0]) < 1000.0
+        self.assertGreater(float(np.mean(z[inner] > 0)), 0.9)
+
+    def test_plateaus_get_their_surface_after_the_sea_level_solve(self):
+        from mapgen import plateaus as PL
+        from mapgen.crust import center_dist
+        plats = [{"name": "p", "center": [0.0, 120.0], "area_km2": 3.0e6, "top_m": [1500.0, 3000.0]},
+                 {"name": "q", "center": [40.0, 150.0], "area_km2": 1.0e6, "top_m": [1000.0, 1500.0],
+                  "islands": True}]
+        za = EL.run(scenario(SUBDUCT, self.LAND, 0.1))["elevation_m"]
+        b = scenario({**SUBDUCT, "plateau": plats}, self.LAND, 0.1)
+        zb = EL.run(b)["elevation_m"]
+        g, ids = b.grid, b.data["plateau_id"]
+        away = (center_dist(g, [0.0, 120.0]) > 2500.0) & (center_dist(g, [40.0, 150.0]) > 1500.0)
+        np.testing.assert_allclose(zb[away], za[away], atol=1e-3)
+        _, short = PL.semi_axes(plats[0])
+        core = np.flatnonzero(ids == 0)
+        core = core[(1.0 - PL.rho(g.xyz[core], plats[0], b.seed, g.radius_km)) * short > PL.MARGIN_KM]
+        self.assertGreater(len(core), 5)
+        self.assertTrue(1500.0 - PL.RELIEF_M <= float(np.median(-zb[core])) <= 3000.0 + PL.RELIEF_M)
+        self.assertLessEqual(zb[ids == 0].max(), PL.HIDDEN_MAX_M)
+        near_q = center_dist(g, [40.0, 150.0]) < 1500.0
+        self.assertGreater(zb[near_q].max(), 0.0, "the island plateau breaks the surface")
+
+    def test_land_added_marks_what_the_patch_makes_land(self):
+        from mapgen.crust import center_dist
+        from mapgen.graph import ocean_mask
+        a = EL.run(scenario(SUBDUCT, self.LAND, 0.05))
+        tect = {**SUBDUCT, "land_patch": [{"name": "fill", "center": [0.0, 120.0], "radius_km": 1500.0}]}
+        b = scenario(tect, self.LAND, 0.05)
+        ob = EL.run(b)
+        g = b.grid
+        land_a = ~ocean_mask(g, a["elevation_m"])
+        land_b = ~ocean_mask(g, ob["elevation_m"])
+        added = ob["land_added"]
+        self.assertGreater(float(added[center_dist(g, [0.0, 120.0]) < 1000.0].mean()), 0.9)
+        self.assertFalse((added & land_a).any(), "land without the patch is not added land")
+        self.assertFalse((added & ~land_b).any(), "added land is land")
+        self.assertFalse(a["land_added"].any(), "no patch, no plateau: nothing added")
diff --git a/tests/test_environment.py b/tests/test_environment.py
new file mode 100644
index 0000000..0f4c5c6
--- /dev/null
+++ b/tests/test_environment.py
@@ -0,0 +1,173 @@
+import unittest
+
+import numpy as np
+
+from mapgen import environment as EN
+from tests.helpers import make_ctx
+
+
+def name(bio, p, tmin):
+    z, _ = EN.holdridge(np.array([bio]), np.array([p]), np.array([tmin]))
+    return EN.HOLDRIDGE_NAMES[z[0]]
+
+
+class HoldridgeTest(unittest.TestCase):
+    def test_count(self):
+        self.assertEqual(len(EN.HOLDRIDGE_NAMES), 38)
+
+    def test_table(self):
+        self.assertEqual(name(26, 3000, 20), "tropical moist forest")
+        self.assertEqual(name(26, 9000, 20), "tropical rain forest")
+        self.assertEqual(name(26, 100, 20), "tropical desert")
+        self.assertEqual(name(0.5, 100, -30), "polar desert")
+        self.assertEqual(name(8, 300, -10), "cool temperate steppe")
+        self.assertEqual(name(15, 1500, -5), "warm temperate moist forest")
+        self.assertEqual(name(15, 1500, 5), "subtropical moist forest")
+        self.assertEqual(name(4, 700, -20), "boreal wet forest")
+        self.assertEqual(name(2, 200, -25), "subpolar moist tundra")
+
+    def test_polar_desert_with_summer_monsoon(self):
+        tag = EN.seasonality(np.array([180.0, 20.0]), np.array([20.0, 180.0]), np.array([80.0, 80.0]),
+                             np.array([-30.0, -30.0]), np.array([100.0, 100.0]), EN.DEFAULTS)
+        self.assertEqual(tag[0], EN.SEAS_W)   # NH: wet June (summer), dry December
+        self.assertEqual(tag[1], EN.SEAS_S)   # NH: dry summer
+        sh = EN.seasonality(np.array([20.0]), np.array([180.0]), np.array([-80.0]), np.array([-30.0]),
+                            np.array([100.0]), EN.DEFAULTS)
+        self.assertEqual(sh[0], EN.SEAS_W)    # SH summer is December
+
+    def test_tropical_monsoon(self):
+        tag = EN.seasonality(np.array([3000.0]), np.array([200.0]), np.array([15.0]), np.array([20.0]),
+                             np.array([1800.0]), EN.DEFAULTS)
+        self.assertEqual(tag[0], EN.SEAS_M)
+
+
+class StageTest(unittest.TestCase):
+    def test_run(self):
+        ctx = make_ctx(3)
+        g = ctx.grid
+        n = g.n
+        land = (np.abs(g.lat) < 40) & (np.abs(g.lon) < 60)
+        z = np.where(land, 300.0, -4000.0)
+        z = np.where(land & (np.abs(g.lon) < 4), 4000.0, z)
+        zeros = np.zeros(n)
+        ctx.data.update({
+            "elevation_eroded_m": z, "continental": land, "age_class": np.where(land, 1, 0).astype(np.int8),
+            "d_over_km": np.full(n, np.inf), "T_mean": np.where(np.abs(g.lat) > 30, -5.0, 20.0),
+            "T_min": np.full(n, -10.0), "T_range": np.full(n, 20.0), "P_ann": np.full(n, 800.0),
+            "P_jun": np.full(n, 800.0), "P_dec": np.full(n, 800.0), "biotemp": np.full(n, 10.0),
+            "dist_ocean_km": np.where(land, 1000.0, 0.0), "river": np.zeros(n, bool),
+            "strahler": np.zeros(n, np.int8), "discharge_km3_yr": zeros, "lake": np.zeros(n, bool),
+            "salt_flat": np.zeros(n, bool)})
+        out = EN.run(ctx)
+        self.assertEqual(out["landform"][g.cell_index(0.0, 0.0)], EN.LF_MOUNTAINS)
+        self.assertEqual(out["landform"][g.cell_index(0.0, 150.0)], EN.LF_OCEAN)
+        self.assertEqual(out["ground"][g.cell_index(35.0, 40.0)], EN.GR_PERMAFROST)
+        for k in ("holdridge", "hold_region", "seasonality", "landform", "relief_m", "lithology", "ground",
+                  "coal_potential", "iron_potential"):
+            self.assertEqual(len(out[k]), n)
+
+
+class OceanMaskEnvTest(unittest.TestCase):
+    def test_inland_basin_not_ocean_landform(self):
+        ctx = make_ctx(3)
+        g = ctx.grid
+        n = g.n
+        land = (np.abs(g.lat) < 40) & (np.abs(g.lon) < 60)
+        z = np.where(land, 300.0, -4000.0)
+        basin = (np.abs(g.lat) < 6) & (np.abs(g.lon) < 6)
+        z[basin] = -30.0
+        zeros = np.zeros(n)
+        ctx.data.update({
+            "elevation_eroded_m": z, "ocean": ~land, "continental": land, "age_class": np.where(land, 1, 0).astype(np.int8),
+            "d_over_km": np.full(n, np.inf), "T_mean": np.full(n, 20.0), "T_min": np.full(n, 10.0),
+            "T_range": np.full(n, 10.0), "P_ann": np.full(n, 800.0), "P_jun": np.full(n, 800.0),
+            "P_dec": np.full(n, 800.0), "biotemp": np.full(n, 20.0), "dist_ocean_km": np.where(land, 1000.0, 0.0),
+            "river": np.zeros(n, bool), "strahler": np.zeros(n, np.int8), "discharge_km3_yr": zeros,
+            "lake": np.zeros(n, bool), "salt_flat": np.zeros(n, bool)})
+        out = EN.run(ctx)
+        self.assertFalse(np.any(out["landform"][basin] == EN.LF_OCEAN))
+
+
+class ReliefResolutionTest(unittest.TestCase):
+    def test_relief_window_is_in_km(self):
+        from tests.helpers import small_grid
+        med = []
+        for res in (2, 3):
+            g = small_grid(res)
+            z = 1000.0 * np.sin(10.0 * g.xyz[:, 0]) + 800.0 * np.cos(9.0 * g.xyz[:, 1])
+            _, relief = EN.landform(g, z, np.ones(g.n, np.int8), np.ones(g.n, np.int8), np.full(g.n, np.inf),
+                                    np.full(g.n, 800.0), np.zeros(g.n, bool), relief_km=700.0)
+            med.append(np.median(relief))
+        self.assertAlmostEqual(med[0] / med[1], 1.0, delta=0.3)
+
+
+class WetlandRuleTest(unittest.TestCase):
+    def test_wetlands_need_real_rivers_on_real_flats(self):
+        from types import SimpleNamespace
+        from mapgen import graph as G
+        g = small_grid_env()
+        n = g.n
+        z = np.where(g.lat > -30, 100.0 + 0.0001 * g.lat, -3000.0)          # a very flat humid plain
+        common = dict(lit=np.full(n, EN.LI_GRANITE, np.int8), t_mean=np.full(n, 15.0), t_min=np.full(n, 5.0),
+                      p_ann=np.full(n, 1200.0), dist_ocean=np.full(n, 2000.0), river=np.zeros(n, bool),
+                      strahler=np.zeros(n, np.int8), lake=np.zeros(n, bool), salt_flat=np.zeros(n, bool))
+        def wet(q):
+            gr = EN.ground(g, z, common["lit"], common["t_mean"], common["t_min"], common["p_ann"], common["dist_ocean"],
+                           common["river"], common["strahler"], np.full(n, q), common["lake"], common["salt_flat"],
+                           EN.DEFAULTS, g.lat <= -30)
+            return np.mean(gr[g.lat > 0] == EN.GR_WETLAND)
+        self.assertLess(wet(3.0), 0.05)          # small streams don't make a wetland
+        self.assertGreater(wet(8.0), 0.9)        # a big river across a flat does
+
+
+def small_grid_env():
+    from tests.helpers import small_grid
+    return small_grid(3)
+
+
+class DepositsTest(unittest.TestCase):
+    def test_shares_rules_and_main(self):
+        from mapgen import minerals as MN
+        from mapgen import environment as EN
+        from mapgen.crust import CRATON, POST_OROGEN
+        g = make_ctx(4).grid
+        n = g.n
+        land = np.abs(g.lat) < 60
+        mountain = land & (g.lon > 0) & (g.lon < 60)
+        f = {"land": land, "z": np.where(mountain, 3500.0, 200.0), "lit": np.where(mountain, EN.LI_METAMORPHIC, EN.LI_SANDSTONE).astype(np.int8),
+             "age": np.where(mountain, POST_OROGEN, CRATON), "lf": np.where(mountain, EN.LF_MOUNTAINS, EN.LF_PLAIN),
+             "relief": np.where(mountain, 2500.0, 50.0), "t_mean": np.full(n, 15.0), "p_ann": np.full(n, 900.0),
+             "d_coll": np.full(n, np.inf), "salt": np.zeros(n, bool), "endo": np.zeros(n, bool),
+             "dist_ocean": np.full(n, 500.0), "coal": land & ~mountain}
+        bits, main = MN.place(g, f, 7)
+        self.assertEqual(len(MN.DEPOSITS), 32)
+        self.assertFalse(bits[~land].any())
+        for name in ("ironstone", "oil_gas", "iron_high", "tungsten"):
+            has = (bits & np.uint32(MN.BIT[name])) > 0
+            share = dict((d[0], d[2]) for d in MN.DEPOSITS)[name]
+            self.assertLessEqual(has.sum(), round(share * land.sum()) + 1, name)
+            self.assertTrue(has.any(), name)
+        for name in MN.DEEP[:2]:                          # the good stuff: in the mountains only here
+            has = (bits & np.uint32(MN.BIT[name])) > 0
+            self.assertTrue(np.all(mountain[has]), name)
+        some = bits > 0
+        self.assertTrue(np.all(main[some] > 0) and np.all(main[~some] == 0))
+        k = np.flatnonzero(some)[0]
+        self.assertTrue(bits[k] & np.uint32(1 << (int(main[k]) - 1)))
+
+    def test_deposits_are_sprinkled_not_one_hotspot(self):
+        from mapgen import minerals as MN
+        g = make_ctx(4).grid
+        land = np.abs(g.lat) < 60
+        hot = land & (g.lat > 0) & (g.lat < 30) & (g.lon > 0) & (g.lon < 40)      # a world-class district
+        score = np.where(hot, 3.0, np.where(land, 1.0, 0.0))
+        prov = MN.provinces(g, 7)
+        sel = MN._rank_select(score, 0.02, land, prov)
+        self.assertEqual(sel.sum(), round(0.02 * land.sum()))
+        self.assertFalse(sel[~land].any())
+        self.assertGreater((sel & hot).sum(), 0.3 * sel.sum())                 # the district keeps the most
+        self.assertGreater((sel & ~hot).sum(), 0.3 * sel.sum())                # the rest is sprinkled
+        self.assertGreater(len(np.unique(prov[sel & ~hot])), 0.3 * len(np.unique(prov[land & ~hot])))
+        self.assertFalse(MN._rank_select(np.where(hot, 1.0, 0.0), 0.02, land, prov)[~hot].any())   # rules still rule
+        poor = MN._rank_select(np.where(hot, 10.0, np.where(land, 1.0, 0.0)), 0.02, land, prov)
+        self.assertFalse(poor[~hot].any())                                     # far poorer ground: no local mines
diff --git a/tests/test_eras.py b/tests/test_eras.py
new file mode 100644
index 0000000..96b7c82
--- /dev/null
+++ b/tests/test_eras.py
@@ -0,0 +1,303 @@
+import shutil
+import importlib.util
+import tempfile
+import unittest
+from pathlib import Path
+
+import numpy as np
+
+from mapgen import eras as ER, events as EV, pipeline as P
+from mapgen.graph import components, distance_to
+from mapgen.sphere import gc_dist_km, latlon_to_xyz, xyz_to_latlon
+from tests.test_render import small_world
+
+R = 12742.0
+ERAS = """
+[[event]]
+name = "cut"
+kind = "disintegrate"
+center = {cut}
+radius_km = 3000.0
+depth_m = 3000.0
+[[event]]
+name = "aura"
+kind = "zone"
+shape = "landmass"
+seed = {seed}
+reach_km = [500.0, 1000.0]
+edge_km = 5.0
+fields = {{ gravity_g = 0.35, pressure_bar = 2.0, o2_fraction = 0.35, fire_reactivity = {fire} }}
+[[event]]
+name = "rise"
+kind = "volcano"
+center = {vol}
+radius_km = 900.0
+peak_m = 1500.0
+peak_mode = "absolute"
+[[event]]
+name = "ghost"
+kind = "zone"
+shape = "landmass"
+seed = {seed}
+reach_km = [100.0, 200.0]
+edge_km = 5.0
+fields = {{ fire_reactivity = 0.2 }}
+[eras]
+order = ["before", "glow", "after", "risen", "haunt"]
+default = "after"
+[eras.before]
+label = "Before"
+events = []
+[eras.glow]
+label = "Glow"
+events = ["aura"]
+[eras.after]
+label = "After"
+events = ["cut"]
+[eras.risen]
+label = "Risen"
+events = ["rise"]
+years = 8000
+[eras.haunt]
+label = "Haunt"
+events = ["ghost"]
+"""
+
+
+def cli():
+    spec = importlib.util.spec_from_file_location("mapgen_cli", Path(ER.__file__).parents[1] / "mapgen.py")
+    m = importlib.util.module_from_spec(spec)
+    spec.loader.exec_module(m)                          # the script (mapgen/ is the package)
+    return m
+
+
+def latlon(p):
+    lat, lon = xyz_to_latlon(np.asarray(p, dtype=np.float64))
+    return [round(float(lat), 3), round(float(lon), 3)]
+
+
+class ErasTest(unittest.TestCase):
+    @classmethod
+    def setUpClass(cls):
+        cls.root = Path(tempfile.mkdtemp())
+        small_world(cls.root)
+        cls.base = P.build(cls.root, 2, log=lambda m: None)
+        g, d = cls.base.grid, cls.base.data
+        land = ~np.asarray(d["ocean"])
+        lab = components(g, land)
+        big = lab == np.bincount(lab[land]).argmax()
+        mid = g.xyz[big].mean(axis=0)
+        cls.cut = latlon(g.xyz[big][np.argmax(g.xyz[big] @ (mid / np.linalg.norm(mid)))])
+        far = gc_dist_km(g.xyz, latlon_to_xyz(*cls.cut), R)
+        cls.seed = latlon(g.xyz[np.argmax(np.where(big, far, -1.0))])     # stays land after the cut
+        cls.vol_i = int(np.argmax(distance_to(g, land)))                   # the loneliest sea
+        cls.vol = latlon(g.xyz[cls.vol_i])
+        cls.g0 = float(cls.base.cfg["planet"]["gravity_g"])
+        cls.write_eras()
+        for n in ("glow", "after", "risen", "haunt"):
+            ER.build_era(cls.root, 2, n, log=lambda m: None, base=cls.base)
+
+    def test_low_memory_era_writes_identical_files(self):
+        # oracle: the normal-mode era built in setUpClass
+        tmp = Path(tempfile.mkdtemp())
+        try:
+            shutil.copytree(self.root, tmp, dirs_exist_ok=True)
+            shutil.rmtree(ER.era_dir(tmp, 2, "glow"))
+            ER.build_era(tmp, 2, "glow", log=lambda m: None, low_memory=True)
+            a, b = ER.era_dir(self.root, 2, "glow"), ER.era_dir(tmp, 2, "glow")
+            fa = sorted(p.relative_to(a) for p in a.rglob("*") if p.is_file())
+            self.assertEqual(fa, sorted(p.relative_to(b) for p in b.rglob("*") if p.is_file()))
+            for f in fa:
+                self.assertEqual((a / f).read_bytes(), (b / f).read_bytes(), str(f))
+        finally:
+            shutil.rmtree(tmp)
+
+    @classmethod
+    def write_eras(cls, fire=0.5):
+        (cls.root / "config" / "eras.toml").write_text(ERAS.format(cut=cls.cut, seed=cls.seed, vol=cls.vol, fire=fire))
+
+    def load(self, name=None):
+        d = self.root / "out" / "r2" if name is None else ER.era_dir(self.root, 2, name)
+        with np.load(d / "cells.npz") as z:
+            return {k: z[k] for k in z.files}
+
+    def test_outside_the_mask_is_the_base(self):
+        b = self.load()
+        for n in ("glow", "after", "risen"):
+            e = self.load(n)
+            m = e["era_mask"]
+            self.assertTrue(m.any() and not m.all(), n)
+            for k, v in b.items():
+                if v.shape[:1] == m.shape:
+                    np.testing.assert_array_equal(e[k][~m], v[~m], err_msg=f"{n}: {k}")
+
+    def test_deposits_are_the_base_geology(self):
+        b = self.load()                                      # an event moves ground; it neither makes nor unmakes ore
+        for n in ("glow", "after", "risen"):
+            e = self.load(n)
+            for k in ("deposits", "deposit_main", "iron_potential"):
+                np.testing.assert_array_equal(e[k], b[k], err_msg=f"{n}: {k}")
+
+    def test_disintegration_only_lowers_and_reaches_depth(self):
+        b, e = self.load(), self.load("after")
+        zb, ze = b["elevation_eroded_m"].astype(np.float64), e["elevation_eroded_m"].astype(np.float64)
+        self.assertTrue(np.all(ze <= zb + 1e-3))
+        d = gc_dist_km(b["g_xyz"], latlon_to_xyz(*self.cut), R)
+        _, z_ref = EV.disintegrate(b["g_xyz"], zb, ~b.get("open_water", b["ocean"]), {"center": self.cut, "radius_km": 3000.0,
+                                                                "depth_m": 3000.0}, R)
+        inner = d < 1500.0
+        self.assertTrue(inner.any())
+        self.assertTrue(np.all(ze[inner] <= z_ref - 3000.0 * np.sqrt(0.75) + 1e-3))
+        self.assertTrue((e["ocean"] | e["lake"])[inner].all(), "the bowl fills with water")
+
+    def test_zone_fields_hold_inside_the_landmass(self):
+        b, e = self.load(), self.load("glow")
+        land, w = e["zone_aura_land"], e["zone_aura"]
+        self.assertTrue(land.any())
+        self.assertTrue(np.all(w[land] == 1.0))
+        np.testing.assert_allclose(e["gravity_g"][land], 0.35, rtol=1e-6)
+        np.testing.assert_allclose(e["fire_reactivity"][land], 0.5)
+        np.testing.assert_allclose(e["po2_bar"], e["o2_fraction"] * e["pressure_bar"], rtol=1e-5)
+        np.testing.assert_array_equal(e["gravity_g"][w == 0], b["gravity_g"][w == 0])
+
+    def test_zone_follows_the_landmass_before_the_era(self):
+        b, e = self.load(), self.load("after")
+        d = gc_dist_km(b["g_xyz"], latlon_to_xyz(*self.cut), R)
+        lost = ~b["open_water"] & e["open_water"] & (d < 1500.0)
+        self.assertTrue(lost.any())
+        self.assertTrue(e["zone_aura_land"][lost].all(), "the old coast, not the bitten one")
+        np.testing.assert_allclose(e["gravity_g"][lost], 0.35, rtol=1e-6)
+
+    def test_plants_in_a_zone(self):
+        from mapgen import environment as EN
+        e = self.load("glow")
+        land = e["zone_aura_land"] & ~e["ocean"]
+        zone, _ = EN.holdridge(e["biotemp"], e["P_ann"] * np.sqrt(2.0), e["T_min"])
+        np.testing.assert_array_equal(e["holdridge"][land], zone[land])
+        np.testing.assert_allclose(e["plant_height_x"], self.g0 / e["gravity_g"], rtol=1e-5)
+
+    def test_era_years_set_its_erosion(self):
+        e = self.load("risen")
+        self.assertGreater(float(e["elevation_eroded_m"][self.vol_i]), 1400.0, "8,000 years barely wear a cone")
+        self.assertNotEqual(ER.fingerprint("k", [("a", [], 8000)]), ER.fingerprint("k", [("a", [], 9000)]))
+
+    def test_steps_list_each_eras_own_events_and_years(self):
+        from mapgen import config as C
+        st = ER.steps(C.load(self.root)[1], "risen")
+        self.assertEqual([(n, [e["name"] for e in ev], y) for n, ev, y in st],
+                         [("before", [], None), ("glow", ["aura"], None), ("after", ["cut"], None),
+                          ("risen", ["rise"], 8000)])
+
+    def test_a_later_zone_follows_the_land_of_its_own_era(self):
+        b, e = self.load(), self.load("haunt")
+        d = gc_dist_km(b["g_xyz"], latlon_to_xyz(*self.cut), R)
+        lost = ~b["open_water"] & self.load("after")["open_water"] & (d < 1500.0)
+        self.assertTrue(lost.any())
+        self.assertFalse(e["zone_ghost_land"][lost].any(), "ghost came after the cut: the bitten coast")
+        self.assertTrue(e["zone_aura_land"][lost].all(), "aura still follows the land before the cut")
+
+    def test_a_new_label_needs_no_rebuild(self):
+        import json
+        f = ER.era_dir(self.root, 2, "glow") / "cells.npz"
+        before = f.stat().st_mtime_ns
+        p = self.root / "config" / "eras.toml"
+        p.write_text(p.read_text().replace('label = "Glow"', 'label = "The Glow"'))
+        try:
+            logs = []
+            ER.build_era(self.root, 2, "glow", log=logs.append)
+            self.assertIn("era glow: cached", logs)
+            meta = json.loads((ER.era_dir(self.root, 2, "glow") / "cells_meta.json").read_text())
+            self.assertEqual(meta["era"]["label"], "The Glow")
+            self.assertEqual(f.stat().st_mtime_ns, before)
+        finally:
+            self.write_eras()
+            ER.build_era(self.root, 2, "glow", log=lambda m: None)
+
+    def test_a_renamed_era_whose_folder_moved_needs_no_rebuild(self):
+        import json
+        old, new = ER.era_dir(self.root, 2, "glow"), ER.era_dir(self.root, 2, "shine")
+        before = (old / "cells.npz").stat().st_mtime_ns
+        p = self.root / "config" / "eras.toml"
+        p.write_text(p.read_text().replace('"glow"', '"shine"').replace("[eras.glow]", "[eras.shine]"))
+        old.rename(new)
+        try:
+            logs = []
+            ER.build_era(self.root, 2, "shine", log=logs.append)
+            self.assertIn("era shine: cached", logs)
+            self.assertEqual(json.loads((new / "cells_meta.json").read_text())["era"]["name"], "shine")
+            self.assertEqual((new / "cells.npz").stat().st_mtime_ns, before)
+        finally:
+            new.rename(old)
+            self.write_eras()
+            ER.build_era(self.root, 2, "glow", log=lambda m: None)
+
+    def test_one_lake_id_never_names_two_lakes(self):
+        b, e = self.load(), self.load("after")
+        m, ids = e["era_mask"], e["lake_id"]
+        for k in np.unique(ids[m & (ids >= 0)]).tolist():
+            if (ids[~m] == k).any():
+                np.testing.assert_array_equal(ids == k, b["lake_id"] == k, err_msg=f"lake {k}")
+
+    def test_zone_events_never_change_height(self):
+        b, e = self.load(), self.load("glow")
+        for k in ("elevation_m", "elevation_eroded_m", "z_surface_m"):
+            np.testing.assert_array_equal(e[k], b[k], err_msg=k)
+
+    def test_a_volcano_raises_new_land(self):
+        e = self.load("risen")
+        self.assertFalse(e["ocean"][self.vol_i])
+        self.assertGreater(float(e["elevation_eroded_m"][self.vol_i]), 0.0)
+
+    def test_era_is_cached_until_its_events_change(self):
+        logs = []
+        ER.build_era(self.root, 2, "glow", log=logs.append)
+        self.assertIn("era glow: cached", logs)
+        first = self.load("glow")
+        self.write_eras(fire=0.6)
+        try:
+            logs, base_logs = [], []
+            ER.build_era(self.root, 2, "glow", log=logs.append)
+            self.assertNotIn("era glow: cached", logs)
+            P.build(self.root, 2, log=base_logs.append)
+            self.assertTrue(base_logs and all(line.endswith("cached") for line in base_logs),
+                            "tuning an event never rebuilds the base")
+        finally:
+            self.write_eras()
+            ER.build_era(self.root, 2, "glow", log=lambda m: None)
+        again = self.load("glow")                             # repeatable: the same inputs, the same arrays
+        for k in first:
+            np.testing.assert_array_equal(again[k], first[k], err_msg=k)
+
+    def test_build_command_builds_eras_and_removes_stale_ones(self):
+        stale = ER.era_dir(self.root, 2, "gone")
+        stale.mkdir(parents=True, exist_ok=True)
+        self.assertEqual(cli().main(["build", "--res", "2"], root=self.root), 0)
+        self.assertFalse(stale.exists())
+        for n in ("glow", "after", "risen"):
+            self.assertTrue((ER.era_dir(self.root, 2, n) / "cells.npz").exists(), n)
+        self.assertFalse(ER.era_dir(self.root, 2, "before").exists(), "the base era is the base build")
+        self.assertEqual(cli().main(["era", "glow", "--res", "2"], root=self.root), 0)
+        self.assertEqual(cli().main(["era", "nowhen", "--res", "2"], root=self.root), 2)
+
+
+class MergeLakeIdsTest(unittest.TestCase):
+    def test_kept_lakes_keep_ids_changed_and_new_lakes_get_fresh_ones(self):
+        base = np.array([0, 0, 1, 1, -1, 2, -1], np.int32)
+        new = np.array([5, 5, 0, 0, 0, -1, 1], np.int32)   # 5 = base 0 as it was; 0 = base 1 grown; 1 new; 2 gone
+        out = ER.merge_lake_ids(base, new, np.ones(7, bool))
+        np.testing.assert_array_equal(out, [0, 0, 3, 3, 3, -1, 4])
+        self.assertEqual(out.dtype, np.int32)
+
+    def test_outside_the_mask_stays_the_base(self):
+        base = np.array([0, 0, -1, 1])
+        new = np.array([3, 3, 7, 7])
+        out = ER.merge_lake_ids(base, new, np.array([False, False, True, True]))
+        np.testing.assert_array_equal(out, [0, 0, 2, 2])
+
+
+class SameTest(unittest.TestCase):
+    def test_nan_equals_nan_in_float_arrays_only(self):
+        self.assertTrue(ER._same(np.array([np.nan, 1.0]), np.array([np.nan, 1.0])))
+        self.assertFalse(ER._same(np.array([np.nan, 1.0]), np.array([np.nan, 2.0])))
+        self.assertTrue(ER._same(np.array([1, 2]), np.array([1, 2])))
+        self.assertTrue(ER._same(np.array([True]), np.array([True])))
diff --git a/tests/test_erosion.py b/tests/test_erosion.py
new file mode 100644
index 0000000..2adefc0
--- /dev/null
+++ b/tests/test_erosion.py
@@ -0,0 +1,111 @@
+import unittest
+
+import numpy as np
+
+from mapgen import erosion as ER
+from mapgen.crust import CRATON
+from tests.helpers import make_ctx
+
+
+def mountain_ctx(k=None):
+    cfg = {"build": {"land_fraction": 0.05}}
+    if k is not None:
+        cfg["erosion"] = {"k": k}
+    ctx = make_ctx(3, cfg=cfg)
+    g = ctx.grid
+    d = g.radius_km * np.arccos(np.clip(g.xyz @ g.xyz[g.cell_index(0.0, 0.0)], -1, 1))
+    z = np.where(d < 6000, 6000 * np.clip(1 - d / 6000, 0, None) ** 1.2 + 50, -4000.0)
+    ctx.data.update({"elevation_m": z.astype(np.float32), "age_class": np.full(g.n, CRATON, np.int8)})
+    return ctx, z
+
+
+class ErosionTest(unittest.TestCase):
+    def test_erode_lowers_land(self):
+        ctx, z0 = mountain_ctx()
+        kmult = np.full(ctx.grid.n, 1.5)
+        z1 = ER.erode(ctx.grid, z0, kmult, ER.DEFAULTS)
+        land = z0 > 0
+        self.assertLess(z1[land].mean(), z0[land].mean())
+        self.assertLessEqual(z1.max(), z0.max() + 1e-6)
+
+    def test_stage_keeps_land_fraction(self):
+        ctx, _ = mountain_ctx()
+        z = ER.run(ctx)["elevation_eroded_m"]
+        g = ctx.grid
+        self.assertTrue(np.all(np.isfinite(z)))
+        self.assertAlmostEqual(g.area_km2[z > 0].sum() / g.area_km2.sum(), 0.05, delta=0.01)
+
+    def test_stable_with_huge_k(self):
+        ctx, _ = mountain_ctx(k=1e3)
+        z = ER.run(ctx)["elevation_eroded_m"]
+        self.assertTrue(np.all(np.isfinite(z)))
+        self.assertGreaterEqual(z.min(), -11000)
+
+    def test_stage_outputs_connected_ocean(self):
+        ctx, z0 = mountain_ctx()
+        g = ctx.grid
+        c = g.xyz[g.cell_index(0.0, 0.0)]
+        d = g.radius_km * np.arccos(np.clip(g.xyz @ c, -1, 1))
+        z0 = np.where(d < 700, -500.0, z0)                     # deep interior pit in the mountain block
+        ctx.data["elevation_m"] = z0.astype(np.float32)
+        out = ER.run(ctx)
+        self.assertIn("ocean", out)
+        self.assertFalse(out["ocean"][d < 300].any())
+        self.assertTrue(out["ocean"][d > 7000].all())
+
+    def test_stage_adds_the_elevation_stages_added_land(self):
+        ctx, z0 = mountain_ctx()
+        ctx.cfg["build"]["land_fraction"] = 0.02
+        g = ctx.grid
+        mtn = z0 > 0
+        ctx.data.update({"continental": mtn, "continental_base": mtn})   # lifted shelf: no new crust at all
+        added = mtn & (g.lon >= 0)
+        ctx.data["land_added"] = added
+        extra = g.area_km2[added].sum() / g.area_km2.sum()
+        z = ER.run(ctx)["elevation_eroded_m"]
+        self.assertAlmostEqual(g.area_km2[z > 0].sum() / g.area_km2.sum(), 0.02 + extra, delta=0.005)
+
+class BaseLevelTest(unittest.TestCase):
+    def test_erosion_keeps_coastline(self):
+        ctx = make_ctx(3)
+        g = ctx.grid
+        block = (np.abs(g.lat) < 30) & (np.abs(g.lon) < 50)
+        z0 = np.where(block, 800.0 + 2500.0 * np.exp(-((g.lon + 45) / 6.0) ** 2), -6000.0)   # coastal range by a deep sea
+        z1 = ER.erode(g, z0, np.full(g.n, 1.5), ER.DEFAULTS)
+        np.testing.assert_array_equal(z1 > 0, z0 > 0)
+        self.assertTrue(np.all(z1[~block] == z0[~block]))
+
+
+class HillslopeResolutionTest(unittest.TestCase):
+    def test_hillslope_smoothing_is_resolution_independent(self):
+        from tests.helpers import small_grid
+        drops = []
+        for res in (3, 4):
+            g = small_grid(res)
+            c = g.xyz[g.cell_index(20.0, 0.0)]
+            d = g.radius_km * np.arccos(np.clip(g.xyz @ c, -1, 1))
+            z0 = np.where(g.lat > -40, 500.0 + 3000.0 * np.exp(-(d / 600.0) ** 2), -4000.0)
+            P = dict(ER.DEFAULTS, k=0.0)                       # isolate hillslope smoothing
+            z1 = ER.erode(g, z0, np.ones(g.n), P)
+            drops.append(z0.max() - z1[d < 150].max())
+        self.assertAlmostEqual(drops[0] / max(drops[1], 1e-9), 1.0, delta=0.35)
+
+
+class SedimentFillTest(unittest.TestCase):
+    def test_shallow_basins_fill_deep_ones_remain(self):
+        g = make_ctx(3).grid
+        z0 = np.where(g.lat > -20, 200.0 + 2.0 * (g.lat + 20), -3000.0)
+        depth = {}
+        for name, (lat, dz) in {"shallow": (20.0, 60.0), "deep": (50.0, 2500.0)}.items():
+            c = g.xyz[g.cell_index(lat, 0.0)]
+            d = g.radius_km * np.arccos(np.clip(g.xyz @ c, -1, 1))
+            z0 = np.where(d < 700, z0 - dz * (1 - d / 700), z0)
+            depth[name] = d < 250
+        from mapgen.graph import priority_flood
+        def fill_after(frac):
+            z1 = ER.erode(g, z0, np.full(g.n, 1.0), dict(ER.DEFAULTS, k=0.0, sediment_fill=frac))
+            return priority_flood(g, z1, z1 <= -1000) - z1
+        self.assertGreater(fill_after(0.0)[depth["shallow"]].max(), 20.0)     # without sediment it stays a basin
+        fill = fill_after(ER.DEFAULTS["sediment_fill"])
+        self.assertLess(fill[depth["shallow"]].max(), 0.2 * 60.0)             # noise-scale basins mostly filled
+        self.assertGreater(fill[depth["deep"]].max(), 20.0)
diff --git a/tests/test_events.py b/tests/test_events.py
new file mode 100644
index 0000000..1b3d51a
--- /dev/null
+++ b/tests/test_events.py
@@ -0,0 +1,121 @@
+import unittest
+
+import numpy as np
+
+from mapgen import events as EV
+from mapgen.pipeline import StageError
+from mapgen.sphere import east_north, gc_dist_km, great_circle_point, latlon_to_xyz
+from tests.helpers import small_grid
+
+R = 12742.0
+CUT = {"name": "cut", "kind": "disintegrate", "center": [0.0, 0.0], "radius_km": 3000.0, "depth_m": 3000.0}
+
+
+class DisintegrateTest(unittest.TestCase):
+    def test_bowl_only_lowers_and_reaches_its_depth(self):
+        g = small_grid(3)
+        d = gc_dist_km(g.xyz, latlon_to_xyz(0.0, 0.0), R)
+        z = np.where(d < 4000, 800.0, -4000.0)
+        z2, z_ref = EV.disintegrate(g.xyz, z, z > 0, CUT, R)
+        self.assertAlmostEqual(z_ref, 800.0)
+        self.assertTrue(np.all(z2 <= z))
+        np.testing.assert_array_equal(z2[d >= 3000], z[d >= 3000])
+        c = int(np.argmin(d))
+        self.assertAlmostEqual(z2[c], 800.0 - 3000.0 * np.sqrt(1 - (d[c] / 3000.0) ** 2), places=6)
+
+
+    def test_disintegration_over_open_sea_has_no_rim_height(self):
+        g = small_grid(3)
+        z2, z_ref = EV.disintegrate(g.xyz, np.full(g.n, -1000.0), np.zeros(g.n, bool), CUT, R)
+        self.assertIsNone(z_ref)
+        self.assertLess(float(z2[g.cell_index(0.0, 0.0)]), -2990.0, "the bowl hangs from 0 m")
+
+    def test_heights_warn_when_a_cut_has_no_rim_land(self):
+        from mapgen import eras as ER
+        g = small_grid(3)
+        logs = []
+        ER._heights(g, np.full(g.n, -1000.0), [dict(CUT)], 5.0e6, logs.append, "x")
+        self.assertTrue(any("no land in its rim" in m for m in logs), logs)
+
+
+class VolcanoTest(unittest.TestCase):
+    def test_shapes_and_modes(self):
+        g = small_grid(3)
+        d = gc_dist_km(g.xyz, latlon_to_xyz(0.0, 90.0), R)
+        c = int(np.argmin(d))
+        z = np.full(g.n, -3000.0)
+        above = EV.volcano(g.xyz, z, {"center": [0.0, 90.0], "radius_km": 900.0, "peak_m": 2000.0}, R)
+        self.assertTrue(np.all(above >= z))
+        self.assertAlmostEqual(above[c], -3000.0 + 2000.0 * (1 - d[c] / 900.0) ** 1.5)
+        absolute = EV.volcano(g.xyz, z, {"center": [0.0, 90.0], "radius_km": 900.0, "peak_m": 1500.0,
+                                         "peak_mode": "absolute"}, R)
+        self.assertGreater(absolute[c], 0.0, "a volcano above 0 m makes land")
+        np.testing.assert_array_equal(absolute[d >= 900.0], z[d >= 900.0])
+        cal = EV.volcano(g.xyz, np.zeros(g.n), {"center": [0.0, 90.0], "radius_km": 2000.0, "peak_m": 2000.0,
+                                                "shape": "caldera"}, R)
+        ring = (d > 300.0) & (d < 700.0)
+        self.assertLess(cal[c], cal[ring].max(), "a caldera's middle sits below its rim")
+
+
+class ZoneWeightTest(unittest.TestCase):
+    def test_sharp_circle_edge_ramps_over_edge_km(self):
+        ev = {"name": "z", "shape": "circle", "center": [10.0, 20.0], "radius_km": 800.0, "edge_km": 5.0}
+        p0 = latlon_to_xyz(10.0, 20.0)
+        e, n = east_north(p0[None])
+        s = np.arange(780.0, 810.0, 0.5)
+        pts = np.array([great_circle_point(p0, e[0], si, R) for si in s])
+        w = EV.zone_weight(pts, ev, R, 1296)
+        self.assertTrue(np.all(w[s <= 794.5] == 1.0))
+        self.assertTrue(np.all(w[s >= 800.5] == 0.0))
+        ramp = s[(w > 0) & (w < 1)]
+        self.assertLessEqual(ramp.max() - ramp.min(), 5.0)
+
+    def test_landmass_reach_varies_between_its_bounds(self):
+        g = small_grid(3)
+        land = (np.abs(g.lat) < 20) & (np.abs(g.lon) < 30)
+        ev = {"name": "aura", "shape": "landmass", "seed": [0.0, 0.0], "reach_km": [500.0, 1000.0], "edge_km": 5.0}
+        lm = EV.landmass(g, land, ev["seed"], ev["name"])
+        np.testing.assert_array_equal(lm, land)
+        w = EV.zone_weight(g.xyz, ev, R, 1296, g.xyz[lm])
+        self.assertTrue(np.all(w[land] == 1.0))
+        from scipy.spatial import cKDTree
+        chord, _ = cKDTree(g.xyz[lm]).query(g.xyz)
+        d = 2 * np.arcsin(chord / 2) * R
+        self.assertTrue(np.all(w[d > 1000.0] == 0.0))
+        self.assertTrue(np.all(w[(d > 0) & (d < 495.0)] == 1.0))
+
+    def test_landmass_seed_in_the_sea_names_the_event(self):
+        g = small_grid(3)
+        land = (np.abs(g.lat) < 20) & (np.abs(g.lon) < 30)
+        with self.assertRaises(StageError) as cm:
+            EV.landmass(g, land, [0.0, 150.0], "impact-aftermath")
+        self.assertIn("impact-aftermath", str(cm.exception))
+        self.assertIn("nearest land", str(cm.exception))
+
+
+class ApplyZoneTest(unittest.TestCase):
+    EV = {"name": "aura", "fields": {"gravity_g": 0.35, "pressure_bar": 2.0, "o2_fraction": 0.35,
+                                     "fire_reactivity": 0.5}}
+
+    def base(self, n=4):
+        return {"gravity_g": np.full(n, 1.05, np.float32), "o2_fraction": np.array([0.21, 0.3, 0.1, 0.21]),
+                "pressure_bar": np.array([1.0, 1.0, np.exp(-1), 1.0]), "fire_reactivity": np.ones(n),
+                "po2_bar": np.array([0.21, 0.3, 0.1 * np.exp(-1), 0.21])}
+
+    def test_zone_values_are_absolute_over_any_base(self):
+        f = self.base()
+        w = np.array([1.0, 1.0, 1.0, 0.0])
+        z = np.array([0.0, -50.0, 8000.0, 0.0])
+        out = EV.apply_zone(f, w, self.EV, z, 8000.0)
+        np.testing.assert_allclose(out["o2_fraction"][:3], 0.35)
+        np.testing.assert_allclose(out["gravity_g"][:3], 0.35, rtol=1e-6)
+        np.testing.assert_allclose(out["pressure_bar"][:3], [2.0, 2.0, 2.0 * np.exp(-1)])
+        np.testing.assert_allclose(out["po2_bar"], out["o2_fraction"] * out["pressure_bar"])
+        self.assertEqual(out["gravity_g"].dtype, np.float32)
+        for k in out:
+            self.assertEqual(out[k][3], f[k][3], f"{k}: untouched outside the zone")
+
+    def test_half_weight_blends(self):
+        out = EV.apply_zone(self.base(), np.full(4, 0.5), self.EV, np.zeros(4), 8000.0)
+        self.assertAlmostEqual(float(out["fire_reactivity"][0]), 0.75)
+        self.assertAlmostEqual(float(out["pressure_bar"][0]), 1.5)
diff --git a/tests/test_geo.py b/tests/test_geo.py
new file mode 100644
index 0000000..d438299
--- /dev/null
+++ b/tests/test_geo.py
@@ -0,0 +1,84 @@
+import unittest
+from types import SimpleNamespace
+
+import numpy as np
+
+from mapgen import geo as GE
+from mapgen.plates import CONV
+from tests.helpers import small_grid
+
+
+def _area(ring):
+    a = np.asarray(ring, dtype=float)
+    return 0.5 * float(np.sum(a[:-1, 0] * a[1:, 1] - a[1:, 0] * a[:-1, 1]))
+
+
+class GeoTest(unittest.TestCase):
+    def test_split_antimeridian(self):
+        parts = GE.split_antimeridian([[178.0, 0.0], [179.5, 1.0], [-179.5, 1.0], [-178.0, 0.0]])
+        self.assertEqual(len(parts), 2)
+        self.assertEqual(parts[0][-1], [179.5, 1.0])
+
+    def test_river_split_at_antimeridian(self):
+        g = SimpleNamespace(n=4, lat=np.zeros(4), lon=np.array([170.0, 178.0, -178.0, -170.0]))
+        recv = np.array([1, 2, 3, 3])
+        feats = GE.river_lines(g, recv, np.array([True, True, True, False]), np.array([1, 1, 2, 0], np.int8),
+                               np.array([1.0, 2.0, 3.0, 0.0]))
+        self.assertEqual(len(feats), 2)
+        for f in feats:
+            xs = [c[0] for c in f["geometry"]["coordinates"]]
+            self.assertLess(max(xs) - min(xs), 180)
+
+    def test_contour_square(self):
+        m = np.zeros((10, 20), bool)
+        m[3:7, 5:9] = True
+        rings = GE.contours(m)
+        self.assertEqual(len(rings), 1)
+        r = np.array(rings[0])
+        self.assertEqual(tuple(r[0]), tuple(r[-1]))
+        self.assertTrue(4 <= r[:, 0].min() and r[:, 0].max() <= 9 and 2 <= r[:, 1].min() and r[:, 1].max() <= 7)
+
+    def test_contour_ring_split(self):
+        m = np.zeros((10, 20), bool)
+        m[3:7, :3] = True
+        m[3:7, 17:] = True       # one island wrapping the antimeridian
+        feats = GE.contour_features(m, "land")
+        self.assertEqual(len(feats), 1)
+        geom = feats[0]["geometry"]
+        self.assertEqual(geom["type"], "MultiPolygon")           # fillable, split at ±180
+        for poly in geom["coordinates"]:
+            ring = poly[0]
+            self.assertEqual(ring[0], ring[-1])
+            xs = [c[0] for c in ring]
+            self.assertTrue(-180 <= min(xs) and max(xs) <= 180 and max(xs) - min(xs) < 180)
+            self.assertGreater(_area(ring), 0)                   # counter-clockwise exterior
+
+    def test_polygon_orientation_and_holes(self):
+        m = np.zeros((20, 40), bool)
+        m[4:16, 8:24] = True
+        m[8:12, 13:19] = False    # a lake-shaped hole
+        feats = GE.contour_features(m, "land")
+        self.assertEqual(len(feats), 1)
+        geom = feats[0]["geometry"]
+        self.assertEqual(geom["type"], "Polygon")
+        self.assertEqual(len(geom["coordinates"]), 2)
+        self.assertGreater(_area(geom["coordinates"][0]), 0)      # exterior CCW
+        self.assertLess(_area(geom["coordinates"][1]), 0)         # hole CW
+
+    def test_river_features_have_single_order(self):
+        g = SimpleNamespace(n=6, lat=np.arange(6.0), lon=np.zeros(6))
+        recv = np.array([1, 2, 3, 4, 5, 5])
+        order = np.array([1, 1, 2, 2, 3, 0], np.int8)
+        feats = GE.river_lines(g, recv, np.array([True] * 5 + [False]), order, np.arange(6.0))
+        self.assertEqual([f["properties"]["order"] for f in feats], [1, 2, 3])
+        for f in feats:
+            self.assertGreaterEqual(len(f["geometry"]["coordinates"]), 2)
+
+    def test_plate_boundaries(self):
+        g = small_grid(1)
+        plate = (g.lon > 0).astype(np.int16)
+        bt = np.full(g.n, CONV, np.int8)
+        feats = GE.boundary_features(g, plate, bt)
+        self.assertEqual(len(feats), 1)
+        self.assertEqual(feats[0]["properties"]["type"], "convergent")
+        self.assertGreater(len(feats[0]["geometry"]["coordinates"]), 5)
diff --git a/tests/test_graph.py b/tests/test_graph.py
new file mode 100644
index 0000000..ca324cc
--- /dev/null
+++ b/tests/test_graph.py
@@ -0,0 +1,362 @@
+import unittest
+
+import numpy as np
+
+from mapgen import graph as G
+from mapgen.sphere import latlon_to_xyz, east_north
+from tests.helpers import small_grid
+
+
+class GraphTest(unittest.TestCase):
+    def setUp(self):
+        self.g = small_grid(2)
+
+    def test_mean_max_min_diffuse(self):
+        g = self.g
+        f = np.zeros(g.n)
+        f[0] = 6.0
+        m = G.nbr_mean(g, f)
+        self.assertAlmostEqual(m[g.nbr_idx[g.nbr_ptr[0]]], 6.0 / g.counts[g.nbr_idx[g.nbr_ptr[0]]])
+        self.assertEqual(G.nbr_max(g, f)[g.nbr_idx[g.nbr_ptr[0]]], 6.0)
+        d = G.diffuse(g, f, 10)
+        self.assertLess(d.max(), 6.0)
+        self.assertGreater(np.count_nonzero(d > 1e-6), 20)
+
+    def test_gradient_of_linear_field(self):
+        g = self.g
+        f = g.xyz[:, 2] * g.radius_km  # height ∝ z → gradient points north near the equator
+        gr = G.gradient(g, f)
+        eq = np.abs(g.lat) < 20
+        _, n = east_north(g.xyz[eq])
+        cos_to_north = np.sum(gr[eq] * n, axis=1) / np.linalg.norm(gr[eq], axis=1)
+        self.assertGreater(np.median(cos_to_north), 0.99)
+        self.assertAlmostEqual(float(np.median(np.linalg.norm(gr[eq], axis=1))), 1.0, delta=0.1)
+
+    def test_distance_matches_great_circle(self):
+        g = self.g
+        i = g.cell_index(0.0, 0.0)
+        d = G.distance_to(g, np.arange(g.n) == i)
+        gc = g.radius_km * np.arccos(np.clip(g.xyz @ g.xyz[i], -1, 1))
+        far = gc > 3000
+        ratio = d[far] / gc[far]
+        self.assertTrue(np.all(ratio >= 0.999) and np.median(ratio) < 1.15)
+        self.assertTrue(np.all(np.isinf(G.distance_to(g, np.zeros(g.n, bool)))))
+
+    def test_nearest_source_labels(self):
+        g = self.g
+        a, b = g.cell_index(0.0, -90.0), g.cell_index(0.0, 90.0)
+        _, src = G.nearest_source(g, [a, b])
+        self.assertEqual(src[g.cell_index(0.0, -60.0)], a)
+        self.assertEqual(src[g.cell_index(0.0, 60.0)], b)
+
+    def test_priority_flood_fills_basin_and_drains(self):
+        g = self.g
+        ocean = g.lat < -30
+        z = np.where(ocean, -1000.0, 500.0 + 10 * g.lat)
+        pit = g.cell_index(40.0, 0.0)
+        z[pit] = -50.0  # land pit below sea level, not connected to ocean
+        zf = G.priority_flood(g, z, ocean)
+        self.assertGreater(zf[pit], z[pit])
+        recv, slope, dist = G.steepest_receivers(g, zf)
+        recv[ocean] = np.flatnonzero(ocean)
+        levels = G.receiver_levels(recv)
+        self.assertEqual(sum(len(l) for l in levels), g.n)
+        self.assertTrue(np.all(ocean[levels[0]]))
+
+    def test_priority_flood_needs_sink(self):
+        with self.assertRaisesRegex(ValueError, "no sink"):
+            G.priority_flood(self.g, np.ones(self.g.n), np.zeros(self.g.n, bool))
+
+    def test_accumulate_conserves(self):
+        g = self.g
+        ocean = g.lat < -30
+        z = np.where(ocean, -1000.0, 1000.0 + 20 * g.lat)
+        zf = G.priority_flood(g, z, ocean)
+        recv, _, _ = G.steepest_receivers(g, zf)
+        recv[ocean] = np.flatnonzero(ocean)
+        lv = G.receiver_levels(recv)
+        w = np.where(ocean, 0.0, 1.0)
+        acc = G.accumulate(recv, lv, w)
+        self.assertAlmostEqual(acc[lv[0]].sum(), w.sum())
+
+    def test_cycle_detected(self):
+        with self.assertRaisesRegex(ValueError, "cycle"):
+            G.receiver_levels(np.array([1, 0, 2]))
+
+    def test_components(self):
+        g = self.g
+        m = (np.abs(g.lat) < 10) & (np.abs(g.lon) < 20) | (np.abs(g.lat - 50) < 8) & (np.abs(g.lon) < 20)
+        lab = G.components(g, m)
+        self.assertEqual(len(np.unique(lab[m])), 2)
+        self.assertTrue(np.all(lab[~m] == -1))
+
+
+class SmoothKmTest(unittest.TestCase):
+    def test_constant_preserved_and_spike_decays(self):
+        g = small_grid(3)
+        np.testing.assert_allclose(G.smooth_km(g, np.full(g.n, 5.0), 1000.0), 5.0, rtol=1e-4)
+        i = g.cell_index(0.0, 0.0)
+        s = G.smooth_km(g, (np.arange(g.n) == i).astype(float), 1000.0)
+        d = g.radius_km * np.arccos(np.clip(g.xyz @ g.xyz[i], -1, 1))
+        near, far = s[(d > 800) & (d < 1200)].mean(), s[(d > 2800) & (d < 3200)].mean()
+        self.assertTrue(near > far > 0)
+
+    def test_resolution_independent(self):
+        vals = []
+        for res in (2, 3):
+            g = small_grid(res)
+            d = g.radius_km * np.arccos(np.clip(g.xyz @ g.xyz[g.cell_index(0.0, 0.0)], -1, 1))
+            s = G.smooth_km(g, (d < 2000).astype(float), 1500.0)
+            vals.append(s[g.cell_index(0.0, 30.0)])      # ~6700 km away
+        self.assertAlmostEqual(vals[0], vals[1], delta=0.25 * max(vals))
+        self.assertGreater(min(vals), 0.005)
+
+
+class OceanMaskTest(unittest.TestCase):
+    def test_inland_depression_is_not_ocean(self):
+        g = small_grid(3)
+        z = np.where(g.lat < 0, -3000.0, 500.0)
+        c = g.xyz[g.cell_index(40.0, 0.0)]
+        d = g.radius_km * np.arccos(np.clip(g.xyz @ c, -1, 1))
+        z[d < 600] = -50.0                                   # interior basin below sea level
+        ocean = G.ocean_mask(g, z, 1.0e6)
+        self.assertTrue(ocean[g.lat < -5].all())
+        self.assertFalse(ocean[d < 600].any())
+
+
+class SmoothKmRobustTest(unittest.TestCase):
+    def test_zero_and_tiny_fields(self):
+        g = small_grid(3)
+        np.testing.assert_array_equal(G.smooth_km(g, np.zeros(g.n), 25.0), 0.0)
+        f = np.where(g.lat > 0, 1e-9, 0.0)
+        s = G.smooth_km(g, f, 25.0)
+        self.assertTrue(np.all(np.isfinite(s)) and s.max() <= 1e-9 * (1 + 1e-6))
+
+    def test_femto_scale_field(self):
+        g = small_grid(3)
+        rng = np.random.default_rng(0)
+        f = np.where(rng.random(g.n) < 0.06, rng.random(g.n) * 8e-14, 0.0)
+        s = G.smooth_km(g, f, 25.0)
+        self.assertTrue(np.all(np.isfinite(s)))
+        self.assertAlmostEqual(float(s.sum() / f.sum()), 1.0, delta=0.05)
+
+
+def _flood_ref(g, z, sink_mask, eps=0.01):
+    """The original pure-Python priority flood (oracle for the compiled one)."""
+    import heapq
+    has_open = np.bincount(g.src, weights=(~sink_mask)[g.dst].astype(np.float64), minlength=g.n) > 0
+    zf, done = np.asarray(z, dtype=np.float64).tolist(), sink_mask.tolist()
+    ptr, idx = g.nbr_ptr.tolist(), g.nbr_idx.tolist()
+    heap = [(zf[i], i) for i in np.flatnonzero(sink_mask & has_open).tolist()]
+    heapq.heapify(heap)
+    while heap:
+        zc, c = heapq.heappop(heap)
+        for k in range(ptr[c], ptr[c + 1]):
+            n = idx[k]
+            if not done[n]:
+                done[n] = True
+                zf[n] = max(zf[n], zc + eps)
+                heapq.heappush(heap, (zf[n], n))
+    return np.array(zf)
+
+
+def _steepest_ref(g, z):
+    """The original lexsort version (oracle)."""
+    slope = (z[g.src] - z[g.dst]) / g.edge_km
+    first = np.lexsort((-slope, g.src))[g.nbr_ptr[:-1]]
+    s = slope[first]
+    down = s > 0
+    return (np.where(down, g.dst[first], np.arange(g.n)), np.where(down, s, 0.0),
+            np.where(down, g.edge_km[first], np.inf))
+
+
+def _accumulate_ref(recv, levels, w):
+    acc = np.asarray(w, dtype=np.float64).copy()
+    for lv in reversed(levels[1:]):
+        acc += np.bincount(recv[lv], weights=acc[lv], minlength=len(acc))
+    return acc
+
+
+class FastPathsTest(unittest.TestCase):
+    """The speed-ups (numba flood, sort-free receivers, per-level accumulate) give bit-identical results."""
+
+    def fields(self, g):
+        rng = np.random.default_rng(7)
+        rough = rng.normal(0, 300, g.n)
+        flat = np.round(rng.normal(0, 2, g.n))            # many exact ties: tie order matters
+        pits = np.where(rng.random(g.n) < 0.2, -50.0, rough)
+        return {"rough": rough, "flat": flat, "pits": pits}
+
+    def test_flood_matches_reference(self):
+        g = small_grid(3)
+        for name, z in self.fields(g).items():
+            sink = z < np.quantile(z, 0.1)
+            for eps in (0.01, 0.0):
+                with self.subTest(name=name, eps=eps):
+                    want = _flood_ref(g, z, sink, eps)
+                    got = G.priority_flood(g, z, sink, eps)
+                    self.assertTrue(np.array_equal(got, want))
+                    py = G._flood_py(z, sink, g.nbr_ptr, g.nbr_idx,
+                                     np.flatnonzero(sink & (np.bincount(g.src, weights=(~sink)[g.dst].astype(float),
+                                                                       minlength=g.n) > 0)), eps)
+                    self.assertTrue(np.array_equal(py, want))
+
+    def test_steepest_receivers_match_reference(self):
+        g = small_grid(3)
+        for name, z in self.fields(g).items():
+            with self.subTest(name=name):
+                for got, want in zip(G.steepest_receivers(g, z), _steepest_ref(g, z)):
+                    self.assertTrue(np.array_equal(got, want))
+
+    def test_accumulate_matches_reference(self):
+        g = small_grid(3)
+        rng = np.random.default_rng(3)
+        for name, z in self.fields(g).items():
+            with self.subTest(name=name):
+                zf = G.priority_flood(g, z, z < np.quantile(z, 0.1))
+                recv, _, _ = G.steepest_receivers(g, zf)
+                lv = G.receiver_levels(recv)
+                w = rng.random(g.n) * 1e3 * np.where(rng.random(g.n) < 0.1, -0.0, 1.0)   # with negative zeros
+                got, want = G.accumulate(recv, lv, w), _accumulate_ref(recv, lv, w)
+                self.assertTrue(np.array_equal(got, want))
+                self.assertTrue(np.array_equal(np.signbit(got), np.signbit(want)))
+
+
+class SweepTest(unittest.TestCase):
+    def test_sweep_matches_full_bincount(self):
+        from mapgen import hydrology as HY
+
+        def sweep_ref(recv, levels, water, outlets, cap):
+            acc = np.asarray(water, dtype=np.float64).copy()
+            loss = np.zeros(len(acc))
+            is_out = np.zeros(len(acc), bool)
+            is_out[outlets] = True
+            cap_cell = np.zeros(len(acc))
+            cap_cell[outlets] = cap
+            for lv in reversed(levels[1:]):
+                push = acc[lv].copy()
+                o = is_out[lv]
+                if o.any():
+                    cells = lv[o]
+                    lost = np.minimum(acc[cells], cap_cell[cells])
+                    loss[cells] = lost
+                    push[o] = acc[cells] - lost
+                acc += np.bincount(recv[lv], weights=push, minlength=len(acc))
+            return acc, loss
+
+        g = small_grid(3)
+        rng = np.random.default_rng(11)
+        z = rng.normal(0, 300, g.n)
+        zf = G.priority_flood(g, z, z < np.quantile(z, 0.1))
+        recv, _, _ = G.steepest_receivers(g, zf)
+        lv = G.receiver_levels(recv)
+        water = rng.random(g.n) * np.where(rng.random(g.n) < 0.1, -0.0, 1.0)     # with negative zeros
+        outlets = rng.choice(g.n, 200, replace=False)
+        cap = rng.random(200) * 2
+        for got, want in zip(HY._sweep(recv, lv, water, outlets, cap), sweep_ref(recv, lv, water, outlets, cap)):
+            self.assertTrue(np.array_equal(got, want))
+            self.assertTrue(np.array_equal(np.signbit(got), np.signbit(want)))
+
+
+class LeavesTest(unittest.TestCase):
+    def test_leaves_all_matches_the_walk(self):
+        from mapgen import hydrology as HY
+        g = small_grid(3)
+        rng = np.random.default_rng(4)
+        for seed in range(3):
+            z = rng.normal(0, 300, g.n)
+            ocean = z < np.quantile(z, 0.2)
+            zf = G.priority_flood(g, z, ocean)
+            lab = G.components(g, ~ocean & (zf - z > 1.0))
+            recv, _, _ = G.steepest_receivers(g, zf)
+            recv = np.where(ocean, np.arange(g.n), recv)
+            lv = G.receiver_levels(recv)
+            xs = np.flatnonzero(lab >= 0)
+            for limit in (100000, 3):
+                want = np.array([HY._leaves(recv, lab, x, limit) for x in xs])
+                got = HY._leaves_all(recv, lab, lv, limit)[xs]
+                self.assertGreater(want.sum(), 0)
+                self.assertTrue(np.array_equal(got, want), (seed, limit))
+
+
+class BicgstabJacobiTest(unittest.TestCase):
+    """The fused solver walks scipy's iterates exactly: same answers bit for bit, same exit codes."""
+    def systems(self):
+        from scipy import sparse
+        rng = np.random.default_rng(0)
+        for n in (2000, 9000):
+            i = np.repeat(np.arange(n), 6)
+            j = (i + rng.integers(-50, 50, len(i))) % n
+            L = sparse.csr_matrix((rng.random(len(i)), (i, j)), shape=(n, n))
+            S = L + L.T
+            yield (sparse.diags(np.asarray(S.sum(1)).ravel()) - S).tocsr() * 40 + sparse.identity(n, format="csr")
+            yield (sparse.identity(n, format="csr") * (1 + np.asarray(L.sum(1)).ravel().max() * 0.6) - L).tocsr()
+
+    def setUp(self):
+        self.min_n = G.JIT_MIN_N
+        G.JIT_MIN_N = 0                                      # the compiled path even on small test systems
+
+    def tearDown(self):
+        G.JIT_MIN_N = self.min_n
+
+    def test_matches_scipy_bicgstab(self):
+        from scipy.sparse import linalg as splinalg
+        rng = np.random.default_rng(1)
+        for k, A in enumerate(self.systems()):
+            for rtol, maxiter in ((1e-6, 5000), (1e-9, 5000), (1e-12, 7)):
+                b = rng.normal(size=A.shape[0])
+                inv = 1.0 / A.diagonal()
+                M = splinalg.LinearOperator(A.shape, matvec=lambda x: inv * x)
+                want = splinalg.bicgstab(A, b, x0=b * 0.5, rtol=rtol, maxiter=maxiter, M=M)
+                got = G.bicgstab_jacobi(A, b, b * 0.5, inv, rtol, maxiter)
+                with self.subTest(k=k, rtol=rtol, maxiter=maxiter):
+                    self.assertEqual(got[1], want[1])
+                    self.assertTrue(np.array_equal(got[0], want[0]))
+
+    def test_zero_right_hand_side_and_zero_start(self):
+        A = next(self.systems())
+        inv = 1.0 / A.diagonal()
+        x, info = G.bicgstab_jacobi(A, np.zeros(A.shape[0]), np.zeros(A.shape[0]), inv, 1e-6, 100)
+        self.assertEqual(info, 0)
+        self.assertFalse(x.any())
+
+
+class PmapTest(unittest.TestCase):
+    def test_order_and_same_floats_as_serial(self):
+        import os
+        from unittest import mock
+        g = small_grid(2)
+        fields = [np.random.default_rng(i).normal(size=g.n) for i in range(4)]
+        serial = [G.smooth_km(g, f, 900.0) for f in fields]
+        with mock.patch.dict(os.environ, {"WORLDGEN_THREADS": "4"}):
+            self.assertEqual(G.workers(), 4)
+            got = G.pmap(lambda f: G.smooth_km(g, f, 900.0), fields)
+        for a, b in zip(serial, got):
+            self.assertTrue(np.array_equal(a, b))
+        with mock.patch.dict(os.environ, {"WORLDGEN_THREADS": "3"}):
+            self.assertEqual(G.pmap(lambda k: k * k, range(9)), [k * k for k in range(9)])
+        with mock.patch.dict(os.environ, {"WORLDGEN_THREADS": "x"}):
+            self.assertEqual(G.workers(), 3)
+
+
+class ComponentsTest(unittest.TestCase):
+    def test_same_labels_as_scipy(self):
+        from scipy import sparse
+        from scipy.sparse import csgraph
+        g = small_grid(3)
+        rng = np.random.default_rng(11)
+
+        def scipy_labels(mask):                                # the previous implementation, as the oracle
+            e = mask[g.src] & mask[g.dst]
+            m = sparse.csr_matrix((np.ones(int(e.sum())), (g.src[e], g.dst[e])), shape=(g.n, g.n))
+            _, lab = csgraph.connected_components(m, directed=False)
+            return np.where(mask, lab, -1)
+        for p in (0.0, 0.2, 0.45, 0.6, 0.9, 1.0):
+            for _ in range(3):
+                mask = rng.random(g.n) < p
+                want, got = scipy_labels(mask), G.components(g, mask)
+                self.assertEqual(got.dtype, want.dtype)
+                self.assertTrue(np.array_equal(got, want), p)
+        smooth = G.smooth_km(g, rng.normal(size=g.n), 2000.0) > 0      # big blobs, many cells each
+        self.assertTrue(np.array_equal(G.components(g, smooth), scipy_labels(smooth)))
+        self.assertGreater(len(np.unique(G.components(g, smooth))), 2)
diff --git a/tests/test_grid.py b/tests/test_grid.py
new file mode 100644
index 0000000..a6b0848
--- /dev/null
+++ b/tests/test_grid.py
@@ -0,0 +1,83 @@
+import unittest
+
+import numpy as np
+
+from mapgen.grid import Grid
+from tests.helpers import small_grid
+
+
+class GridTest(unittest.TestCase):
+    def test_counts_and_pentagons(self):
+        g = small_grid(2)
+        self.assertEqual(g.n, 5882)
+        self.assertTrue(np.all(np.diff(g.ids.astype(np.int64)) > 0))
+        self.assertEqual(int(np.sum(g.counts == 5)), 12)
+        self.assertEqual(int(np.sum(g.counts == 6)), g.n - 12)
+
+    def test_neighbours_symmetric(self):
+        g = small_grid(2)
+        pairs = set(zip(g.src.tolist(), g.dst.tolist()))
+        self.assertTrue(all((b, a) in pairs for a, b in pairs))
+
+    def test_area_sums_to_sphere(self):
+        g = small_grid(2)
+        self.assertAlmostEqual(g.area_km2.sum() / (4 * np.pi * 12742.0**2), 1.0, places=6)
+
+    def test_spacing_and_tangents(self):
+        g = small_grid(2)
+        self.assertTrue(500 < g.spacing_km < 800)  # res 2 ≈ 632 km
+        np.testing.assert_allclose(np.linalg.norm(g.edge_tangents, axis=1), 1.0)
+        np.testing.assert_allclose(np.sum(g.edge_tangents * g.xyz[g.src], axis=1), 0.0, atol=1e-12)
+
+    def test_roundtrip_arrays_and_lookup(self):
+        g = small_grid(2)
+        h = Grid.from_arrays(g.to_arrays(), g.res, g.radius_km)
+        np.testing.assert_array_equal(h.nbr_idx, g.nbr_idx)
+        i = g.cell_index(10.0, 20.0)
+        self.assertLess(abs(g.lat[i] - 10.0), 5.0)
+
+
+class CellParentsTest(unittest.TestCase):
+    def test_same_as_h3_cell_to_parent(self):
+        import h3.api.basic_int as h3
+        import numpy as np
+        from mapgen.grid import cell_parents
+        rng = np.random.default_rng(5)
+        for res in (1, 3, 6, 9):
+            cells = [h3.latlng_to_cell(float(la), float(lo), res + int(k))
+                     for la, lo, k in zip(rng.uniform(-90, 90, 300), rng.uniform(-180, 180, 300), rng.integers(0, 4, 300))]
+            pent = h3.get_pentagons(max(res - 2, 0))[:2] + [h3.latlng_to_cell(10.0, 20.0, max(res - 2, 0))]
+            cells += [c for p in pent for c in h3.cell_to_children(p, res)]                      # incl. pentagons
+            for pr in range(0, res + 1):
+                want = np.array([h3.cell_to_parent(c, pr) for c in cells], np.uint64)
+                self.assertTrue(np.array_equal(cell_parents(cells, pr), want), (res, pr))
+        with self.assertRaises(ValueError):
+            cell_parents([h3.latlng_to_cell(0.0, 0.0, 2)], 3)
+
+
+class EdgeBlocksTest(unittest.TestCase):
+    def test_blocked_edge_fields_same_as_whole(self):
+        from unittest import mock
+        from mapgen import grid as GR
+        from mapgen import plates as PL
+        from mapgen.sphere import tangent_dir
+        g = small_grid(3)
+        xyz, src, dst = g.xyz, g.src, g.dst
+        want_t = tangent_dir(xyz[src], xyz[dst])                      # the whole-array formulas, written out
+        want_km = g.radius_km * np.arccos(np.clip(np.sum(xyz[src] * xyz[dst], axis=1), -1.0, 1.0))
+        vel = np.random.default_rng(2).normal(size=(g.n, 3))
+        rel = vel[dst] - vel[src]
+        along = np.sum(rel * want_t, axis=1)
+        want_tang = np.linalg.norm(rel - along[:, None] * want_t, axis=1)
+        with mock.patch.object(GR, "EDGE_BLOCK", 997):                 # many uneven blocks
+            g2 = small_grid(3)
+            self.assertTrue(np.array_equal(g2.edge_tangents, want_t))
+            self.assertTrue(np.array_equal(g2.edge_km, want_km))
+            conv, tang = PL.edge_convergence(g2, vel)
+        self.assertTrue(np.array_equal(conv, -along))
+        self.assertTrue(np.array_equal(tang, want_tang))
+        self.assertGreater(len(GR.edge_blocks(len(dst), 997)), 3)
+        with mock.patch.object(GR, "EDGE_BLOCK", 997):
+            for idx in (src, dst):
+                self.assertTrue(np.array_equal(GR.rowdot_at(vel, idx, want_t), np.sum(vel[idx] * want_t, axis=1)))
+                self.assertTrue(np.array_equal(GR.rowdot_at(vel, idx, want_t), np.sum(want_t * vel[idx], axis=1)))
diff --git a/tests/test_hydrology.py b/tests/test_hydrology.py
new file mode 100644
index 0000000..3e9314a
--- /dev/null
+++ b/tests/test_hydrology.py
@@ -0,0 +1,133 @@
+import unittest
+
+import numpy as np
+
+from mapgen import graph as G
+from mapgen import hydrology as HY
+from tests.helpers import make_ctx
+
+
+def basin_ctx(p_basin, pet_basin):
+    ctx = make_ctx(3)
+    g = ctx.grid
+    ocean = g.lat < -20
+    z = np.where(ocean, -3000.0, 200.0 + 30.0 * (g.lat + 20))          # land rises northward
+    c = g.xyz[g.cell_index(40.0, 0.0)]
+    d = g.radius_km * np.arccos(np.clip(g.xyz @ c, -1, 1))
+    z = np.where(d < 800, z - 1500 * (1 - d / 800), z)                  # closed basin
+    ctx.data.update({"elevation_eroded_m": z.astype(np.float32), "P_ann": np.full(g.n, p_basin),
+                     "PET": np.full(g.n, pet_basin)})
+    return ctx, d
+
+
+class HydrologyTest(unittest.TestCase):
+    def test_strahler(self):
+        recv = np.array([2, 2, 6, 5, 5, 6, 6])
+        levels = G.receiver_levels(recv)
+        o = HY.strahler(recv, levels, np.ones(7, bool))
+        np.testing.assert_array_equal(o, [1, 1, 2, 1, 1, 2, 3])
+
+    def test_arid_basin_is_endorheic_with_salt(self):
+        ctx, d = basin_ctx(100.0, 1500.0)
+        out = HY.run(ctx)
+        inside = d < 300
+        self.assertTrue(out["endorheic"][inside].all())
+        self.assertTrue(out["salt_flat"][inside].any() or out["lake"][inside].any())
+
+    def test_humid_basin_overflows(self):
+        ctx, d = basin_ctx(2500.0, 400.0)
+        out = HY.run(ctx)
+        inside = d < 300
+        self.assertTrue(out["lake"][inside].all())
+        self.assertFalse(out["endorheic"][inside].any())
+
+    def test_mass_balance_and_rivers(self):
+        ctx, _ = basin_ctx(2500.0, 400.0)
+        out = HY.run(ctx)
+        g = ctx.grid
+        roots = out["recv"] == np.arange(g.n)
+        total = np.sum(out["runoff_mm"] * g.area_km2) * 1e-6
+        self.assertAlmostEqual(out["discharge_km3_yr"][roots].sum(), total, delta=total * 1e-6)
+        self.assertTrue(out["river"].any())
+        self.assertTrue(np.all(out["strahler"][out["river"]] >= 1))
+        land = ctx.data["elevation_eroded_m"] > 0
+        nonendo = land & ~out["endorheic"] & ~roots
+        self.assertTrue(np.all(out["z_filled_m"][out["recv"][nonendo]] < out["z_filled_m"][nonendo]))
+
+
+class InlandDepressionTest(unittest.TestCase):
+    def test_humid_depression_below_sea_level_becomes_lake(self):
+        ctx, d = basin_ctx(2500.0, 400.0)
+        g = ctx.grid
+        z = ctx.data["elevation_eroded_m"].astype(np.float64)
+        z[d < 300] = -40.0
+        ctx.data["elevation_eroded_m"] = z.astype(np.float32)
+        ctx.data["ocean"] = G.ocean_mask(g, z, 1.0e6)
+        out = HY.run(ctx)
+        self.assertTrue(out["lake"][d < 300].all())
+
+
+def two_pit_ctx(p, pet):
+    ctx = make_ctx(3)
+    g = ctx.grid
+    ocean = g.lat < -20
+    z = np.where(ocean, -3000.0, 200.0 + 30.0 * (g.lat + 20))
+    for lon in (-3.5, 3.5):
+        c = g.xyz[g.cell_index(40.0, lon)]
+        d = g.radius_km * np.arccos(np.clip(g.xyz @ c, -1, 1))
+        z = np.where(d < 900, z - (1500 + 200 * (lon > 0)) * (1 - d / 900), z)   # two pits, one depression
+    ctx.data.update({"elevation_eroded_m": z.astype(np.float32), "P_ann": np.full(g.n, p), "PET": np.full(g.n, pet)})
+    return ctx, z
+
+
+class TerminalTest(unittest.TestCase):
+    def test_every_land_root_is_lake_or_salt(self):
+        ctx, z = two_pit_ctx(700.0, 1400.0)
+        out = HY.run(ctx)
+        roots = (out["recv"] == np.arange(len(z))) & (z > 0) | (out["recv"] == np.arange(len(z))) & out["endorheic"]
+        self.assertTrue(out["endorheic"].any())
+        self.assertTrue(np.all(out["lake"][roots] | out["salt_flat"][roots]))
+        self.assertEqual(int((roots & out["endorheic"]).sum()), 1)          # one terminal for the depression
+
+    def test_mass_balance_with_lake_evaporation(self):
+        ctx, z = two_pit_ctx(1000.0, 1200.0)
+        out = HY.run(ctx)
+        g = ctx.grid
+        roots = out["recv"] == np.arange(g.n)
+        total = np.sum(out["runoff_mm"] * g.area_km2) * 1e-6
+        self.assertGreater(out["lake_loss_km3_yr"].sum(), 0)
+        self.assertAlmostEqual(out["discharge_km3_yr"][roots].sum() + out["lake_loss_km3_yr"].sum(), total,
+                               delta=total * 1e-6)
+
+
+class LakeBedTest(unittest.TestCase):
+    def test_lakes_carved_below_level_and_rerun_is_stable(self):
+        ctx, d = basin_ctx(2500.0, 400.0)
+        z0 = np.asarray(ctx.data["elevation_eroded_m"], dtype=np.float64)
+        out = HY.run(ctx)
+        lake, lev, z = out["lake"], out["lake_level_m"], out["elevation_eroded_m"].astype(np.float64)
+        self.assertTrue(lake.any())
+        self.assertTrue(np.all(np.isnan(lev[~lake])))
+        self.assertTrue(np.all(z[lake] <= lev[lake] - HY.DEFAULTS["lake_min_m"] + 1e-3))
+        np.testing.assert_array_equal(z[~lake], z0[~lake].astype(np.float32))   # only lake beds change
+        ctx.data.update({"elevation_eroded_m": out["elevation_eroded_m"]})     # again (eras rerun hydrology)
+        again = HY.run(ctx)
+        np.testing.assert_array_equal(again["lake"], lake)
+        np.testing.assert_allclose(again["lake_level_m"][lake], lev[lake])
+        np.testing.assert_array_equal(again["elevation_eroded_m"], out["elevation_eroded_m"])
+
+    def test_bigger_lakes_deeper(self):
+        from tests.helpers import make_ctx as mk
+        g = mk(3).grid
+        P = {**HY.DEFAULTS, "lake_depth_exp": 0.3}
+        z = np.zeros(g.n)
+        lake = np.zeros(g.n, bool)
+        small, big = g.cell_index(10.0, 0.0), np.flatnonzero(g.xyz @ g.xyz[g.cell_index(-30.0, 90.0)] > 0.97)
+        lake[small] = lake[big] = True
+        lid = np.where(lake, 0, -1)
+        lid[big] = 1
+        bed = HY.lake_beds(g, z, np.where(lake, 0.0, np.nan), lake, lid, np.zeros(g.n, bool), np.zeros(g.n, bool),
+                           np.full(g.n, 1000.0), 1, P)
+        self.assertLess(bed[big].min(), -P["lake_min_m"])
+        self.assertLessEqual(bed[small], -P["lake_min_m"])
+        self.assertTrue(np.all(bed[~lake] == 0))
diff --git a/tests/test_ice_fields.py b/tests/test_ice_fields.py
new file mode 100644
index 0000000..08fda18
--- /dev/null
+++ b/tests/test_ice_fields.py
@@ -0,0 +1,98 @@
+import unittest
+
+import numpy as np
+
+from mapgen import fields as FI
+from mapgen import ice as IC
+from tests.helpers import make_ctx
+
+
+class IceTest(unittest.TestCase):
+    def test_sheet_glacier_sea_ice(self):
+        ctx = make_ctx(3)
+        g = ctx.grid
+        polar_land = g.lat > 70
+        z = np.where(polar_land, 500.0, -3000.0)
+        peak = g.cell_index(0.0, 0.0)
+        z[peak] = 6000.0
+        t_mean = np.where(g.lat > 60, -25.0, 15.0)
+        t_mean[peak] = -9.0          # too warm for a sheet, summer below zero → glacier
+        t_sum = t_mean + 8
+        t_win = t_mean - 8
+        ctx.data.update({"elevation_eroded_m": z, "T_mean": t_mean, "T_jun": t_sum, "T_dec": t_win,
+                         "P_ann": np.full(g.n, 300.0)})
+        out = IC.run(ctx)
+        self.assertEqual(out["ice"][g.cell_index(85.0, 0.0)], IC.ICE_SHEET)
+        self.assertEqual(out["ice"][peak], IC.ICE_GLACIER)
+        self.assertEqual(out["ice"][g.cell_index(65.0, 0.0)], IC.ICE_SEA_PERENNIAL)
+        self.assertEqual(out["ice"][g.cell_index(0.0, 150.0)], IC.ICE_NONE)
+        i = g.cell_index(85.0, 0.0)
+        self.assertGreater(out["z_surface_m"][i], z[i] + 200)
+
+
+class FieldsTest(unittest.TestCase):
+    def test_o2_and_gravity(self):
+        ctx = make_ctx(1)
+        n = ctx.grid.n
+        zs = np.zeros(n)
+        zs[1] = 8000.0
+        m_o2 = np.zeros(n)
+        m_o2[2] = 1.0
+        m_g = np.zeros(n)
+        m_g[3] = -1.0
+        ctx.data.update({"z_surface_m": zs, "m_o2_zones": m_o2, "m_gravity_zones": m_g})
+        out = FI.run(ctx)
+        self.assertAlmostEqual(out["po2_bar"][0], 0.21)
+        self.assertAlmostEqual(out["po2_bar"][1], 0.21 * np.exp(-1), places=6)
+        self.assertAlmostEqual(out["po2_bar"][2], 0.21 * 1.5)
+        self.assertAlmostEqual(out["gravity_g"][0], 1.05)
+        self.assertAlmostEqual(out["gravity_g"][3], 1.05 * 0.3)
+
+    def test_pressure_o2_fraction_and_fire(self):
+        ctx = make_ctx(1)
+        n = ctx.grid.n
+        zs = np.zeros(n)
+        zs[1] = 8000.0
+        m_o2 = np.zeros(n)
+        m_o2[2] = 1.0
+        ctx.data.update({"z_surface_m": zs, "m_o2_zones": m_o2, "m_gravity_zones": np.zeros(n)})
+        out = FI.run(ctx)
+        self.assertAlmostEqual(out["pressure_bar"][0], 1.0)
+        self.assertAlmostEqual(out["pressure_bar"][1], np.exp(-1))
+        self.assertAlmostEqual(out["o2_fraction"][2], 0.21 * 1.5)
+        np.testing.assert_allclose(out["po2_bar"], out["o2_fraction"] * out["pressure_bar"])
+        self.assertTrue(np.all(out["fire_reactivity"] == 1.0))
+        np.testing.assert_allclose(out["plant_height_x"], ctx.cfg["planet"]["gravity_g"] / out["gravity_g"])
+        np.testing.assert_allclose(FI.pressure(2.0, np.array([-100.0, 8000.0]), 8000.0), [2.0, 2.0 * np.exp(-1)])
+
+class IceOceanMaskTest(unittest.TestCase):
+    def test_cold_inland_basin_is_not_sea_ice(self):
+        ctx = make_ctx(3)
+        g = ctx.grid
+        land = g.lat > 50
+        z = np.where(land, 300.0, -3000.0)
+        basin = (g.lat > 60) & (g.lat < 68) & (np.abs(g.lon) < 20)
+        z[basin] = -20.0
+        t = np.where(g.lat > 50, -5.0, 10.0)
+        ctx.data.update({"elevation_eroded_m": z, "ocean": ~land, "T_mean": t, "T_jun": t + 8, "T_dec": t - 8,
+                         "P_ann": np.full(g.n, 300.0)})
+        out = IC.run(ctx)
+        self.assertFalse(np.isin(out["ice"][basin], [IC.ICE_SEA_SEASONAL, IC.ICE_SEA_PERENNIAL]).any())
+
+
+class SummerMeltRuleTest(unittest.TestCase):
+    def test_ice_sheets_need_summers_below_freezing(self):
+        ctx = make_ctx(3)
+        g = ctx.grid
+        land = g.lat > 55
+        z = np.where(land, 300.0, -3000.0)
+        t_mean = np.where(g.lat > 75, -25.0, -12.0)           # 55–75°: cold on average, but summers melt
+        t_sum = np.where(g.lat > 75, -8.0, 4.0)
+        peak = g.cell_index(0.0, 0.0)
+        z[peak], t_mean[peak], t_sum[peak] = 6000.0, -9.0, -1.0
+        ctx.data.update({"elevation_eroded_m": z, "ocean": ~(land | (np.arange(g.n) == peak)), "T_mean": t_mean,
+                         "T_jun": t_sum, "T_dec": t_mean - 10, "P_ann": np.full(g.n, 300.0)})
+        out = IC.run(ctx)
+        self.assertEqual(out["ice"][g.cell_index(65.0, 0.0)], IC.ICE_NONE)        # tundra, not ice
+        self.assertEqual(out["ice"][g.cell_index(85.0, 0.0)], IC.ICE_SHEET)
+        self.assertEqual(out["ice"][peak], IC.ICE_GLACIER)                          # small cold spot = glacier
diff --git a/tests/test_low_memory.py b/tests/test_low_memory.py
new file mode 100644
index 0000000..c5ed333
--- /dev/null
+++ b/tests/test_low_memory.py
@@ -0,0 +1,83 @@
+import shutil
+import tempfile
+import tracemalloc
+import unittest
+from pathlib import Path
+
+import numpy as np
+
+from mapgen import pipeline as P
+from mapgen import render as RN
+from mapgen.testing import small_world
+
+
+def _files(root: Path) -> dict:
+    return {p.relative_to(root).as_posix(): p.read_bytes()
+            for p in sorted(root.rglob("*")) if p.is_file() and "cache" not in p.parts}
+
+
+class LowMemoryTest(unittest.TestCase):
+    def test_low_memory_build_writes_identical_files(self):
+        # oracle: the normal-mode build of the same world; every output byte must match
+        outs = []
+        for low in (False, True):
+            tmp = Path(tempfile.mkdtemp())
+            try:
+                small_world(tmp)
+                P.build(tmp, 2, log=lambda m: None, low_memory=low)
+                outs.append({**_files(tmp / "out"), **_files(tmp / "previews")})
+            finally:
+                shutil.rmtree(tmp)
+        self.assertEqual(sorted(outs[0]), sorted(outs[1]))
+        for k in outs[0]:
+            self.assertEqual(outs[0][k], outs[1][k], k)
+
+    def test_chunked_relief_matches_and_uses_less_memory(self):
+        rng = np.random.default_rng(0)
+        H, W = 512, 1024
+        z = rng.normal(0, 2000, (H, W))
+        hs = rng.random((H, W))
+        zone = rng.integers(0, 38, (H, W)).astype(np.int16)
+        ground = rng.integers(0, 9, (H, W)).astype(np.int8)
+        ice = rng.integers(0, 5, (H, W)).astype(np.int8)
+        lake = rng.random((H, W)) < 0.05
+        peaks = []
+        outs = []
+        for low in (False, True):
+            tracemalloc.start()
+            outs.append(RN.relief_rgb(z, hs, zone, ground, ice, lake, low_memory=low))
+            peaks.append(tracemalloc.get_traced_memory()[1])
+            tracemalloc.stop()
+        np.testing.assert_array_equal(outs[0], outs[1])
+        self.assertLess(peaks[1], 0.4 * peaks[0])                     # measured ≈ 0.28
+
+
+class ProjectionMemoryTest(unittest.TestCase):
+    def test_sampling_converts_one_channel_at_a_time(self):
+        from mapgen import projections as PJ
+        rng = np.random.default_rng(1)
+        img = rng.integers(0, 256, (512, 1024, 3), dtype=np.uint8)
+        lat = rng.uniform(-90, 90, (64, 64))
+        lon = rng.uniform(-180, 180, (64, 64))
+        whole = img.astype(np.float64)                        # oracle: the whole-image float conversion
+        want = np.stack([PJ.sample_equirect(whole[..., k], lat, lon) for k in range(3)], axis=-1)
+        del whole
+        tracemalloc.start()
+        got = PJ._sample_rgb(img, lat, lon)
+        peak = tracemalloc.get_traced_memory()[1]
+        tracemalloc.stop()
+        np.testing.assert_array_equal(got, want)
+        self.assertLess(peak, 0.5 * img.size * 8)               # no full-size float64 copy
+
+
+class HillshadeTest(unittest.TestCase):
+    def test_chunked_hillshade_matches_and_uses_less_memory(self):
+        z = np.random.default_rng(2).normal(0, 1500, (700, 1400))        # 700 rows: chunk edges fall mid-raster
+        peaks, outs = [], []
+        for low in (False, True):
+            tracemalloc.start()
+            outs.append(RN.hillshade(z, 6371.0, low_memory=low))
+            peaks.append(tracemalloc.get_traced_memory()[1])
+            tracemalloc.stop()
+        np.testing.assert_array_equal(outs[0], outs[1])
+        self.assertLess(peaks[1], 0.5 * peaks[0])
diff --git a/tests/test_ocean.py b/tests/test_ocean.py
new file mode 100644
index 0000000..81ad8e1
--- /dev/null
+++ b/tests/test_ocean.py
@@ -0,0 +1,258 @@
+import unittest
+
+import numpy as np
+
+from mapgen import ocean as OC
+from mapgen.sphere import east_north
+from tests.helpers import small_grid
+
+DAY = 31.149
+KM_PER_DEG = 12742.0 * np.pi / 180.0
+
+
+def basin(g):
+    """Ocean box 15–45° N between two meridional coasts at ±60° lon (everything else land)."""
+    return (g.lat > 15) & (g.lat < 45) & (np.abs(g.lon) < 60)
+
+
+def gyre_wind(g):
+    """Trades (from the east) at 15° |lat|, westerlies at 45° |lat|: −8 … +8 m/s east."""
+    e, _ = east_north(g.xyz)
+    U = -8.0 * np.cos(np.radians((np.abs(g.lat) - 15.0) / 30.0 * 180.0))
+    return U[:, None] * e
+
+
+def params(**over):
+    P = dict(OC.DEFAULTS)
+    P["friction_days"] = 2.0          # r/β ≈ 750 km: resolved on the res-3 test grid (≈240 km spacing)
+    P.update(over)
+    return P
+
+
+def northward(g, u):
+    _, n = east_north(g.xyz)
+    return np.sum(u * n, axis=1)
+
+
+class GyreTest(unittest.TestCase):
+    @classmethod
+    def setUpClass(cls):
+        cls.g = small_grid(3)
+        cls.ocean = basin(cls.g)
+        cls.u = OC.currents(cls.g, cls.ocean, gyre_wind(cls.g), params(), DAY)
+
+    def test_clockwise_with_western_intensification(self):
+        g, v = self.g, northward(self.g, self.u)
+        band = self.ocean & (g.lat > 25) & (g.lat < 35)
+        west_km = (g.lon + 60.0) * KM_PER_DEG * np.cos(np.radians(30.0))
+        east_km = (60.0 - g.lon) * KM_PER_DEG * np.cos(np.radians(30.0))
+        west, east = band & (west_km < 1500), band & (east_km < 1500)
+        self.assertGreater(v[west].max(), 0.0)                          # northward along the west coast
+        self.assertLess(v[east].min(), 0.0)                             # southward in the east
+        self.assertGreater(v[west].max(), 3.0 * -v[east].min())         # the western boundary current is the fast one
+
+    def test_land_is_still(self):
+        g = self.g
+        psi = OC.streamfunction(g, self.ocean, gyre_wind(g), params(), DAY)
+        sp = np.linalg.norm(OC.velocity(g, psi), axis=1)
+        inland = ~self.ocean
+        for _ in range(2):                                              # two cells away from any sea
+            inland = inland & ~np.bincount(g.src, weights=(~inland[g.dst]).astype(float), minlength=g.n).astype(bool)
+        self.assertLess(sp[inland].max(), 0.01 * sp[self.ocean].max())
+
+    def test_zero_on_land_in_output(self):
+        self.assertTrue(np.all(self.u[~self.ocean] == 0.0))
+
+    def test_southern_hemisphere_is_anticlockwise(self):
+        g = self.g
+        ocean = (g.lat < -15) & (g.lat > -45) & (np.abs(g.lon) < 60)
+        v = northward(g, OC.currents(g, ocean, gyre_wind(g), params(), DAY))
+        band = ocean & (g.lat < -25) & (g.lat > -35)
+        west = band & (g.lon < -45)
+        self.assertLess(v[west].min(), 0.0)                             # southward along the west coast
+        self.assertGreater(-v[west].min(), 3.0 * max(v[band & (g.lon > 45)].max(), 1e-9))
+
+
+class ChannelTest(unittest.TestCase):
+    def test_westerlies_drive_eastward_flow_downwind_without_rotation(self):
+        g = small_grid(3)
+        ocean = (g.lat > 30) & (g.lat < 60)
+        e, _ = east_north(g.xyz)
+        wind = 8.0 * e
+        mid = (g.lat > 40) & (g.lat < 50)
+        u = OC.currents(g, ocean, wind, params(), 1e7)                  # ~no rotation: flow is downwind
+        ue, un = np.sum(u * e, axis=1), northward(g, u)
+        self.assertGreater(ue[mid].mean(), 0.0)
+        self.assertLess(np.abs(un[mid]).mean(), 0.05 * ue[mid].mean())
+        u = OC.currents(g, ocean, wind, params(), DAY)
+        self.assertGreater(np.sum(u * e, axis=1)[mid].mean(), 0.0)
+
+
+class DayLengthTest(unittest.TestCase):
+    def test_slower_rotation_widens_boundary_current(self):
+        g = small_grid(3)
+        ocean = basin(g)
+        band = ocean & (g.lat > 27) & (g.lat < 33)
+        widths = []
+        for day in (DAY, 4 * DAY):
+            v = northward(g, OC.currents(g, ocean, gyre_wind(g), params(friction_days=4.0), day))
+            half = band & (v > 0.5 * v[band].max())
+            widths.append(g.lon[half].max() - g.lon[band].min())
+        self.assertGreater(widths[1], widths[0])
+
+
+class EdgeCaseTest(unittest.TestCase):
+    def test_no_ocean_gives_zeros(self):
+        g = small_grid(2)
+        u = OC.currents(g, np.zeros(g.n, bool), gyre_wind(g), params(), DAY)
+        self.assertTrue(np.all(u == 0.0))
+
+    def test_all_ocean_globe_is_finite_at_the_poles(self):
+        g = small_grid(3)
+        e, _ = east_north(g.xyz)
+        wind = (-6.0 * np.cos(np.radians(3 * g.lat)))[:, None] * e
+        u = OC.currents(g, np.ones(g.n, bool), wind, params(), DAY)
+        self.assertTrue(np.all(np.isfinite(u)))
+        self.assertLess(np.linalg.norm(u, axis=1).max(), 20.0)
+
+
+class CoarseTest(unittest.TestCase):
+    def test_coarse_fallback_keeps_the_gyre(self):
+        g = small_grid(3)
+        ocean = basin(g)
+        v = northward(g, OC.currents(g, ocean, gyre_wind(g), params(friction_days=1.0, direct_max_cells=1000), DAY))
+        band = ocean & (g.lat > 25) & (g.lat < 35)
+        self.assertTrue(np.all(np.isfinite(v)))
+        self.assertGreater(v[band & (g.lon < -40)].mean(), 0.0)
+        self.assertLess(v[band & (g.lon > -20)].mean(), 0.0)
+
+
+class SSTTest(unittest.TestCase):
+    def test_still_water_keeps_equilibrium(self):
+        g = small_grid(3)
+        ocean = basin(g)
+        T_eq = 30.0 - 0.5 * np.abs(g.lat)
+        T = OC.sst(g, ocean, np.zeros((g.n, 3)), T_eq, params(kappa_m2s=0.0))
+        np.testing.assert_allclose(T, T_eq, atol=1e-6)
+
+    def test_boundary_current_warms_the_west_side(self):
+        g = small_grid(3)
+        ocean = basin(g)
+        u = OC.currents(g, ocean, gyre_wind(g), params(), DAY)
+        T_eq = 30.0 - 0.5 * np.abs(g.lat)
+        a = OC.sst(g, ocean, u, T_eq, params()) - T_eq
+        band = ocean & (g.lat > 32) & (g.lat < 40)
+        self.assertGreater(a[band & (g.lon < -45)].mean(), a[band & (g.lon > 45)].mean() + 0.2)
+        self.assertGreater(a[band & (g.lon < -45)].mean(), 0.0)
+
+
+class UpwellingTest(unittest.TestCase):
+    def setUp(self):
+        self.g = small_grid(3)
+        _, self.n = east_north(self.g.xyz)
+        self.P = params(upwell_coast_km=0.0)
+
+    def coast(self, ocean, lo, hi):
+        g = self.g
+        touches_land = np.bincount(g.src, weights=(~ocean[g.dst]).astype(float), minlength=g.n) > 0
+        return ocean & touches_land & (g.lat > lo) & (g.lat < hi)
+
+    def test_west_coast_upwells(self):
+        g = self.g
+        ocean = ~((g.lon > 0) & (g.lon < 90) & (g.lat > 10) & (g.lat < 50))   # continent east of the sea
+        w = OC.upwelling(g, ocean, -6.0 * self.n, self.P, DAY)                 # equatorward wind (north)
+        self.assertGreater(w[self.coast(ocean, 20, 40) & (g.lon < 0)].mean(), 0.0)
+
+    def test_east_coast_downwells(self):
+        g = self.g
+        ocean = ~((g.lon > -90) & (g.lon < 0) & (g.lat > 10) & (g.lat < 50))  # continent west of the sea
+        w = OC.upwelling(g, ocean, -6.0 * self.n, self.P, DAY)
+        self.assertLess(w[self.coast(ocean, 20, 40) & (g.lon > 0)].mean(), 0.0)
+
+    def test_southern_west_coast_upwells_under_equatorward_wind(self):
+        g = self.g
+        ocean = ~((g.lon > 0) & (g.lon < 90) & (g.lat < -10) & (g.lat > -50))
+        w = OC.upwelling(g, ocean, 6.0 * self.n, self.P, DAY)                  # equatorward in the south = north
+        self.assertGreater(w[self.coast(ocean, -40, -20) & (g.lon < 0)].mean(), 0.0)
+
+    def test_trades_upwell_at_the_equator(self):
+        g = self.g
+        e, _ = east_north(g.xyz)
+        w = OC.upwelling(g, np.ones(g.n, bool), -6.0 * e, self.P, DAY)
+        self.assertGreater(w[np.abs(g.lat) < 3].mean(), 0.0)
+        self.assertLess(w[(np.abs(g.lat) > 8) & (np.abs(g.lat) < 20)].mean(), w[np.abs(g.lat) < 3].mean())
+
+    def test_zero_on_land(self):
+        g = self.g
+        ocean = g.lat < 0
+        self.assertTrue(np.all(OC.upwelling(g, ocean, -6.0 * self.n, params(), DAY)[~ocean] == 0.0))
+
+
+class ProductivityTest(unittest.TestCase):
+    def test_upwelling_coast_beats_gyre_centre(self):
+        g = small_grid(3)
+        ocean = ~((g.lon > 0) & (g.lon < 90) & (g.lat > 10) & (g.lat < 50))
+        _, n = east_north(g.xyz)
+        w = OC.upwelling(g, ocean, -6.0 * n, params(), DAY)
+        p = OC.productivity(g, ocean, w, np.full(g.n, -4000.0), np.zeros(g.n), np.full(g.n, 20.0), params())
+        touches_land = np.bincount(g.src, weights=(~ocean[g.dst]).astype(float), minlength=g.n) > 0
+        coast = ocean & touches_land & (g.lat > 25) & (g.lat < 35) & (g.lon < 0)
+        centre = g.cell_index(30.0, -60.0)
+        self.assertGreater(p[coast].max(), 0.1)
+        self.assertGreater(p[coast].max(), 10.0 * p[centre])
+
+    def test_shelf_and_light(self):
+        g = small_grid(3)
+        ocean = g.lat < 0
+        zero = np.zeros(g.n)
+        shelf = OC.productivity(g, ocean, zero, np.full(g.n, -100.0), zero, np.full(g.n, 25.0), params())
+        deep = OC.productivity(g, ocean, zero, np.full(g.n, -3000.0), zero, np.full(g.n, 25.0), params())
+        i = g.cell_index(-10.0, 0.0)
+        self.assertAlmostEqual(shelf[i], 0.4 * (0.3 + 0.7 * np.cos(np.radians(g.lat[i]))), places=6)  # a_s·light
+        self.assertEqual(deep[i], 0.0)
+        self.assertGreater(shelf[i], shelf[g.cell_index(-80.0, 0.0)])
+
+    def test_sharp_sst_front_is_productive(self):
+        g = small_grid(3)
+        ocean = g.lat < 0
+        zero, deep = np.zeros(g.n), np.full(g.n, -4000.0)
+        front = 15.0 + 8.0 * np.tanh((g.lat + 35.0) / 2.0)             # 16 °C across ≈ 4° (≈ 900 km): ≈ 2 °C/100 km
+        flat = 25.0 + 0.1 * g.lat                                       # background pole-ward cooling, ≈ 0.05 °C/100 km
+        pf = OC.productivity(g, ocean, zero, deep, zero, front, params())
+        pb = OC.productivity(g, ocean, zero, deep, zero, flat, params())
+        at, far = g.cell_index(-35.0, 0.0), g.cell_index(-60.0, 0.0)
+        self.assertGreater(pf[at], 0.15)
+        self.assertLess(pf[far], 0.02)
+        self.assertLess(pb.max(), 1e-9)                                 # gentle background gradients add nothing
+
+    def test_land_sea_contrast_is_not_a_front(self):
+        g = small_grid(3)
+        ocean = g.lat < 0
+        zero, deep = np.zeros(g.n), np.full(g.n, -4000.0)
+        sst_c = np.where(ocean, 20.0, -10.0)                            # land temperatures differ wildly
+        self.assertLess(OC.productivity(g, ocean, zero, deep, zero, sst_c, params()).max(), 1e-9)
+
+    def test_range_and_land(self):
+        g = small_grid(3)
+        ocean = g.lat < 0
+        p = OC.productivity(g, ocean, np.full(g.n, 1e4), np.zeros(g.n), np.full(g.n, 40.0), np.zeros(g.n), params())
+        self.assertTrue(np.all((p >= 0) & (p <= 1)))
+        self.assertTrue(np.all(p[~ocean] == 0))
+
+
+class GaugeTest(unittest.TestCase):
+    def test_no_spurious_vortex_at_the_gauge(self):
+        # cell 0 (79° N, 38° E on H3 grids) at sea in a world with a continent: pinning ψ there must not leave a point
+        # vortex (the solve's compatibility residual) — the speed around cell 0 stays within the ocean's p99
+        g = small_grid(3)
+        land = (np.abs(g.lat) < 40) & (np.abs(g.lon) < 50)
+        ocean = ~land
+        self.assertTrue(ocean[0])
+        e, _ = east_north(g.xyz)
+        wind = (-8.0 * np.cos(np.radians(3.0 * g.lat)))[:, None] * e
+        sp = np.linalg.norm(OC.currents(g, ocean, wind, params(), DAY), axis=1)
+        near = np.zeros(g.n, bool)
+        near[0] = True
+        for _ in range(2):
+            near = near | (np.bincount(g.src, weights=near[g.dst].astype(float), minlength=g.n) > 0)
+        self.assertLessEqual(sp[near].max(), np.percentile(sp[ocean], 99))
diff --git a/tests/test_pipeline.py b/tests/test_pipeline.py
new file mode 100644
index 0000000..adbe3c4
--- /dev/null
+++ b/tests/test_pipeline.py
@@ -0,0 +1,185 @@
+import shutil
+import tempfile
+import unittest
+from pathlib import Path
+
+import numpy as np
+
+from mapgen import pipeline as P
+from mapgen.testing import FIXTURE_TOML
+from tests.helpers import ROOT, make_ctx
+
+TECT = """
+[[plate]]
+id = "a"
+seed = [0.0, 0.0]
+kind = "continental"
+motion = [90.0, 3.0]
+[[plate]]
+id = "b"
+seed = [0.0, 90.0]
+kind = "oceanic"
+motion = [270.0, 3.0]
+"""
+
+
+class PipelineTest(unittest.TestCase):
+    def setUp(self):
+        self.tmp = Path(tempfile.mkdtemp())
+        (self.tmp / "config").mkdir()
+        shutil.copy(FIXTURE_TOML, self.tmp / "config" / "world.toml")
+        (self.tmp / "config" / "tectonics.toml").write_text(TECT)
+        self.log = []
+
+    def tearDown(self):
+        shutil.rmtree(self.tmp)
+
+    def _build(self, **kw):
+        self.log = []
+        return P.build(self.tmp, 1, stop="grid", log=self.log.append, **kw)
+
+    def test_cache_hit(self):
+        ctx = self._build()
+        self.assertEqual(ctx.grid.n, 842)
+        self.assertFalse(any("cached" in m for m in self.log))
+        ctx = self._build()
+        self.assertTrue(any("grid: cached" in m for m in self.log))
+        self.assertEqual(ctx.grid.n, 842)
+
+    def test_config_change_invalidates(self):
+        self._build()
+        p = self.tmp / "config" / "world.toml"
+        p.write_text(p.read_text().replace("seed = 1296", "seed = 1297"))
+        self._build()
+        self.assertFalse(any("cached" in m for m in self.log))
+
+    def test_from_forces_rerun(self):
+        self._build()
+        self._build(start="grid")
+        self.assertFalse(any("cached" in m for m in self.log))
+
+    def test_need_reports_missing(self):
+        ctx = make_ctx(1)
+        with self.assertRaisesRegex(P.StageError, "sk_land"):
+            ctx.need("sk_land")
+
+
+class OceanExportTest(unittest.TestCase):
+    def test_cells_and_fields_carry_the_ocean(self):
+        import json
+        from mapgen.testing import built_world
+        root = built_world()
+        outs = sorted((root / "out").glob("r*/cells.npz"))
+        self.assertTrue(outs)
+        for cells in outs + sorted((root / "out").glob("r*/eras/*/cells.npz")):   # every era re-runs climate
+            z = np.load(cells)
+            for k in ("current", "current_speed", "sst", "upwelling", "productivity"):
+                self.assertIn(k, z.files, f"{cells}: {k}")
+        meta = json.loads((outs[0].parent / "fields.json").read_text())
+        for k in ("sst", "productivity", "current_speed", "upwelling"):
+            self.assertIn(k, meta["continuous"])
+        layers = [L["id"] for L in json.loads((outs[0].parent / "viewer" / "layers.json").read_text())]
+        self.assertIn("currents", layers)
+
+
+class NpzMapsTest(unittest.TestCase):
+    def test_same_arrays_as_np_load_writable_and_file_untouched(self):
+        import hashlib
+        import tempfile
+        import numpy as np
+        from pathlib import Path
+        from mapgen import pipeline as P
+        rng = np.random.default_rng(0)
+        arrays = {"f": rng.normal(size=1000), "i": rng.integers(0, 9, (40, 3)).astype(np.int16),
+                  "b": rng.random(77) < 0.5, "F": np.asfortranarray(rng.random((30, 4))), "s": np.array(2.5),
+                  "e": np.zeros((0, 3), np.float32), "_key": np.array("abc"), "u": np.arange(5, dtype=np.uint64),
+                  "big": rng.normal(size=(3_000_000,)), "nc": rng.random((50, 8))[:, ::3], "x" * 200: np.arange(7)}
+        with tempfile.TemporaryDirectory() as t:
+            for save in (np.savez, np.savez_compressed, P.save_npz_aligned):
+                f = Path(t) / f"{save.__name__}.npz"
+                save(f, **arrays)
+                before = hashlib.sha256(f.read_bytes()).hexdigest()
+                got = P.npz_maps(f)
+                if save is P.save_npz_aligned:                  # every member mapped in place, at aligned addresses
+                    for k, v in got.items():
+                        if v.size:
+                            self.assertEqual(v.ctypes.data % P.ALIGN, 0, k)
+                            b = v
+                            while isinstance(b, (np.ndarray, memoryview)):
+                                b = b.base if isinstance(b, np.ndarray) else b.obj
+                            self.assertIsInstance(b, __import__("mmap").mmap, k)
+                with np.load(f) as z:
+                    self.assertEqual(sorted(got), sorted(z.files))
+                    for k in z.files:
+                        self.assertEqual(got[k].dtype, z[k].dtype, k)
+                        self.assertEqual(got[k].shape, z[k].shape, k)
+                        self.assertTrue(np.array_equal(got[k], z[k]), k)
+                        self.assertTrue(np.array_equal(got[k], arrays[k]), k)
+                for k, v in got.items():                        # whatever is mapped is aligned (else: read)
+                    b = v
+                    while isinstance(b, (np.ndarray, memoryview)):
+                        b = b.base if isinstance(b, np.ndarray) else b.obj
+                    if isinstance(b, __import__("mmap").mmap):
+                        self.assertEqual(v.ctypes.data % P.ALIGN, 0, (save.__name__, k))
+                got["f"][:] = -1.0                       # private copy-on-write: allowed, file unchanged
+                got["F"][0, 0] = 9.0
+                self.assertEqual(float(got["F"][0, 0]), 9.0)
+                self.assertEqual(hashlib.sha256(f.read_bytes()).hexdigest(), before, save.__name__)
+                self.assertTrue(np.array_equal(np.load(f)["f"], arrays["f"]))
+
+    def test_low_memory_rebuild_from_cache_identical(self):
+        import shutil
+        import tempfile
+        from pathlib import Path
+        from mapgen import pipeline as P
+        from mapgen.testing import small_world
+        outs = []
+        tmp = Path(tempfile.mkdtemp())
+        try:
+            small_world(tmp)
+            for low in (False, True, True):                # normal; low (fresh cache hit); low again
+                ctx = P.build(tmp, 2, log=lambda m: None, low_memory=low)
+                outs.append({k: v for k, v in ctx.data.items()})
+            for k, v in outs[0].items():
+                for o in outs[1:]:
+                    self.assertTrue(np.array_equal(np.asarray(o[k]), np.asarray(v), equal_nan=np.asarray(v).dtype.kind == "f"), k)
+            self.assertFalse(list((tmp / "out" / "cache" / "r2").glob("*.tmp-*")), "no temporary files left")
+            import mmap
+            ctx = P.build(tmp, 2, start="climate", log=lambda m: None, low_memory=True)   # freshly computed
+            for k in ("T_jun", "P_ann", "z_surface_m", "plate"):
+                b = ctx.data[k]
+                while isinstance(b, (np.ndarray, memoryview)):
+                    b = b.base if isinstance(b, np.ndarray) else b.obj
+                self.assertIsInstance(b, mmap.mmap, f"{k}: kept on disk, not in RAM")
+                self.assertTrue(np.array_equal(ctx.data[k], outs[0][k], equal_nan=True), k)
+        finally:
+            shutil.rmtree(tmp)
+
+
+class AlignedSumTest(unittest.TestCase):
+    def test_misaligned_data_sums_differently(self):
+        """Why the cache is aligned: numpy sums misaligned float64 data in another order (else drop the alignment)."""
+        import numpy as np
+        x = np.random.default_rng(0).normal(size=2_000_000)
+        buf = bytearray(8 * len(x) + 8)
+        a = np.frombuffer(memoryview(buf)[4:4 + 8 * len(x)], np.float64)
+        np.copyto(a, x)
+        self.assertFalse(a.flags.aligned)
+        self.assertTrue(np.array_equal(a, x))
+        if a.sum() == x.sum():
+            self.skipTest("this numpy sums misaligned data the same way")
+        self.assertNotEqual(a.sum(), x.sum())
+
+
+class BlasSpinTest(unittest.TestCase):
+    def test_import_sets_a_short_blas_spin_unless_given(self):
+        import os
+        import subprocess
+        import sys
+        code = "import os, mapgen, numpy; print(os.environ['OPENBLAS_THREAD_TIMEOUT'])"
+        env = {k: v for k, v in os.environ.items() if k != "OPENBLAS_THREAD_TIMEOUT"}
+        out = subprocess.run([sys.executable, "-c", code], capture_output=True, text=True, env=env, cwd=os.getcwd())
+        self.assertEqual(out.stdout.strip(), "4")
+        out = subprocess.run([sys.executable, "-c", code], capture_output=True, text=True,
+                             env={**env, "OPENBLAS_THREAD_TIMEOUT": "9"}, cwd=os.getcwd())
+        self.assertEqual(out.stdout.strip(), "9")
diff --git a/tests/test_plateaus.py b/tests/test_plateaus.py
new file mode 100644
index 0000000..00da0ed
--- /dev/null
+++ b/tests/test_plateaus.py
@@ -0,0 +1,78 @@
+# map/tests/test_plateaus.py
+import json
+import unittest
+
+import numpy as np
+
+from mapgen import plateaus as PL
+from mapgen.sphere import east_north, great_circle_point, latlon_to_xyz
+from tests.helpers import small_grid
+
+R = 12742.0
+P1 = {"name": "east-flank", "center": [-39.9, -18.0], "area_km2": 3.0e6, "elongation": 1.9, "azimuth_deg": 30.0,
+      "top_m": [1500.0, 3000.0]}
+P2 = {"name": "plateau-02", "center": [-0.4, 176.6], "area_km2": 1.0e6, "top_m": [1000.0, 1500.0], "islands": True}
+
+
+def at(lat, lon):
+    return latlon_to_xyz(lat, lon)[None]
+
+
+class PlateauGeometryTest(unittest.TestCase):
+    def test_centre_inside_and_area_close_to_config(self):
+        g = small_grid(4)
+        for p in (P1, P2):
+            self.assertLess(PL.rho(at(*p["center"]), p, 1296, R)[0], 1e-6)
+            area = g.area_km2[PL.rho(g.xyz, p, 1296, R) < 1.0].sum()
+            self.assertAlmostEqual(area / p["area_km2"], 1.0, delta=0.2, msg=p["name"])
+
+    def test_long_axis_follows_the_azimuth(self):
+        a, b = PL.semi_axes(P1)
+        self.assertAlmostEqual(a / b, 1.9)
+        self.assertAlmostEqual(np.pi * a * b, P1["area_km2"], delta=1.0)
+        c = latlon_to_xyz(*P1["center"])
+        e, n = east_north(c[None])
+        az = np.radians(30.0)
+        along = np.sin(az) * e[0] + np.cos(az) * n[0]
+        across = np.cos(az) * e[0] - np.sin(az) * n[0]
+        self.assertLess(PL.rho(great_circle_point(c, along, 0.8 * a, R)[None], P1, 1296, R)[0], 1.0)
+        self.assertGreater(PL.rho(great_circle_point(c, across, 0.8 * a, R)[None], P1, 1296, R)[0], 1.0)
+
+    def test_features_are_deterministic_inside_and_json_ready(self):
+        f = PL.features(P2, 1296, R)
+        self.assertEqual(f, PL.features(P2, 1296, R))
+        json.dumps(f)
+        self.assertTrue(6 <= len(f["cones"]) <= 20)
+        self.assertTrue(1 <= len(f["calderas"]) <= 3)
+        self.assertTrue(3 <= sum(c["island"] for c in f["cones"]) <= 6)
+        for c in f["cones"] + f["calderas"]:
+            self.assertLess(PL.rho(at(c["lat"], c["lon"]), P2, 1296, R)[0], 0.86)
+        self.assertFalse(any(c["island"] for c in PL.features(P1, 1296, R)["cones"]))
+        self.assertEqual(PL.features({**P1, "vent": 0.0}, 1296, R)["cones"], [])
+
+    def test_surface_depth_follows_top_m(self):
+        g = small_grid(4)
+        inner = g.xyz[PL.rho(g.xyz, P1, 1296, R) < 0.5]
+        z = PL.surface(inner, P1, {"cones": [], "calderas": []}, 1296, R, np.full(len(inner), -5900.0))
+        self.assertTrue(np.all(z <= -P1["top_m"][0] + PL.RELIEF_M + 1e-6))
+        self.assertTrue(np.all(z >= -P1["top_m"][1] - PL.RELIEF_M - 1e-6))
+
+    def test_apply_hidden_clamp_islands_and_nothing_else(self):
+        g = small_grid(4)
+        plats = [P1, P2]
+        ids = PL.cell_ids(g.xyz, plats, 1296, R)
+        self.assertGreater((ids == 0).sum(), 100)
+        z = PL.apply(g, np.full(g.n, -6000.0), plats, ids, 1296)
+        self.assertLessEqual(z[ids == 0].max(), PL.HIDDEN_MAX_M)
+        self.assertGreater(z.max(), 0.0, "island plateau: small volcanic islands")
+        changed_outside = (ids < 0) & (z != -6000.0)
+        self.assertLessEqual(changed_outside.sum(), 6, "only islands' nearest cells may lie outside the outline")
+
+    def test_outline_across_the_antimeridian_and_near_the_pole(self):
+        p = {"name": "plateau-11", "center": [77.0, 179.0], "area_km2": 0.9e6, "top_m": [1000.0, 1500.0]}
+        self.assertLess(PL.rho(at(77.0, 179.8), p, 1296, R)[0], 0.3)
+        self.assertLess(PL.rho(at(77.0, -179.8), p, 1296, R)[0], 0.3)
+        q = {"name": "polar", "center": [88.0, 0.0], "area_km2": 1.0e6, "top_m": [1000.0, 1500.0]}
+        self.assertLess(PL.rho(at(89.9, 90.0), q, 1296, R)[0], 1.0)
+        self.assertLess(PL.rho(at(89.9, -90.0), q, 1296, R)[0], 1.0)
+        self.assertTrue(np.all(np.isfinite(PL.rho(small_grid(3).xyz, q, 1296, R))))
diff --git a/tests/test_plates.py b/tests/test_plates.py
new file mode 100644
index 0000000..9b71bb9
--- /dev/null
+++ b/tests/test_plates.py
@@ -0,0 +1,73 @@
+import unittest
+
+import numpy as np
+
+from mapgen import plates as PL
+from mapgen.pipeline import StageError
+from tests.helpers import make_ctx
+
+
+def two_plates(az_a, az_b, speed=5.0, kinds=("oceanic", "oceanic")):
+    return {"plate": [
+        {"id": "a", "seed": [0.0, -40.0], "kind": kinds[0], "motion": [az_a, speed]},
+        {"id": "b", "seed": [0.0, 40.0], "kind": kinds[1], "motion": [az_b, speed]},
+    ]}
+
+
+def ctx_for(tect, land=None, res=2):
+    ctx = make_ctx(res, tect=tect)
+    n = ctx.grid.n
+    ctx.data.update({"sk_land": np.zeros(n) if land is None else land, "m_land_hint": np.zeros(n)})
+    return ctx
+
+
+class PlatesTest(unittest.TestCase):
+    def test_convergent_then_divergent(self):
+        for az_a, az_b, want in ((90.0, 270.0, PL.CONV), (270.0, 90.0, PL.DIV)):
+            ctx = ctx_for(two_plates(az_a, az_b))
+            out = PL.run(ctx)
+            g = ctx.grid
+            near = (np.abs(g.lon) < 12) & (np.abs(g.lat) < 30) & (out["bnd_type"] > 0)
+            self.assertGreater(near.sum(), 3)
+            self.assertTrue(np.mean(out["bnd_type"][near] == want) > 0.8, (az_a, want))
+
+    def test_every_cell_assigned_and_seeds_own_cells(self):
+        ctx = ctx_for(two_plates(90.0, 270.0))
+        out = PL.run(ctx)
+        g = ctx.grid
+        self.assertEqual(set(np.unique(out["plate"])), {0, 1})
+        self.assertEqual(out["plate"][g.cell_index(0.0, -40.0)], 0)
+        self.assertEqual(out["plate"][g.cell_index(0.0, 40.0)], 1)
+
+    def test_continent_stays_on_its_plate(self):
+        tect = {"plate": [
+            {"id": "c", "seed": [0.0, 0.0], "kind": "continental", "motion": [0.0, 1.0]},
+            {"id": "o", "seed": [0.0, 60.0], "kind": "oceanic", "motion": [0.0, 1.0]},
+        ]}
+        ctx = make_ctx(2, tect=tect)
+        g = ctx.grid
+        land = ((np.abs(g.lat) < 20) & (np.abs(g.lon) < 40)).astype(float)
+        ctx.data.update({"sk_land": land, "m_land_hint": np.zeros(g.n)})
+        out = PL.run(ctx)
+        self.assertTrue(np.all(out["plate"][land > 0.5] == 0))
+
+    def test_duplicate_seed_cell_errors(self):
+        tect = two_plates(90.0, 270.0)
+        tect["plate"][1]["seed"] = [0.01, -40.01]
+        with self.assertRaisesRegex(StageError, "same cell"):
+            PL.run(ctx_for(tect))
+
+    def test_zero_speed_plate(self):
+        out = PL.run(ctx_for(two_plates(90.0, 270.0, speed=0.0)))
+        self.assertTrue(np.all(np.isfinite(out["vel"])))
+        self.assertTrue(np.allclose(out["vel"], 0.0))
+
+
+class MeanderTest(unittest.TestCase):
+    def test_boundaries_meander(self):
+        ctx = ctx_for(two_plates(90.0, 270.0, kinds=("continental", "continental")), res=4)
+        g = ctx.grid
+        ctx.data["sk_land"] = ((np.abs(g.lat) < 40) & (np.abs(g.lon) < 70)).astype(float)
+        out = PL.run(ctx)
+        b = (out["bnd_type"] > 0) & (np.abs(g.lat) < 35) & (np.abs(g.lon) < 40)
+        self.assertGreater(g.lon[b].std(), 3.0)
diff --git a/tests/test_projections.py b/tests/test_projections.py
new file mode 100644
index 0000000..f1a5f90
--- /dev/null
+++ b/tests/test_projections.py
@@ -0,0 +1,55 @@
+import unittest
+
+import numpy as np
+
+from mapgen import projections as PJ
+
+
+class ProjectionTest(unittest.TestCase):
+    def test_equal_earth_roundtrip(self):
+        lat = np.array([0.0, 30.0, -60.0, 89.0, 10.0])
+        lon = np.array([0.0, 120.0, -170.0, 45.0, -179.0])
+        x, y = PJ.equal_earth_forward(lat, lon)
+        la, lo = PJ.equal_earth_inverse(x, y)
+        np.testing.assert_allclose(la, lat, atol=1e-6)
+        np.testing.assert_allclose(lo, lon, atol=1e-6)
+
+    def test_mollweide_roundtrip(self):
+        lat = np.array([0.0, 45.0, -80.0, 89.9])
+        lon = np.array([0.0, -100.0, 170.0, 20.0])
+        x, y = PJ.mollweide_forward(lat, lon)
+        la, lo = PJ.mollweide_inverse(x, y)
+        np.testing.assert_allclose(la, lat, atol=1e-6)
+        np.testing.assert_allclose(lo, lon, atol=1e-6)
+
+    def test_equal_area(self):
+        # equal lat/lon boxes at the equator and at 60° should differ in area by cos(lat) — as on the sphere
+        def box_area(f, lat0):
+            la = np.array([lat0, lat0, lat0 + 1, lat0 + 1]); lo = np.array([0.0, 1.0, 1.0, 0.0])
+            x, y = f(la, lo)
+            return 0.5 * abs(np.dot(x, np.roll(y, -1)) - np.dot(y, np.roll(x, -1)))
+        for f in (PJ.equal_earth_forward, PJ.mollweide_forward):
+            ratio = box_area(f, 60.0) / box_area(f, 0.0)
+            self.assertAlmostEqual(ratio, (np.sin(np.radians(61)) - np.sin(np.radians(60))) / np.sin(np.radians(1)),
+                                   delta=0.01)
+
+    def test_orthographic_centre_and_limb(self):
+        lat, lon, ok = PJ.orthographic_inverse(np.array([0.0, 0.99, 1.5]), np.array([0.0, 0.0, 0.0]), 40.0, 10.0)
+        self.assertAlmostEqual(lat[0], 40.0, places=9)
+        self.assertAlmostEqual(lon[0], 10.0, places=9)
+        self.assertTrue(ok[0] and ok[1] and not ok[2])
+
+    def test_reproject_image(self):
+        img = np.zeros((90, 180, 3), np.uint8)
+        img[:, :90] = (255, 0, 0)            # western hemisphere red
+        img[:, 90:] = (0, 0, 255)            # eastern blue
+        out = PJ.reproject(img, "equal_earth", 400)
+        h, w = out.shape[:2]
+        self.assertTrue(1.9 < w / h < 2.2)
+        self.assertEqual(tuple(out[h // 2, w // 4]), (255, 0, 0))
+        self.assertEqual(tuple(out[h // 2, 3 * w // 4]), (0, 0, 255))
+        self.assertEqual(tuple(out[2, 2]), PJ.BACKGROUND)      # outside the map outline
+        globe = PJ.globe(img, 0.0, -90.0, 200)
+        self.assertEqual(globe.shape[:2], (200, 200))
+        self.assertGreater(int(globe[100, 100, 0]), 150)        # centred on the red hemisphere
+        self.assertEqual(tuple(globe[2, 2]), PJ.BACKGROUND)
diff --git a/tests/test_render.py b/tests/test_render.py
new file mode 100644
index 0000000..4c25cb0
--- /dev/null
+++ b/tests/test_render.py
@@ -0,0 +1,239 @@
+import json
+import re
+import shutil
+import tempfile
+import unittest
+from pathlib import Path
+
+import numpy as np
+from PIL import Image
+
+from mapgen import pipeline as P
+from mapgen import render as RN
+from tests.helpers import small_grid
+from mapgen.sphere import east_north
+from mapgen.testing import FIXTURE_TOML
+
+from mapgen.testing import TECT, small_world  # noqa: F401  (other tests import them from here)
+
+
+class RenderTest(unittest.TestCase):
+    def test_encode_roundtrip(self):
+        v = np.array([-11000.0, 0.0, 8848.3])
+        np.testing.assert_allclose(RN.decode(RN.encode(v, 0.5, -12000.0), 0.5, -12000.0), v, atol=0.25)
+
+    def test_full_pipeline_small(self):
+        tmp = Path(tempfile.mkdtemp())
+        try:
+            small_world(tmp)
+            ctx = P.build(tmp, 2, log=lambda m: None)
+            out = tmp / "out" / "r2"
+            meta = json.loads((out / "fields.json").read_text())
+            im = Image.open(out / "raster" / meta["continuous"]["elevation"]["file"])
+            self.assertEqual(im.size, (256, 128))
+            raw = np.array(im)
+            e = meta["continuous"]["elevation"]
+            z = RN.decode(raw, e["scale"], e["offset"])
+            self.assertTrue(-11500 < z.min() < 0 < z.max() < 13000)
+            cells = np.load(out / "cells.npz")
+            self.assertEqual(len(cells["g_ids"]), ctx.grid.n)
+            self.assertIn("holdridge", cells.files)
+            self.assertTrue((tmp / "previews" / "r2" / "contact_sheet.png").exists())
+            self.assertTrue((out / "geo" / "coast.geojson").exists())
+            for f in ("proj_equal_earth.png", "proj_mollweide.png", "globes_sheet.png", "globe_north_pole.png"):
+                self.assertTrue((tmp / "previews" / "r2" / f).exists(), f)
+            vidx = json.loads((out / "viewer" / "layers.json").read_text())
+            self.assertEqual([L["id"] for L in vidx], ["relief", "biomes", "elevation", "temperature", "rainfall",
+                                                       "seasonality", "landform", "ground", "ice", "deposits", "plates", "o2",
+                                                       "gravity", "pressure", "fire", "seabed", "minerals",
+                                                       "bottom_temp", "sediment", "vent_potential", "currents", "sst", "productivity"])
+            self.assertEqual(Image.open(out / "viewer" / "biomes.jpg").size, (256, 128))
+            self.assertEqual(len(vidx[1]["legend"]["items"]), 38)
+            cm = json.loads((out / "cells_meta.json").read_text())
+            self.assertEqual(cm["res"], 2)
+            self.assertEqual(len(cm["legends"]["holdridge"]), 38)
+        finally:
+            shutil.rmtree(tmp)
+
+    def test_deterministic(self):
+        outs = []
+        for _ in range(2):
+            tmp = Path(tempfile.mkdtemp())
+            try:
+                small_world(tmp)
+                outs.append(P.build(tmp, 2, stop="erosion", log=lambda m: None).data["elevation_eroded_m"])
+            finally:
+                shutil.rmtree(tmp)
+        np.testing.assert_array_equal(outs[0], outs[1])
+
+    def test_viewer_has_the_sea_floor_and_zone_layers(self):
+        from mapgen.testing import built_world
+        root = built_world()
+        ids = [L["id"] for L in json.loads((root / "out" / "r2" / "viewer" / "layers.json").read_text())]
+        for k in ("seabed", "minerals", "bottom_temp", "sediment", "vent_potential", "pressure", "fire"):
+            self.assertIn(k, ids)
+
+
+class DevWidthTest(unittest.TestCase):
+    def test_dev_builds_use_dev_raster_width(self):
+        tmp = Path(tempfile.mkdtemp())
+        try:
+            small_world(tmp)
+            w = tmp / "config" / "world.toml"
+            w.write_text(w.read_text().replace("dev_raster_width = 256", "dev_raster_width = 128"))
+            P.build(tmp, 2, log=lambda m: None)                      # res 2 != res_final → dev width
+            meta = json.loads((tmp / "out" / "r2" / "fields.json").read_text())
+            self.assertEqual(meta["width"], 128)
+        finally:
+            shutil.rmtree(tmp)
+
+
+class RiverDrawTest(unittest.TestCase):
+    def test_rivers_are_drawn_and_seam_skipped(self):
+        from types import SimpleNamespace
+        g = SimpleNamespace(n=4, lat=np.array([0.0, 0.0, 0.0, 0.0]), lon=np.array([-90.0, 0.0, 90.0, 179.0]))
+        recv = np.array([1, 2, 2, 0])
+        img = np.zeros((20, 40, 3), np.uint8)
+        out = RN.draw_rivers(img, g, recv, np.array([True, True, False, True]), np.array([2, 2, 0, 1], np.int8))
+        row = out[9:11].max(axis=0)
+        self.assertTrue(np.all(row[11:29, 2] > 100))           # blue along lon −90..90
+        self.assertEqual(int(out[:, 32:, 2].sum()), 0)          # the 179 → −90 segment wraps the seam: not drawn
+
+
+class RegionalWidthTest(unittest.TestCase):
+    def test_finer_than_final_uses_full_width(self):
+        tmp = Path(tempfile.mkdtemp())
+        try:
+            small_world(tmp)
+            w = tmp / "config" / "world.toml"
+            w.write_text(w.read_text().replace("res_final = 5", "res_final = 1").replace("res_dev = 4", "res_dev = 1")
+                         .replace("dev_raster_width = 256", "dev_raster_width = 128"))
+            P.build(tmp, 2, log=lambda m: None)                  # res 2 > res_final: regional/finer build
+            meta = json.loads((tmp / "out" / "r2" / "fields.json").read_text())
+            self.assertEqual(meta["width"], 256)
+        finally:
+            shutil.rmtree(tmp)
+
+
+class PixelCoastTest(unittest.TestCase):
+    def test_pixel_land_rule(self):
+        ocean_k = np.array([[[True, True, True], [False, False, False], [True, False, False], [True, False, False]]])
+        z = np.array([[50.0, -30.0, 10.0, -10.0]])
+        np.testing.assert_array_equal(RN.pixel_land(ocean_k, z), [[False, True, True, False]])
+
+
+class ReliefRiverOrderTest(unittest.TestCase):
+    def test_first_order_streams_hidden(self):
+        from types import SimpleNamespace
+        g = SimpleNamespace(n=3, lat=np.zeros(3), lon=np.array([-90.0, 0.0, 90.0]))
+        img = np.zeros((20, 40, 3), np.uint8)
+        out = RN.draw_rivers(img, g, np.array([1, 2, 2]), np.array([True, True, False]), np.array([1, 1, 0], np.int8))
+        self.assertEqual(int(out.sum()), 0)
+
+
+class RasterRangeTest(unittest.TestCase):
+    def test_continuous_rasters_hold_every_configurable_value(self):
+        from mapgen import render as RD
+        from mapgen.config import ZONE_FIELDS
+        top = lambda name: RD.CONTINUOUS[name][1] * 65535 + RD.CONTINUOUS[name][2]
+        g_lo = ZONE_FIELDS["gravity_g"][0]
+        need = {"pressure": ZONE_FIELDS["pressure_bar"][1], "o2_fraction": ZONE_FIELDS["o2_fraction"][1],
+                "fire": ZONE_FIELDS["fire_reactivity"][1], "gravity": ZONE_FIELDS["gravity_g"][1],
+                "po2": ZONE_FIELDS["o2_fraction"][1] * ZONE_FIELDS["pressure_bar"][1], "plant_height": 5.0 / g_lo}
+        for name, v in need.items():
+            self.assertGreaterEqual(top(name), v, name)
+
+
+class CurrentArrowsTest(unittest.TestCase):
+    def test_eastward_current_draws_horizontal_arrows_at_sea_only(self):
+        g = small_grid(2)
+        e, _ = east_north(g.xyz)
+        ocean = g.lat < 0
+        base = np.zeros((128, 256, 3), np.uint8)
+        out = RN.draw_currents(base.copy(), g, 0.5 * e, ocean)
+        ys, xs = np.nonzero(out.any(axis=2))
+        self.assertGreater(len(ys), 50)
+        self.assertTrue(np.all(ys >= 60))                          # north half (land) untouched (row 64 = equator)
+        still = RN.draw_currents(base.copy(), g, np.zeros((g.n, 3)), ocean)
+        self.assertFalse(still.any())                              # no current, no arrows
+
+    def test_new_layers_registered(self):
+        for name in ("sst", "productivity", "current_speed", "upwelling"):
+            self.assertIn(name, RN.CONTINUOUS)
+            self.assertIn(name, RN.RAMPS)
+        ids = [v[0] for v in RN.VIEWER_LAYERS]
+        for vid in ("currents", "sst", "productivity"):
+            self.assertIn(vid, ids)
+
+
+class ByRowsTest(unittest.TestCase):
+    """Rasters are coloured and encoded CHUNK rows at a time: the same bytes as in one piece."""
+    def test_ramp_and_encode_same_as_whole(self):
+        from unittest import mock
+        import numpy as np
+        from mapgen import render as R
+        v = np.random.default_rng(3).normal(0, 2000, size=(300, 37))
+        v[5, 5], v[7, 7] = -np.inf, np.inf
+        stops = [[0, 0, 0], [10, 200, 30], [255, 255, 255]]
+        with mock.patch.object(R, "CHUNK", 10**6):
+            want_r, want_e = R._ramp(v, -3000.0, 4000.0, stops), R.encode(v, 0.5, -12000.0)
+        with mock.patch.object(R, "CHUNK", 64):
+            got_r, got_e = R._ramp(v, -3000.0, 4000.0, stops), R.encode(v, 0.5, -12000.0)
+        self.assertEqual(got_r.shape, (300, 37, 3))
+        self.assertEqual(got_r.dtype, np.uint8)
+        self.assertEqual(got_e.dtype, np.uint16)
+        self.assertTrue(np.array_equal(got_r, want_r))
+        self.assertTrue(np.array_equal(got_e, want_e))
+        # hand-derived: midpoint of a 0..1 ramp between 0 and 200 → 100
+        self.assertEqual(R._ramp(np.full((200, 2), 0.5), 0.0, 1.0, [[0, 0, 0], [200, 200, 200]])[150, 1, 0], 100)
+        self.assertEqual(int(R.encode(np.full((200, 2), 3.0), 0.5, -12000.0)[199, 0]), 24006)
+
+
+class CompiledRenderTest(unittest.TestCase):
+    def test_ramp_same_bytes_as_numpy_formula(self):
+        import numpy as np
+        from mapgen import render as R
+        rng = np.random.default_rng(7)
+        v = np.concatenate([rng.normal(0, 3000, 199_997), [np.inf, -np.inf, 0.0, -0.0, 4000.0, -3000.0, 500.0]])
+        v = v.reshape(-1, 7)
+        for lo, hi, stops in ((-3000, 4000, [[0, 0, 0], [10, 200, 30], [255, 255, 255]]),
+                              (-6500.0, 0.0, [[11, 43, 90], [30, 90, 150], [143, 198, 224]]),
+                              (0.1, 0.35, [[1, 2, 3], [250, 9, 77], [3, 255, 100], [255, 255, 255]]),
+                              (0, 1, [[0, 0, 0], [255, 255, 255]])):
+            st = np.asarray(stops, dtype=np.float64)
+            t = np.clip((v - lo) / (hi - lo), 0, 1) * (len(st) - 1)          # the numpy statements, written out
+            i = np.minimum(t.astype(np.int64), len(st) - 2)
+            f = (t - i)[..., None]
+            want = (st[i] * (1 - f) + st[i + 1] * f).astype(np.uint8)
+            got = R._ramp(v, lo, hi, stops)
+            self.assertEqual(got.dtype, np.uint8)
+            self.assertTrue(np.array_equal(got, want), (lo, hi))
+
+    def test_sample_cont_same_as_numpy_sum(self):
+        import numpy as np
+        from mapgen import render as R
+        rng = np.random.default_rng(8)
+        v = rng.normal(size=5000) * 10.0 ** rng.integers(-6, 6, 5000)
+        for k in (1, 3, 5):
+            idx = rng.integers(0, 5000, (300, 170, k)).astype(np.int32)
+            w = rng.random((300, 170, k)).astype(np.float32)
+            want = np.sum(v[idx] * w, axis=-1)
+            self.assertTrue(np.array_equal(R.sample_cont(v, idx, w), want), k)
+
+    def test_writer_finishes_everything_and_raises_errors(self):
+        import threading
+        from mapgen import render as R
+        done = []
+        w = R._Writer(threads=2, depth=2)
+        for i in range(9):
+            w(lambda i=i: done.append(i))
+        w.close()
+        self.assertEqual(sorted(done), list(range(9)))
+        w = R._Writer()
+        w(lambda: (_ for _ in ()).throw(OSError("disk full")))
+        with self.assertRaises(OSError):
+            w.close()
+        with self.assertRaises(OSError):                 # surfaced at the latest by the next call over depth
+            w = R._Writer(threads=1, depth=1)
+            w(lambda: (_ for _ in ()).throw(OSError("disk full")))
+            w(lambda: None)
diff --git a/tests/test_seabed.py b/tests/test_seabed.py
new file mode 100644
index 0000000..607aeec
--- /dev/null
+++ b/tests/test_seabed.py
@@ -0,0 +1,94 @@
+import json
+import unittest
+
+import numpy as np
+
+from mapgen import seabed as SB
+from mapgen.graph import distance_to
+from tests.helpers import make_ctx
+
+
+def sea_ctx(tect=None):
+    ctx = make_ctx(3, tect=tect or {"plate": []})
+    g = ctx.grid
+    land = (np.abs(g.lat) < 30) & (np.abs(g.lon) < 40)
+    ridge = np.abs(g.lon - 120) < 1.0
+    d_div = distance_to(g, ridge)
+    age = np.where(land, 0.0, np.minimum(d_div / 30.0, 200.0))
+    ctx.data.update({"elevation_eroded_m": np.where(land, 500.0, -4500.0).astype(np.float32), "ocean": ~land,
+                     "continental": land, "ocean_age_myr": age.astype(np.float32), "d_div_km": d_div,
+                     "d_over_km": np.full(g.n, np.inf), "d_sub_km": np.full(g.n, np.inf),
+                     "T_mean": 25.0 - 0.5 * np.abs(g.lat), "vel": np.zeros((g.n, 3))})
+    return ctx
+
+
+class SeabedTest(unittest.TestCase):
+    def test_land_has_no_sea_floor(self):
+        ctx = sea_ctx()
+        out = SB.run(ctx)
+        land = ~ctx.data["ocean"]
+        for k in ("vent_potential", "seabed_type", "seabed_mineral", "bottom_temp_c", "sediment_m"):
+            self.assertTrue(np.all(np.asarray(out[k])[land] == 0), k)
+
+    def test_vents_on_the_ridge_not_far_from_it(self):
+        ctx = sea_ctx()
+        out = SB.run(ctx)
+        d = ctx.data["d_div_km"]
+        self.assertGreaterEqual(out["vent_potential"][d == 0].min(), 0.5)
+        self.assertLess(out["vent_potential"][ctx.data["ocean"] & (d > 3000)].max(), 0.05)
+
+    def test_sediment_thickens_with_age_away_from_land(self):
+        ctx = sea_ctx()
+        out = SB.run(ctx)
+        g = ctx.grid
+        young = ctx.data["ocean"] & (ctx.data["d_div_km"] < 300) & (np.abs(g.lon) > 90)
+        old = ctx.data["ocean"] & (ctx.data["ocean_age_myr"] > 100) & (np.abs(g.lon) > 90)
+        self.assertLess(out["sediment_m"][young].mean(), out["sediment_m"][old].mean())
+
+    def test_bottom_temperature(self):
+        ctx = sea_ctx()
+        g = ctx.grid
+        eq, polar = g.cell_index(0.0, 170.0), g.cell_index(-80.0, 170.0)
+        z = np.asarray(ctx.data["elevation_eroded_m"]).copy()
+        shallow = g.cell_index(-20.0, 170.0)
+        z[shallow] = -100.0
+        ctx.data["elevation_eroded_m"] = z
+        out = SB.run(ctx)
+        self.assertAlmostEqual(out["bottom_temp_c"][eq], 1.0 + 3.0 * np.cos(np.radians(g.lat[eq])) ** 2, places=4)
+        self.assertLess(out["bottom_temp_c"][polar], 1.3)
+        deep = 1.0 + 3.0 * np.cos(np.radians(g.lat[shallow])) ** 2
+        t = 25.0 - 0.5 * abs(g.lat[shallow])
+        self.assertAlmostEqual(out["bottom_temp_c"][shallow], deep + 0.875 * (t - deep), places=3)
+
+    def test_plateau_volcanic_field_makes_volcanic_sea_floor(self):
+        plat = {"name": "p", "center": [0.0, -150.0], "area_km2": 3.0e6, "top_m": [1500.0, 3000.0]}
+        ctx = sea_ctx({"plate": [], "plateau": [plat]})
+        out = SB.run(ctx)
+        from mapgen.crust import center_dist
+        near = center_dist(ctx.grid, [0.0, -150.0]) < 800.0
+        self.assertTrue(np.isin(out["seabed_type"][near], [SB.SB_VOLCANIC, SB.SB_VENTS]).any())
+        self.assertGreater(out["vent_potential"][near].max(), 0.3)
+
+    def test_deterministic(self):
+        a, b = SB.run(sea_ctx()), SB.run(sea_ctx())
+        for k in a:
+            np.testing.assert_array_equal(a[k], b[k])
+
+
+class SeabedRenderTest(unittest.TestCase):
+    def test_rasters_legends_and_plateaus_in_the_build(self):
+        from mapgen.testing import built_world
+        out = built_world() / "out" / "r2"
+        fields = json.loads((out / "fields.json").read_text())
+        for name in ("pressure", "o2_fraction", "fire", "vent_potential", "bottom_temp", "sediment",
+                     "plant_height"):
+            self.assertIn(name, fields["continuous"])
+            self.assertTrue((out / "raster" / f"{name}.png").exists(), name)
+        for name in ("seabed_type", "seabed_mineral"):
+            self.assertIn(name, fields["categorical"])
+        meta = json.loads((out / "cells_meta.json").read_text())
+        self.assertEqual(meta["legends"]["seabed_type"], SB.SEABED_NAMES)
+        self.assertEqual(meta["plateaus"], [])
+        with np.load(out / "cells.npz") as z:
+            for k in ("vent_potential", "seabed_type", "pressure_bar", "fire_reactivity", "plateau_id"):
+                self.assertIn(k, z.files)
diff --git a/tests/test_sketch.py b/tests/test_sketch.py
new file mode 100644
index 0000000..f827943
--- /dev/null
+++ b/tests/test_sketch.py
@@ -0,0 +1,151 @@
+import shutil
+import tempfile
+import unittest
+from pathlib import Path
+
+import numpy as np
+from PIL import Image
+
+from mapgen import sketch as SK
+from mapgen.pipeline import StageError
+from tests.helpers import make_ctx
+
+
+class SampleTest(unittest.TestCase):
+    def test_bilinear_centres_and_wrap(self):
+        img = np.arange(32, dtype=np.float64).reshape(4, 8)
+        # pixel (1, 2) centre: lon = 2.5/8*360-180 = -67.5, lat = 90-1.5/4*180 = 22.5
+        self.assertAlmostEqual(SK.sample_equirect(img, np.array([22.5]), np.array([-67.5]))[0], img[1, 2])
+        a = SK.sample_equirect(img, np.array([22.5]), np.array([179.999]))[0]
+        b = SK.sample_equirect(img, np.array([22.5]), np.array([-179.999]))[0]
+        self.assertAlmostEqual(a, b, places=2)  # wraps across the antimeridian
+        self.assertAlmostEqual(SK.sample_equirect(img, np.array([90.0]), np.array([-157.5]))[0], img[0, 0])
+
+    def test_mask_8bit_and_any_size(self):
+        tmp = Path(tempfile.mkdtemp())
+        try:
+            Image.fromarray(np.full((19, 37), 255, np.uint8), "L").save(tmp / "a.png")
+            Image.fromarray(np.full((10, 20), 32768, np.uint16)).save(tmp / "b.png")
+            self.assertGreater(SK.load_mask(tmp / "a.png").min(), 0.99)
+            self.assertAlmostEqual(float(np.abs(SK.load_mask(tmp / "b.png")).max()), 0.0)
+        finally:
+            shutil.rmtree(tmp)
+
+
+class StageTest(unittest.TestCase):
+    def setUp(self):
+        self.tmp = Path(tempfile.mkdtemp())
+
+    def tearDown(self):
+        shutil.rmtree(self.tmp)
+
+    def test_requires_import(self):
+        ctx = make_ctx(1, root=self.tmp)
+        with self.assertRaisesRegex(StageError, "new-world"):
+            SK.run(ctx)
+
+    def test_outputs(self):
+        (self.tmp / "sketch").mkdir()
+        land = np.zeros((100, 200), np.uint8)
+        land[:, :100] = 255  # western hemisphere is land
+        for n in SK.SKETCH:
+            Image.fromarray(land if n == "land" else np.zeros_like(land), "L").save(self.tmp / "sketch" / f"{n}.png")
+        ctx = make_ctx(2, root=self.tmp)
+        out = SK.run(ctx)
+        g = ctx.grid
+        self.assertGreater(out["sk_land"][g.lon < -10].mean(), 0.95)
+        self.assertLess(out["sk_land"][g.lon > 10].mean(), 0.05)
+        for m in SK.MASKS:
+            self.assertTrue(np.all(out[f"m_{m}"] == 0))
+
+
+class MaskFormatsTest(unittest.TestCase):
+    def setUp(self):
+        self.tmp = Path(tempfile.mkdtemp())
+
+    def tearDown(self):
+        shutil.rmtree(self.tmp)
+
+    def test_gimp_formats(self):
+        Image.fromarray(np.full((8, 16), 128, np.uint8), "L").save(self.tmp / "grey.png")
+        pal = Image.new("P", (16, 8), 0)
+        pal.putpalette([128, 128, 128, 255, 255, 255] + [0] * 762)
+        pal.save(self.tmp / "pal_grey.png")
+        pal.paste(1, (0, 0, 16, 8))
+        pal.save(self.tmp / "pal_white.png")
+        Image.new("1", (16, 8), 1).save(self.tmp / "bit_white.png")
+        Image.new("RGBA", (16, 8), (0, 0, 0, 0)).save(self.tmp / "clear.png")
+        self.assertEqual(float(np.abs(SK.load_mask(self.tmp / "grey.png")).max()), 0.0)
+        self.assertEqual(float(np.abs(SK.load_mask(self.tmp / "pal_grey.png")).max()), 0.0)
+        self.assertAlmostEqual(float(SK.load_mask(self.tmp / "pal_white.png").min()), 1.0)
+        self.assertAlmostEqual(float(SK.load_mask(self.tmp / "bit_white.png").min()), 1.0)
+        self.assertEqual(float(np.abs(SK.load_mask(self.tmp / "clear.png")).max()), 0.0)
+
+    def test_corrupt_mask_is_stage_error(self):
+        (self.tmp / "bad.png").write_bytes(b"not a png")
+        with self.assertRaisesRegex(StageError, "bad.png"):
+            SK.load_mask(self.tmp / "bad.png")
+
+
+class SketchWarpTest(unittest.TestCase):
+    def setUp(self):
+        self.tmp = Path(tempfile.mkdtemp())
+        (self.tmp / "sketch").mkdir()
+        lat = 90 - (np.arange(100) + 0.5) * 1.8
+        lon = (np.arange(200) + 0.5) * 1.8 - 180
+        LA, LO = np.meshgrid(lat, lon, indexing="ij")
+        land = ((np.abs(LA) < 35) & (np.abs(LO) < 60)).astype(np.uint8) * 255
+        for n in SK.SKETCH:
+            Image.fromarray(land if n == "land" else np.zeros_like(land), "L").save(self.tmp / "sketch" / f"{n}.png")
+
+    def tearDown(self):
+        shutil.rmtree(self.tmp)
+
+    def _land(self, **over):
+        ctx = make_ctx(3, root=self.tmp, cfg={"sketch": over} if over else None)
+        return SK.run(ctx)["sk_land"] > 0.5, ctx.grid
+
+    def test_warp_off_matches_drawing(self):
+        warped, g = self._land(warp_km=0.0, detail_warp_km=0.0)
+        direct = SK.sample_equirect(SK.load_png01(self.tmp / "sketch" / "land.png"), g.lat, g.lon) > 0.5
+        np.testing.assert_array_equal(warped, direct)
+
+    def test_default_warp_deviates_but_keeps_continent(self):
+        straight, g = self._land(warp_km=0.0, detail_warp_km=0.0)
+        warped, _ = self._land()
+        iou = (warped & straight).sum() / (warped | straight).sum()
+        self.assertTrue(0.35 < iou < 0.85, iou)
+        again, _ = self._land()
+        np.testing.assert_array_equal(warped, again)
+
+
+class SketchMoveTest(unittest.TestCase):
+    def setUp(self):
+        self.tmp = Path(tempfile.mkdtemp())
+        (self.tmp / "sketch").mkdir()
+        lat = 90 - (np.arange(200) + 0.5) * 0.9
+        lon = (np.arange(400) + 0.5) * 0.9 - 180
+        LA, LO = np.meshgrid(lat, lon, indexing="ij")
+        land = (((LA - 0) ** 2 + (LO - 0) ** 2 < 15 ** 2) | ((LA + 40) ** 2 + (LO - 100) ** 2 < 10 ** 2))
+        for n in SK.SKETCH:
+            Image.fromarray((land * 255).astype(np.uint8) if n == "land" else np.zeros(land.shape, np.uint8), "L") \
+                .save(self.tmp / "sketch" / f"{n}.png")
+
+    def tearDown(self):
+        shutil.rmtree(self.tmp)
+
+    def _land(self, moves):
+        ctx = make_ctx(4, root=self.tmp, cfg={"sketch": {"warp_km": 0.0, "detail_warp_km": 0.0, "moves": moves}})
+        return SK.run(ctx)["sk_land"] > 0.5, ctx.grid
+
+    def test_move_rotates_one_continent_and_keeps_its_area(self):
+        before, g = self._land([])
+        after, _ = self._land([{"at": [0.0, 0.0], "to": [30.0, 0.0]}])
+        near = lambda la, lo, r: g.radius_km * np.arccos(np.clip(g.xyz @ g.xyz[g.cell_index(la, lo)], -1, 1)) < r
+        self.assertTrue(before[near(0.0, 0.0, 800)].all() and not after[near(0.0, 0.0, 800)].any())
+        self.assertTrue(after[near(30.0, 0.0, 800)].all())
+        other = near(-40.0, 100.0, 600)
+        np.testing.assert_array_equal(before[other], after[other])            # the other island stays
+        a = g.area_km2
+        moved_area = a[after & ~other].sum() / a[before & ~other].sum()
+        self.assertAlmostEqual(moved_area, 1.0, delta=0.05)                   # rotation preserves area
diff --git a/tests/test_sphere_noise.py b/tests/test_sphere_noise.py
new file mode 100644
index 0000000..a3c325e
--- /dev/null
+++ b/tests/test_sphere_noise.py
@@ -0,0 +1,86 @@
+import unittest
+
+import numpy as np
+
+from mapgen import noise as N
+from mapgen import sphere as S
+
+
+class SphereTest(unittest.TestCase):
+    def test_roundtrip(self):
+        lat = np.array([0.0, 45.0, -60.0, 89.0])
+        lon = np.array([0.0, 120.0, -170.0, 10.0])
+        la, lo = S.xyz_to_latlon(S.latlon_to_xyz(lat, lon))
+        np.testing.assert_allclose(la, lat, atol=1e-9)
+        np.testing.assert_allclose(lo, lon, atol=1e-9)
+
+    def test_east_north_at_pole(self):
+        p = S.latlon_to_xyz(np.array([90.0, -90.0, 0.0]), np.array([0.0, 0.0, 0.0]))
+        e, n = S.east_north(p)
+        self.assertTrue(np.all(np.isfinite(e)) and np.all(np.isfinite(n)))
+        np.testing.assert_allclose(np.linalg.norm(e, axis=1), 1.0)
+        np.testing.assert_allclose(np.sum(e * p, axis=1), 0.0, atol=1e-12)
+        # at the equator/prime meridian east = +y, north = +z
+        np.testing.assert_allclose(e[2], [0, 1, 0], atol=1e-12)
+        np.testing.assert_allclose(n[2], [0, 0, 1], atol=1e-12)
+
+    def test_motion_velocity_matches_request(self):
+        R = 12742.0
+        om = S.motion_to_omega(10.0, 20.0, 90.0, 5.0, R)  # due east, 5 cm/yr
+        p = S.latlon_to_xyz(np.array([10.0]), np.array([20.0]))
+        v = S.velocity(p, om, R)[0]
+        e, n = S.east_north(p)
+        self.assertAlmostEqual(float(np.dot(v, e[0])), 0.05, places=6)
+        self.assertAlmostEqual(float(np.dot(v, n[0])), 0.0, places=6)
+
+    def test_zero_speed_is_zero_velocity(self):
+        om = S.motion_to_omega(0.0, 0.0, 45.0, 0.0, 12742.0)
+        self.assertTrue(np.allclose(om, 0.0))
+
+    def test_rotate_and_azimuth(self):
+        p = S.latlon_to_xyz(np.array([0.0]), np.array([0.0]))
+        e, n = S.east_north(p)
+        r = S.rotate_about(p, e, np.pi / 2)  # CCW seen from outside: east -> north
+        np.testing.assert_allclose(r, n, atol=1e-12)
+        q = S.latlon_to_xyz(np.array([0.0, 10.0]), np.array([10.0, 0.0]))
+        np.testing.assert_allclose(S.azimuth_deg(p[0], q), [90.0, 0.0], atol=1e-9)
+
+    def test_great_circle_point(self):
+        p0 = S.latlon_to_xyz(np.array([0.0]), np.array([0.0]))[0]
+        e, _ = S.east_north(p0[None])
+        R = 12742.0
+        q = S.great_circle_point(p0, e[0], np.pi * R / 2, R)  # quarter turn east
+        np.testing.assert_allclose(q, [0, 1, 0], atol=1e-12)
+
+
+class NoiseTest(unittest.TestCase):
+    def test_deterministic_and_bounded(self):
+        rng = np.random.default_rng(0)
+        p = rng.normal(size=(5000, 3))
+        p /= np.linalg.norm(p, axis=1, keepdims=True)
+        a = N.fbm(p, 42)
+        b = N.fbm(p, 42)
+        c = N.fbm(p, 43)
+        np.testing.assert_array_equal(a, b)
+        self.assertFalse(np.allclose(a, c))
+        self.assertTrue(np.all(np.abs(a) <= 1.0))
+        self.assertGreater(a.std(), 0.05)
+        r = N.ridged(p, 1)
+        self.assertTrue(np.all((r >= 0) & (r <= 1)))
+
+    def test_continuity(self):
+        p = np.array([[0.3, 0.4, 0.5]])
+        d = N.fbm(p + 1e-6, 5) - N.fbm(p, 5)
+        self.assertLess(abs(d[0]), 1e-3)
+
+
+class CompiledNoiseTest(unittest.TestCase):
+    def test_compiled_value_noise_matches_numpy(self):
+        from mapgen import noise as N
+        if N._noise_jit is None:
+            self.skipTest("numba not installed")
+        rng = np.random.default_rng(0)
+        for scale in (1.0, 37.0, 1e4, 1e7):
+            for seed in (0, 4242, 123456789, -5, 2 ** 40):
+                p = rng.normal(0, scale, (2000, 3))
+                self.assertTrue(np.array_equal(N.value_noise(p, seed), N._value_noise(p, seed)), (scale, seed))
diff --git a/tests/test_viewer_export.py b/tests/test_viewer_export.py
new file mode 100644
index 0000000..b34f2f9
--- /dev/null
+++ b/tests/test_viewer_export.py
@@ -0,0 +1,33 @@
+import json
+import shutil
+import tempfile
+import unittest
+from pathlib import Path
+
+import numpy as np
+from PIL import Image
+
+from mapgen import viewer_export as VE
+
+
+class ViewerExportTest(unittest.TestCase):
+    def test_write_all(self):
+        tmp = Path(tempfile.mkdtemp())
+        try:
+            rgb = np.zeros((20, 40, 3), np.uint8)
+            rgb[:, :20] = (200, 30, 30)
+            VE.write_all(tmp / "viewer", [
+                {"id": "a", "name": "A", "rgb": rgb, "legend": None},
+                {"id": "b", "name": "B", "rgb": rgb,
+                 "legend": VE.categorical_legend(["x", "y"], np.array([[255, 0, 0], [0, 0, 255]]))},
+                {"id": "c", "name": "C", "rgb": rgb,
+                 "legend": VE.continuous_legend("m", -10.0, 10.0, [[0, 0, 0], [128, 128, 128], [255, 255, 255]])},
+            ])
+            idx = json.loads((tmp / "viewer" / "layers.json").read_text())
+            self.assertEqual([L["id"] for L in idx], ["a", "b", "c"])
+            self.assertEqual(idx[1]["legend"]["items"][1], {"name": "y", "color": "#0000ff"})
+            self.assertEqual(idx[2]["legend"]["stops"], [[-10.0, "#000000"], [0.0, "#808080"], [10.0, "#ffffff"]])
+            im = Image.open(tmp / "viewer" / "a.jpg")
+            self.assertEqual((im.size, im.format), ((40, 20), "JPEG"))
+        finally:
+            shutil.rmtree(tmp)
diff --git a/tests/test_zones.py b/tests/test_zones.py
new file mode 100644
index 0000000..f5b73c1
--- /dev/null
+++ b/tests/test_zones.py
@@ -0,0 +1,66 @@
+# map/tests/test_zones.py
+import tempfile
+import unittest
+from pathlib import Path
+
+import numpy as np
+from PIL import Image
+
+from mapgen import sketch as SK, zones as ZN
+from mapgen.sphere import east_north, great_circle_point, latlon_to_xyz
+from tests.helpers import make_ctx
+
+R = 12742.0
+R1 = {"name": "R1", "field": "o2", "center": [-5.6, -118.8], "radius_km": 3200.0, "v": 1.0}
+L1 = {"name": "L1", "field": "gravity", "center": [-26.3, -133.7], "radius_km": 2200.0, "v": -0.6}
+
+
+def ray(zone, az_deg, step_km=5.0):
+    p0 = latlon_to_xyz(*zone["center"])
+    e, n = east_north(p0[None])
+    a = np.radians(az_deg)
+    t = np.sin(a) * e[0] + np.cos(a) * n[0]
+    s = np.arange(0.0, 1.4 * zone["radius_km"], step_km)
+    return s, np.array([great_circle_point(p0, t, si, R) for si in s])
+
+
+class ZoneProfileTest(unittest.TestCase):
+    def test_peak_at_the_centre_and_nothing_beyond_the_warped_radius(self):
+        for z in (R1, L1):
+            s, pts = ray(z, 40.0)
+            f = ZN.contribution(pts, z, 1296, R)
+            self.assertAlmostEqual(f[0], z["v"])
+            self.assertTrue(np.all(f[s > z["radius_km"] / (1 - ZN.WARP) + 1] == 0.0))
+
+    def test_steepest_change_is_hard_to_notice_over_a_days_walk(self):
+        """≤ 0.23 O₂ points and ≤ 0.014 g per 30 km, for the strongest O₂ and low-gravity zones."""
+        for z, per_unit, limit in ((R1, 0.21 * 0.5 * 100.0, 0.23), (L1, 1.05 * 0.7, 0.014)):
+            worst = 0.0
+            for az in range(0, 360, 5):
+                _, pts = ray(z, az)
+                f = ZN.contribution(pts, z, 1296, R) * per_unit
+                worst = max(worst, float(np.max(np.abs(np.diff(f)))) / 5.0 * 30.0)
+            self.assertLessEqual(worst, limit, z["name"])
+
+    def test_outline_is_not_a_circle(self):
+        at = [ZN.contribution(ray(R1, az)[1][[400]], R1, 1296, R)[0] for az in range(0, 360, 30)]   # 2,000 km out
+        self.assertGreater(max(at) - min(at), 0.02)
+
+
+class SketchZonesTest(unittest.TestCase):
+    def test_config_zones_add_to_painted_masks_and_clip(self):
+        root = Path(tempfile.mkdtemp())
+        (root / "sketch").mkdir()
+        (root / "masks").mkdir()
+        for n in SK.SKETCH:
+            Image.fromarray(np.zeros((10, 20), np.uint8), "L").save(root / "sketch" / f"{n}.png")
+        Image.fromarray(np.full((10, 20), 191, np.uint8), "L").save(root / "masks" / "o2_zones.png")   # ≈ +0.5
+        zones = [{**R1, "center": [0.0, 0.0]}, {**L1, "center": [0.0, 0.0]}]
+        ctx = make_ctx(2, tect={"plate": [], "zone": zones}, root=root)
+        out = SK.run(ctx)
+        g = ctx.grid
+        c, far = g.cell_index(0.0, 0.0), g.cell_index(0.0, 180.0)
+        self.assertEqual(out["m_o2_zones"][c], 1.0, "mask + zone, clipped")
+        self.assertAlmostEqual(out["m_o2_zones"][far], 0.498, places=2)
+        self.assertAlmostEqual(out["m_gravity_zones"][c], -0.6, delta=0.05)
+        self.assertEqual(out["m_gravity_zones"][far], 0.0)