"""
Validation framework for satellite-derived water extraction and
erosion/accretion analysis, using the JRC Global Surface Water
dataset as reference.
"""

import numpy as np


JRC_WATER_URL = (
    "https://storage.googleapis.com/global-surface-water/"
    "downloads/occurrence/occurrence_90w_60n_30s_0e.tif"
)


def compute_confusion_matrix(extracted_mask, reference_mask):
    """
    Compute pixel-wise confusion matrix metrics.

    Parameters
    ----------
    extracted_mask : np.ndarray (bool)
        Water mask from plugin.
    reference_mask : np.ndarray (bool)
        Reference water mask (e.g., JRC).

    Returns
    -------
    dict with TP, TN, FP, FN, OA, precision, recall, F1.
    """
    extracted = extracted_mask.astype('bool')
    reference = reference_mask.astype('bool')

    TP = np.sum(extracted & reference)
    TN = np.sum(~extracted & ~reference)
    FP = np.sum(extracted & ~reference)
    FN = np.sum(~extracted & reference)

    total = TP + TN + FP + FN
    OA = (TP + TN) / total if total > 0 else 0.0
    precision = TP / (TP + FP) if (TP + FP) > 0 else 0.0
    recall = TP / (TP + FN) if (TP + FN) > 0 else 0.0
    F1 = (
        2 * precision * recall / (precision + recall)
        if (precision + recall) > 0
        else 0.0
    )
    kappa_num = (
        (TP + TN) * total
        - (TP + FP) * (TP + FN)
        - (FN + TN) * (FP + TN)
    )
    kappa_den = (
        total ** 2
        - (TP + FP) * (TP + FN)
        - (FN + TN) * (FP + TN)
    )
    kappa = kappa_num / kappa_den if kappa_den > 0 else 0.0

    return {
        'TP': int(TP),
        'TN': int(TN),
        'FP': int(FP),
        'FN': int(FN),
        'OA': float(OA),
        'precision': float(precision),
        'recall': float(recall),
        'F1': float(F1),
        'kappa': float(kappa),
    }


def sample_agreement_points(plugin_gdf, reference_gdf,
                            n_samples=100, seed=42):
    """
    Generate random sample points and compare classifications.

    Parameters
    ----------
    plugin_gdf, reference_gdf : GeoDataFrame
    n_samples : int
    seed : int

    Returns
    -------
    dict with agreement_pct, kappa, and per-class metrics.
    """
    rng = np.random.default_rng(seed)

    bounds = plugin_gdf.total_bounds
    xs = rng.uniform(bounds[0], bounds[2], n_samples)
    ys = rng.uniform(bounds[1], bounds[3], n_samples)

    from shapely.geometry import Point
    import geopandas as gpd

    samples = gpd.GeoDataFrame(
        {'geometry': [Point(x, y) for x, y in zip(xs, ys)]},
        crs=plugin_gdf.crs
    )

    plugin_vals = _spatial_join_value(samples, plugin_gdf)
    ref_vals = _spatial_join_value(samples, reference_gdf)

    valid = (plugin_vals != -1) & (ref_vals != -1)
    if valid.sum() == 0:
        return {'error': 'No valid sample points'}

    plugin_vals = plugin_vals[valid]
    ref_vals = ref_vals[valid]

    agreement = (plugin_vals == ref_vals).sum() / len(plugin_vals)
    return {'agreement': float(agreement), 'n_samples': int(valid.sum())}


def _spatial_join_value(samples_gdf, polygon_gdf):
    """Check if each sample point intersects water polygons."""
    values = np.full(len(samples_gdf), -1, dtype='int32')
    for i, pt in enumerate(samples_gdf.geometry):
        hits = polygon_gdf[polygon_gdf.intersects(pt)]
        values[i] = 1 if not hits.empty else 0
    return values


def validate_erosion_accretion(
    plugin_erosion_gdf, plugin_accretion_gdf,
    reference_erosion_gdf=None,
    reference_accretion_gdf=None,
    buffer_m=30,
):
    """
    Validate erosion/accretion zones by spatial overlap analysis.

    Uses buffer-based overlap or direct comparison with reference
    erosion/accretion polygons (e.g., from manual interpretation
    or higher-resolution data).

    Parameters
    ----------
    plugin_erosion_gdf, plugin_accretion_gdf : GeoDataFrame
    reference_erosion_gdf, reference_accretion_gdf : GeoDataFrame or None
    buffer_m : float
        Buffer distance for overlap tolerance.

    Returns
    -------
    dict of validation metrics.
    """
    results = {}

    for label, plugin, reference in [
        ('erosion', plugin_erosion_gdf, reference_erosion_gdf),
        ('accretion', plugin_accretion_gdf, reference_accretion_gdf),
    ]:
        if reference is None or reference.empty:
            results[label] = {'validated': False,
                              'reason': 'No reference data'}
            continue
        if plugin is None or plugin.empty:
            results[label] = {'validated': False,
                              'reason': 'No plugin output'}
            continue

        try:
            plugin_union = plugin.geometry.unary_union
            ref_union = reference.geometry.unary_union

            if buffer_m > 0:
                plugin_union = plugin_union.buffer(buffer_m)

            overlap = plugin_union.intersection(ref_union)
            overlap_area = overlap.area if not overlap.is_empty else 0
            plugin_area = (
                plugin_union.area if not plugin_union.is_empty
                else 1
            )
            ref_area = ref_union.area if not ref_union.is_empty else 1

            overlap_pct_plugin = overlap_area / plugin_area * 100 \
                if plugin_area > 0 else 0
            overlap_pct_ref = overlap_area / ref_area * 100 \
                if ref_area > 0 else 0

            results[label] = {
                'validated': True,
                'overlap_area_m2': overlap_area,
                'plugin_area_m2': plugin_area,
                'reference_area_m2': ref_area,
                'overlap_pct_of_plugin': float(overlap_pct_plugin),
                'overlap_pct_of_reference': float(overlap_pct_ref),
            }
        except Exception as exc:
            results[label] = {'validated': False,
                              'error': str(exc)}

    return results
