# -*- coding: utf-8 -*-
"""Core Google Earth Engine helpers for the RUPP GEE LULC Toolkit.

Complete version:
- LULC classes: Tree → Tree/Forest, Agri → Cropland.
- Google Satellite Embedding (yearly) support with year spinbox.
- Gamma=2.0 applied after min/max stretch (configurable).
- Reset function to restore UI defaults on close.
- Example handlers for embedding toggle and build button.
"""

import json
import urllib.request
import os
import tempfile

from qgis.PyQt.QtCore import QSettings
from qgis.core import (
    QgsProject, QgsRasterLayer, QgsVectorLayer, QgsGeometry,
    QgsCoordinateReferenceSystem, QgsCoordinateTransform,
    QgsPalettedRasterRenderer,
    QgsContrastEnhancement,
    QgsMultiBandColorRenderer,
)
from qgis.PyQt.QtGui import QColor

# ---------------------------------------------------------------- classes
# (code, name, hex color)  - updated names
# Level 1 (classification level). Codes 1-6 keep their old meaning;
# 7 Shrubland is new, following the Cambodia LU2003 categories.
LULC_CLASSES = [
    (1, "Forest",      "#1a7431"),
    (2, "Agriculture", "#f2d43d"),
    (3, "Barren",      "#b08a5a"),
    (4, "Built-up",    "#d7191c"),
    (5, "Water",       "#2b7cb5"),
    (6, "Grassland",   "#a3d977"),
    (7, "Shrubland",   "#8f9e3d"),
    (8, "Other",       "#9e9e9e"),
]

# Indices computed from the 6 common optical bands
# (blue, green, red, nir, swir1, swir2) - works for S2 and Landsat 8/9.
OPTICAL_INDICES = [
    # vegetation
    "NDVI", "GNDVI", "EVI", "EVI2", "SAVI", "MSAVI", "OSAVI",
    "ARVI", "VARI", "GCVI",
    # water / moisture
    "NDWI", "MNDWI", "NDMI", "AWEI_NSH", "AWEI_SH", "WRI",
    # built-up / bare / burn
    "NDBI", "UI", "IBI", "BSI", "NBR",
    # tasseled cap
    "TCB", "TCG", "TCW",
]

# Sentinel-2 only (red-edge based)
S2_ONLY_INDICES = ["NDRE1", "NDRE2", "CIRE"]

SAR_INDICES = ["VH_VV", "RVI", "RFDI"]

# renamed band lists (order matters for temporal reducers)
S2_BAND_NAMES = ["blue", "green", "red", "re1", "re2", "re3",
                 "nir", "nirn", "swir1", "swir2"]
LS_BAND_NAMES = ["blue", "green", "red", "nir", "swir1", "swir2"]

# ESA WorldCover class labels (for sampling hints)
WC_LABELS = {
    10: "Tree cover", 20: "Shrubland", 30: "Grassland", 40: "Cropland",
    50: "Built-up", 60: "Bare/sparse", 70: "Snow/ice", 80: "Water",
    90: "Herb. wetland", 95: "Mangroves", 100: "Moss/lichen",
}

WGS84 = QgsCoordinateReferenceSystem("EPSG:4326")


# ---------------------------------------------------------------- init
def ee_project_id():
    return QSettings().value("rupp_gee_lulc/project_id", "", type=str)


def set_ee_project_id(pid):
    QSettings().setValue("rupp_gee_lulc/project_id", pid)


def init_ee():
    """Initialize Earth Engine. Authenticates on first use."""
    import ee
    pid = ee_project_id() or None
    try:
        ee.Initialize(project=pid)
    except (ImportError, ModuleNotFoundError, RuntimeError):  # GEE/Import
        ee.Authenticate()
        ee.Initialize(project=pid)
    return ee


# ---------------------------------------------------------------- retry
def _retry(fn, attempts=4, base=2.0):
    """Retry transient GEE/network errors (503, 429, timeouts) with
    exponential backoff, like the browser Code Editor does silently."""
    import time
    last = None
    for k in range(attempts):
        try:
            return fn()
        except Exception as e:
            last = e
            msg = str(e)
            transient = any(t in msg for t in (
                "503", "502", "500", "429", "Service Unavailable",
                "Too Many Requests", "timed out", "timeout",
                "Connection reset", "RemoteDisconnected"))
            if not transient or k == attempts - 1:
                raise
            time.sleep(base * (2 ** k))
    raise last


# ---------------------------------------------------------------- masking
def _mask_s2(img):
    scl = img.select("SCL")
    bad = scl.eq(3).Or(scl.eq(8)).Or(scl.eq(9)).Or(scl.eq(10)).Or(scl.eq(11))
    out = (img.updateMask(bad.Not())
              .select(["B2", "B3", "B4", "B5", "B6", "B7",
                       "B8", "B8A", "B11", "B12"],
                      ["blue", "green", "red", "re1", "re2", "re3",
                       "nir", "nirn", "swir1", "swir2"])
              .divide(10000.0))
    return out.copyProperties(img, ["system:time_start"])


def _mask_landsat(img):
    qa = img.select("QA_PIXEL")
    mask = qa.bitwiseAnd(0b11111).eq(0)   # fill/dilated/cirrus/cloud/shadow
    opt = (img.select(["SR_B2", "SR_B3", "SR_B4", "SR_B5", "SR_B6", "SR_B7"],
                      ["blue", "green", "red", "nir", "swir1", "swir2"])
              .multiply(0.0000275).add(-0.2))
    return opt.updateMask(mask).copyProperties(img, ["system:time_start"])


# ---------------------------------------------------------------- collections
def optical_collection(sensor, aoi, start, end, cloud):
    """sensor: 'S2', 'L8', 'L9', 'L89'."""
    import ee
    if sensor == "S2":
        col = (ee.ImageCollection("COPERNICUS/S2_SR_HARMONIZED")
               .filterBounds(aoi).filterDate(start, end)
               .filter(ee.Filter.lte("CLOUDY_PIXEL_PERCENTAGE", cloud))
               .map(_mask_s2))
        return col, 10
    ids = {"L8": ["LANDSAT/LC08/C02/T1_L2"],
           "L9": ["LANDSAT/LC09/C02/T1_L2"],
           "L89": ["LANDSAT/LC08/C02/T1_L2", "LANDSAT/LC09/C02/T1_L2"]}[sensor]
    col = ee.ImageCollection(ids[0])
    for i in ids[1:]:
        col = col.merge(ee.ImageCollection(i))
    col = (col.filterBounds(aoi).filterDate(start, end)
              .filter(ee.Filter.lte("CLOUD_COVER", cloud))
              .map(_mask_landsat))
    return col, 30


