"""
Raster Create-Study-Area overlay primitive.

Direct analogue of the ArcGIS Pro Spatial Analyst "Combine" tool that
the EMDS-for-ArcGIS-Pro raster Create-Study-Area used. QGIS ships the
same algorithm via GRASS as **`grass:r.cross`** (or `grass7:r.cross`
on older builds): it produces a unique-combination output raster + a
category table holding each input's value per output integer code.

This module wraps that call and enriches the bare category table with
the user-selected RAT columns from each input (joined on each input's
own pixel value), so the final output RAT has one row per unique cell
combination with `(per input × carried bands) + (per input × joined
RAT columns)`. The result is then written via RatWriter so engines
downstream can consume the RAT as a pandas DataFrame.

The primitive is callable end-to-end from the raster Create-Study-Area
dialog, and is designed to run inside a QgsTask-backed background
worker so large rasters don't freeze the UI.

Testability: `processing.run` and `RatReader` / `RatWriter` are
injectable via the `RasterOverlay` constructor. Unit tests stub them
to verify the algorithm id and parameter shape passed to r.cross,
plus the RAT enrichment join, without QGIS.
"""

from __future__ import annotations

from dataclasses import dataclass
from pathlib import Path
from typing import Any, Callable, Dict, List, Optional, Tuple

import pandas as pd

from .models import LayerKind, LayerRef, RasterOverlaySpec, Study
from .rat_io import RAT_VALUE_COLUMN, RatReader, RatWriter


# Algorithm id used by QGIS Processing for GRASS r.cross. Newer QGIS
# builds expose `grass:r.cross`; the older `grass7:` namespace is
# kept as a fallback.
GRASS_R_CROSS_ID_PRIMARY = "grass:r.cross"
GRASS_R_CROSS_ID_FALLBACK = "grass7:r.cross"


class RasterOverlayError(RuntimeError):
    """Raised when the raster overlay step cannot complete.

    Wraps GRASS / Processing failures, missing inputs, or RAT-join
    inconsistencies. Caller (the raster Create-Study-Area dialog or a
    QgsTask wrapper) renders the message to the UI.
    """


@dataclass
class RasterOverlayResult:
    """What the overlay primitive produced.

    output_raster_path: filesystem path to the unique-combination raster
    rat:                final enriched RAT as a pandas DataFrame
    cell_count_total:   number of pixels considered (None if the
                        Processing output didn't report it)
    unique_combinations: number of rows in the RAT
    """

    output_raster_path: str
    rat: pd.DataFrame
    cell_count_total: Optional[int]
    unique_combinations: int


# Type aliases for the injection seams. Kept narrow so tests don't
# need to import qgis just to construct stand-ins.
ProcessingRun = Callable[[str, Dict[str, Any]], Dict[str, Any]]


def _default_processing_run(algorithm_id: str, parameters: Dict[str, Any]) -> Dict[str, Any]:
    """Production default: dispatch through QGIS Processing."""
    from qgis import processing  # type: ignore
    return processing.run(algorithm_id, parameters)


