# -*- coding: utf-8 -*-
"""
Plugin QGIS : Sen Hydro - Delimitation de Bassins Versants
------------------------------------------------
Auteur       : Adiouma FALL - Geo Senegal 2026
Base sur     : l'API HTTP publique et gratuite "Global Watersheds"
               (https://mghydro.com/watersheds), qui applique la
               methode du package Python open-source "delineator"
               (Matthew Heberger, MIT License) sur les donnees
               MERIT-Hydro/HydroSHEDS.
               https://github.com/mheberger/delineator

Interface inspiree du site https://mghydro.com/watersheds : options de
calcul (precision, simplification, beautification), personnalisation du
style d'affichage (couleurs, transparence, epaisseur), et rapport du
bassin versant (surface, perimetre, longueur du reseau), le tout
directement integre dans QGIS.

Note technique :
    Depuis la version 1.10, le calcul n'est plus effectue localement
    (package Python "delineator" installe dans un environnement dedie) :
    le plugin appelle directement, en HTTP asynchrone (QgsNetworkAccessManager),
    l'API "Global Watersheds" (mghydro.com/app/watershed_api et
    upstream_rivers_api). Cela supprime toute dependance Python locale
    (numpy, geopandas, rasterio...) et le mecanisme d'installation
    (venv/vendor) qui l'accompagnait, au prix d'une connexion internet
    desormais necessaire a chaque delimitation.

Limitation connue :
    Le mode "Aval - Tracer le trajet d'ecoulement" (flowpath_api sur
    mghydro.com) n'est pas encore integre a Sen Hydro. Seule la
    delimitation amont du bassin versant est disponible ici.
"""

import json
import os
import tempfile
import time

from qgis.PyQt.QtCore import Qt, QUrl, QTimer
from qgis.PyQt.QtGui import QIcon, QColor
from qgis.PyQt.QtNetwork import QNetworkRequest
from qgis.PyQt.QtWidgets import (
    QAction,
    QMessageBox,
    QFileDialog,
    QApplication,
)
from qgis.core import (
    QgsProject,
    QgsVectorLayer,
    QgsCoordinateReferenceSystem,
    QgsCoordinateTransform,
    QgsSingleSymbolRenderer,
    QgsFillSymbol,
    QgsLineSymbol,
    QgsMarkerSymbol,
    QgsGeometry,
    QgsDistanceArea,
    QgsUnitTypes,
    QgsMessageLog,
    Qgis,
    QgsCategorizedSymbolRenderer,
    QgsRendererCategory,
    QgsWkbTypes,
    QgsFeature,
    QgsNetworkAccessManager,
)
from qgis.gui import QgsMapToolEmitPoint

from .dockwidget import SenHydroDockWidget


GLOBAL_WATERSHEDS_API_BASE = "https://mghydro.com/app"
API_REQUEST_TIMEOUT_MS = 60000


LAYER_LABELS = {
    "watershed": "Bassin versant",
    "rivers": "Reseau hydrographique",
    "outlet": "Exutoire",
}


