"""
Calcul de la densité de population maximale dans l'empreinte du GRC (Ground Risk
Class) à partir d'un fichier KML de mission et de la grille de données carroyées
INSEE Filosofi (200 m).

Source des données : Insee, "Données au carreau de 200m (y compris données
imputées)" — Filosofi 2021.
Page : https://www.insee.fr/fr/statistiques/8735162
Fichier (shapefile, ~240 Mo, couvre France métropolitaine + Martinique + La
Réunion en un seul zip) :
https://www.insee.fr/fr/statistiques/fichier/8735162/Filosofi2021_carreaux_200m_shp.zip
Documentation des champs : https://www.insee.fr/fr/statistiques/fichier/8735106/documentation_donnees-carroyees_filosofi2021.pdf

Champs utilisés :
- "ind"        : nombre d'individus dans le carreau (200 m x 200 m)
- "idcar_200m" : identifiant du carreau
- projection France métropolitaine : Lambert 93 (EPSG:2154)

⚠️ IMPORTANT — Limite de l'environnement de développement :
Le téléchargement (insee.fr) est bloqué par l'allowlist réseau du bac à sable où
ce module a été écrit et ne peut donc pas être testé ici avec les vraies données.
La logique géométrique (parsing KML, intersection, pondération par recouvrement)
a été validée avec des données synthétiques (voir test_insee_density.py). Sur le
poste/serveur de déploiement (accès internet normal), une préparation explicite
(cf. `prepare_insee_data()` / bouton "Préparer les données INSEE" dans l'appli)
télécharge le fichier (~240 Mo) puis le convertit UNE FOIS en GeoPackage indexé
spatialement (index R-tree) : c'est cette étape de conversion qui rend les
calculs suivants rapides. Sans elle, le shapefile brut (2,3 millions de
carreaux, sans index) doit être scanné intégralement à CHAQUE calcul, ce qui
peut prendre de nombreuses minutes et donne l'impression que l'appli est bloquée.
"""
import json
import os
import re
import threading
import time
import zipfile
import xml.etree.ElementTree as ET

import requests
from shapely.geometry import Polygon, MultiPolygon
from shapely.ops import unary_union
import shapely

BASE_DIR = os.path.dirname(os.path.abspath(__file__))
CACHE_DIR = os.path.join(BASE_DIR, "data", "insee_cache")
os.makedirs(CACHE_DIR, exist_ok=True)

INSEE_ZIP_URL = "https://www.insee.fr/fr/statistiques/fichier/8735162/Filosofi2021_carreaux_200m_shp.zip"
INSEE_ZIP_PATH = os.path.join(CACHE_DIR, "Filosofi2021_carreaux_200m_shp.zip")
INSEE_EXTRACT_DIR = os.path.join(CACHE_DIR, "extracted")
INSEE_INDEXED_PATH = os.path.join(CACHE_DIR, "filosofi_200m_indexed.gpkg")
STATUS_PATH = os.path.join(CACHE_DIR, "status.json")

# Jeu de données utilisé pour tous les calculs de densité (source unique de
# vérité, réutilisée dans l'interface et en pied de page de la fiche).
# Millésime le plus récent publié par l'INSEE (diffusion février 2026).
INSEE_MILLESIME = "Filosofi 2021"
INSEE_DATASET_LABEL = "INSEE Filosofi 2021 — carreaux 200 m (diffusion février 2026)"
INSEE_DATASET_URL = "https://www.insee.fr/fr/statistiques/8735162"
INSEE_FOOTER_NOTE = "Densités de population estimées d'après les données carroyées INSEE Filosofi 2021 (carreaux 200 m)."

LAMBERT93 = "EPSG:2154"
WGS84 = "EPSG:4326"

KML_NS = {"kml": "http://www.opengis.net/kml/2.2"}

_lock = threading.Lock()


class InseeDataUnavailable(Exception):
    """Levée quand les données carroyées INSEE ne sont ni en cache, ni
    téléchargeables (pas de réseau, allowlist, erreur serveur...)."""


class InseeDataNotReady(Exception):
    """Levée quand les données ne sont pas encore préparées (téléchargement /
    indexation jamais lancés, ou en cours). Contient le statut courant."""
    def __init__(self, status):
        self.status = status
        super().__init__(status.get("message", "Données INSEE non prêtes."))


# ---------------------------------------------------------------------------
# 0. Statut de préparation (persisté sur disque, lisible par un autre process)
# ---------------------------------------------------------------------------