class RasterOverlay:
    """Drive the raster Create-Study-Area overlay end-to-end.

    Construct with optional injected dependencies for testability:

        overlay = RasterOverlay()
        result  = overlay.run(study)

    Or under test:

        overlay = RasterOverlay(
            processing_run=fake_processing_run,
            rat_reader=fake_rat_reader,
            rat_writer=fake_rat_writer,
        )
    """

    def __init__(
        self,
        *,
        processing_run: ProcessingRun = _default_processing_run,
        rat_reader: Optional[RatReader] = None,
        rat_writer: Optional[RatWriter] = None,
        grass_algorithm_ids: Tuple[str, ...] = (
            GRASS_R_CROSS_ID_PRIMARY,
            GRASS_R_CROSS_ID_FALLBACK,
        ),
        layer_resolver: Optional[Callable[[str], Any]] = None,
    ):
        self._processing_run = processing_run
        self._rat_reader = rat_reader or RatReader()
        self._rat_writer = rat_writer or RatWriter()
        self._grass_algorithm_ids = grass_algorithm_ids
        self._layer_resolver = layer_resolver or _default_layer_resolver

    # ---- Public entry point --------------------------------------------

    def run(self, study: Study) -> RasterOverlayResult:
        """Execute the overlay for a fully-populated raster Study.

        Pre-conditions:
          - study.kind == LayerKind.RASTER
          - study.input_layer_refs has at least 1 raster LayerRef
            (single-input studies use a passthrough path — Phase 3.5e)
          - study.raster_overlay is populated (snap_to, output_path…)
        """
        self._validate_study(study)
        spec = study.raster_overlay  # validated non-None above
        ordered_refs = self._ordered_input_refs(study)

        # Resolve the live QgsRasterLayer for each named input so we
        # can grab their RATs and pass them to r.cross. The dialog
        # populated `study.input_layer_refs[*].name`; those names must
        # match layers loaded in the current QGIS project.
        resolved_layers = [
            (ref, self._layer_resolver(ref.name)) for ref in ordered_refs
        ]
        for ref, layer in resolved_layers:
            if layer is None:
                raise RasterOverlayError(
                    f"Input raster '{ref.name}' is not loaded in the "
                    f"current QGIS project. Add it back or rebuild the "
                    f"study area."
                )

        # Phase 3.5e: degenerate single-raster path skips r.cross
        # entirely. ArcGIS Pro's Combine also handles single-input
        # silently; we emit the input as the output (copying sidecars)
        # and pass its RAT through filtered to the user's carry-forward
        # selection.
        if len(resolved_layers) == 1:
            return self._run_single_raster_passthrough(resolved_layers[0], spec)

        # Phase 3.5d: pre-align misaligned inputs to the snap-to grid
        # via gdal:warpreproject before handing them to GRASS r.cross.
        # r.cross is stricter than ArcGIS Pro's Combine and will fail
        # or mis-align silently when CRS / extent / cell size differ.
        resolved_layers = self._prealign_inputs(resolved_layers, spec)

        # 1. Run r.cross to produce the unique-combination output raster.
        cross_result = self._run_grass_r_cross(
            resolved_layers,
            spec,
        )
        output_raster_path = self._extract_output_raster_path(cross_result, spec)

        # 2. r.cross also emits a category table mapping each output
        #    code -> tuple of input values. Pull that into a DataFrame.
        #    The shape is: one row per unique code, columns are
        #    Value (the output code) plus one per input layer holding
        #    that input's pixel value.
        category_df = self._read_r_cross_category_table(
            cross_result,
            output_raster_path,
            input_names=[ref.name for ref in ordered_refs],
        )

        # 3. Enrich: for each input, left-join its selected RAT columns
        #    onto the category table using the input's value column as
        #    the join key. This is the load-bearing step that gives
        #    EMDS engines the RAT columns they expect to see.
        enriched_rat = self._enrich_with_input_rats(
            category_df=category_df,
            resolved_layers=resolved_layers,
        )

        # 4. Write the enriched RAT alongside the output raster (native
        #    QgsRasterAttributeTable if available; sidecar .vat.dbf
        #    otherwise — see RatWriter for the fallback policy).
        self._rat_writer.write(output_raster_path, band=1, df=enriched_rat)

        return RasterOverlayResult(
            output_raster_path=output_raster_path,
            rat=enriched_rat,
            cell_count_total=cross_result.get("CELL_COUNT_TOTAL")
                if isinstance(cross_result, dict) else None,
            unique_combinations=len(enriched_rat.index),
        )

    # ---- Phase 3.5e: single-raster degenerate path -------------------

    def _run_single_raster_passthrough(
        self,
        ref_layer_tuple: "Tuple[LayerRef, Any]",
        spec: RasterOverlaySpec,
    ) -> RasterOverlayResult:
        """Single-input raster study: copy input to output, build canonical RAT.

        No unique-combination step is needed because there's nothing to
        combine. We still emit the output in the canonical multi-input
        RAT shape (Value + <input_name> + <input_name>__<col>...) so
        downstream consumers (RasterStudyInputLoader, BandMapper,
        engines) don't care which path produced it.
        """
        from .raster_output_writer import RasterOutputWriter
        ref, layer = ref_layer_tuple

        # 1. Copy the input raster to the output path (sidecars + all).
        source_path = self._get_layer_source_path(layer)
        if not source_path:
            raise RasterOverlayError(
                f"Single-raster input '{ref.name}' has no filesystem "
                "source path; cannot copy to output. Save the raster "
                "to disk first."
            )
        RasterOutputWriter._copy_raster_with_sidecars(
            source_path, spec.output_path
        )

        # 2. Build the output RAT from the input's RAT.
        input_rat = self._rat_reader.read(layer, band=1)
        if input_rat is None or input_rat.empty:
            if ref.rat_columns_to_carry:
                raise RasterOverlayError(
                    f"Input raster '{ref.name}' has no readable RAT "
                    f"but the study asked to carry RAT columns "
                    f"{ref.rat_columns_to_carry}. Either provide a "
                    f"RAT for the input or unselect those columns."
                )
            # No RAT available; emit an empty canonical shell so the
            # writer / readers see the expected columns even if no
            # rows exist yet.
            output_rat = pd.DataFrame({
                RAT_VALUE_COLUMN: pd.Series(dtype="Int64"),
                ref.name: pd.Series(dtype="Int64"),
            })
        else:
            cols_lower = {col.lower(): col for col in input_rat.columns}
            value_col = cols_lower.get(RAT_VALUE_COLUMN.lower())
            if value_col is None:
                raise RasterOverlayError(
                    f"Input raster '{ref.name}' RAT is missing the "
                    f"'{RAT_VALUE_COLUMN}' join key column."
                )
            missing = [
                c for c in ref.rat_columns_to_carry
                if c not in input_rat.columns
            ]
            if missing:
                raise RasterOverlayError(
                    f"Input raster '{ref.name}' RAT is missing "
                    f"requested columns: {missing}. Available: "
                    f"{list(input_rat.columns)}."
                )
            output_rat = pd.DataFrame()
            output_rat[RAT_VALUE_COLUMN] = input_rat[value_col]
            # Per-input value column = Value (output code == input
            # value in the degenerate case). Matches the multi-input
            # shape convention so BandMapper.resolve_source_column
            # works unchanged.
            output_rat[ref.name] = input_rat[value_col]
            for col in ref.rat_columns_to_carry:
                output_rat[f"{ref.name}__{col}"] = input_rat[col]

        # 3. Write the RAT alongside the copied raster.
        self._rat_writer.write(spec.output_path, band=1, df=output_rat)

        return RasterOverlayResult(
            output_raster_path=spec.output_path,
            rat=output_rat,
            cell_count_total=None,
            unique_combinations=len(output_rat.index),
        )

    @staticmethod
    def _get_layer_source_path(layer: Any) -> str:
        """Pull the on-disk source path out of a layer (best-effort)."""
        source_attr = getattr(layer, "source", None)
        if callable(source_attr):
            try:
                return str(source_attr())
            except Exception:
                return ""
        return ""

    # ---- Phase 3.5d: snap-to + resampling pre-alignment --------------

    # Numeric tolerance for cell-size + extent comparison. Real-world
    # geo data carries floating-point fuzz; an exact equality test
    # would warp every input even when grids match. 1e-7 in projected
    # units handles the typical 1mm-of-floating-point-noise case.
    _ALIGNMENT_TOLERANCE = 1e-7

    # QGIS data types that count as "integer" for the resampling
    # guard. Mirrors the Qgis.DataType enum values; we compare by
    # name to avoid pinning to a specific QGIS version's enum int.
    _INTEGER_DATA_TYPE_NAMES = frozenset({
        "Byte", "UInt16", "Int16", "UInt32", "Int32", "CInt16", "CInt32",
    })

    # Resampling methods that destroy categorical meaning. Reject for
    # integer-typed inputs unless the user explicitly opts in (future
    # `allow_categorical_resampling` flag — not exposed today).
    _NON_CATEGORICAL_RESAMPLINGS = frozenset({"bilinear", "cubic", "cubicspline", "lanczos"})

    # GDAL/Processing names for each high-level resampling choice.
    # The dialog stores names from RasterOverlaySpec; gdal:warpreproject
    # takes integer codes via the RESAMPLING parameter. The integer
    # codes here match QGIS's gdal:warpreproject documentation.
    _GDAL_RESAMPLING_CODES = {
        "nearest": 0,
        "bilinear": 1,
        "cubic": 2,
        "cubicspline": 3,
        "lanczos": 4,
        "average": 5,
        "mode": 6,
        "max": 7,
        "min": 8,
        "med": 9,
        "q1": 10,
        "q3": 11,
    }

    def _prealign_inputs(
        self,
        resolved_layers: "List[Tuple[LayerRef, Any]]",
        spec: RasterOverlaySpec,
    ) -> "List[Tuple[LayerRef, Any]]":
        """Warp misaligned inputs to the snap-to grid before r.cross.

        The snap-to layer is always first in `resolved_layers` (see
        `_ordered_input_refs`); every other input is compared against
        it and warped via gdal:warpreproject when CRS / extent / cell
        size differ. Returns a new resolved-layers list with the
        warped layers swapped in; the snap-to entry is unchanged.

        When grid properties can't be read off a layer (test contexts
        with duck-typed stand-ins), the alignment check defaults to
        "already aligned" — the path is intentionally narrow to the
        production case where QgsRasterLayer has the full API.
        """
        if not resolved_layers:
            return resolved_layers
        snap_ref, snap_layer = resolved_layers[0]
        snap_props = self._extract_grid_properties(snap_layer)

        out: "List[Tuple[LayerRef, Any]]" = [(snap_ref, snap_layer)]
        for ref, layer in resolved_layers[1:]:
            layer_props = self._extract_grid_properties(layer)
            if self._grids_aligned(snap_props, layer_props):
                out.append((ref, layer))
                continue
            # Misaligned — categorical guard then warp.
            self._guard_resampling_for_data_type(
                spec.resampling, layer, ref.name,
            )
            warped_layer = self._warp_to_snap(
                layer, snap_props, spec, ref.name,
            )
            out.append((ref, warped_layer))
        return out

    @staticmethod
    def _extract_grid_properties(layer: Any) -> Optional[Dict[str, Any]]:
        """Return CRS / extent / cell-size for a layer, or None if not exposed.

        Production: QgsRasterLayer has crs(), extent(),
        rasterUnitsPerPixelX/Y(). Tests with duck-typed fakes can stub
        these. When the layer doesn't expose the full API we return
        None and the alignment check defaults to "aligned" (no warp).
        """
        try:
            crs_attr = getattr(layer, "crs", None)
            extent_attr = getattr(layer, "extent", None)
            x_size_attr = getattr(layer, "rasterUnitsPerPixelX", None)
            y_size_attr = getattr(layer, "rasterUnitsPerPixelY", None)
            if not all((crs_attr, extent_attr, x_size_attr, y_size_attr)):
                return None
            crs = crs_attr()
            extent = extent_attr()
            return {
                "crs_authid": (
                    crs.authid() if hasattr(crs, "authid") else str(crs)
                ),
                "x_min": float(extent.xMinimum()),
                "y_min": float(extent.yMinimum()),
                "x_max": float(extent.xMaximum()),
                "y_max": float(extent.yMaximum()),
                "x_size": float(x_size_attr()),
                "y_size": float(y_size_attr()),
            }
        except Exception:
            return None

    @classmethod
    def _grids_aligned(
        cls,
        a: Optional[Dict[str, Any]],
        b: Optional[Dict[str, Any]],
    ) -> bool:
        """True if both grids share CRS + extent + cell size (within tolerance)."""
        if a is None or b is None:
            return True  # can't compare; defensive default
        if a.get("crs_authid") != b.get("crs_authid"):
            return False
        for key in ("x_min", "y_min", "x_max", "y_max", "x_size", "y_size"):
            if abs(float(a[key]) - float(b[key])) > cls._ALIGNMENT_TOLERANCE:
                return False
        return True

    @classmethod
    def _guard_resampling_for_data_type(
        cls,
        resampling: str,
        layer: Any,
        layer_name: str,
    ) -> None:
        """Reject bilinear/cubic for integer-typed (categorical) rasters.

        Those resamplings invent half-pixel values that destroy
        categorical meaning. For typical EMDS raster studies
        (land-cover codes, hazard classes) integer + nearest/mode is
        the correct combination.

        Detection is best-effort via QgsRasterDataProvider.dataType(1);
        unknown / non-QGIS layers default to "allow" so the guard
        doesn't false-positive in test contexts.
        """
        if resampling not in cls._NON_CATEGORICAL_RESAMPLINGS:
            return
        data_type_name = cls._read_data_type_name(layer)
        if data_type_name is None:
            return  # can't tell; defensive default
        if data_type_name in cls._INTEGER_DATA_TYPE_NAMES:
            raise RasterOverlayError(
                f"Input raster '{layer_name}' is an integer-typed "
                f"(categorical) raster ({data_type_name}) but the "
                f"raster overlay spec requested '{resampling}' "
                f"resampling. Bilinear / cubic / cubicspline / "
                f"lanczos invent half-pixel values that destroy "
                f"categorical meaning. Use 'nearest' (the default) "
                f"or 'mode' for categorical rasters."
            )

    @staticmethod
    def _read_data_type_name(layer: Any) -> Optional[str]:
        """Resolve the band-1 GDAL data-type name, or None if unavailable."""
        try:
            provider_attr = getattr(layer, "dataProvider", None)
            if not callable(provider_attr):
                return None
            provider = provider_attr()
            data_type_attr = getattr(provider, "dataType", None)
            if not callable(data_type_attr):
                return None
            raw = data_type_attr(1)
            # raw may be a Qgis.DataType enum or int. Try str() first.
            name = getattr(raw, "name", None)
            if name:
                return str(name)
            return str(raw)
        except Exception:
            return None

    def _warp_to_snap(
        self,
        layer: Any,
        snap_props: Optional[Dict[str, Any]],
        spec: RasterOverlaySpec,
        layer_name: str,
    ) -> Any:
        """Dispatch gdal:warpreproject to align `layer` onto the snap-to grid."""
        if snap_props is None:
            # No snap-to grid info; can't construct a meaningful warp
            # call. Return the original layer rather than failing —
            # r.cross will surface its own error if alignment matters.
            return layer
        resampling_code = self._GDAL_RESAMPLING_CODES.get(
            spec.resampling, 0,  # default to nearest
        )
        params = {
            "INPUT": layer,
            "TARGET_CRS": snap_props["crs_authid"],
            "TARGET_RESOLUTION": float(snap_props["x_size"]),
            "TARGET_EXTENT": (
                f"{snap_props['x_min']},{snap_props['x_max']},"
                f"{snap_props['y_min']},{snap_props['y_max']}"
                f" [{snap_props['crs_authid']}]"
            ),
            "RESAMPLING": resampling_code,
            "OUTPUT": "TEMPORARY_OUTPUT",
        }
        try:
            result = self._processing_run("gdal:warpreproject", params)
        except Exception as exc:
            raise RasterOverlayError(
                f"gdal:warpreproject failed while aligning input "
                f"'{layer_name}' to the snap-to grid: {exc}"
            ) from exc
        output_path = ""
        if isinstance(result, dict):
            for key in ("OUTPUT", "output", "Output"):
                value = result.get(key)
                if isinstance(value, str) and value:
                    output_path = value
                    break
        if not output_path:
            # Tests may stub processing.run to return a layer object
            # directly under OUTPUT. Pass that through.
            if isinstance(result, dict):
                for key in ("OUTPUT", "output", "Output"):
                    value = result.get(key)
                    if value is not None and not isinstance(value, str):
                        return value
            raise RasterOverlayError(
                f"gdal:warpreproject returned no usable OUTPUT for "
                f"input '{layer_name}'."
            )
        return self._open_raster_path(output_path)

    # ---- Validation ----------------------------------------------------

    @staticmethod
    def _validate_study(study: Study) -> None:
        if study.kind != LayerKind.RASTER:
            raise RasterOverlayError(
                f"RasterOverlay.run requires study.kind == RASTER; "
                f"got {study.kind.value}."
            )
        refs = study.input_layer_refs
        # Phase 3.5e: single-input raster studies are valid (the run()
        # method dispatches them to a passthrough path that skips
        # r.cross). Zero inputs are still an error.
        if len(refs) < 1:
            raise RasterOverlayError(
                "RasterOverlay.run requires at least 1 raster input layer; "
                "got 0."
            )
        non_raster = [r for r in refs if r.kind != LayerKind.RASTER]
        if non_raster:
            raise RasterOverlayError(
                "RasterOverlay.run found non-raster inputs in "
                f"input_layer_refs: {[r.name for r in non_raster]}. "
                "Study homogeneity violation — non-raster inputs belong "
                "in study.display_layers, not input_layer_refs."
            )
        if study.raster_overlay is None:
            raise RasterOverlayError(
                "RasterOverlay.run requires study.raster_overlay to be "
                "populated (snap_to / resampling / output_path)."
            )
        for ref in refs:
            if not ref.bands_to_carry and not ref.rat_columns_to_carry:
                raise RasterOverlayError(
                    f"Input raster '{ref.name}' carries neither a band "
                    "value nor any RAT columns into the output. Each "
                    "input must contribute at least one column."
                )

    @staticmethod
    def _ordered_input_refs(study: Study) -> List[LayerRef]:
        """Return inputs with snap_to first (it defines the output grid).

        r.cross treats inputs symmetrically but the dialog asked the
        user which input defines the snap-to grid; preserving that as
        the first element keeps the column ordering stable across
        re-runs and makes diff-debugging the RAT readable.
        """
        snap_to = study.raster_overlay.snap_to if study.raster_overlay else ""
        if not snap_to:
            return list(study.input_layer_refs)
        snap = [r for r in study.input_layer_refs if r.name == snap_to]
        rest = [r for r in study.input_layer_refs if r.name != snap_to]
        return snap + rest

    # ---- GRASS r.cross dispatch ---------------------------------------

    def _run_grass_r_cross(
        self,
        resolved_layers: List[Tuple[LayerRef, Any]],
        spec: RasterOverlaySpec,
    ) -> Dict[str, Any]:
        """Try the primary algorithm id, fall back to grass7: namespace.

        Parameter shape mirrors what processing.algorithmHelp("grass:r.cross")
        documents:
          - input: list of input raster layers
          - -z: ignore cells where any input is NULL (we pass False so
                NULL combinations get their own code, matching ArcGIS
                Combine semantics)
          - output: filesystem path for the unique-combination raster
          - GRASS_REGION_PARAMETER, GRASS_REGION_CELLSIZE_PARAMETER:
                set from the snap_to layer in the dialog
        """
        input_layers = [layer for _, layer in resolved_layers]
        parameters: Dict[str, Any] = {
            "input": input_layers,
            "-z": False,
            "output": spec.output_path,
            # GRASS-region overrides — left empty here means "use the
            # union extent of inputs at the snap_to cell size". The
            # dialog can override by populating raster_overlay with
            # explicit extent/cellsize. We leave that wiring to the
            # dialog so this primitive stays narrowly scoped to the
            # overlay step itself.
        }

        last_error: Optional[Exception] = None
        for algorithm_id in self._grass_algorithm_ids:
            try:
                return self._processing_run(algorithm_id, parameters)
            except Exception as exc:  # noqa: BLE001
                last_error = exc
                continue
        raise RasterOverlayError(
            f"GRASS r.cross failed under all known algorithm ids "
            f"({list(self._grass_algorithm_ids)}). "
            f"Confirm the GRASS provider is enabled in QGIS Processing. "
            f"Last error: {last_error}"
        ) from last_error

    @staticmethod
    def _extract_output_raster_path(
        cross_result: Dict[str, Any],
        spec: RasterOverlaySpec,
    ) -> str:
        """Pull the output path from the processing result, falling back to spec."""
        if isinstance(cross_result, dict):
            for key in ("output", "OUTPUT", "Output"):
                value = cross_result.get(key)
                if isinstance(value, str) and value:
                    return value
        return spec.output_path

    # ---- Category-table parsing + RAT enrichment ----------------------

    def _read_r_cross_category_table(
        self,
        cross_result: Dict[str, Any],
        output_raster_path: str,
        input_names: List[str],
    ) -> pd.DataFrame:
        """Pull r.cross's per-code category table into a DataFrame.

        GRASS r.cross writes the category mapping into the raster's
        category metadata. QGIS surfaces this via the output raster's
        RAT once the layer is loaded. We re-use RatReader for that.

        Schema after this step:
            Value | <input_1_name> | <input_2_name> | …
        where Value is the output code and each input column holds
        that input's pixel value for that code.

        If the Processing result already carries a parsed mapping
        (some builds do via a `CATEGORIES` key), use that instead
        of touching the raster.
        """
        # Fast path: some Processing builds emit a pre-parsed mapping.
        if isinstance(cross_result, dict):
            categories = cross_result.get("CATEGORIES") or cross_result.get("categories")
            if isinstance(categories, list) and categories:
                return self._categories_list_to_df(categories, input_names)

        # Standard path: read the RAT of the output raster.
        raster_layer = self._layer_resolver(output_raster_path) or self._open_raster_path(output_raster_path)
        df = self._rat_reader.read(raster_layer, band=1)
        if df is None or df.empty:
            raise RasterOverlayError(
                f"r.cross output raster '{output_raster_path}' has no "
                f"readable category table. Cannot enrich the RAT."
            )
        return self._normalize_category_df(df, input_names)

    @staticmethod
    def _categories_list_to_df(
        categories: List[Any],
        input_names: List[str],
    ) -> pd.DataFrame:
        """Convert a Processing CATEGORIES list into the normalized df.

        Two known shapes:
          1. list of dicts: [{"value": 1, "inputs": [11, 0]}, …]
          2. list of tuples: [(1, [11, 0]), …]
        """
        rows = []
        for entry in categories:
            if isinstance(entry, dict):
                value = entry.get("value")
                input_values = entry.get("inputs") or entry.get("values") or []
            elif isinstance(entry, (list, tuple)) and len(entry) == 2:
                value, input_values = entry
            else:
                continue
            row: Dict[str, Any] = {RAT_VALUE_COLUMN: value}
            for idx, name in enumerate(input_names):
                row[name] = input_values[idx] if idx < len(input_values) else None
            rows.append(row)
        return pd.DataFrame(rows, columns=[RAT_VALUE_COLUMN] + list(input_names))

    @staticmethod
    def _normalize_category_df(
        df: pd.DataFrame,
        input_names: List[str],
    ) -> pd.DataFrame:
        """Ensure the category df has Value + one column per input."""
        # Case-insensitive lookup for the value column; rename if needed.
        cols_lower = {col.lower(): col for col in df.columns}
        if RAT_VALUE_COLUMN.lower() in cols_lower:
            actual = cols_lower[RAT_VALUE_COLUMN.lower()]
            if actual != RAT_VALUE_COLUMN:
                df = df.rename(columns={actual: RAT_VALUE_COLUMN})
        else:
            # r.cross emits the category as a string "label"; we'd
            # need to parse it. Refuse here with a clear message.
            raise RasterOverlayError(
                f"Output RAT is missing a '{RAT_VALUE_COLUMN}' column. "
                f"Got columns: {list(df.columns)}."
            )

        # Ensure one column per input layer; r.cross typically names
        # them with a numeric or input-name convention. If the column
        # count matches, rename them to our canonical input names.
        non_value_cols = [c for c in df.columns if c != RAT_VALUE_COLUMN]
        if len(non_value_cols) == len(input_names) and non_value_cols != input_names:
            df = df.rename(columns=dict(zip(non_value_cols, input_names)))
        return df

    def _enrich_with_input_rats(
        self,
        category_df: pd.DataFrame,
        resolved_layers: List[Tuple[LayerRef, Any]],
    ) -> pd.DataFrame:
        """Left-join each input's selected RAT columns onto the category table.

        For each input layer:
          1. Read its RAT (may return None if the input is single-value).
          2. Filter the RAT to (Value, *ref.rat_columns_to_carry).
          3. Prefix every joined column with the input name so multi-input
             RATs that happen to share column names don't collide.
          4. Left-join on `category_df[input.name] == rat[Value]`.

        Inputs with no RAT and only bands_to_carry are passed through:
        their input column already holds the band value, which IS the
        carried column. No-op for them.
        """
        enriched = category_df.copy()

        for ref, layer in resolved_layers:
            if not ref.rat_columns_to_carry:
                continue

            input_rat = self._rat_reader.read(layer, band=1)
            if input_rat is None or input_rat.empty:
                # User asked for RAT columns but the input has no RAT;
                # this is a configuration error caught here rather than
                # producing silently-missing columns in the output.
                raise RasterOverlayError(
                    f"Input raster '{ref.name}' has no readable RAT but "
                    f"the study area asked to carry RAT columns "
                    f"{ref.rat_columns_to_carry}. Either provide a RAT "
                    f"for the input or unselect those columns."
                )

            cols_lower = {col.lower(): col for col in input_rat.columns}
            if RAT_VALUE_COLUMN.lower() not in cols_lower:
                raise RasterOverlayError(
                    f"Input raster '{ref.name}' RAT is missing the "
                    f"'{RAT_VALUE_COLUMN}' join key column."
                )
            value_col = cols_lower[RAT_VALUE_COLUMN.lower()]

            missing = [
                c for c in ref.rat_columns_to_carry
                if c not in input_rat.columns
            ]
            if missing:
                raise RasterOverlayError(
                    f"Input raster '{ref.name}' RAT is missing requested "
                    f"columns: {missing}. Available: {list(input_rat.columns)}."
                )

            slice_df = input_rat[[value_col] + list(ref.rat_columns_to_carry)].copy()
            slice_df = slice_df.rename(
                columns={
                    value_col: ref.name,  # join key matches enriched col
                    **{col: f"{ref.name}__{col}" for col in ref.rat_columns_to_carry},
                }
            )
            enriched = enriched.merge(slice_df, on=ref.name, how="left")

        return enriched

    # ---- Layer resolution helpers --------------------------------------

    @staticmethod
    def _open_raster_path(path: str) -> Any:
        """Open a raster file as a QgsRasterLayer (fallback when the
        layer isn't loaded in the project)."""
        try:
            from qgis.core import QgsRasterLayer  # type: ignore
        except ImportError as exc:
            raise RasterOverlayError(
                f"Cannot open raster '{path}': QGIS not available ({exc})."
            )
        layer = QgsRasterLayer(path, Path(path).stem)
        if not layer.isValid():
            raise RasterOverlayError(
                f"Could not open output raster '{path}'."
            )
        return layer


def _default_layer_resolver(name_or_path: str) -> Optional[Any]:
    """Look up a layer by name in the live QGIS project, or None.

    Resolver is injectable to keep RasterOverlay testable without QGIS.
    """
    try:
        from qgis.core import QgsProject  # type: ignore
    except ImportError:
        return None
    layers = QgsProject.instance().mapLayersByName(name_or_path)
    return layers[0] if layers else None