def s1_collection(aoi, start, end):
    import ee
    col = (ee.ImageCollection("COPERNICUS/S1_GRD")
           .filterBounds(aoi).filterDate(start, end)
           .filter(ee.Filter.eq("instrumentMode", "IW"))
           .filter(ee.Filter.listContains(
               "transmitterReceiverPolarisation", "VV"))
           .filter(ee.Filter.listContains(
               "transmitterReceiverPolarisation", "VH"))
           .select(["VV", "VH"]))
    return col, 10


def composite(col, method):
    return {"Median": col.median, "Mean": col.mean,
            "Min": col.min, "Max": col.max}[method]()


# ---------------------------------------------------------------- indices
# Tasseled cap coefficients: order blue, green, red, nir, swir1, swir2
_TC_LANDSAT = {   # Baig et al. (2014), Landsat 8 OLI
    "TCB": [0.3029, 0.2786, 0.4733, 0.5599, 0.5080, 0.1872],
    "TCG": [-0.2941, -0.2430, -0.5424, 0.7276, 0.0713, -0.1608],
    "TCW": [0.1511, 0.1973, 0.3283, 0.3407, -0.7117, -0.4559],
}
_TC_S2 = {        # Shi & Xu (2019), Sentinel-2 MSI
    "TCB": [0.3510, 0.3813, 0.3437, 0.7196, 0.2396, 0.1949],
    "TCG": [-0.3599, -0.3533, -0.4734, 0.6633, 0.0087, -0.2856],
    "TCW": [0.2578, 0.2305, 0.0883, 0.1071, -0.7611, -0.5308],
}


def add_optical_indices(img, names, sensor="S2"):
    """Add the requested indices as bands. sensor: 'S2' or 'L*'."""

    def b(n):
        return img.select(n)

    def nd(a, c):
        return lambda: img.normalizedDifference([a, c])

    def ex(expr, extra=None):
        d = {"B": b("blue"), "G": b("green"), "R": b("red"),
             "N": b("nir"), "S1": b("swir1"), "S2b": b("swir2")}
        if extra:
            d.update(extra)
        return lambda: img.expression(expr, d)

    def tc(key):
        c = (_TC_S2 if sensor == "S2" else _TC_LANDSAT)[key]
        def make():
            return (b("blue").multiply(c[0])
                    .add(b("green").multiply(c[1]))
                    .add(b("red").multiply(c[2]))
                    .add(b("nir").multiply(c[3]))
                    .add(b("swir1").multiply(c[4]))
                    .add(b("swir2").multiply(c[5])))
        return make

    F = {
        "NDVI": nd("nir", "red"),
        "GNDVI": nd("nir", "green"),
        "NDWI": nd("green", "nir"),
        "MNDWI": nd("green", "swir1"),
        "NDMI": nd("nir", "swir1"),
        "NBR": nd("nir", "swir2"),
        "NDBI": nd("swir1", "nir"),
        "UI": nd("swir2", "nir"),
        "EVI": ex("2.5*((N-R)/(N+6*R-7.5*B+1))"),
        "EVI2": ex("2.5*((N-R)/(N+2.4*R+1))"),
        "SAVI": ex("1.5*((N-R)/(N+R+0.5))"),
        "MSAVI": ex("(2*N+1-sqrt((2*N+1)*(2*N+1)-8*(N-R)))/2"),
        "OSAVI": ex("(N-R)/(N+R+0.16)"),
        "ARVI": ex("(N-(2*R-B))/(N+(2*R-B))"),
        "VARI": ex("(G-R)/(G+R-B)"),
        "GCVI": ex("N/G-1"),
        "AWEI_NSH": ex("4*(G-S1)-(0.25*N+2.75*S2b)"),
        "AWEI_SH": ex("B+2.5*G-1.5*(N+S1)-0.25*S2b"),
        "WRI": ex("(G+R)/(N+S1)"),
        "BSI": ex("((S1+R)-(N+B))/((S1+R)+(N+B))"),
        "IBI": ex(
            "(2*S1/(S1+N)-(N/(N+R)+G/(G+S1)))/"
            "(2*S1/(S1+N)+(N/(N+R)+G/(G+S1)))"),
        "TCB": tc("TCB"), "TCG": tc("TCG"), "TCW": tc("TCW"),
        # Sentinel-2 red-edge only
        "NDRE1": nd("nirn", "re1"),
        "NDRE2": nd("nirn", "re2"),
        "CIRE": lambda: b("nirn").divide(b("re1")).subtract(1),
    }
    out = img
    for n in names:
        if n in F:
            out = out.addBands(F[n]().rename(n).toFloat())
    return out


def add_sar_indices(img, names):
    import ee
    d = {}
    d["VH_VV"] = img.select("VH").subtract(img.select("VV"))   # dB ratio
    d["RFDI"] = (img.select("VV").subtract(img.select("VH"))
                 .divide(img.select("VV").add(img.select("VH"))))
    vv_lin = ee.Image(10).pow(img.select("VV").divide(10))
    vh_lin = ee.Image(10).pow(img.select("VH").divide(10))
    d["RVI"] = vh_lin.multiply(4).divide(vv_lin.add(vh_lin))
    d["RFDI"] = (img.select("VV").subtract(img.select("VH"))
                 .divide(img.select("VV").add(img.select("VH"))))
    out = img
    for n in names:
        if n in d:
            out = out.addBands(d[n].rename(n).toFloat())
    return out


# ---------------------------------------------------------------- WorldCover
def worldcover_image():
    import ee
    return ee.ImageCollection("ESA/WorldCover/v200").first().select("Map")


# ---------------------------------------------------------------- GAUL admin
_gaul_cache = {}


def gaul_cambodia():
    import ee
    return (ee.FeatureCollection("FAO/GAUL/2015/level0")
            .filter(ee.Filter.eq("ADM0_NAME", "Cambodia")))


def gaul_provinces():
    import ee
    return (ee.FeatureCollection("FAO/GAUL/2015/level1")
            .filter(ee.Filter.eq("ADM0_NAME", "Cambodia")))


def gaul_districts(province):
    import ee
    return (ee.FeatureCollection("FAO/GAUL/2015/level2")
            .filter(ee.Filter.eq("ADM0_NAME", "Cambodia"))
            .filter(ee.Filter.eq("ADM1_NAME", province)))


def list_province_names():
    if "prov" not in _gaul_cache:
        _gaul_cache["prov"] = sorted(
            gaul_provinces().aggregate_array("ADM1_NAME").getInfo())
    return _gaul_cache["prov"]


