import rasterio
import numpy as np
from rasterio.warp import calculate_default_transform, reproject, Resampling
from pathlib import Path

def read_raster_safe(path: Path, band: int = 1) -> tuple[np.ndarray, dict]:
    """Safe raster reading with nodata handling."""
    with rasterio.open(path) as src:
        data = src.read(band)
        profile = src.profile.copy()
        nodata = src.nodata
        if nodata is not None:
            data = np.where(data == nodata, np.nan, data.astype(float))
    return data, profile

def reproject_match(src_path: Path, ref_path: Path,
                    resampling: str = 'nearest') -> tuple[np.ndarray, dict]:
    """
    Reproject raster to the CRS and grid of a reference raster.
    IMPORTANT: classifications → use 'nearest' (not 'bilinear').
    Reference: Foody (2002) RSE 80(1):185
    """
    from rasterio.enums import Resampling as RS
    method = getattr(RS, resampling)
    with rasterio.open(ref_path) as ref:
        ref_crs, ref_transform = ref.crs, ref.transform
        ref_height, ref_width = ref.height, ref.width
    with rasterio.open(src_path) as src:
        transform, width, height = calculate_default_transform(
            src.crs, ref_crs, ref_width, ref_height,
            left=ref.bounds.left, bottom=ref.bounds.bottom,
            right=ref.bounds.right, top=ref.bounds.top)
        kwargs = src.meta.copy()
        kwargs.update({'crs': ref_crs, 'transform': transform,
                       'width': ref_width, 'height': ref_height})
        data = np.empty((src.count, ref_height, ref_width), dtype=src.dtypes[0])
        reproject(source=rasterio.band(src, list(range(1, src.count + 1))),
                  destination=data, src_transform=src.transform,
                  src_crs=src.crs, dst_transform=transform,
                  dst_crs=ref_crs, resampling=method)
    return data, kwargs

CRS_RECOMMENDATIONS = {
    'global_analysis':    'ESRI:54012',   # Eckert IV (equal-area)
    'south_america':      'EPSG:4087',    # World Equidistant Cylindrical
    'brazil_national':    'EPSG:4674',    # SIRGAS 2000
    'landsat_default':    'EPSG:4326',    # WGS84 geographic
    'sentinel_utm_z22s':  'EPSG:32722',   # UTM Zone 22S
}
# ⚠️ ALWAYS verify CRS before spatial operations

DN (Digital Number)
  ↓ radiometric calibration
Radiance (W/m²/sr/μm)
  ↓ divide by solar irradiance + cos(SZA)
TOA Reflectance (%)
  ↓ atmospheric correction (6S, Sen2Cor, LEDAPS, LaSRC)
Surface Reflectance (%)
  ↓ topographic correction (optional)
Normalised Reflectance

# Example: TOA reflectance for Landsat 8 OLI
def dn_to_toa_reflectance(dn: np.ndarray, M_rho: float, A_rho: float,
                           sun_elevation_deg: float) -> np.ndarray:
    """
    Convert DN → TOA Reflectance (Landsat Collection 2).
    M_rho, A_rho: coefficients from MTL file (REFLECTANCE_MULT/ADD_BAND_x)
    Reference: USGS Landsat Collection 2 Science Product Guide
    """
    toa = M_rho * dn + A_rho
    sun_elevation_rad = np.deg2rad(sun_elevation_deg)
    return toa / np.sin(sun_elevation_rad)

def compute_index(band_a: np.ndarray, band_b: np.ndarray,
                  formula: str = 'ndvi') -> np.ndarray:
    """Compute normalised spectral indices."""
    eps = 1e-8
    if formula == 'ndvi':     # (NIR - Red) / (NIR + Red)
        return (band_a - band_b) / (band_a + band_b + eps)
    elif formula == 'ndwi':   # (Green - NIR) / (Green + NIR)  Gao 1996
        return (band_b - band_a) / (band_b + band_a + eps)
    elif formula == 'nbr':    # (NIR - SWIR2) / (NIR + SWIR2)
        return (band_a - band_b) / (band_a + band_b + eps)