def _write_status(state, message="", progress=None):
    payload = {"state": state, "message": message, "progress": progress, "updated_at": time.time()}
    tmp = STATUS_PATH + ".tmp"
    with open(tmp, "w", encoding="utf-8") as f:
        json.dump(payload, f)
    os.replace(tmp, STATUS_PATH)
    return payload


def get_status():
    """État de préparation des données : not_started / downloading /
    converting / ready / error."""
    if os.path.exists(INSEE_INDEXED_PATH):
        return {"state": "ready", "message": "Données INSEE prêtes.", "progress": 100}
    if os.path.exists(STATUS_PATH):
        try:
            with open(STATUS_PATH, encoding="utf-8") as f:
                return json.load(f)
        except Exception:
            pass
    return {"state": "not_started", "message": "Données INSEE non préparées.", "progress": 0}


def prepare_insee_data(force=False):
    """Télécharge (si besoin) puis convertit le shapefile national en
    GeoPackage indexé spatialement. Bloquant — à appeler depuis un thread
    d'arrière-plan (voir app.py) plutôt que dans la requête HTTP de l'utilisateur."""
    if not _lock.acquire(blocking=False):
        return get_status()  # préparation déjà en cours dans un autre thread
    try:
        if os.path.exists(INSEE_INDEXED_PATH) and not force:
            return _write_status("ready", "Données INSEE déjà prêtes.", 100)

        try:
            _write_status("downloading", "Téléchargement du fichier INSEE (~240 Mo)...", 5)
            shp_path = _ensure_shapefile_downloaded()

            _write_status("converting", "Indexation spatiale (conversion en GeoPackage)...", 60)
            import geopandas as gpd
            gdf = gpd.read_file(shp_path)
            os.makedirs(os.path.dirname(INSEE_INDEXED_PATH), exist_ok=True)
            tmp_gpkg = INSEE_INDEXED_PATH + ".tmp"
            gdf.to_file(tmp_gpkg, driver="GPKG", layer="carreaux_200m")
            os.replace(tmp_gpkg, INSEE_INDEXED_PATH)

            return _write_status("ready", "Données INSEE prêtes.", 100)
        except Exception as e:
            return _write_status("error", f"Échec de préparation des données INSEE : {e}", None)
    finally:
        _lock.release()


def _looks_like_zip(path):
    """Vrai si le fichier commence par la signature d'une archive ZIP (PK...)."""
    try:
        with open(path, "rb") as f:
            return f.read(4)[:2] == b"PK"
    except Exception:
        return False


def _ensure_shapefile_downloaded():
    shp_path = _find_cached_shapefile()
    if shp_path:
        return shp_path

    # Un éventuel zip déjà présent mais invalide (ancien téléchargement
    # interrompu, page d'erreur HTML enregistrée à la place du zip...) est
    # supprimé pour être re-téléchargé proprement.
    if os.path.exists(INSEE_ZIP_PATH) and not _looks_like_zip(INSEE_ZIP_PATH):
        os.remove(INSEE_ZIP_PATH)

    if not os.path.exists(INSEE_ZIP_PATH):
        try:
            _download(INSEE_ZIP_URL, INSEE_ZIP_PATH)
        except Exception as e:
            raise InseeDataUnavailable(
                f"Impossible de télécharger les données INSEE ({INSEE_ZIP_URL}) : {e}. "
                "Vérifiez la connexion internet du serveur, ou téléchargez le fichier "
                "manuellement (voir README) et placez-le dans data/insee_cache/."
            )

    # Validation : le serveur INSEE peut renvoyer une page HTML (erreur, mur de
    # consentement...) au lieu du zip -> on le détecte avant de tenter l'extraction.
    if not _looks_like_zip(INSEE_ZIP_PATH):
        size = os.path.getsize(INSEE_ZIP_PATH) if os.path.exists(INSEE_ZIP_PATH) else 0
        try:
            os.remove(INSEE_ZIP_PATH)
        except OSError:
            pass
        raise InseeDataUnavailable(
            "Le fichier téléchargé depuis l'INSEE n'est pas une archive ZIP valide "
            f"(taille reçue : {size} octets ; il s'agit probablement d'une page d'erreur). "
            "Réessayez ; si le problème persiste, téléchargez manuellement "
            "« Carreau 200m – Shapefile » depuis https://www.insee.fr/fr/statistiques/8735162 "
            "et placez le .zip dans data/insee_cache/ (nom attendu : "
            f"{os.path.basename(INSEE_ZIP_PATH)}), puis relancez la préparation."
        )

    os.makedirs(INSEE_EXTRACT_DIR, exist_ok=True)
    try:
        with zipfile.ZipFile(INSEE_ZIP_PATH) as zf:
            zf.extractall(INSEE_EXTRACT_DIR)
    except zipfile.BadZipFile:
        os.remove(INSEE_ZIP_PATH)
        raise InseeDataUnavailable(
            "L'archive INSEE téléchargée est corrompue (téléchargement incomplet ?). "
            "Elle a été supprimée ; relancez la préparation pour re-télécharger."
        )

    shp_path = _find_cached_shapefile()
    if not shp_path:
        raise InseeDataUnavailable(
            "Le fichier INSEE a été téléchargé mais aucun shapefile métropole n'a été trouvé "
            f"dans {INSEE_EXTRACT_DIR}. Vérifiez le contenu de l'archive."
        )
    return shp_path