def list_district_names(province):
    key = "dist_" + province
    if key not in _gaul_cache:
        _gaul_cache[key] = sorted(
            gaul_districts(province).aggregate_array("ADM2_NAME").getInfo())
    return _gaul_cache[key]


# ---------------------------------------------------------------- geometry
def to_wgs84_geometry(geom, source_crs, max_vertices=8000):
    """QgsGeometry -> simplified WGS84 QgsGeometry (local, no network)."""
    g = QgsGeometry(geom)
    if source_crs != WGS84:
        ct = QgsCoordinateTransform(source_crs, WGS84, QgsProject.instance())
        g.transform(ct)
    tol = 0.0001   # ~11 m
    while (g.constGet() is not None
           and g.constGet().nCoordinates() > max_vertices and tol < 0.05):
        g = g.simplify(tol)
        tol *= 2
    return g


def qgis_geom_to_ee(geom, source_crs, max_vertices=8000):
    """QgsGeometry -> ee.Geometry (WGS84). Simplifies very detailed
    boundaries (e.g. shapefiles) so the GEE request stays small."""
    import ee
    g = to_wgs84_geometry(geom, source_crs, max_vertices)
    return ee.Geometry(json.loads(g.asJson(6)))


def _aoi_layer_from_geojson(gj, name):
    """Write a UNIQUE temp GeoJSON (avoids Windows file locks) and load it
    styled as a red boundary. Returns the layer or None."""
    import os
    import re
    import uuid
    if gj.get("type") == "GeometryCollection":
        # keep only the polygon parts (dissolve can return mixed parts)
        coords = []
        for part in gj.get("geometries", []):
            if part["type"] == "Polygon":
                coords.append(part["coordinates"])
            elif part["type"] == "MultiPolygon":
                coords.extend(part["coordinates"])
        if not coords:
            return None
        gj = {"type": "MultiPolygon", "coordinates": coords}
    fc = {"type": "FeatureCollection",
          "features": [{"type": "Feature", "properties": {"name": name},
                        "geometry": gj}]}
    safe = re.sub(r"[^A-Za-z0-9_-]+", "_", name)[:40] or "aoi"
    path = os.path.join(tempfile.gettempdir(), "aoi_{}_{}.geojson".format(
        safe, uuid.uuid4().hex[:8]))
    with open(path, "w", encoding="utf-8") as f:
        json.dump(fc, f)
    proj = QgsProject.instance()
    for lyr0 in list(proj.mapLayers().values()):
        if lyr0.name() == "AOI - " + name:        # replace older copies
            proj.removeMapLayer(lyr0.id())
    lyr = QgsVectorLayer(path, "AOI - " + name, "ogr")
    if not lyr.isValid():
        return None
    try:
        r = lyr.renderer()
        if r is not None and r.symbol() is not None:
            sym = r.symbol()
            sym.symbolLayer(0).setBrushStyle(0)              # no fill
            sym.symbolLayer(0).setStrokeWidth(0.7)
            sym.symbolLayer(0).setStrokeColor(QColor("#e31a1c"))
    except (AttributeError, TypeError, RuntimeError):  # QGIS geometry
        pass
    proj.addMapLayer(lyr)
    return lyr


def local_geometry_to_layer(wgs84_geom, name):
    """AOI boundary layer built from a LOCAL geometry - no GEE call."""
    gj = json.loads(wgs84_geom.asJson(6))
    return _aoi_layer_from_geojson(gj, name)


def ee_geometry_to_layer(ee_geom, name, simplify_m=100):
    """AOI boundary layer fetched from GEE, with retry on 5xx errors."""
    def _fetch():
        try:
            return ee_geom.dissolve(simplify_m).simplify(simplify_m).getInfo()
        except (AttributeError, TypeError, RuntimeError):  # QGIS geometry
            return ee_geom.simplify(simplify_m).getInfo()
    gj = _retry(_fetch)
    return _aoi_layer_from_geojson(gj, name)


# ---------------------------------------------------------------- layers
def add_ee_image_layer(image, vis, name):
    """(Keep for other visualizations – not used for stack loading)"""
    import ee
    mapid = _retry(lambda: ee.Image(image).getMapId(vis))
    url = mapid["tile_fetcher"].url_format
    uri = "type=xyz&zmin=0&zmax=20&url=" + url
    lyr = QgsRasterLayer(uri, name, "wms")
    if lyr.isValid():
        QgsProject.instance().addMapLayer(lyr)
    return lyr


def add_google_satellite():
    url = "https://mt1.google.com/vt/lyrs%3Ds%26x%3D{x}%26y%3D{y}%26z%3D{z}"
    uri = "type=xyz&zmin=0&zmax=21&url=" + url
    for lyr in QgsProject.instance().mapLayers().values():
        if lyr.name() == "Google Satellite":
            return lyr
    lyr = QgsRasterLayer(uri, "Google Satellite", "wms")
    if lyr.isValid():
        QgsProject.instance().addMapLayer(lyr)
        root = QgsProject.instance().layerTreeRoot()
        node = root.findLayer(lyr.id())
        if node:
            clone = node.clone()
            root.insertChildNode(-1, clone)
            node.parent().removeChildNode(node)
    return lyr


def _safe_download(url, path):
    """Download URL to path using only HTTP/HTTPS handlers (Bandit B310 safe)."""
    parsed = urllib.parse.urlparse(url)
    if parsed.scheme not in ("https", "http"):
        raise ValueError("Invalid URL scheme: {}".format(parsed.scheme))
    # Build opener with only HTTP/HTTPS handlers - blocks file:// etc.
    opener = urllib.request.build_opener(
        urllib.request.HTTPHandler,
        urllib.request.HTTPSHandler,
        urllib.request.HTTPErrorProcessor,
        urllib.request.HTTPDefaultErrorHandler,
        urllib.request.HTTPRedirectHandler,
    )
    with opener.open(url) as response:
        with open(path, "wb") as out_file:
            out_file.write(response.read())
    return path


def download_geotiff(image, region, scale, path, bands=None):
    img = image
    if bands:
        img = image.select(bands)
    def _do():
        url = img.getDownloadURL({
            "region": region, "scale": scale,
            "format": "GEO_TIFF", "filePerBand": False})
        return _safe_download(url, path)
    return _retry(_do)


# Water recolored to match the Water bodies vector layer
WATER_VECTOR_COLOR = "#96c8e6"