class WatershedDelineatorPlugin:
    """
    Plugin principal.
    Cree par Adiouma FALL - Geo Senegal 2026.
    """

    COMMUNES_LAYER_NAME = "Limites_communes"
    AQUIFER_LAYER_NAME = "unites_aquiferes_corrige"
    SOIL_LAYER_NAME = "Morpho_Pedo"

    def __init__(self, iface):
        self.iface = iface
        self.canvas = iface.mapCanvas()
        self.action = None
        self.map_tool = None
        self.dock = None
        self._senegal_boundary_4326 = None
        self._senegal_boundary_simplified_4326 = None
        self._senegal_engine = None
        self._aquifer_layer_cache = None
        self._soil_layer_cache = None
        self._pending_replies = []

    def initGui(self):
        icon_path = os.path.join(os.path.dirname(__file__), "icon.png")
        self.action = QAction(
            QIcon(icon_path),
            "Sen Hydro - Delimiter un bassin versant (par Adiouma FALL - Geo Senegal)",
            self.iface.mainWindow(),
        )
        self.action.setCheckable(True)
        self.action.triggered.connect(self.toggle_panel)
        self.iface.addToolBarIcon(self.action)
        self.iface.addPluginToMenu("&Sen Hydro", self.action)

        self.map_tool = QgsMapToolEmitPoint(self.canvas)
        self.map_tool.canvasClicked.connect(self.on_map_click)

        self.dock = SenHydroDockWidget(self.iface.mainWindow())
        self.iface.addDockWidget(Qt.RightDockWidgetArea, self.dock)
        self.dock.hide()
        self.dock.btn_activate.toggled.connect(self.toggle_tool)
        self.dock.reportExportRequested.connect(self.export_report)
        self.dock.visibilityChanged.connect(self._on_dock_visibility_changed)

    def unload(self):
        self.iface.removeToolBarIcon(self.action)
        self.iface.removePluginMenu("&Sen Hydro", self.action)
        if self.map_tool is not None:
            self.canvas.unsetMapTool(self.map_tool)
        if self.dock is not None:
            self.iface.removeDockWidget(self.dock)
        self._abort_pending_replies()

    def _abort_pending_replies(self):
        """Annule toute requete HTTP en cours vers l'API Global Watersheds."""
        for reply in self._pending_replies:
            try:
                reply.abort()
            except Exception as e:
                # L'annulation d'une requete deja terminee/detruite peut
                # legitimement echouer ; on journalise plutot que d'ignorer
                # silencieusement l'erreur (evite le pattern try/except/pass).
                QgsMessageLog.logMessage(
                    "Sen Hydro: echec de l'annulation d'une requete en attente: %s" % e,
                    "Sen Hydro", Qgis.Warning,
                )
        self._pending_replies = []

    def toggle_panel(self, checked):
        if self.dock is None:
            return
        self.dock.setVisible(checked)

    def _on_dock_visibility_changed(self, visible):
        if self.action is not None:
            self.action.setChecked(visible)
        if not visible and self.dock is not None:
            self.dock.btn_activate.setChecked(False)

    def toggle_tool(self, checked):
        if checked:
            self.canvas.setMapTool(self.map_tool)
        else:
            self.canvas.unsetMapTool(self.map_tool)

    def _build_api_url(self, endpoint, lat, lon, config):
        """
        Construit l'URL de l'API Global Watersheds (mghydro.com) pour
        l'endpoint donne ("watershed_api" ou "upstream_rivers_api").
        """
        params = {
            "lat": "%.6f" % lat,
            "lng": "%.6f" % lon,
            "precision": "high" if config.get("high_res", True) else "low",
            "simplify": "true" if config.get("simplify", False) else "false",
            "beautify": "true" if config.get("smooth", False) else "false",
        }
        query = "&".join("%s=%s" % (k, v) for k, v in params.items())
        return "%s/%s?%s" % (GLOBAL_WATERSHEDS_API_BASE, endpoint, query)

    def _start_api_request(self, endpoint, lat, lon, config, on_finished):
        """
        Lance une requete GET asynchrone vers l'API Global Watersheds, sans
        bloquer l'interface de QGIS (contrairement a un appel HTTP
        synchrone). `on_finished(reply)` est appele une fois la reponse
        recue (succes ou erreur). Une echeance de securite interrompt la
        requete si le serveur ne repond pas.
        """
        url = self._build_api_url(endpoint, lat, lon, config)
        request = QNetworkRequest(QUrl(url))
        reply = QgsNetworkAccessManager.instance().get(request)
        self._pending_replies.append(reply)

        timer = QTimer()
        timer.setSingleShot(True)

        def _cleanup_and_call():
            timer.stop()
            if reply in self._pending_replies:
                self._pending_replies.remove(reply)
            on_finished(reply)

        timer.timeout.connect(reply.abort)
        reply.finished.connect(_cleanup_and_call)
        timer.start(API_REQUEST_TIMEOUT_MS)
        return reply

    @staticmethod
    def _reply_http_status(reply):
        return reply.attribute(QNetworkRequest.HttpStatusCodeAttribute)

    def _http_error_message(self, reply):
        """Traduit une reponse HTTP en erreur en un message utilisateur clair."""
        status = self._reply_http_status(reply)
        if status is None:
            return (
                "Impossible de contacter le service de delimitation "
                "(mghydro.com). Verifiez votre connexion internet."
            )
        if status == 400:
            return "La requete envoyee au service de delimitation etait invalide."
        if status == 404:
            return None  # cas gere specifiquement par l'appelant (pas d'erreur bloquante)
        if status in (502, 503):
            return "Le service de delimitation est temporairement sature. Reessayez dans quelques instants."
        if status == 500:
            return (
                "Une erreur est survenue sur le serveur de delimitation. "
                "Reessayez plus tard, ou contactez l'auteur (matt@mghydro.com) "
                "si le probleme persiste."
            )
        return "Le service de delimitation a repondu avec le code d'erreur %s." % status

    def on_map_click(self, point, button):
        source_crs = self.canvas.mapSettings().destinationCrs()
        dest_crs = QgsCoordinateReferenceSystem("EPSG:4326")
        transform = QgsCoordinateTransform(source_crs, dest_crs, QgsProject.instance())
        point_wgs84 = transform.transform(point)
        lon, lat = point_wgs84.x(), point_wgs84.y()

        config = self.dock.get_config() if self.dock is not None else {}
        self._abort_pending_replies()

        tmp_dir = tempfile.mkdtemp(prefix="watershed_")
        start_time = time.monotonic()

        status_bar = self.iface.mainWindow().statusBar()
        status_bar.showMessage("Sen Hydro : calcul du bassin versant en cours...")
        QApplication.setOverrideCursor(Qt.WaitCursor)

        self._start_api_request(
            "watershed_api", lat, lon, config,
            lambda reply: self._on_watershed_reply(reply, lat, lon, config, tmp_dir, start_time),
        )

    def _on_watershed_reply(self, reply, lat, lon, config, tmp_dir, start_time):
        status_bar = self.iface.mainWindow().statusBar()

        if reply.error() or self._reply_http_status(reply) != 200:
            QApplication.restoreOverrideCursor()
            status_bar.clearMessage()
            status = self._reply_http_status(reply)
            if status == 404:
                QMessageBox.warning(
                    self.iface.mainWindow(),
                    "Aucun resultat",
                    "Aucun bassin versant n'a pu etre calcule pour ce point "
                    "(par exemple, point situe en mer ou hors du domaine de "
                    "donnees MERIT-Hydro). Essayez un autre emplacement, de "
                    "preference sur ou pres d'un cours d'eau.",
                )
            else:
                QMessageBox.critical(
                    self.iface.mainWindow(), "Erreur de delimitation", self._http_error_message(reply),
                )
            reply.deleteLater()
            return

        body = bytes(reply.readAll())
        reply.deleteLater()
        try:
            geojson = json.loads(body)
        except Exception as e:
            QApplication.restoreOverrideCursor()
            status_bar.clearMessage()
            QMessageBox.critical(
                self.iface.mainWindow(),
                "Erreur de lecture des resultats",
                "La reponse du service de delimitation n'a pas pu etre lue :\n%s" % e,
            )
            return

        features = geojson.get("features") or []
        if not features:
            QApplication.restoreOverrideCursor()
            status_bar.clearMessage()
            QMessageBox.warning(
                self.iface.mainWindow(),
                "Aucun resultat",
                "Aucune couche n'a ete generee pour ce point. Essayez un autre "
                "emplacement, de preference sur ou pres d'un cours d'eau.",
            )
            return

        props = features[0].get("properties") or {}
        try:
            area_km2 = float(props.get("area_km2"))
        except (TypeError, ValueError):
            area_km2 = None
        outlet_lat = props.get("outlet_lat", lat)
        outlet_lng = props.get("outlet_lng", lon)

        watershed_path = os.path.join(tmp_dir, "watershed.geojson")
        with open(watershed_path, "wb") as f:
            f.write(body)

        outlet_path = os.path.join(tmp_dir, "outlet.geojson")
        outlet_geojson = {
            "type": "FeatureCollection",
            "features": [
                {
                    "type": "Feature",
                    "properties": {"type": "requested"},
                    "geometry": {"type": "Point", "coordinates": [lon, lat]},
                },
                {
                    "type": "Feature",
                    "properties": {"type": "snapped"},
                    "geometry": {"type": "Point", "coordinates": [outlet_lng, outlet_lat]},
                },
            ],
        }
        with open(outlet_path, "w", encoding="utf-8") as f:
            json.dump(outlet_geojson, f)

        paths = {"watershed": watershed_path, "outlet": outlet_path}
        stats = {}
        if area_km2 is not None:
            stats["area_km2"] = area_km2

        show_rivers = self.dock.chk_show_rivers.isChecked() if self.dock is not None else True
        if show_rivers:
            status_bar.showMessage("Sen Hydro : recuperation du reseau hydrographique...")
            self._start_api_request(
                "upstream_rivers_api", lat, lon, config,
                lambda reply: self._on_rivers_reply(reply, lat, lon, config, tmp_dir, start_time, paths, stats),
            )
        else:
            self._log_duration(start_time, lat, lon)
            QApplication.restoreOverrideCursor()
            status_bar.clearMessage()
            self._finalize_results(paths, stats, lat, lon)

    def _on_rivers_reply(self, reply, lat, lon, config, tmp_dir, start_time, paths, stats):
        status_bar = self.iface.mainWindow().statusBar()
        status = self._reply_http_status(reply)

        if not reply.error() and status == 200:
            try:
                body = bytes(reply.readAll())
                rivers_geojson = json.loads(body)
                if rivers_geojson.get("features"):
                    rivers_path = os.path.join(tmp_dir, "rivers.geojson")
                    with open(rivers_path, "wb") as f:
                        f.write(body)
                    paths["rivers"] = rivers_path
            except Exception as e:
                QgsMessageLog.logMessage(
                    "Sen Hydro: reponse upstream_rivers_api illisible: %s" % e,
                    "Sen Hydro", Qgis.Warning,
                )
        elif status != 404:
            # 404 = pas de reseau hydrographique resolu pour ce bassin (cas
            # frequent pour un tres petit bassin) : ce n'est pas une erreur
            # bloquante, le bassin versant reste affiche sans les rivieres.
            QgsMessageLog.logMessage(
                "Sen Hydro: echec de upstream_rivers_api (%s): %s"
                % (status, reply.errorString()),
                "Sen Hydro", Qgis.Warning,
            )

        reply.deleteLater()
        self._log_duration(start_time, lat, lon)
        QApplication.restoreOverrideCursor()
        status_bar.clearMessage()
        self._finalize_results(paths, stats, lat, lon)

    def _log_duration(self, start_time, lat, lon):
        QgsMessageLog.logMessage(
            "Sen Hydro: calcul termine en %.1f s pour (%.5f, %.5f)"
            % (time.monotonic() - start_time, lat, lon),
            "Sen Hydro", Qgis.Info,
        )

    def _finalize_results(self, paths, stats, lat, lon):
        style = self.dock.get_style() if self.dock is not None else None
        show_watershed = self.dock.chk_show_watershed.isChecked() if self.dock is not None else True
        outlet_both = self.dock.chk_outlet_both.isChecked() if self.dock is not None else True
        clip_senegal = self.dock.chk_clip_senegal.isChecked() if self.dock is not None else False

        loaded_layers = {}
        for key, path in paths.items():
            if key == "watershed" and not show_watershed:
                continue
            layer = self._add_layer(path, LAYER_LABELS.get(key, key))
            if layer is not None:
                if key == "outlet" and not outlet_both:
                    self._keep_only_snapped_outlet(layer)
                loaded_layers[key] = layer
                if style is not None:
                    self._apply_style(layer, key, style)

        clip_warning = None
        if clip_senegal:
            senegal_geom = self._get_senegal_boundary()
            if senegal_geom is None:
                clip_warning = (
                    "La couche '%s' est introuvable (ou vide) dans le projet : "
                    "le decoupage aux limites du Senegal a ete ignore." % self.COMMUNES_LAYER_NAME
                )
            else:
                proj = QgsProject.instance()
                for key in list(loaded_layers.keys()):
                    old_layer = loaded_layers[key]
                    clipped_layer = self._clip_layer_to_senegal(old_layer)
                    if clipped_layer is not old_layer:
                        proj.removeMapLayer(old_layer.id())
                        proj.addMapLayer(clipped_layer)
                        if style is not None:
                            self._apply_style(clipped_layer, key, style)
                        loaded_layers[key] = clipped_layer

        # L'API ne fournit ni perimetre, ni longueur/ordre du reseau
        # hydrographique : on les calcule cote client sur les couches
        # effectivement affichees (apres decoupage Senegal eventuel).
        if "watershed" in loaded_layers:
            stats = self._recompute_watershed_stats(loaded_layers["watershed"], stats)
        if "rivers" in loaded_layers:
            stats = self._recompute_river_stats(loaded_layers["rivers"], stats)

        cross_aquifers = (
            self.dock.chk_cross_aquifers.isChecked() if self.dock is not None else False
        )
        aquifer_warning = None
        soil_warning = None
        if cross_aquifers and "watershed" in loaded_layers:
            aquifer_layer = self._get_reference_layer(
                "unites_aquiferes.gpkg", self.AQUIFER_LAYER_NAME, "_aquifer_layer_cache"
            )
            if aquifer_layer is None:
                stats = dict(stats)
                stats["aquifer_layer_missing"] = True
                aquifer_warning = (
                    "La couche des unites aquiferes ('%s') est introuvable "
                    "(ni embarquee dans le plugin, ni dans le projet) : le "
                    "croisement avec les nappes phreatiques a ete ignore."
                    % self.AQUIFER_LAYER_NAME
                )
            else:
                stats = self._cross_with_reference_layer(
                    loaded_layers["watershed"], aquifer_layer, "nom",
                    "Unites aquiferes du bassin versant", stats, "aquifer",
                    (180, 90, 30),
                )

            soil_layer = self._get_reference_layer(
                "morpho_pedo.gpkg", self.SOIL_LAYER_NAME, "_soil_layer_cache"
            )
            if soil_layer is None:
                stats = dict(stats)
                stats["soil_layer_missing"] = True
                soil_warning = (
                    "La couche des types de sol ('%s') est introuvable (ni "
                    "embarquee dans le plugin, ni dans le projet) : le "
                    "croisement avec la pedologie a ete ignore."
                    % self.SOIL_LAYER_NAME
                )
            else:
                stats = self._cross_with_reference_layer(
                    loaded_layers["watershed"], soil_layer, "MSDNOM",
                    "Types de sol du bassin versant", stats, "soil",
                    (100, 140, 40),
                )

        if self.dock is not None:
            self.dock.show_report(stats, lat, lon)

        if clip_warning:
            QMessageBox.information(self.iface.mainWindow(), "Decoupage Senegal", clip_warning)
        if aquifer_warning:
            QMessageBox.information(self.iface.mainWindow(), "Nappes phreatiques", aquifer_warning)
        if soil_warning:
            QMessageBox.information(self.iface.mainWindow(), "Types de sol", soil_warning)

        self.canvas.refresh()

    def _keep_only_snapped_outlet(self, layer):
        """
        Ne conserve que le point exutoire 'snapped' (ajuste a la riviere),
        qui est celui reellement utilise pour le calcul du bassin versant.
        """
        type_idx = layer.fields().indexOf("type")
        if type_idx < 0:
            return
        layer.startEditing()
        to_delete = [
            f.id() for f in layer.getFeatures()
            if str(f.attribute("type")).lower() != "snapped"
        ]
        for fid in to_delete:
            layer.deleteFeature(fid)
        layer.commitChanges()
        layer.updateExtents()

    def _add_layer(self, path, layer_name):
        layer = QgsVectorLayer(path, layer_name, "ogr")
        if layer.isValid():
            QgsProject.instance().addMapLayer(layer)
            return layer
        QMessageBox.warning(
            self.iface.mainWindow(),
            "Couche invalide",
            "Impossible de charger la couche : %s" % layer_name,
        )
        return None

    def _apply_style(self, layer, key, style):
        opacity = style.get("opacity", 0.35)
        width = style.get("line_width", 0.6)

        try:
            if key == "watershed":
                color = QColor(style["watershed_color"])
                symbol = QgsFillSymbol.createSimple({})
                symbol.setColor(color)
                symbol.setOpacity(opacity)
                if symbol.symbolLayerCount() > 0:
                    symbol.symbolLayer(0).setStrokeColor(color)
                    symbol.symbolLayer(0).setStrokeWidth(width)
                layer.setRenderer(QgsSingleSymbolRenderer(symbol))
            elif key == "rivers":
                color = QColor(style["river_color"])
                sorder_idx = layer.fields().indexOf("sorder")
                if sorder_idx >= 0:
                    orders = sorted({
                        f.attribute("sorder") for f in layer.getFeatures()
                        if f.attribute("sorder") is not None
                    })
                    categories = []
                    max_order = max(orders) if orders else 1
                    for order in orders:
                        # Plus l'ordre de Strahler est eleve, plus le trait est
                        # epais : les cours d'eau principaux ressortent mieux.
                        order_width = width * (0.6 + 1.4 * (order / max(max_order, 1)))
                        symbol = QgsLineSymbol.createSimple({})
                        symbol.setColor(color)
                        symbol.setWidth(order_width)
                        categories.append(
                            QgsRendererCategory(
                                order, symbol, "Ordre %s" % order
                            )
                        )
                    layer.setRenderer(
                        QgsCategorizedSymbolRenderer("sorder", categories)
                    )
                else:
                    symbol = QgsLineSymbol.createSimple({})
                    symbol.setColor(color)
                    symbol.setWidth(width)
                    layer.setRenderer(QgsSingleSymbolRenderer(symbol))
            elif key == "outlet":
                color = QColor(style["watershed_color"])
                type_idx = layer.fields().indexOf("type")
                if type_idx >= 0:
                    # Point "demande" (clic) : marqueur creux (contour seul)
                    requested_symbol = QgsMarkerSymbol.createSimple({})
                    requested_symbol.setColor(QColor(255, 255, 255, 0))
                    if requested_symbol.symbolLayerCount() > 0:
                        requested_symbol.symbolLayer(0).setStrokeColor(color)
                        requested_symbol.symbolLayer(0).setStrokeWidth(0.8)
                    requested_symbol.setSize(3)

                    # Point "ajuste a la riviere" : marqueur plein
                    snapped_symbol = QgsMarkerSymbol.createSimple({})
                    snapped_symbol.setColor(color)
                    snapped_symbol.setSize(3)

                    categories = [
                        QgsRendererCategory(
                            "requested", requested_symbol, "Point demande (clic)"
                        ),
                        QgsRendererCategory(
                            "snapped", snapped_symbol, "Point ajuste a la riviere"
                        ),
                    ]
                    layer.setRenderer(
                        QgsCategorizedSymbolRenderer("type", categories)
                    )
                else:
                    symbol = QgsMarkerSymbol.createSimple({})
                    symbol.setColor(color)
                    symbol.setSize(3)
                    layer.setRenderer(QgsSingleSymbolRenderer(symbol))
        except Exception as style_err:
            # Le style est une amelioration cosmetique : en cas d'echec,
            # la couche reste affichee avec le style par defaut de QGIS.
            # On journalise plutot que d'ignorer silencieusement l'erreur.
            QgsMessageLog.logMessage(
                "Sen Hydro: echec de l'application du style pour '%s': %s"
                % (key, style_err),
                "Sen Hydro",
                Qgis.Warning,
            )

        layer.triggerRepaint()

    def _get_senegal_boundary(self):
        """
        Retourne (et met en cache) la geometrie de la frontiere du Senegal,
        en EPSG:4326 (meme CRS que les couches produites par 'delineator').

        Deux sources sont essayees, dans cet ordre :
        1. Le fichier embarque dans le plugin (data/senegal_boundary.gpkg),
           genere une fois pour toutes a partir des limites communales
           officielles. Fonctionne sans configuration, dans n'importe quel
           projet QGIS.
        2. A defaut (fichier absent/corrompu), une couche nommee
           'Limites_communes' presente dans le projet ouvert, dissoute a la
           volee (comportement historique, conserve en secours).
        """
        if self._senegal_boundary_4326 is not None:
            return self._senegal_boundary_4326

        # --- Source 1 : fichier embarque dans le plugin ---
        bundled_path = os.path.join(
            os.path.dirname(__file__), "data", "senegal_boundary.gpkg"
        )
        if os.path.isfile(bundled_path):
            bundled_layer = QgsVectorLayer(bundled_path, "senegal_boundary", "ogr")
            if bundled_layer.isValid() and bundled_layer.featureCount() > 0:
                feat = next(bundled_layer.getFeatures(), None)
                if feat is not None and feat.geometry() is not None and not feat.geometry().isEmpty():
                    self._senegal_boundary_4326 = QgsGeometry(feat.geometry())
                    return self._senegal_boundary_4326

        # --- Source 2 (secours) : couche 'Limites_communes' du projet ---
        communes_layer = None
        for layer in QgsProject.instance().mapLayers().values():
            if layer.name() == self.COMMUNES_LAYER_NAME:
                communes_layer = layer
                break

        if communes_layer is None:
            return None

        geoms = [
            f.geometry() for f in communes_layer.getFeatures()
            if f.geometry() is not None and not f.geometry().isEmpty()
        ]
        if not geoms:
            return None

        union_geom = QgsGeometry.unaryUnion(geoms)
        if union_geom is None or union_geom.isEmpty():
            return None

        transform = QgsCoordinateTransform(
            communes_layer.crs(),
            QgsCoordinateReferenceSystem("EPSG:4326"),
            QgsProject.instance(),
        )
        union_geom.transform(transform)

        self._senegal_boundary_4326 = union_geom
        return union_geom

    def _get_senegal_clip_engine(self):
        """
        Retourne (et met en cache) un couple (geometrie Senegal simplifiee,
        moteur GEOS "prepared") utilise pour des tests d'intersection tres
        rapides, meme avec des milliers d'entites. La frontiere brute du
        Senegal comporte ~55 000 sommets ; une intersection geometrique
        classique contre une geometrie aussi detaillee, repetee pour chaque
        entite (des centaines de troncons de riviere), est lente. La
        simplification (tolerance ~50 m, negligeable a l'echelle d'un
        bassin versant) et la preparation GEOS accelerent cela d'un facteur
        5 a 10.
        """
        if self._senegal_engine is not None:
            return self._senegal_boundary_simplified_4326, self._senegal_engine

        senegal_geom = self._get_senegal_boundary()
        if senegal_geom is None:
            return None, None

        simplified = senegal_geom.simplify(0.0005)  # ~50 m a cette latitude
        if simplified is None or simplified.isEmpty():
            simplified = senegal_geom

        engine = QgsGeometry.createGeometryEngine(simplified.constGet())
        engine.prepareGeometry()

        self._senegal_boundary_simplified_4326 = simplified
        self._senegal_engine = engine
        return simplified, engine

    def _clip_layer_to_senegal(self, layer):
        """
        Retourne une NOUVELLE couche memoire contenant uniquement la partie
        des entites de `layer` situee sur le territoire senegalais.

        On reconstruit une couche memoire plutot que d'editer `layer` en
        place (startEditing/changeGeometry/deleteFeature) : sur des couches
        de plusieurs centaines/milliers d'entites (reseau hydrographique),
        le systeme d'edition/undo de QGIS est beaucoup plus lent que de
        simplement construire une nouvelle couche avec les geometries deja
        decoupees.
        """
        simplified_senegal, engine = self._get_senegal_clip_engine()
        if engine is None:
            return layer

        geom_type_str = QgsWkbTypes.displayString(layer.wkbType())
        mem_layer = QgsVectorLayer(
            "%s?crs=%s" % (geom_type_str, layer.crs().authid()),
            layer.name(),
            "memory",
        )
        mem_provider = mem_layer.dataProvider()
        mem_provider.addAttributes(layer.fields())
        mem_layer.updateFields()

        new_features = []
        for feat in layer.getFeatures():
            geom = feat.geometry()
            if geom is None or geom.isEmpty():
                continue
            if not engine.intersects(geom.constGet()):
                continue
            clipped = simplified_senegal.intersection(geom)
            if clipped is None or clipped.isEmpty():
                continue
            new_feat = QgsFeature(mem_layer.fields())
            new_feat.setGeometry(clipped)
            new_feat.setAttributes(feat.attributes())
            new_features.append(new_feat)

        mem_provider.addFeatures(new_features)
        mem_layer.updateExtents()
        return mem_layer

    def _make_distance_area(self):
        da = QgsDistanceArea()
        da.setEllipsoid("WGS84")
        da.setSourceCrs(
            QgsCoordinateReferenceSystem("EPSG:4326"), QgsProject.instance().transformContext()
        )
        return da

    def _recompute_watershed_stats(self, layer, stats):
        da = self._make_distance_area()
        total_area_km2 = 0.0
        total_perimeter_km = 0.0
        for feat in layer.getFeatures():
            geom = feat.geometry()
            if geom is None or geom.isEmpty():
                continue
            total_area_km2 += da.measureArea(geom) / 1_000_000.0
            total_perimeter_km += da.measurePerimeter(geom) / 1000.0
        stats = dict(stats)
        if layer.featureCount() > 0:
            stats["area_km2"] = total_area_km2
            stats["perimeter_km"] = total_perimeter_km
        else:
            stats.pop("area_km2", None)
            stats.pop("perimeter_km", None)
            stats["clipped_outside_senegal"] = True
        return stats

    def _recompute_river_stats(self, layer, stats):
        da = self._make_distance_area()
        total_length_km = 0.0
        sorder_idx = layer.fields().indexOf("sorder")
        max_sorder = None
        for feat in layer.getFeatures():
            geom = feat.geometry()
            if geom is not None and not geom.isEmpty():
                total_length_km += da.measureLength(geom) / 1000.0
            if sorder_idx >= 0:
                value = feat.attribute("sorder")
                if value is not None:
                    max_sorder = value if max_sorder is None else max(max_sorder, value)
        stats = dict(stats)
        stats["river_reach_count"] = layer.featureCount()
        if layer.featureCount() > 0:
            stats["river_length_km"] = total_length_km
        else:
            stats.pop("river_length_km", None)
        if max_sorder is not None:
            stats["max_stream_order"] = int(max_sorder)
        else:
            stats.pop("max_stream_order", None)
        return stats

    def _get_reference_layer(self, bundled_filename, project_layer_name, cache_attr_name):
        """
        Retourne (et met en cache) une couche de reference vectorielle :
        1. En priorite, le fichier embarque dans le plugin
           (data/<bundled_filename>) : fonctionne sans configuration, dans
           n'importe quel projet QGIS.
        2. A defaut, une couche du meme nom presente dans le projet ouvert.
        Retourne None si aucune des deux sources n'est disponible.
        """
        cached = getattr(self, cache_attr_name)
        if cached is not None:
            return cached

        layer = None
        bundled_path = os.path.join(
            os.path.dirname(__file__), "data", bundled_filename
        )
        if os.path.isfile(bundled_path):
            candidate = QgsVectorLayer(bundled_path, project_layer_name, "ogr")
            if candidate.isValid() and candidate.featureCount() > 0:
                layer = candidate

        if layer is None:
            for lyr in QgsProject.instance().mapLayers().values():
                if lyr.name() == project_layer_name:
                    layer = lyr
                    break

        if layer is not None:
            setattr(self, cache_attr_name, layer)
        return layer

    def _cross_with_reference_layer(
        self, watershed_layer, ref_layer, name_field, result_layer_name,
        stats, stats_prefix, fill_color,
    ):
        """
        Croise le bassin versant delimite avec une couche de reference
        (unites aquiferes, types de sol, etc.), decoupee par le contour du
        bassin versant. Ajoute une couche memoire d'intersection au projet
        et complete les statistiques du rapport avec, pour le prefixe
        `stats_prefix` ("aquifer" ou "soil") :
          - <prefix>_units       : liste triee des libelles rencontres
          - <prefix>_unit_areas_km2 : dict libelle -> surface (km2)
          - <prefix>_area_km2    : surface totale (km2)
        """
        watershed_geoms = [
            f.geometry() for f in watershed_layer.getFeatures()
            if f.geometry() is not None and not f.geometry().isEmpty()
        ]
        if not watershed_geoms:
            return stats

        watershed_union = QgsGeometry.unaryUnion(watershed_geoms)
        if watershed_union is None or watershed_union.isEmpty():
            return stats

        transform = QgsCoordinateTransform(
            watershed_layer.crs(), ref_layer.crs(), QgsProject.instance()
        )
        watershed_union_proj = QgsGeometry(watershed_union)
        watershed_union_proj.transform(transform)

        geom_type_str = QgsWkbTypes.displayString(ref_layer.wkbType())
        mem_layer = QgsVectorLayer(
            "%s?crs=%s" % (geom_type_str, ref_layer.crs().authid()),
            result_layer_name,
            "memory",
        )
        mem_provider = mem_layer.dataProvider()
        mem_provider.addAttributes(ref_layer.fields())
        mem_layer.updateFields()

        da = QgsDistanceArea()
        da.setEllipsoid("WGS84")
        da.setSourceCrs(ref_layer.crs(), QgsProject.instance().transformContext())

        name_idx = ref_layer.fields().indexOf(name_field)

        new_features = []
        unit_areas = {}
        for feat in ref_layer.getFeatures():
            geom = feat.geometry()
            if geom is None or geom.isEmpty():
                continue
            if not geom.intersects(watershed_union_proj):
                continue
            clipped = geom.intersection(watershed_union_proj)
            if clipped is None or clipped.isEmpty():
                continue
            area_km2 = da.measureArea(clipped) / 1_000_000.0
            new_feat = QgsFeature(mem_layer.fields())
            new_feat.setGeometry(clipped)
            new_feat.setAttributes(feat.attributes())
            new_features.append(new_feat)
            label = feat.attribute(name_field) if name_idx >= 0 else None
            if not label:
                label = "unite %s" % feat.id()
            unit_areas[label] = unit_areas.get(label, 0.0) + area_km2

        mem_provider.addFeatures(new_features)
        mem_layer.updateExtents()

        if new_features:
            proj = QgsProject.instance()
            proj.addMapLayer(mem_layer)
            symbol = QgsFillSymbol.createSimple({})
            symbol.setColor(QColor(*fill_color))
            symbol.setOpacity(0.45)
            mem_layer.setRenderer(QgsSingleSymbolRenderer(symbol))
            mem_layer.triggerRepaint()

        stats = dict(stats)
        if unit_areas:
            stats["%s_units" % stats_prefix] = sorted(unit_areas.keys())
            stats["%s_unit_areas_km2" % stats_prefix] = unit_areas
            stats["%s_area_km2" % stats_prefix] = sum(unit_areas.values())
        return stats

    def export_report(self):
        if self.dock is None:
            return
        text = self.dock.get_report_text()
        if not text:
            return
        path, _ = QFileDialog.getSaveFileName(
            self.iface.mainWindow(),
            "Exporter le rapport",
            "rapport_bassin_versant.txt",
            "Fichiers texte (*.txt)",
        )
        if not path:
            return
        with open(path, "w", encoding="utf-8") as f:
            f.write(text)
        QMessageBox.information(
            self.iface.mainWindow(),
            "Rapport exporte",
            "Le rapport a ete enregistre :\n%s" % path,
        )
