#!/usr/bin/env python3
"""Construit les GeoJSON de la carte à partir de l'extraction OSM et des calques
uMap du collectif « Le Nuage était sous nos pieds ».

Usage : build.py [zone]   (france par défaut)
Entrées : data/<zone>/dc-osm.geojsonseq (scripts/extract_osm.sh), data/nuage/*.json (France),
          data/ref/ne_50m_admin_0_countries.geojson (Natural Earth)
Sorties : web/data/<zone>/osm.geojson, nuage.geojson (France), stats.json
"""
import argparse
import json
import math
import re
import shutil
from datetime import date
from collections import Counter, defaultdict
from pathlib import Path

ROOT = Path(__file__).resolve().parent.parent
DATA = ROOT / "data"
COUNTRIES = DATA / "ref" / "ne_50m_admin_0_countries.geojson"
GRID = 0.05  # taille de maille (degrés) de l'index spatial

# Distance max (m) pour apparier un point du Nuage à un objet OSM
MATCH_RADIUS = 400
# Niveaux supposés quand building:levels est absent
DEFAULT_LEVELS = 2
# Part bâtie d'une emprise industrial=data_center sans bâtiment cartographié
LANDUSE_BUILT_RATIO = 0.4
POWER_KEYS = ("data_center:power", "power_supply:capacity", "power_rating",
              "power:input", "capacity:powered", "capacity:power_supply")

LAYER_STATUS = {
    "Data centers existant": "existant",
    "Les data centers en construction": "construction",
    "Projet de data center validés, en phase d'instruction": "instruction",
    "Projet de data center annoncé": "annonce",
    "Sites \"clef en main\"": "site_cle_en_main",
    "Les oppositions aux data centers": "opposition",
    "Recours juridique": "recours",
}


def is_dc(tags):
    return (tags.get("telecom") == "data_center"
            or tags.get("building") == "data_center"
            or tags.get("industrial") == "data_center"
            or tags.get("man_made") == "data_center"
            or "data_center" in tags)


def ring_area(ring):
    """Aire (m²) d'un anneau lon/lat, projection équirectangulaire locale."""
    if len(ring) < 4:
        return 0.0
    lat0 = math.radians(sum(p[1] for p in ring) / len(ring))
    k = 6371008.8 * math.pi / 180
    pts = [(p[0] * k * math.cos(lat0), p[1] * k) for p in ring]
    return abs(sum(x1 * y2 - x2 * y1 for (x1, y1), (x2, y2) in zip(pts, pts[1:]))) / 2


def polygons(geom):
    if geom["type"] == "Polygon":
        return [geom["coordinates"]]
    if geom["type"] == "MultiPolygon":
        return geom["coordinates"]
    return []


def area(geom):
    return sum(ring_area(p[0]) - sum(ring_area(h) for h in p[1:]) for p in polygons(geom))


def centroid(geom):
    if geom["type"] == "Point":
        return tuple(geom["coordinates"][:2])
    rings = [p[0] for p in polygons(geom)]
    pts = [pt for r in rings for pt in r[:-1]]
    return (sum(p[0] for p in pts) / len(pts), sum(p[1] for p in pts) / len(pts))


def point_in_ring(x, y, ring):
    inside = False
    for (x1, y1), (x2, y2) in zip(ring, ring[1:]):
        if (y1 > y) != (y2 > y) and x < (x2 - x1) * (y - y1) / (y2 - y1) + x1:
            inside = not inside
    return inside


def point_in_geom(pt, geom):
    return any(point_in_ring(pt[0], pt[1], p[0]) and
               not any(point_in_ring(pt[0], pt[1], h) for h in p[1:])
               for p in polygons(geom))


def bbox(geom):
    pts = [pt for p in polygons(geom) for pt in p[0]]
    return (min(p[0] for p in pts), min(p[1] for p in pts),
            max(p[0] for p in pts), max(p[1] for p in pts))