def _lulc_display_classes():
    out = []
    for c, n, h in LULC_CLASSES:
        if n.lower().startswith("water"):
            h = WATER_VECTOR_COLOR
        out.append((c, n, h))
    return out


def apply_lulc_style(raster_layer, scheme="lulc"):
    if scheme == "lulc":
        cls = _lulc_display_classes()
    else:
        cls = scheme_classes(scheme)
    classes = [QgsPalettedRasterRenderer.Class(c, QColor(h), n)
               for c, n, h in cls]
    r = QgsPalettedRasterRenderer(raster_layer.dataProvider(), 1, classes)
    raster_layer.setRenderer(r)
    raster_layer.triggerRepaint()


def apply_lu2003_style(raster_layer):
    """Paletted renderer with the 40 LU2003 classes."""
    from qgis.core import QgsPalettedRasterRenderer
    classes = [QgsPalettedRasterRenderer.Class(
        c, QColor(lu2003_color(c)), "{} {}".format(LU2003[c][1],
                                                   LU2003[c][2]))
        for c in sorted(LU2003)]
    r = QgsPalettedRasterRenderer(raster_layer.dataProvider(), 1, classes)
    raster_layer.setRenderer(r)
    raster_layer.triggerRepaint()


def lulc_palette():
    # use the SAME display colours as the final styled layer, so the GEE
    # preview tiles and the downloaded "LULC Classes" layer match (water is
    # recoloured to the vector water colour in both).
    return [h.lstrip("#") for _, _, h in _lulc_display_classes()]


# ---------------------------------------------------------------- base bands
S2_BASE = ["blue", "green", "red", "re1", "re2", "re3",
           "nir", "nirn", "swir1", "swir2"]
LS_BASE = ["blue", "green", "red", "nir", "swir1", "swir2"]

# Preset GEE assets for AOI selection (Sarath's account)
AOI_ASSETS = [
    "users/yamsarath168/2026/KH_Khet",
    "users/yamsarath168/2026/KH_Srok",
    "users/yamsarath168/2026/KH_Khum",
    "users/yamsarath168/Kh_Admin/Kh_Country_area",
]


# ---------------------------------------------------------------- temporal
def add_percentile_bands(stack, col, base_bands):
    """P20 (dry-season minimum) and P80 (wet-season peak) composites."""
    import ee
    p20 = (col.reduce(ee.Reducer.percentile([20]))
           .rename([b + "_p20" for b in base_bands]))
    p80 = (col.reduce(ee.Reducer.percentile([80]))
           .rename([b + "_p80" for b in base_bands]))
    return stack.addBands(p20.toFloat()).addBands(p80.toFloat())


def add_std_bands(stack, col, base_bands):
    """Per-band temporal standard deviation (variability)."""
    import ee
    std = (col.reduce(ee.Reducer.stdDev())
           .rename([b + "_std" for b in base_bands]))
    return stack.addBands(std.toFloat())


def add_temporal_ndvi(stack, col):
    """NDVI mean / std / p20 / p80 — paddy rice phenology signature."""
    import ee

    def _ndvi(i):
        return (i.normalizedDifference(["nir", "red"]).rename("NDVI")
                .copyProperties(i, ["system:time_start"]))
    nc = col.map(_ndvi)
    out = stack
    out = out.addBands(nc.reduce(ee.Reducer.mean())
                       .rename("NDVI_mean").toFloat())
    out = out.addBands(nc.reduce(ee.Reducer.stdDev())
                       .rename("NDVI_std").toFloat())
    out = out.addBands(nc.reduce(ee.Reducer.percentile([20]))
                       .rename("NDVI_p20").toFloat())
    out = out.addBands(nc.reduce(ee.Reducer.percentile([80]))
                       .rename("NDVI_p80").toFloat())
    return out


def landsat_lst(aoi, start, end):
    """Landsat 8+9 median LST in Kelvin (ST_B10 USGS scaling)."""
    import ee

    def _mask(img):
        qa = img.select("QA_PIXEL")
        return img.updateMask(
            qa.bitwiseAnd(1 << 3).eq(0).And(qa.bitwiseAnd(1 << 4).eq(0)))
    col = (ee.ImageCollection("LANDSAT/LC08/C02/T1_L2")
           .merge(ee.ImageCollection("LANDSAT/LC09/C02/T1_L2"))
           .filterBounds(aoi).filterDate(start, end).map(_mask))
    return (col.select("ST_B10").median()
            .multiply(0.00341802).add(149).rename("LST").toFloat())


def terrain_bands():
    """NASADEM elevation + slope."""
    import ee
    dem = ee.Image("NASA/NASADEM_HGT/001").select("elevation")
    return dem.rename("elevation").addBands(
        ee.Terrain.slope(dem).rename("slope")).toFloat()


# ---------------------------------------------------------------- embeddings
def embedding_image(aoi, year):
    """Google Satellite Embedding annual mosaic (64 bands A00..A63).
    Available from 2017 onward."""
    import ee
    col = (ee.ImageCollection("GOOGLE/SATELLITE_EMBEDDING/V1/ANNUAL")
           .filterBounds(aoi)
           .filterDate("{}-01-01".format(year), "{}-01-01".format(year + 1)))
    return col.mosaic(), col.size()


# ---------------------------------------------------------------- GEE assets
def asset_fc(path):
    import ee
    return ee.FeatureCollection(path)


def asset_field_names(path):
    fc = asset_fc(path)
    names = fc.first().propertyNames().getInfo() or []
    return sorted(n for n in names if not n.startswith("system:"))


# Cambodia admin code fields per boundary asset (unique keys)
ADMIN_CODE_FIELDS = ["KHUM_CODE", "SROK_CODE", "KHET_CODE",
                     "COM_CODE", "DIS_CODE", "PRO_CODE", "GID", "OBJECTID"]
# parent-name fields used to disambiguate duplicate names
ADMIN_PARENT_FIELDS = ["KHET_NAM_1", "KHET_NAME", "PRO_NAME", "KHET_NAM",
                       "SROK_NAM_1", "SROK_NAME", "DIS_NAME"]


def asset_code_field(path):
    """Best unique-code field present in the asset, else None."""
    names = set(asset_field_names(path))
    for c in ADMIN_CODE_FIELDS:
        if c in names:
            return c
    return None


def asset_parent_field(path, exclude=None):
    names = set(asset_field_names(path))
    for c in ADMIN_PARENT_FIELDS:
        if c in names and c != exclude:
            return c
    return None


def asset_field_values(path, field):
    fc = asset_fc(path)
    vals = fc.aggregate_array(field).distinct().sort().getInfo() or []
    return [str(v) for v in vals]


