1 files changed, 117 insertions, 0 deletions
diff --git a/mapgen/newworld.py b/mapgen/newworld.py
new file mode 100644
index 0000000..80109b0
--- /dev/null
+++ b/mapgen/newworld.py
@@ -0,0 +1,117 @@
+"""`mapgen.py --world DIR new-world`: a random starting world (sketch/*.png continents, config/tectonics.toml plates),
+so a build runs without drawing anything. Used for the moons; edit or redraw afterwards."""
+from __future__ import annotations
+
+from pathlib import Path
+
+import numpy as np
+from PIL import Image
+from scipy import ndimage
+
+from .noise import fbm, ridged
+from .sphere import latlon_to_xyz
+
+W, H = 2000, 1000
+
+
+def _spread(rng, n, min_deg, avoid=(), max_lat=70.0, tries=4000):
+    pts = list(avoid)
+    out = []
+    for _ in range(tries):
+        if len(out) == n:
+            break
+        lat = np.degrees(np.arcsin(rng.uniform(np.sin(np.radians(-max_lat)), np.sin(np.radians(max_lat)))))
+        lon = rng.uniform(-180, 180)
+        p = latlon_to_xyz(lat, lon)
+        if all(np.degrees(np.arccos(np.clip(p @ latlon_to_xyz(*q), -1, 1))) >= min_deg for q in pts):
+            pts.append((lat, lon))
+            out.append((float(lat), float(lon)))
+    return out
+
+
+def _centroid(mask, LAT, LON):
+    w = np.cos(np.radians(LAT)) * mask
+    p = (latlon_to_xyz(LAT, LON) * w[..., None]).sum(axis=(0, 1))
+    p /= np.linalg.norm(p)
+    return float(np.degrees(np.arcsin(p[2]))), float(np.degrees(np.arctan2(p[1], p[0])))
+
+
+def make(root: Path, seed: int, continents: int = 6, land_share: float = 0.27, toward=None, ridge: float = 0.86,
+         force: bool = False) -> None:
+    root = Path(root)
+    sk, cfg = root / "sketch", root / "config"
+    targets = [sk / "land.png", cfg / "tectonics.toml"]
+    if not force and any(p.exists() for p in targets):
+        raise SystemExit("new-world: sketch/ or config/tectonics.toml already exist (use --force to overwrite)")
+    rng = np.random.default_rng(seed)
+    w, h = W // 2, H // 2
+    LAT, LON = np.meshgrid(90.0 - (np.arange(h) + 0.5) * 180.0 / h, -180.0 + (np.arange(w) + 0.5) * 360.0 / w,
+                           indexing="ij")
+    xyz = latlon_to_xyz(LAT, LON).reshape(-1, 3)
+
+    centres = _spread(rng, continents, 35.0, max_lat=55.0)
+    field = np.zeros(len(xyz))
+    for lat, lon in centres:
+        r = np.radians(rng.uniform(18.0, 34.0))
+        d = np.arccos(np.clip(xyz @ latlon_to_xyz(lat, lon), -1, 1))
+        field = np.maximum(field, np.clip(1.0 - d / r, 0.0, None) * rng.uniform(0.8, 1.2))
+    field += 0.55 * fbm(xyz, seed + 7, 6, 2.2)
+    if toward is not None:
+        field += 1.2 * (xyz @ latlon_to_xyz(*toward))
+    area = np.cos(np.radians(LAT)).ravel()
+    order = np.argsort(-field)
+    cut = field[order][np.searchsorted(np.cumsum(area[order]) / area.sum(), land_share)]
+    land = (field > cut).reshape(h, w)
+
+    mount = (ridged(xyz, seed + 11, 5, 3.0) > ridge).reshape(h, w) & ndimage.binary_erosion(land, iterations=6)
+
+    def save(name, m):
+        Image.fromarray((m * 255).astype(np.uint8), "L").resize((W, H), Image.BILINEAR).point(
+            lambda v: 255 if v > 127 else 0).save(sk / f"{name}.png")
+
+    sk.mkdir(parents=True, exist_ok=True)
+    cfg.mkdir(parents=True, exist_ok=True)
+    save("land", land)
+    save("mountains", mount)
+    for name in ("desert", "rainforest", "trench"):
+        save(name, np.zeros_like(land))
+
+    lab, n = ndimage.label(land)
+    for a, b in zip(lab[:, 0], lab[:, -1]):
+        if a and b and a != b:
+            lab[lab == b] = a
+    sizes = sorted(((int((lab == i).sum()), i) for i in np.unique(lab) if i), reverse=True)
+    plates, seeds = [], []
+    for k, (px, i) in enumerate(sizes[:continents]):
+        m = lab == i
+        if px > 0.035 * m.size:                                    # a big landmass: two plates, a collision belt between
+            lat, lon = _centroid(m, LAT, LON)
+            ang = rng.uniform(0, np.pi)
+            side = (np.sin(ang) * (LAT - lat) + np.cos(ang) * ((LON - lon + 180) % 360 - 180) * np.cos(np.radians(lat))) > 0
+            parts = [m & side, m & ~side]
+        else:
+            parts = [m]
+        for j, part in enumerate(parts):
+            if part.sum() < 50:
+                continue
+            lat, lon = _centroid(part, LAT, LON)
+            seeds.append((lat, lon))
+            plates.append((f"continent-{k + 1}{'ab'[j] if len(parts) > 1 else ''}", lat, lon, "continental",
+                           rng.uniform(0, 360), rng.uniform(2.0, 5.0)))
+    for k, (lat, lon) in enumerate(_spread(rng, 14 - len(plates) // 2, 28.0, seeds, max_lat=85.0)):
+        plates.append((f"ocean-{k + 1}", lat, lon, "oceanic", rng.uniform(0, 360), rng.uniform(3.0, 8.0)))
+
+    ocean = np.argwhere(~ndimage.binary_dilation(land, iterations=10))
+    def sea_point():
+        y, x = ocean[rng.integers(len(ocean))]
+        return float(LAT[y, x]), float(LON[y, x])
+    f = lambda v: f"[{v[0]:.1f}, {v[1]:.1f}]"
+    t = [f"# Generated by `mapgen.py new-world --seed {seed}`. Edit freely (see README.md → Tectonics).",
+         "# seed = [lat, lon] where the plate grows from; motion = [azimuth° clockwise from north, speed cm/yr].", ""]
+    for pid, lat, lon, kind, az, sp in plates:
+        t += ["[[plate]]", f'id = "{pid}"', f"seed = {f((lat, lon))}", f'kind = "{kind}"', f"motion = [{az:.0f}.0, {sp:.1f}]", ""]
+    if len(ocean):
+        t += ["[[hotspot]]", 'name = "island-chain"', f"center = {f(sea_point())}", "length_km = 600.0", ""]
+    (cfg / "tectonics.toml").write_text("\n".join(t))
+
+    print(f"new-world: {len(plates)} plates, sketch/ and config/tectonics.toml written (seed {seed})")