def ensure_insee_indexed():
    """Retourne le chemin du GeoPackage indexé, prêt pour des lectures bbox
    rapides. Lève InseeDataNotReady si la préparation n'a pas encore été
    lancée / est en cours / a échoué (le message précise l'état)."""
    status = get_status()
    if status["state"] == "ready" and os.path.exists(INSEE_INDEXED_PATH):
        return INSEE_INDEXED_PATH
    raise InseeDataNotReady(status)


def _find_cached_shapefile():
    if not os.path.isdir(INSEE_EXTRACT_DIR):
        return None
    candidates = []
    for root, _dirs, files in os.walk(INSEE_EXTRACT_DIR):
        for fn in files:
            if fn.lower().endswith(".shp"):
                candidates.append(os.path.join(root, fn))
    if not candidates:
        return None
    # Écarte les fichiers Martinique/Réunion si identifiables par le nom
    metro = [c for c in candidates if not re.search(r"(martinique|reunion|971|972|974)", c, re.I)]
    return (metro or candidates)[0]


def _download(url, dest_path, chunk_size=1024 * 1024):
    # En-têtes "navigateur" : sans User-Agent, insee.fr peut renvoyer une page
    # d'erreur/consentement au lieu du fichier. On suit aussi les redirections.
    headers = {
        "User-Agent": "Mozilla/5.0 (Windows NT 10.0; Win64; x64) "
                      "AppleWebKit/537.36 (KHTML, like Gecko) Chrome/124.0 Safari/537.36",
        "Accept": "application/zip,application/octet-stream,*/*",
    }
    with requests.get(url, stream=True, timeout=120, headers=headers, allow_redirects=True) as r:
        r.raise_for_status()
        tmp = dest_path + ".part"
        with open(tmp, "wb") as f:
            for chunk in r.iter_content(chunk_size=chunk_size):
                if chunk:
                    f.write(chunk)
        os.replace(tmp, dest_path)


# ---------------------------------------------------------------------------
# 2. Parsing du KML de mission
# ---------------------------------------------------------------------------

def _parse_coords(text):
    coords = []
    for token in text.strip().split():
        parts = token.split(",")
        lon, lat = float(parts[0]), float(parts[1])
        coords.append((lon, lat))
    return coords


def extract_polygons_from_kml(kml_path):
    """Retourne une liste de (nom_placemark, shapely.Polygon en WGS84)."""
    tree = ET.parse(kml_path)
    root = tree.getroot()
    results = []
    for placemark in root.iter("{http://www.opengis.net/kml/2.2}Placemark"):
        name_el = placemark.find("kml:name", KML_NS)
        name = name_el.text.strip() if name_el is not None and name_el.text else ""
        for poly_el in placemark.iter("{http://www.opengis.net/kml/2.2}Polygon"):
            outer = poly_el.find(".//kml:outerBoundaryIs/kml:LinearRing/kml:coordinates", KML_NS)
            if outer is None or not outer.text:
                continue
            ring = _parse_coords(outer.text)
            if len(ring) < 3:
                continue
            holes = []
            for inner in poly_el.findall(".//kml:innerBoundaryIs/kml:LinearRing/kml:coordinates", KML_NS):
                if inner.text:
                    holes.append(_parse_coords(inner.text))
            results.append((name, Polygon(ring, holes)))
    return results


BUFFER_KEYWORDS = ("buffer", "tampon", "grc", "adjacent", "espace adjacent",
                   "zone tampon", "ground risk")