def asset_value_rows(path, name_field):
    """Return [(label, name, code), ...] where label disambiguates
    duplicate names using the parent province/district and the unique
    code. Solves the 'Chres x4 shows once / merges all' problem."""
    fc = asset_fc(path)
    code_field = asset_code_field(path)
    parent_field = asset_parent_field(path, exclude=name_field)
    props = [name_field]
    if code_field:
        props.append(code_field)
    if parent_field:
        props.append(parent_field)

    def _fetch():
        fc2 = fc.select(props)
        return fc2.toList(fc2.size()).getInfo()
    feats = _retry(_fetch)

    rows = []
    seen_names = {}
    for f in feats:
        pr = f.get("properties", {})
        name = str(pr.get(name_field, "")).strip()
        if not name:
            continue
        code = str(pr.get(code_field, "")) if code_field else ""
        parent = str(pr.get(parent_field, "")) if parent_field else ""
        seen_names[name] = seen_names.get(name, 0) + 1
        rows.append([name, code, parent])

    out = []
    for name, code, parent in rows:
        label = name
        if seen_names[name] > 1:            # duplicate -> disambiguate
            if parent:
                label = "{} ({})".format(name, parent)
            if code:
                label = "{}  [{}]".format(label, code)
        elif code:
            label = "{}  [{}]".format(name, code)
        out.append((label, name, code))
    out.sort(key=lambda r: r[0])
    return out, code_field


# ---------------------------------------------------------------- classifiers
CLASSIFIERS = ["Random Forest", "Gradient Tree Boost", "CART", "SVM (RBF)"]


def make_classifier(name, trees=300):
    import ee
    if name == "Random Forest":
        return ee.Classifier.smileRandomForest(trees)
    if name == "Gradient Tree Boost":
        return ee.Classifier.smileGradientTreeBoost(
            numberOfTrees=trees, shrinkage=0.05,
            samplingRate=0.7, maxNodes=16)
    if name == "CART":
        return ee.Classifier.smileCart(maxNodes=32)
    return ee.Classifier.libsvm(kernelType="RBF", gamma=0.5, cost=10)


# ---------------------------------------------------------------- local stack
DOWNLOAD_LIMIT_BYTES = 45 * 1024 * 1024   # stay under the ~48 MB GEE limit


def fit_download_scale(aoi, n_bands, native_scale):
    """Coarsest-necessary scale (m) so a multiband GeoTIFF fits the GEE
    direct-download limit. Returns native scale when it already fits."""
    import math
    area = aoi.area(1).getInfo()             # m^2
    needed = math.sqrt(max(area, 1.0) * n_bands * 4.0 / DOWNLOAD_LIMIT_BYTES)
    return max(int(native_scale), int(math.ceil(needed)))


def apply_rgb_style(layer, band_names, kind, gamma=2.0):
    """Multiband color rendering with min/max stretch + gamma (brightness
    filter). User can change R/G/B bands in Layer Properties > Symbology."""

    def bidx(n, fallback=1):
        return band_names.index(n) + 1 if n in band_names else fallback

    if kind == "sar":
        r, g, b, mn, mx = bidx("VV"), bidx("VH", 2), bidx("VV"), -25.0, 0.0
    elif kind == "embedding":
        r, g, b, mn, mx = (bidx("A01"), bidx("A16", 2), bidx("A09", 3),
                           -0.3, 0.3)
    else:
        r, g, b, mn, mx = (bidx("red"), bidx("green", 2), bidx("blue", 3),
                           0.02, 0.35)
    prov = layer.dataProvider()
    ren = QgsMultiBandColorRenderer(prov, r, g, b)
    for setter, band in ((ren.setRedContrastEnhancement, r),
                         (ren.setGreenContrastEnhancement, g),
                         (ren.setBlueContrastEnhancement, b)):
        ce = QgsContrastEnhancement(prov.dataType(band))
        ce.setContrastEnhancementAlgorithm(
            QgsContrastEnhancement.ContrastEnhancementAlgorithm.StretchToMinimumMaximum)
        ce.setMinimumValue(mn)
        ce.setMaximumValue(mx)
        setter(ce)
    layer.setRenderer(ren)
    try:
        # gamma lives on the brightness/contrast filter in QGIS
        layer.brightnessFilter().setGamma(float(gamma))
    except (AttributeError, TypeError, RuntimeError):  # QGIS geometry
        pass
    layer.triggerRepaint()


# ---------------------------------------------------------------- Drive export
DRIVE_FOLDER = "GEE_folder"


def export_to_drive(image, region, scale, name, folder=DRIVE_FOLDER):
    """Start an Export.image.toDrive task. Returns the task."""
    import ee
    import re
    safe = re.sub(r"[^A-Za-z0-9_-]+", "_", name)[:95] or "gee_export"
    task = ee.batch.Export.image.toDrive(
        image=image, description=safe, folder=folder,
        fileNamePrefix=safe, region=region, scale=scale,
        maxPixels=1e13, fileFormat="GeoTIFF")
    task.start()
    return safe


# ---------------------------------------------------------------- band preset
# Curated ~30-band set (Phnom Penh 2025 paper, REDUCED_BANDS style)
IMPORTANT_BANDS = [
    # median optical
    "blue", "green", "red", "re1", "nir", "swir1", "swir2",
    # percentile composites (dry/wet season)
    "blue_p20", "green_p80", "red_p20", "red_p80",
    "swir1_p20", "swir1_p80", "swir2_p20", "swir2_p80",
    # temporal variability
    "swir1_std", "swir2_std", "NDVI_std",
    # temporal NDVI
    "NDVI", "NDVI_mean", "NDVI_p80",
    # vegetation + red-edge indices
    "EVI", "NDRE1", "NDRE2", "CIRE",
    # water + urban indices
    "AWEI_NSH", "NDWI", "NDBI", "BSI",
    # SAR (fusion only)
    "VV", "VH", "RFDI",
    # extras
    "LST", "elevation", "slope",
]
# indices that must be ticked for the preset
IMPORTANT_INDEX_CHECKS = {"NDVI", "EVI", "NDWI", "AWEI_NSH", "NDBI", "BSI",
                          "NDRE1", "NDRE2", "CIRE", "RFDI"}


def select_important_bands(stack):
    """Keep only the curated bands that exist in the stack (server-side)."""
    import ee
    names = stack.bandNames().filter(
        ee.Filter.inList("item", ee.List(IMPORTANT_BANDS)))
    return stack.select(names)