class Grid:
    """Index spatial minimal : polygones rangés par mailles de leur emprise."""

    def __init__(self, items, size=GRID):
        self.size = size
        self.cells = defaultdict(list)
        for it in items:
            x0, y0, x1, y1 = it["bb"]
            for i in range(int(math.floor(x0 / size)), int(math.floor(x1 / size)) + 1):
                for j in range(int(math.floor(y0 / size)), int(math.floor(y1 / size)) + 1):
                    self.cells[(i, j)].append(it)

    def at(self, pt):
        return self.cells.get((int(math.floor(pt[0] / self.size)), int(math.floor(pt[1] / self.size))), [])

    def near(self, pt, radius_deg):
        i0, j0 = int(math.floor(pt[0] / self.size)), int(math.floor(pt[1] / self.size))
        r = int(math.ceil(radius_deg / self.size))
        seen = set()
        for i in range(i0 - r, i0 + r + 1):
            for j in range(j0 - r, j0 + r + 1):
                for it in self.cells.get((i, j), []):
                    if id(it) not in seen:
                        seen.add(id(it))
                        yield it


def dist_m(a, b):
    lat = math.radians((a[1] + b[1]) / 2)
    dx = (a[0] - b[0]) * math.cos(lat) * 111320
    dy = (a[1] - b[1]) * 110540
    return math.hypot(dx, dy)


def parse_mw(text):
    """Valeurs de puissance (MW) citées dans un texte libre."""
    vals = []
    for num, unit in re.findall(r"(\d+(?:[.,]\d+)?)\s*(GW|MW|mégawatts?|kW)\b", text or "", re.I):
        v = float(num.replace(",", "."))
        unit = unit.lower()
        vals.append(v * 1000 if unit == "gw" else v / 1000 if unit == "kw" else v)
    return vals


def parse_tag_power(tags):
    for k in POWER_KEYS:
        if k in tags:
            vals = parse_mw(tags[k]) or parse_mw(tags[k] + " MW")
            if vals:
                return vals[0], k
    return None, None


def field(desc, *labels):
    for label in labels:
        m = re.search(rf"^\s*{label}s?\s*:?\s*:\s*(.+)$", desc, re.I | re.M)
        if m:
            return m.group(1).strip(" .")
    return None


def load_osm(zone):
    feats = []
    with open(DATA / zone / "dc-osm.geojsonseq") as f:
        for line in f:
            line = line.strip("\x1e\n ")
            if not line:
                continue
            ft = json.loads(line)
            tags = ft["properties"]
            osm_type = tags.pop("@type")
            osm_id = tags.pop("@id")
            meta = {k[1:]: tags.pop(k) for k in [k for k in tags if k.startswith("@")]}
            if not is_dc(tags):
                continue
            g = ft["geometry"]
            if g["type"] != "Point" and not polygons(g):
                continue
            feats.append({"type": "Feature", "geometry": g,
                          "id_osm": f"{osm_type}/{osm_id}", "tags": tags, "meta": meta,
                          "c": centroid(g), "area": area(g) if g["type"] != "Point" else 0})
    return feats


def dedupe(feats):
    """Marque les doublons : point dans un polygone, polygone englobant d'autres."""
    polys = [f for f in feats if f["area"] > 0]
    for f in polys:
        f["bb"] = bbox(f["geometry"])
    grid = Grid(polys)

    def containers(pt, skip):
        return [p for p in grid.at(pt) if p is not skip
                and p["bb"][0] <= pt[0] <= p["bb"][2] and p["bb"][1] <= pt[1] <= p["bb"][3]
                and point_in_geom(pt, p["geometry"])]

    for f in feats:
        f["role"] = "principal"
    for f in feats:
        if f["area"] == 0:
            if containers(f["c"], f):
                f["role"] = "doublon_point"
        else:
            for outer in containers(f["c"], f):
                if outer["area"] > f["area"]:
                    outer["role"] = "emprise"


def estimate(f):
    tags = f["tags"]
    mw, key = parse_tag_power(tags)
    floor = None
    if f["area"] > 0:
        try:
            levels = max(1, float(tags.get("building:levels", "").replace(",", ".")))
            levels_src = "tag"
        except ValueError:
            levels, levels_src = DEFAULT_LEVELS, "défaut"
        if "building" in tags:
            floor = f["area"] * levels
        else:
            floor = f["area"] * LANDUSE_BUILT_RATIO * levels
        f["levels"] = levels
        f["levels_src"] = levels_src
    f["mw_tag"] = mw
    f["mw_key"] = key
    f["floor"] = floor