def _outermost_polygon(polygons):
    """Sélectionne le polygone le plus EXTÉRIEUR d'un jeu de zones emboîtées
    (géographie de vol ⊂ contingence ⊂ buffer). Le buffer est, par
    construction, le tracé le plus éloigné du polygone d'origine : on prend
    donc en priorité un polygone qui contient (ou contient presque) tous les
    autres ; à défaut, celui d'aire maximale.

    `polygons` = liste de (nom, shapely.Polygon)."""
    geoms = [p for _, p in polygons]
    if len(geoms) == 1:
        return geoms[0]

    # Candidat "englobant" : celui qui recouvre le plus les autres.
    best = None
    best_cover = -1.0
    others_union_area = unary_union(geoms).area
    for g in geoms:
        # part de l'union totale couverte par ce seul polygone
        cover = g.area / others_union_area if others_union_area > 0 else 0.0
        if cover > best_cover:
            best_cover = cover
            best = g
    # Si le meilleur candidat couvre l'essentiel (>= 95%) de l'union, c'est le
    # contour extérieur : on le renvoie. Sinon (zones disjointes, ex: plusieurs
    # parcelles distinctes), on renvoie l'union de toutes les zones.
    if best_cover >= 0.95:
        return best
    return unary_union(geoms)


def find_grc_zone(kml_path, keywords=BUFFER_KEYWORDS):
    """Zone de référence pour la densité MAXIMALE au sol (section B).

    Le KML fourni a déjà été enrichi par l'outil de préparation des vols : au
    polygone d'origine (géographie de vol) se sont ajoutés la contingence puis
    le buffer. C'est la zone la plus EXTÉRIEURE (le buffer) qu'il faut retenir.

    Stratégie : si un placemark porte un nom évoquant le buffer/tampon, on le
    prend (le nom exact varie selon le logiciel, d'où une liste de mots-clés) ;
    sinon on sélectionne géométriquement le polygone le plus extérieur."""
    polygons = extract_polygons_from_kml(kml_path)
    if not polygons:
        raise ValueError(f"Aucun polygone trouvé dans {kml_path}")

    for kw in keywords:
        matches = [p for name, p in polygons if kw.lower() in name.lower()]
        if matches:
            # plusieurs placemarks "buffer" -> on garde le plus extérieur d'entre eux
            return _outermost_polygon([(n, p) for n, p in polygons
                                       if kw.lower() in n.lower()])

    return _outermost_polygon(polygons)


def find_parcelle_zone(kml_path, parcelle_ref):
    """Cherche le(s) polygone(s) dont le nom de placemark contient la référence
    de parcelle donnée (ex: '396'). À défaut, retombe sur find_grc_zone (toute
    la zone de mission)."""
    polygons = extract_polygons_from_kml(kml_path)
    if not polygons:
        raise ValueError(f"Aucun polygone trouvé dans {kml_path}")
    if parcelle_ref:
        named = [(name, p) for name, p in polygons if str(parcelle_ref).lower() in name.lower()]
        if named:
            # zone extérieure (buffer) de cette parcelle précise
            return _outermost_polygon(named)
    return find_grc_zone(kml_path)


def _to_lambert93(geom_wgs84):
    from shapely.ops import transform
    from pyproj import Transformer
    to_l93 = Transformer.from_crs(WGS84, LAMBERT93, always_xy=True).transform
    return transform(to_l93, geom_wgs84)


# ---------------------------------------------------------------------------
# 3. Calcul de densité (intersection pondérée avec la grille INSEE)
# ---------------------------------------------------------------------------

CARREAU_AREA_KM2 = 0.2 * 0.2  # carreaux INSEE Filosofi de 200m x 200m = 0.04 km²


def _read_intersecting_grid(zone_l93):
    import geopandas as gpd

    # Lit uniquement les carreaux dans l'emprise de la zone : rapide (millisecondes
    # à quelques secondes) car le GeoPackage indexé (voir prepare_insee_data) a un
    # index spatial R-tree exploité nativement par GDAL/pyogrio, contrairement au
    # shapefile brut d'origine (2,3M carreaux, lecture séquentielle sans index).
    gpkg_path = ensure_insee_indexed()
    minx, miny, maxx, maxy = zone_l93.bounds
    grid = gpd.read_file(gpkg_path, bbox=(minx, miny, maxx, maxy))
    if grid.crs is not None and grid.crs.to_epsg() != 2154:
        grid = grid.to_crs(epsg=2154)

    if grid.empty:
        raise InseeDataUnavailable(
            "Aucun carreau INSEE trouvé dans l'emprise de la zone (hors couverture "
            "métropole/Martinique/Réunion, ou zone trop petite / mal placée)."
        )
    return grid