# ---------------------------------------------------------------- tiled DL
def download_stack_smart(image, aoi, scale, path, n_bands,
                         limit=DOWNLOAD_LIMIT_BYTES, max_tiles=25,
                         progress=None):
    """Download a multiband GeoTIFF at native scale. If the request is
    bigger than the GEE per-request limit (~48 MB), the AOI bounding box
    is split into tiles, each tile is downloaded, and GDAL merges them
    into one file. Returns the number of tiles used."""
    import ee
    import math

    bbox = aoi.bounds(1)
    area = _retry(lambda: bbox.area(1).getInfo())                       # m^2
    scale_use = int(scale)
    total_bytes = (area / (scale_use * scale_use)) * n_bands * 4.0
    n = max(1, int(math.ceil(math.sqrt(total_bytes / (limit * 0.85)))))
    if n * n > max_tiles:
        # auto-resample the LOCAL COPY so it fits into max_tiles requests;
        # the in-memory stack / Drive export stay at native resolution
        factor = math.sqrt((n * n) / float(max_tiles))
        scale_use = int(math.ceil(scale_use * factor))
        total_bytes = (area / (scale_use * scale_use)) * n_bands * 4.0
        n = max(1, int(math.ceil(math.sqrt(total_bytes / (limit * 0.85)))))

    if n == 1:
        if progress:
            progress(0.1, "Requesting GeoTIFF from GEE ...")
        download_geotiff(image, aoi, scale_use, path)
        if progress:
            progress(1.0, "Download finished.")
        return 1, scale_use

    coords = bbox.coordinates().getInfo()[0]
    xs = [c[0] for c in coords]
    ys = [c[1] for c in coords]
    xmin, xmax, ymin, ymax = min(xs), max(xs), min(ys), max(ys)
    dx = (xmax - xmin) / n
    dy = (ymax - ymin) / n

    parts = []
    k = 0
    for i in range(n):
        for j in range(n):
            k += 1
            if progress:
                progress((k - 1) / float(n * n),
                         "Downloading tile {} / {} ...".format(k, n * n))
            rect = ee.Geometry.Rectangle(
                [xmin + i * dx, ymin + j * dy,
                 xmin + (i + 1) * dx, ymin + (j + 1) * dy])
            part = "{}_part_{}_{}.tif".format(path[:-4], i, j)
            download_geotiff(image, rect, scale_use, part)
            parts.append(part)

    if progress:
        progress(0.97, "Merging {} tiles ...".format(len(parts)))
    from osgeo import gdal
    import os as _os
    vrt_path = path[:-4] + ".vrt"
    vrt = gdal.BuildVRT(vrt_path, parts)
    gdal.Translate(path, vrt,
                   creationOptions=["COMPRESS=LZW", "BIGTIFF=IF_SAFER"])
    vrt = None
    try:
        _os.remove(vrt_path)
        for p in parts:
            _os.remove(p)
    except (OSError, IOError):  # File operation
        pass
    return n * n, scale_use


# ---------------------------------------------------------------- GEE assets
ASSET_STACK_ROOT = "users/yamsarath168/QGIS_stacks"


def load_asset_image(asset_id):
    """Load an exported stack back: returns (image, band_names, scale,
    geometry, kind-guess)."""
    import ee
    img = ee.Image(asset_id)
    names = img.bandNames().getInfo()
    try:
        scale = int(round(img.projection().nominalScale().getInfo()))
        if scale <= 0 or scale > 10000:
            scale = 10
    except (AttributeError, TypeError, RuntimeError):  # QGIS geometry
        scale = 10
    if names and all(n.startswith("A") and n[1:].isdigit() for n in names):
        kind = "embedding"
    elif "VV" in names and "red" not in names:
        kind = "sar"
    else:
        kind = "optical"
    return img, names, scale, img.geometry(), kind


# ---------------------------------------------------------------- Asset export
def ensure_asset_folder(folder):
    """Create the asset folder (and missing parents) if needed."""
    import ee
    parts = folder.strip("/").split("/")
    # walk down from the account root, creating missing folders
    # e.g. users/yamsarath168 always exists; create deeper levels
    for i in range(2, len(parts)):
        path = "/".join(parts[:i + 1])
        try:
            ee.data.getAsset(path)
            continue
        except (OSError, IOError):  # File operation
            pass
        try:
            ee.data.createAsset({"type": "FOLDER"}, path)
        except Exception as e:
            raise RuntimeError(
                "Could not create asset folder '{}': {}".format(path, e))


def export_to_asset(image, region, scale, asset_id):
    """Start an Export.image.toAsset task. Returns the task description.
    Creates the parent folder automatically when it does not exist."""
    import ee
    import re
    parent = asset_id.rsplit("/", 1)[0]
    ensure_asset_folder(parent)
    desc = re.sub(r"[^A-Za-z0-9_-]+", "_",
                  asset_id.rsplit("/", 1)[-1])[:95] or "gee_stack"
    task = ee.batch.Export.image.toAsset(
        image=image, description=desc, assetId=asset_id,
        region=region, scale=scale, maxPixels=1e13)
    task.start()
    return desc


def load_image_asset(asset_id):
    """Load an image asset back: returns (image, band_names, scale, kind)."""
    import ee
    img = ee.Image(asset_id)
    band_names = img.bandNames().getInfo()
    try:
        scale = int(round(
            img.select(0).projection().nominalScale().getInfo()))
    except (AttributeError, TypeError, RuntimeError):  # QGIS geometry
        scale = 10
    if any(b.startswith("A0") or b.startswith("A1") for b in band_names[:3]):
        kind = "embedding"
    elif "red" in band_names:
        kind = "optical"
    elif "VV" in band_names:
        kind = "sar"
    else:
        kind = "optical"
    return img, band_names, scale, kind