# Quick reference by category:
SPECTRAL_INDICES = {
    # Vegetation
    'NDVI':   '(NIR - Red) / (NIR + Red)',          # Tucker 1979
    'EVI':    '2.5*(NIR-Red)/(NIR+6*Red-7.5*Blue+1)', # Liu & Huete 1995
    'SAVI':   '(NIR-Red)/(NIR+Red+L)*(1+L)',        # Huete 1988, L=0.5
    'MSAVI':  '(2*NIR+1 - sqrt((2*NIR+1)²-8*(NIR-Red)))/2', # Qi 1994
    'NDRE':   '(RedEdge - Red) / (RedEdge + Red)',  # Gitelson 1994
    'LAI':    'from NDVI via empirical model',
    # Water
    'NDWI':   '(Green - NIR) / (Green + NIR)',      # McFeeters 1996
    'MNDWI':  '(Green - SWIR1) / (Green + SWIR1)', # Xu 2006
    'AWEI':   '4*(Green-SWIR1)-(0.25*NIR+2.75*SWIR2)', # Feyisa 2014
    # Urban / Bare Soil
    'NDBI':   '(SWIR1 - NIR) / (SWIR1 + NIR)',     # Zha 2003
    'BSI':    '((SWIR1+Red)-(NIR+Blue))/((SWIR1+Red)+(NIR+Blue))', # Rikimaru 2002
    # Fire
    'NBR':    '(NIR - SWIR2) / (NIR + SWIR2)',      # Key & Benson 2006
    'dNBR':   'pre_NBR - post_NBR',
    'RBR':    'dNBR / (pre_NBR + 1.001)',           # Parks 2014
    # SAR
    'VH/VV':  'VH / VV (log scale)',
    'RVI':    '4*VH / (VV + VH)',                   # Kim 2012
}

SAR Image (SLC / GRD)
  ↓ Radiometric calibration (sigma0 / beta0 / gamma0)
  ↓ Speckle filtering (Lee / Refined Lee / Gamma-MAP)
  ↓ Terrain correction (RTC — Range-Doppler / Ellipsoid)
  ↓ Convert to dB: 10 * log10(sigma0)
  ↓ Geocoding (target EPSG)
ARD Product (Analysis-Ready SAR)

from scipy.ndimage import uniform_filter
import numpy as np

def lee_filter(img: np.ndarray, size: int = 7) -> np.ndarray:
    """
    Lee filter for speckle reduction in SAR images.
    Reference: Lee (1980) IEEE Trans. Pattern Anal. Mach. Intell. 2(2):165
    """
    img_mean = uniform_filter(img, size)
    img_sq_mean = uniform_filter(img**2, size)
    img_variance = img_sq_mean - img_mean**2
    overall_variance = np.var(img)
    img_weights = img_variance / (img_variance + overall_variance)
    return img_mean + img_weights * (img - img_mean)

def image_differencing(img_t1: np.ndarray, img_t2: np.ndarray,
                        threshold: float = None) -> np.ndarray:
    """
    Change detection by simple image differencing.
    threshold: if None, uses mean + 2*std (automatic).
    """
    diff = np.abs(img_t2.astype(float) - img_t1.astype(float))
    if threshold is None:
        threshold = diff.mean() + 2 * diff.std()
    return diff > threshold

from scipy.signal import savgol_filter
import numpy as np

def savitzky_golay_smooth(ts: np.ndarray, window: int = 7,
                           polyorder: int = 3) -> np.ndarray:
    """
    Savitzky-Golay smoothing for vegetation time series.
    Preserves peaks and troughs better than moving averages.
    Reference: Chen et al. (2004) RSE 91(3-4):332
    """
    return savgol_filter(ts, window_length=window, polyorder=polyorder)

def harmonic_regression(t: np.ndarray, y: np.ndarray,
                         n_harmonics: int = 3) -> np.ndarray:
    """
    Harmonic regression for seasonality modelling (phenology).
    Reference: Zhu & Woodcock (2014) RSE 144:152
    """
    X = [np.ones_like(t), t]
    for k in range(1, n_harmonics + 1):
        X.append(np.cos(2 * np.pi * k * t / 365))
        X.append(np.sin(2 * np.pi * k * t / 365))
    X = np.column_stack(X)
    coeffs, _, _, _ = np.linalg.lstsq(X, y, rcond=None)
    return X @ coeffs

def temporal_frequency(stack: np.ndarray, class_value: int,
                        nodata: int = 0, axis: int = 0) -> np.ndarray:
    """
    Compute the relative frequency of a class in a time series.
    Reference: Defries et al. (2004); Hansen et al. (2013) Science
    """
    valid_mask = stack != nodata
    class_mask = stack == class_value
    return np.where(
        valid_mask.sum(axis=axis) > 0,
        class_mask.sum(axis=axis) / valid_mask.sum(axis=axis),
        np.nan)

def temporal_stability_index(stack: np.ndarray, nodata: int = 0,
                               axis: int = 0) -> np.ndarray:
    """
    Temporal stability index (dominance of the most frequent class).
    Reference: Pérez-Hoyos et al. (2017) Remote Sensing 9(1):2
    """
    from scipy.stats import mode
    dominant_class, dominant_count = mode(stack, axis=axis)
    valid_count = (stack != nodata).sum(axis=axis)
    return np.where(valid_count > 0, dominant_count / valid_count, np.nan)

# ⚠️ NEVER use random split with spatial data — spatial autocorrelation
# will inflate accuracy. Use block spatial cross-validation.

