"""`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})")