# ---------------------------------------------------------------- LU2003
# Cambodia Land Use 2003 (MLMUPC/JICA) - 40 Level-2 classes.
# lu_code: (level1_code, abbreviation, name)
LU2003 = {
    1:  (4, "U",   "Settlement"),
    2:  (4, "I",   "Infrastructure (airfield, playground, cemetery)"),
    3:  (2, "Ar",  "Paddy field"),
    4:  (2, "Al",  "Receding and floating rice fields"),
    5:  (2, "Au",  "Field crop"),
    6:  (2, "As",  "Swidden agriculture (slash and burn)"),
    7:  (2, "Ao",  "Orchard"),
    8:  (2, "Ap",  "Plantation (rubber plantation)"),
    9:  (2, "Av",  "Village garden crop"),
    10: (2, "Ag",  "Garden crop"),
    11: (2, "Arv", "Paddy field with villages"),
    12: (6, "G",   "Grassland (undifferentiated)"),
    13: (6, "Ga",  "Abandoned field covered by grass"),
    14: (6, "Gf",  "Flooded grassland"),
    15: (6, "Gs",  "Grass savannah"),
    16: (6, "Gm",  "Grass with termite mounds"),
    17: (6, "Ms",  "Marsh and swamp"),
    18: (7, "S",   "Shrubland (undifferentiated)"),
    19: (7, "Sa",  "Abandoned field covered by shrub"),
    20: (7, "Sf",  "Flooded shrub"),
    21: (7, "St",  "Woodland and scattered trees (canopy < 10%)"),
    22: (1, "Fe",  "Evergreen broad leafed forest"),
    23: (1, "Fc",  "Coniferous forest"),
    24: (1, "Fd",  "Deciduous forest"),
    25: (1, "Fdo", "Dry deciduous (open) forest"),
    26: (1, "Fx",  "Mixed evergreen and deciduous forest"),
    27: (1, "Fr",  "Riparian forest"),
    28: (1, "Fs",  "Bamboo and secondary forests"),
    29: (1, "Ff",  "Flooded forest"),
    30: (1, "Fm",  "Mangrove forest"),
    31: (1, "Fmd", "Degraded mangrove forest"),
    32: (1, "Fp",  "Forest plantation"),
    33: (5, "Wl",  "Lakes (> 8 ha)"),
    34: (5, "Wp",  "Lakes (< 8 ha)"),
    35: (5, "Wr",  "Reservoir"),
    36: (5, "Ws",  "Shrimp/fish farming and salt pan"),
    37: (5, "Wo",  "Other water (sea, bay)"),
    38: (3, "B",   "Barren land"),
    39: (3, "Bs",  "Sand bank"),
    40: (3, "Br",  "Rock outcrop"),
}


def l1_name(code):
    for c, n, _ in LULC_CLASSES:
        if c == code:
            return n
    return str(code)


def l2_classes_for(l1_code):
    """[(lu_code, abbr, name)] belonging to one Level-1 class."""
    return [(k, v[1], v[2]) for k, v in sorted(LU2003.items())
            if v[0] == l1_code]


def _shade(hex_color, frac):
    """Blend a hex color toward white by frac (0..1)."""
    h = hex_color.lstrip("#")
    r, g, b = (int(h[i:i + 2], 16) for i in (0, 2, 4))
    r = int(r + (255 - r) * frac)
    g = int(g + (255 - g) * frac)
    b = int(b + (255 - b) * frac)
    return "#{:02x}{:02x}{:02x}".format(r, g, b)


def lu2003_color(lu_code):
    """Level-2 color: shade of the parent Level-1 color."""
    l1 = LU2003[lu_code][0]
    parent = dict((c, h) for c, _, h in LULC_CLASSES)[l1]
    siblings = [k for k in sorted(LU2003) if LU2003[k][0] == l1]
    i = siblings.index(lu_code)
    return _shade(parent, min(0.62, 0.10 * i))


def lu2003_palette():
    """Palette for lu_code values 1..40 (for GEE vis min=1 max=40)."""
    return [lu2003_color(c).lstrip("#") for c in sorted(LU2003)]


# ---------------------------------------------------------------- band names
def set_gdal_band_descriptions(tif_path, band_names):
    """Write real band names INTO the GeoTIFF (GDAL band descriptions).
    This is what makes QGIS show 'NDVI' instead of 'Band 1' everywhere -
    Symbology dropdowns, Identify, and the layer tree."""
    try:
        from osgeo import gdal
        ds = gdal.Open(tif_path, gdal.GA_Update)
        if ds is None:
            return False
        for i, name in enumerate(band_names, start=1):
            if i <= ds.RasterCount:
                ds.GetRasterBand(i).SetDescription(name)
        ds.FlushCache()
        ds = None
        return True
    except (OSError, IOError):  # File operation
        return False


def write_band_qml(tif_path, band_names, kind, gamma=2.0):
    """Write a .qml sidecar so QGIS shows real band names (NDVI, SAVI, ...)
    instead of 'Band 20'. Also sets a sensible multiband RGB default."""
    def bidx(n, fb=1):
        return band_names.index(n) + 1 if n in band_names else fb
    if kind == "sar":
        r, g, b = bidx("VV"), bidx("VH", 2), bidx("VV")
    elif kind == "embedding":
        r, g, b = bidx("A01"), bidx("A16", 2), bidx("A09", 3)
    else:
        r, g, b = bidx("red"), bidx("green", 2), bidx("blue", 3)

    entries = "\n".join(
        '   <ChannelName band="{}">{}</ChannelName>'.format(i, n)
        for i, n in enumerate(band_names, start=1))
    qml = '''<!DOCTYPE qgis>
<qgis version="3.28">
  <pipe>
    <rasterrenderer type="multibandcolor" redBand="{r}" greenBand="{g}" blueBand="{b}" opacity="1">
      <bandNames>
{entries}
      </bandNames>
    </rasterrenderer>
    <brightnesscontrast brightness="0" contrast="0" gamma="{gamma}"/>
  </pipe>
</qgis>'''.format(r=r, g=g, b=b, entries=entries, gamma=gamma)
    with open(tif_path[:-4] + ".qml", "w", encoding="utf-8") as fh:
        fh.write(qml)


def set_raster_band_names(layer, band_names):
    """Rename bands on a loaded raster layer so the Symbology dropdowns and
    Identify tool show real names."""
    try:
        for i, n in enumerate(band_names, start=1):
            layer.setBandName(i, n)  # QGIS >= 3.30
    except (AttributeError, TypeError, RuntimeError):  # Safe fallback
        pass


# ---------------------------------------------------------------- layer order
def raise_points_above(names=("sample_points",)):
    """Keep the sample points layer above imagery/index layers so the points
    are never hidden behind a raster."""
    proj = QgsProject.instance()
    root = proj.layerTreeRoot()
    for want in names:
        for lyr in proj.mapLayers().values():
            if lyr.name() != want:
                continue
            node = root.findLayer(lyr.id())
            if node is None:
                continue
            parent = node.parent() or root
            clone = node.clone()
            parent.insertChildNode(0, clone)
            parent.removeChildNode(node)


def show_only_layers(show_names=(), hide_prefixes=()):
    """Tick the given layers and untick any whose name starts with one of
    hide_prefixes. Used to give the labeling step a clean canvas."""
    from qgis.core import QgsLayerTreeLayer
    proj = QgsProject.instance()
    root = proj.layerTreeRoot()
    for node in root.findLayers():
        if not isinstance(node, QgsLayerTreeLayer):
            continue
        lyr = node.layer()
        if lyr is None:
            continue
        nm = lyr.name()
        if any(nm.startswith(p) for p in hide_prefixes):
            node.setItemVisibilityChecked(False)
        if nm in show_names:
            node.setItemVisibilityChecked(True)