from sklearn.model_selection import KFold
import geopandas as gpd

def spatial_block_cv(gdf: gpd.GeoDataFrame, n_splits: int = 5,
                     block_size_km: float = 50.0):
    """
    Cross-validation with spatial blocks to avoid data leakage.
    Reference: Roberts et al. (2017) Ecography 40(8):913
    
    Alternative: use scikit-spatial-models or spatial-cross-val.
    """
    # Generates a block grid over the dataset extent
    # and creates folds based on blocks (not samples)
    raise NotImplementedError("See implementation by Roberts et al. 2017")

preprocessing:
  - [ ] Proper normalisation for RS (divide by 10000 for Landsat/S2 SR)
  - [ ] Balanced classes or class_weight handling
  - [ ] Spatially consistent augmentation (flips, rotations)
  - [ ] Bands in the correct order and documented

training:
  - [ ] Spatial cross-validation (not random split)
  - [ ] Per-class metrics: IoU, F1, precision, recall
  - [ ] Confusion matrix with % and absolutes
  - [ ] Learning curves plotted

validation:
  - [ ] Independent samples (not used in fit)
  - [ ] Stratified by class + region
  - [ ] Area-adjusted accuracy (Olofsson et al. 2014)

def majority_vote_ensemble(products: list[np.ndarray],
                            nodata: int = 0) -> np.ndarray:
    """
    Simple majority voting among multiple products.
    Reference: Herold et al. (2011) RSE 115(6):1559
    """
    stack = np.stack(products, axis=0)
    from scipy.stats import mode
    result, _ = mode(stack, axis=0)
    return result.squeeze()

def weighted_ensemble(products: list[np.ndarray],
                      weights: list[float]) -> np.ndarray:
    """
    Weighted ensemble by product confidence (one-hot + weighted sum).
    Reference: Fritz et al. (2011) ERL 6(4):044005
    """
    raise NotImplementedError("Requires one-hot encoding of categories")

import numpy as np
from sklearn.metrics import cohen_kappa_score

def overall_accuracy(y_true, y_pred):
    """Proportion of agreeing pixels."""
    return np.mean(y_true == y_pred)

def cohens_kappa(y_true, y_pred):
    """Chance-corrected kappa — Cohen (1960)."""
    return cohen_kappa_score(y_true, y_pred)

def dice_coefficient(mask_a: np.ndarray, mask_b: np.ndarray) -> float:
    """Dice/Sørensen coefficient for binary classes."""
    intersection = np.logical_and(mask_a, mask_b).sum()
    return 2 * intersection / (mask_a.sum() + mask_b.sum() + 1e-8)

def area_adjusted_accuracy(confusion_matrix: np.ndarray,
                            mapped_areas: np.ndarray,
                            sample_counts: np.ndarray) -> dict:
    """
    Area-adjusted accuracy and area estimates (design-based).
    
    Reference: Olofsson et al. (2014) RSE 148:42
    REQUIRED for land use/cover change studies.
    
    Args:
        confusion_matrix: (n_classes × n_classes) sample counts
        mapped_areas: mapped area per class (same units)
        sample_counts: samples collected per stratum
    
    Returns:
        dict with adjusted OA, UA, PA, F1 + confidence intervals
    """
    # Implementation following equations 1-5 of Olofsson et al. (2014)
    raise NotImplementedError("Use `area` package or reference implementation")

# Olofsson 2014 CHECKLIST:
# [ ] Stratified sampling (strata = mapped classes)
# [ ] Sample size per stratum documented
# [ ] Confusion matrix in area proportions (not raw counts)
# [ ] Uncertainty (95% CI) reported for OA, UA, PA
# [ ] Estimated area with standard error and CI

# Scientific literature
semantic_scholar_api: "https://api.semanticscholar.org/graph/v1/paper/search?query={TERM}&fields=title,authors,year,abstract,citationCount"
crossref_api:         "https://api.crossref.org/works?query={TERM}&filter=type:journal-article,from-pub-date:2018"
mdpi_remote_sensing:  "https://www.mdpi.com/search?q={TERM}&journal=remotesensing"
eartharxiv:           "https://eartharxiv.org/search/?q={TERM}"

# Product documentation
mapbiomas_github:     "https://github.com/mapbiomas"
copernicus_land:      "https://land.copernicus.eu/global"
glad_umd:             "https://glad.umd.edu/dataset"
esa_worldcover:       "https://esa-worldcover.org"

# Data and catalogs
earthengine_catalog:  "https://developers.google.com/earth-engine/datasets"
stac_planetary:       "https://planetarycomputer.microsoft.com/catalog"
stac_element84:       "https://earth-search.aws.element84.com/v1"
usgs_earthexplorer:   "https://earthexplorer.usgs.gov"
copernicus_hub:       "https://scihub.copernicus.eu"