"""Routing graph climb values per edge (docs/offroad/PLAN.md R10.1b, Appendix E.7; platform since B13.2: every graph).

Each edge is sampled at even steps of at most SAMPLE_M (30 m) from a to b, both ends included, from the land terrain
archive (Mapterhorn, terrarium, z11, read through tools/geodata/lib/pmtiles_read.py with an LRU of decoded tiles), bilinear between
pixel centres across tile edges. The profile is smoothed with a 3-sample median (ends kept), then:

  up / dn  metres climbed / descended from a to b (0.1 m), counting a climb or descent only once it passes
           STEP_M (3 m) from the last counted point: the rule of src/web/src/features/drawing/profile.ts, so DEM
           noise under 3 m peak to peak never adds up. At b, the change since the last counted point is added when it
           carries on in the direction last counted, so a steady grade isn't cut short at every node. An edge whose
           climb never reaches 3 m counts 0, so a route's total climb is computed from its whole profile (R10.5c), not
           the sum of its edges' `up`;
  gmax     the steepest grade (%, 0.1) between samples at least GRADE_M (100 m) apart; an edge shorter than that uses its
           whole length, and one shorter than a DEM pixel (30 m) has none.

An edge with any sample outside the archive gets none of the three (counted in report.json). Values are derived from
Copernicus GLO-30 data: the graph's credits carry LAND_CREDIT (tools/geodata/lib/terrarium.py).
"""
import collections
import math
import os
import sys

import numpy as np

HERE = os.path.dirname(os.path.abspath(__file__))
sys.path.insert(0, os.path.join(HERE, '..', '..', 'lib'))  # tools/geodata/lib
import terrarium  # noqa: E402

R = 6371008.8
SAMPLE_M = 30.0
STEP_M = 3.0
GRADE_M = 100.0
MIN_GRADE_M = 30.0
ZOOM = terrarium.LAND_MAXZOOM
TILE = 512
CACHE_TILES = 128


class Terrain:
    """Elevation at points from tiles(z, x, y) → elevation grid (TILE × TILE, m) or None where there's no tile."""

    def __init__(self, tiles, zoom=ZOOM, size=TILE, cache=CACHE_TILES):
        self.tiles, self.zoom, self.size, self.cache_size = tiles, zoom, size, cache
        self.cache = collections.OrderedDict()

    def tile(self, x, y):
        k = (x, y)
        if k in self.cache:
            self.cache.move_to_end(k)
            return self.cache[k]
        t = self.tiles(self.zoom, x, y)
        if t is not None and t.shape != (self.size, self.size):
            raise ValueError(f'terrain tile {self.zoom}/{x}/{y} is {t.shape[1]}×{t.shape[0]} px, expected {self.size}')
        self.cache[k] = t
        if len(self.cache) > self.cache_size:
            self.cache.popitem(last=False)
        return t

    def pixels(self, ix, iy):
        """Values of global pixels (int arrays), NaN where the archive has no tile."""
        n = self.size << self.zoom
        ix, iy = np.clip(ix, 0, n - 1), np.clip(iy, 0, n - 1)
        tx, ty = ix // self.size, iy // self.size
        out = np.full(ix.shape, np.nan)
        for x, y in sorted(set(zip(tx.tolist(), ty.tolist()))):
            t = self.tile(x, y)
            if t is not None:
                m = (tx == x) & (ty == y)
                out[m] = t[iy[m] - y * self.size, ix[m] - x * self.size]
        return out

    def at(self, lon, lat):
        """Bilinear elevation (m) at each (lon, lat)."""
        px, py = terrarium.world_pixels(lon, lat, self.zoom, self.size)
        px, py = px - 0.5, py - 0.5
        x0, y0 = np.floor(px).astype(np.int64), np.floor(py).astype(np.int64)
        fx, fy = px - x0, py - y0
        v00, v10 = self.pixels(x0, y0), self.pixels(x0 + 1, y0)
        v01, v11 = self.pixels(x0, y0 + 1), self.pixels(x0 + 1, y0 + 1)
        return (v00 * (1 - fx) + v10 * fx) * (1 - fy) + (v01 * (1 - fx) + v11 * fx) * fy


def pmtiles_terrain(path=terrarium.LAND):
    """A Terrain over a terrarium PMTiles archive (GDAL decodes the WebP tiles; the geotools image has it)."""
    import pmtiles_read
    archive = pmtiles_read.Reader(path)

    def tiles(z, x, y):
        data = archive.read(z, x, y)
        return None if data is None else terrarium.decode(data)

    return Terrain(tiles)


def haversine(a, b):
    lon1, lat1, lon2, lat2 = (np.radians(v) for v in (a[:, 0], a[:, 1], b[:, 0], b[:, 1]))
    h = np.sin((lat2 - lat1) / 2) ** 2 + np.cos(lat1) * np.cos(lat2) * np.sin((lon2 - lon1) / 2) ** 2
    return 2 * R * np.arcsin(np.sqrt(np.minimum(h, 1.0)))


def samples(coords, step=SAMPLE_M):
    """(distance along, lon, lat) at even steps of at most `step` m along a line, both ends included."""
    ll = np.asarray(coords, np.float64)[:, :2]
    cum = np.concatenate([[0.0], np.cumsum(haversine(ll[:-1], ll[1:]))])
    total = cum[-1]
    d = np.linspace(0.0, total, max(1, math.ceil(total / step - 1e-9)) + 1) if total > 0 else np.zeros(1)
    lon = np.interp(d, cum, ll[:, 0])
    lat = np.interp(d, cum, ll[:, 1])
    return d, lon, lat


def median3(e):
    if len(e) < 3:
        return e.copy()
    padded = np.concatenate([e[:1], e, e[-1:]])
    return np.median(np.stack([padded[:-2], padded[1:-1], padded[2:]]), axis=0)


def climb(e, step=STEP_M):
    """(up, dn) counted only once a change passes `step` m from the last counted point (profile.ts), plus the change
    since the last counted point at the end when it carries on in the direction last counted."""
    up = dn = 0.0
    anchor, last = e[0], 0
    for v in e[1:]:
        if v - anchor >= step:
            up, anchor, last = up + v - anchor, v, 1
        elif anchor - v >= step:
            dn, anchor, last = dn + anchor - v, v, -1
    rest = e[-1] - anchor
    if last > 0 and rest > 0:
        up += rest
    elif last < 0 and rest < 0:
        dn -= rest
    return up, dn


def max_grade(d, e, window=GRADE_M):
    """Steepest |Δe| / Δd (%) between samples at least `window` m apart (the whole edge when shorter), or None."""
    total = d[-1]
    if total < MIN_GRADE_M:
        return None
    if total < window:
        return abs(e[-1] - e[0]) / total * 100
    j = np.searchsorted(d, d + window - 1e-9)
    ok = j < len(d)
    i, j = np.nonzero(ok)[0], j[ok]
    return float(np.max(np.abs(e[j] - e[i]) / (d[j] - d[i]))) * 100


def edge_climb(coords, terrain):
    """{up, dn, gmax} for an edge's coordinates from a to b, or {} where the terrain has a hole."""
    d, lon, lat = samples(coords)
    raw = terrain.at(lon, lat)
    if np.isnan(raw).any():
        return {}
    e = median3(raw)
    up, dn = climb(e)
    out = {'up': round(up, 1), 'dn': round(dn, 1)}
    g = max_grade(d, e)
    if g is not None:
        out['gmax'] = round(g, 1)
    return out