def add_worldcover_layer():
    """Add ESA WorldCover 2021 as a GEE tile layer (once)."""
    for lyr in QgsProject.instance().mapLayers().values():
        if "worldcover" in lyr.name().lower():
            return lyr
    try:
        ee = init_ee()
        wc = ee.ImageCollection("ESA/WorldCover/v200").first().select("Map")
        vis = {"min": 10, "max": 100,
               "palette": ["006400", "ffbb22", "ffff4c", "f096ff",
                           "fa0000", "b4b4b4", "f0f0f0", "0064c8",
                           "0096a0", "00cf75", "fae6a0"]}
        return add_ee_image_layer(wc, vis, "ESA WorldCover 2021")
    except (ImportError, ModuleNotFoundError, RuntimeError):  # GEE/Import
        return None


# ---------------------------------------------------------------- unsupervised
UNSUPERVISED = ["K-Means (wekaKMeans)", "X-Means (wekaXMeans)",
                "Cascade K-Means (wekaCascadeKMeans)", "LVQ (wekaLVQ)"]


def make_clusterer(name, n_clusters):
    """Return an untrained GEE clusterer for the given algorithm."""
    import ee
    if name.startswith("X-Means"):
        return ee.Clusterer.wekaXMeans(minClusters=2, maxClusters=n_clusters)
    if name.startswith("Cascade"):
        return ee.Clusterer.wekaCascadeKMeans(
            minClusters=2, maxClusters=n_clusters)
    if name.startswith("LVQ"):
        return ee.Clusterer.wekaLVQ(numClusters=n_clusters)
    return ee.Clusterer.wekaKMeans(nClusters=n_clusters)


def unsupervised_classify(stack, aoi, n_clusters=8, scale=10,
                          n_samples=5000, algorithm="K-Means (wekaKMeans)"):
    """Cluster the stack into n_clusters groups (no training labels needed).
    Returns the clustered image (byte) with values 0..n_clusters-1."""
    training = stack.sample(region=aoi, scale=scale, numPixels=n_samples,
                            seed=42, geometries=False)
    clusterer = make_clusterer(algorithm, n_clusters).train(training)
    return stack.cluster(clusterer).clip(aoi).toByte()


def order_reference_layers(top=None, below_points=(), visible=(),
                           unchecked=(), hide_prefixes=()):
    """Arrange the labeling reference layers in the tree.

    - `top`: layer raised to the very top (e.g. sample_points)
    - `below_points`: layers moved to sit directly beneath `top`
    - `visible`: layer names ticked ON
    - `unchecked`: layer names ticked OFF (present but hidden)
    - `hide_prefixes`: any layer whose name starts with one of these is
      unticked (heavy layers we don't want on the canvas)
    """
    from qgis.core import QgsLayerTreeLayer
    proj = QgsProject.instance()
    root = proj.layerTreeRoot()

    def _find(name):
        for lyr in proj.mapLayers().values():
            if lyr.name() == name:
                node = root.findLayer(lyr.id())
                if node is not None:
                    return lyr, node
        return None, None

    # visibility first
    for node in root.findLayers():
        if not isinstance(node, QgsLayerTreeLayer):
            continue
        lyr = node.layer()
        if lyr is None:
            continue
        nm = lyr.name()
        if any(nm.startswith(p) for p in hide_prefixes):
            node.setItemVisibilityChecked(False)
        if nm in visible:
            node.setItemVisibilityChecked(True)
        if nm in unchecked:
            node.setItemVisibilityChecked(False)

    # move `top` to index 0
    order = []
    if top:
        order.append(top)
    order.extend(below_points)
    # re-insert in the desired order at the top of the root
    for i, name in enumerate(order):
        lyr, node = _find(name)
        if node is None:
            continue
        parent = node.parent() or root
        clone = node.clone()
        parent.insertChildNode(i, clone)
        parent.removeChildNode(node)


# ---------------------------------------------------------------- schemes
# ESA WorldCover 2021 classes with official codes and palette colours.
WORLDCOVER_CLASSES = [
    (10, "Tree cover",            "#006400"),
    (20, "Shrubland",             "#ffbb22"),
    (30, "Grassland",             "#ffff4c"),
    (40, "Cropland",              "#f096ff"),
    (50, "Built-up",              "#fa0000"),
    (60, "Bare / sparse veg.",    "#b4b4b4"),
    (70, "Snow and ice",          "#f0f0f0"),
    (80, "Permanent water",       "#0064c8"),
    (90, "Herbaceous wetland",    "#0096a0"),
    (95, "Mangroves",             "#00cf75"),
    (100, "Moss and lichen",      "#fae6a0"),
]

# A user-editable custom scheme. Persisted to a small JSON next to the
# plugin so it survives restarts. Starts empty; the user adds classes.
import json as _json

_CUSTOM_PATH = os.path.join(os.path.dirname(__file__), "custom_classes.json")


def load_custom_classes():
    """Return the user's custom class list [(code, name, hex), ...]."""
    try:
        with open(_CUSTOM_PATH, "r", encoding="utf-8") as fh:
            data = _json.load(fh)
        out = []
        for row in data:
            out.append((int(row["code"]), str(row["name"]),
                        str(row.get("color", "#888888"))))
        return out
    except (OSError, IOError):  # File operation
        return []


def save_custom_classes(classes):
    """Persist the custom class list. `classes` = [(code, name, hex), ...]."""
    try:
        data = [{"code": int(c), "name": n, "color": h}
                for c, n, h in classes]
        with open(_CUSTOM_PATH, "w", encoding="utf-8") as fh:
            _json.dump(data, fh, indent=2)
        return True
    except (OSError, IOError):  # File operation
        return False


# scheme name -> function returning the class list, so callers stay generic
def scheme_classes(scheme):
    if scheme == "worldcover":
        return list(WORLDCOVER_CLASSES)
    if scheme == "custom":
        return load_custom_classes()
    return list(LULC_CLASSES)          # default LULC L1


def scheme_palette(scheme):
    """min, max, palette list for a paletted raster of this scheme."""
    cls = scheme_classes(scheme)
    if not cls:
        return 1, 1, ["888888"]
    codes = [c for c, _, _ in cls]
    return min(codes), max(codes), [h.lstrip("#") for _, _, h in cls]