def _weighted_density(zone_l93):
    """Densité MOYENNE pondérée par recouvrement (population totale estimée
    dans la zone / surface de la zone), à partir de la grille INSEE 200 m.
    À utiliser pour une densité moyenne sur une zone large (ex: rayon de 5 km
    en section D). Retourne (densite_ppl_km2, surface_km2, population_estimee,
    carreaux_utilises)."""
    grid = _read_intersecting_grid(zone_l93)

    zone_area_km2 = zone_l93.area / 1_000_000
    total_population = 0.0
    used = 0
    for _, row in grid.iterrows():
        inter = row.geometry.intersection(zone_l93)
        if inter.is_empty:
            continue
        overlap_ratio = inter.area / row.geometry.area
        total_population += float(row["ind"]) * overlap_ratio
        used += 1

    densite = total_population / zone_area_km2 if zone_area_km2 > 0 else 0.0
    return densite, zone_area_km2, total_population, used


def _max_cell_density(zone_l93):
    """Densité MAXIMALE parmi les carreaux INSEE 200 m qui touchent la zone
    (chaque carreau a sa propre densité = ind / 0.04 km², indépendamment de sa
    fraction de recouvrement avec la zone). À utiliser pour la "densité de
    population maximale dans l'empreinte du GRC" (section B) : c'est le
    carreau le plus dense qui détermine le pire cas de risque au sol, pas une
    moyenne sur toute l'empreinte. Retourne (densite_max_ppl_km2, surface_km2,
    population_carreau_le_plus_dense, carreaux_utilises)."""
    grid = _read_intersecting_grid(zone_l93)

    zone_area_km2 = zone_l93.area / 1_000_000
    max_density = 0.0
    max_pop = 0.0
    used = 0
    for _, row in grid.iterrows():
        if row.geometry.intersects(zone_l93) and not row.geometry.intersection(zone_l93).is_empty:
            used += 1
            pop = float(row["ind"])
            cell_density = pop / CARREAU_AREA_KM2
            if cell_density > max_density:
                max_density = cell_density
                max_pop = pop

    return max_density, zone_area_km2, max_pop, used


def _grc_categorie(densite):
    if densite < 5:
        return "< 5 ppl/km² (iGRC 2)"
    if densite < 50:
        return "< 50 ppl/km² (iGRC 3)"
    if densite < 500:
        return "< 500 ppl/km² (iGRC 4)"
    if densite < 5000:
        return "< 5000 ppl/km² (iGRC 5)"
    return None  # au-delà du barème du formulaire -> traitement manuel (zone contrôlée ?)


def compute_density_from_kml(kml_path):
    """Densité de population MAXIMALE dans l'empreinte du GRC (section B —
    Risques au sol) : le carreau INSEE 200m le plus dense qui touche la zone
    tampon/contingence du KML détermine le pire cas, conformément au libellé
    du formulaire ("Densité de population maximale..."). Retourne
    {"densite_ppl_km2", "categorie_suggeree", "surface_km2",
    "population_estimee", "carreaux_utilises"}."""
    zone_wgs84 = find_grc_zone(kml_path)
    zone_l93 = _to_lambert93(zone_wgs84)
    densite, surface_km2, pop, used = _max_cell_density(zone_l93)
    return {
        "densite_ppl_km2": round(densite, 2),
        "categorie_suggeree": _grc_categorie(densite),
        "surface_km2": round(surface_km2, 4),
        "population_estimee": round(pop, 1),
        "carreaux_utilises": used,
    }


# ---------------------------------------------------------------------------
# 4. Densité moyenne dans un rayon de 5 km autour du centre du buffer
#    (section D — Exigence de confinement, par parcelle)
# ---------------------------------------------------------------------------

RAYON_CONFINEMENT_M = 5000  # 5 km (~3.11 miles)


def compute_avg_density_5km(kml_path, parcelle_ref=None):
    """Densité moyenne pondérée sur un cercle de 5 km de rayon centré sur le
    centroïde du buffer de la parcelle (ou de toute la zone si `parcelle_ref`
    n'est pas fourni ou ne correspond à aucun placemark nommé). Retourne le
    même format que compute_density_from_kml, plus le centre utilisé."""
    zone_wgs84 = find_parcelle_zone(kml_path, parcelle_ref) if parcelle_ref else find_grc_zone(kml_path)
    zone_l93 = _to_lambert93(zone_wgs84)
    center = zone_l93.centroid
    circle_l93 = center.buffer(RAYON_CONFINEMENT_M, resolution=64)

    densite, surface_km2, pop, used = _weighted_density(circle_l93)
    return {
        "densite_ppl_km2": round(densite, 2),
        "surface_km2": round(surface_km2, 4),
        "population_estimee": round(pop, 1),
        "carreaux_utilises": used,
        "centre_l93": (round(center.x, 1), round(center.y, 1)),
    }