def load_nuage():
    meta = json.load(open(DATA / "nuage" / "map.json"))
    names = {l["id"]: l["properties"]["name"].strip() for l in meta["properties"]["datalayers"]}
    out = []
    for uid, name in names.items():
        status = next((v for k, v in LAYER_STATUS.items() if name.startswith(k)), None)
        if status is None:
            continue
        layer = json.load(open(DATA / "nuage" / f"layer-{uid}.json"))
        for ft in layer["features"]:
            p = ft["properties"]
            desc = p.get("description") or ""
            carac = field(desc, "Caractéristique") or ""
            mws = parse_mw(carac) or parse_mw(desc)
            m2 = re.search(r"(\d[\d\s.]*)\s*m(?:²|2)\b", carac or desc)
            out.append({
                "geometry": ft["geometry"],
                "c": tuple(ft["geometry"]["coordinates"][:2]),
                "status": status,
                "layer": name,
                "name": (p.get("name") or "").strip(),
                "description": desc.strip(),
                "operator": p.get("operator") or field(desc, "Entreprise"),
                "statut_txt": field(desc, "Statut"),
                "date": field(desc, "Date"),
                "osm_ref": p.get("@id") if str(p.get("@id", "")).startswith(("node/", "way/", "relation/")) else None,
                "mw": mws[0] if mws else None,
                "surface_m2": int(re.sub(r"\D", "", m2.group(1))) if m2 else None,
                "sources": re.findall(r"https?://\S+", desc),
            })
    return out


def load_countries():
    feats = json.load(open(COUNTRIES))["features"]
    out = []
    for ft in feats:
        p = ft["properties"]
        iso = p.get("ISO_A2_EH") if p.get("ISO_A2_EH") not in (None, "-99") else p.get("ADM0_A3")
        c = {"geometry": ft["geometry"], "iso": iso, "name": p.get("NAME_FR") or p.get("NAME"),
             "continent": p.get("CONTINENT")}
        c["bb"] = bbox(c["geometry"])
        out.append(c)
    return out, Grid(out, size=1.0)


def locate(pt, countries, grid):
    """Pays contenant le point, sinon le plus proche (côtes simplifiées au 1:50 M)."""
    for c in grid.at(pt):
        if point_in_geom(pt, c["geometry"]):
            return c
    best, bd = None, None
    for c in grid.near(pt, 2.0):
        for poly in polygons(c["geometry"]):
            for q in poly[0][::4]:
                d = dist_m(pt, q)
                if bd is None or d < bd:
                    best, bd = c, d
    return best


def suggest_tags(n, osm):
    """Tags qu'une conflation pourrait ajouter à l'objet OSM (ou à un nouvel objet)."""
    tags = osm["tags"] if osm else {}
    s = {}
    if n["status"] in ("existant", "construction") and not osm:
        s["telecom"] = "data_center"
        if n["status"] == "construction":
            s["building"] = "construction"
            s["construction"] = "data_center"
    if n["name"] and "name" not in tags and n["status"] == "existant":
        s["name"] = n["name"]
    if n["operator"] and "operator" not in tags:
        s["operator"] = n["operator"]
    if n["mw"] and not (osm and osm["mw_tag"]):
        s["data_center:power"] = f"{n['mw']:g} MW"
    year = re.search(r"\b(20\d\d)\b", n["date"] or "")
    if year and n["status"] == "construction" and "opening_date" not in tags:
        s["opening_date"] = year.group(1)
    if n["sources"] and s:
        s["source:data_center:power" if "data_center:power" in s else "source"] = n["sources"][0]
    return s


def main():
    ap = argparse.ArgumentParser()
    ap.add_argument("zone", nargs="?", default="france")
    ap.add_argument("--nuage", action=argparse.BooleanOptionalAction, default=None,
                    help="conflation avec la carte du Nuage (défaut : seulement pour france)")
    args = ap.parse_args()
    zone = args.zone
    with_nuage = args.nuage if args.nuage is not None else zone == "france"
    OUT = ROOT / "web" / "data" / zone
    OUT.mkdir(parents=True, exist_ok=True)

    osm = load_osm(zone)
    dedupe(osm)
    countries, cgrid = load_countries()
    for f in osm:
        estimate(f)
        c = locate(f["c"], countries, cgrid)
        f["country"] = c["iso"] if c else None
        f["country_name"] = c["name"] if c else None
        f["continent"] = c["continent"] if c else None
    by_ref = {f["id_osm"]: f for f in osm}
    principal = [f for f in osm if f["role"] == "principal"]

    nuage = load_nuage() if with_nuage else []
    for f in osm:
        f["bb"] = bbox(f["geometry"]) if f["area"] else (*f["c"], *f["c"])
    match_grid = Grid(osm)
    dc_status = ("existant", "construction", "instruction", "annonce")
    for n in nuage:
        match, d = None, None
        if n["osm_ref"] and n["osm_ref"] in by_ref:
            match, d = by_ref[n["osm_ref"]], 0
        elif n["status"] in dc_status:
            # objets principaux d'abord, emprises de campus en dernier recours
            around = list(match_grid.near(n["c"], MATCH_RADIUS / 70000))
            for pool in ([f for f in around if f["role"] == "principal"],
                         [f for f in around if f["role"] == "emprise"]):
                for f in pool:
                    if f["area"] and point_in_geom(n["c"], f["geometry"]):
                        dd = 0
                    else:
                        dd = dist_m(n["c"], f["c"])
                    if dd <= MATCH_RADIUS and (d is None or dd < d):
                        match, d = f, dd
                if match:
                    break
        n["match"] = match["id_osm"] if match else None
        n["match_dist"] = round(d) if d is not None else None
        n["suggest"] = suggest_tags(n, match) if n["status"] in dc_status else {}
        if match and n["name"]:
            match.setdefault("nuage_names", []).append(n["name"])
        if n["status"] in ("existant", "construction"):
            n["conflation"] = "completer" if match and n["suggest"] else "ok" if match else "absent_osm"
        elif n["status"] in ("instruction", "annonce"):
            n["conflation"] = "projet_non_cartographiable"
        else:
            n["conflation"] = None
        if match and n["mw"]:
            match.setdefault("nuage_mw", n["mw"])

    # plusieurs points du Nuage sur un même objet : ne pas proposer de nom unique
    for n in nuage:
        m = by_ref.get(n["match"])
        if m and len(m.get("nuage_names", [])) > 1 and "name" in n["suggest"]:
            n["suggest"].pop("name")
            n["suggest"]["note"] = "Plusieurs sites du Nuage dans cet objet : " + ", ".join(m["nuage_names"])
            if n["conflation"] == "ok":
                n["conflation"] = "completer"

    # --- GeoJSON OSM
    osm_out = []
    for f in osm:
        t = f["tags"]
        props = {
            "id": f["id_osm"], "name": t.get("name"), "operator": t.get("operator"),
            "role": f["role"], "area_m2": round(f["area"]) or None,
            "floor_m2": round(f["floor"]) if f["floor"] else None,
            "levels": f.get("levels"), "levels_src": f.get("levels_src"),
            "mw_tag": f["mw_tag"], "mw_key": f["mw_key"], "nuage_mw": f.get("nuage_mw"),
            "in_nuage": bool(f.get("nuage_names")),
            "tags": json.dumps(t, ensure_ascii=False),
            "postcode": t.get("addr:postcode"),
            "country": f["country"],
            "timestamp": f["meta"].get("timestamp"),
            "user": f["meta"].get("user"),
        }
        osm_out.append({"type": "Feature", "geometry": f["geometry"], "properties": props})
    json.dump({"type": "FeatureCollection", "features": osm_out},
              open(OUT / "osm.geojson", "w"), ensure_ascii=False, separators=(",", ":"))

    nuage_out = []
    for n in nuage:
        props = {k: n[k] for k in ("status", "layer", "name", "description", "operator",
                                    "statut_txt", "date", "mw", "surface_m2", "match",
                                    "match_dist", "conflation")}
        props["suggest"] = json.dumps(n["suggest"], ensure_ascii=False)
        props["n_suggest"] = len(n["suggest"])
        nuage_out.append({"type": "Feature", "geometry": n["geometry"], "properties": props})
    if with_nuage:
        json.dump({"type": "FeatureCollection", "features": nuage_out},
                  open(OUT / "nuage.geojson", "w"), ensure_ascii=False, separators=(",", ":"))

    # --- Statistiques (la puissance estimée par surface est recalculée côté client)
    floors = sorted(f["floor"] for f in principal if f["floor"])
    ops = Counter((f["tags"].get("operator") or "inconnu").strip() for f in principal)
    dept = Counter((f["tags"].get("addr:postcode") or "")[:2] or "?" for f in principal)
    status_mw = defaultdict(lambda: {"count": 0, "with_mw": 0, "mw": 0.0})
    for n in nuage:
        s = status_mw[n["status"]]
        s["count"] += 1
        if n["mw"]:
            s["with_mw"] += 1
            s["mw"] += n["mw"]
    tag_counter = Counter()
    for f in osm:
        for k in ("telecom", "building", "industrial", "man_made"):
            if f["tags"].get(k) == "data_center":
                tag_counter[f"{k}=data_center"] += 1
    key_fill = {k: sum(1 for f in principal if k in f["tags"])
                for k in ("name", "operator", "building:levels", "start_date", "website",
                          "operator:wikidata", "addr:street") + POWER_KEYS}
    conf = Counter(n["conflation"] for n in nuage if n["conflation"])
    # agrégats par pays : nombre, plancher, puissance taguée ; la puissance estimée
    # dépend de la densité choisie et se recalcule côté client
    countries_stats = defaultdict(lambda: {"name": None, "continent": None, "count": 0,
                                           "floor_m2": 0, "points_only": 0, "mw_tag": 0.0})
    for f in principal:
        cs = countries_stats[f["country"] or "?"]
        cs["name"] = f["country_name"] or "Inconnu"
        cs["continent"] = f["continent"]
        cs["count"] += 1
        if f["floor"]:
            cs["floor_m2"] += round(f["floor"])
        else:
            cs["points_only"] += 1
        if f["mw_tag"]:
            cs["mw_tag"] += f["mw_tag"]
    date_file = DATA / zone / "extract-date"
    stats = {
        "zone": zone,
        "generated": date.today().isoformat(),
        "osm_extract": date_file.read_text().strip() if date_file.exists() else None,
        "has_nuage": with_nuage,
        "countries": dict(sorted(countries_stats.items(), key=lambda kv: -kv[1]["count"])),
        "osm_total": len(osm),
        "osm_principal": len(principal),
        "osm_doublons": sum(1 for f in osm if f["role"] != "principal"),
        "osm_with_area": len(floors),
        "osm_points_only": sum(1 for f in principal if not f["floor"]),
        "floor_total_m2": round(sum(floors)),
        "floor_median_m2": round(floors[len(floors) // 2]) if floors else 0,
        "osm_mw_tagged": sum(f["mw_tag"] for f in principal if f["mw_tag"]),
        "osm_n_mw_tagged": sum(1 for f in principal if f["mw_tag"]),
        "tags": dict(tag_counter),
        "key_fill": key_fill,
        "top_operators": ops.most_common(12),
        "top_depts": dept.most_common(10),
        "nuage_status": status_mw,
        "conflation": dict(conf),
        "matched_mw": sum(f["nuage_mw"] for f in principal if f.get("nuage_mw")),
        "params": {"default_levels": DEFAULT_LEVELS, "landuse_ratio": LANDUSE_BUILT_RATIO,
                   "match_radius_m": MATCH_RADIUS},
    }
    json.dump(stats, open(OUT / "stats.json", "w"), ensure_ascii=False, indent=1)
    changes = DATA / zone / "changes.json"
    if changes.exists():
        shutil.copy(changes, OUT / "changes.json")
    print(json.dumps({k: v for k, v in stats.items() if k not in ("top_operators", "top_depts")},
                     ensure_ascii=False, indent=1))


if __name__ == "__main__":
    main()