#!/usr/bin/env python3
"""
Site Grid — Faz 3b: SoilGrids GeoTIFF'lerinden hücre başına tek satır üret.

KALINLIK AĞIRLIKLI 0-30 cm BİRLEŞTİRME (direktif):
    v = (v0_5*5 + v5_15*10 + v15_30*15) / 30
BASİT ORTALAMA YAPILMAZ. Bu, `src/lib/ne-yetisir/adapters/soil.ts` içindeki
`normalizeSoil()` fonksiyonunun yaptığı işin AYNISIDIR; davranış korunur.

DOKU SINIFI BURADA HESAPLANMAZ.
Sebep: `usdaTextureClass()` ve `faoTextureCategory()` TypeScript tarafında
zaten var ve skorlamayı besliyor. Aynı algoritmayı Python'da ikinci kez yazmak
iki uygulamanın zamanla ayrışması demektir. Bu yüzden ızgara YALNIZ sayısal
değerleri (ph, sand, silt, clay, som_pct) taşır; sınıflandırmayı mevcut TS
fonksiyonları yapar. Tek algoritma, tek doğruluk kaynağı.

BİRİMLER (SoilGrids "mapped units" -> hedef birim, d_factor = 10):
    phh2o  68  -> 6.8  pH
    sand  400  -> 40.0 %
    silt  350  -> 35.0 %
    clay  250  -> 25.0 %
    soc    12  ->  1.2 g/kg  -> som_pct = soc_g_kg * 1.724 / 10
NODATA = 0 (doğrulandı: deniz bbox'ında %70,9 sıfır; pH 0 imkânsız).

ÇÖZÜNÜRLÜK KARARI:
SoilGrids yerel çözünürlüğü 250 m (~0,0023°). Türkiye'yi bu çözünürlükte
satırlamak ~13 milyon satır demek. Izgara 0,01° (~1,1 km) tutuluyor ve her
hücrede içine düşen yerel pikseller ORTALANIYOR. Kaybın büyüklüğü tahmin
edilmiyor, ÖLÇÜLÜYOR: 05_compare_soil.py mevcut REST önbelleğiyle karşılaştırır.
"""

from __future__ import annotations

import argparse
import sys
from pathlib import Path

from config import OUT_DIR, PILOT_BBOX, RAW_DIR, TURKEY_BBOX, BBox

try:
    import numpy as np
    import pandas as pd
    import tifffile
except ImportError:
    print("HATA: numpy/pandas/tifffile gerekli. `pip install -r requirements.txt tifffile`", file=sys.stderr)
    raise SystemExit(4) from None

SOIL_DIR = RAW_DIR / "soil"
PROPERTIES = ["phh2o", "sand", "silt", "clay", "soc"]
# (etiket, kalınlık cm)
DEPTHS = [("0-5cm", 5), ("5-15cm", 10), ("15-30cm", 15)]
TOTAL_THICKNESS = 30
D_FACTOR = 10.0
NODATA = 0

SOIL_GRID_STEP = 0.01  # ~1,1 km


def read_layer(prop: str, depth: str, bbox: BBox) -> tuple[np.ndarray, tuple[float, float, float, float]]:
    """GeoTIFF'i oku; (dizi, (lon0, lat0, dlon, dlat)) döndür."""
    path = SOIL_DIR / f"{prop}_{depth}_{bbox.name}.tif"
    if not path.exists():
        raise SystemExit(f"HATA: {path} yok. Önce 03_download_soil.py çalıştır.")
    with tifffile.TiffFile(path) as tf:
        page = tf.pages[0]
        arr = page.asarray().astype(np.float32)
        tags = {t.name: t.value for t in page.tags}
    scale = tags["ModelPixelScaleTag"]
    tie = tags["ModelTiepointTag"]
    # tie = (i, j, k, x, y, z) -> sol üst köşenin lon/lat'ı
    lon0, lat0 = float(tie[3]), float(tie[4])
    return arr, (lon0, lat0, float(scale[0]), float(scale[1]))


def weighted_depth_mean(bbox: BBox, prop: str) -> tuple[np.ndarray, tuple[float, float, float, float]]:
    """
    Üç derinliği KALINLIK AĞIRLIKLI birleştirir.
    Bir derinlik nodata ise o piksel için ağırlığı düşülür (kısmi profil kabul
    edilir); hiçbiri yoksa NaN kalır. soil.ts'in davranışıyla aynı.
    """
    total = None
    weight = None
    geo = None
    for label, thick in DEPTHS:
        arr, g = read_layer(prop, label, bbox)
        geo = geo or g
        if arr.shape != (total.shape if total is not None else arr.shape):
            pass
        valid = arr != NODATA
        contrib = np.where(valid, arr * thick, 0.0)
        w = np.where(valid, float(thick), 0.0)
        total = contrib if total is None else total + contrib
        weight = w if weight is None else weight + w
    with np.errstate(invalid="ignore", divide="ignore"):
        out = np.where(weight > 0, total / weight, np.nan) / D_FACTOR
    return out, geo  # type: ignore[return-value]


def build(bbox: BBox, band_rows: int = 400) -> pd.DataFrame:
    """
    Izgarayı BANT BANT kurar (bellek).

    NEDEN: pilot (2,8 milyon piksel) tek seferde 796 MB peak RAM kullandı.
    Türkiye ~26,7 milyon piksel = 9,5 kat; tek seferde düzleştirmek birkaç GB
    ister ve makineyi zorlar. Enlem bantlarına bölünce her bandın DataFrame'i
    1/N boyutunda kalıyor, sonuç aynı: 0,01° hücreler bant sınırında bölünmez
    çünkü bant yüksekliği (400 piksel ~0,9°) hücre adımından çok büyük ve
    gruplama en sonda bir kez daha birleştiriliyor.
    """
    layers: dict[str, np.ndarray] = {}
    geo = None
    for prop in PROPERTIES:
        arr, g = weighted_depth_mean(bbox, prop)
        layers[prop] = arr
        geo = geo or g
        print(f"  {prop}: {arr.shape}  geçerli %{100 * np.isfinite(arr).mean():.1f}")

    lon0, lat0, dlon, dlat = geo  # type: ignore[misc]
    ny, nx = layers["phh2o"].shape

    lons = lon0 + (np.arange(nx) + 0.5) * dlon
    lats = lat0 - (np.arange(ny) + 0.5) * dlat  # kuzeyden güneye
    cell_lon = np.round(lons / SOIL_GRID_STEP).astype(np.int64)
    cell_lat_all = np.round(lats / SOIL_GRID_STEP).astype(np.int64)

    parts: list[pd.DataFrame] = []
    n_bands = int(np.ceil(ny / band_rows))
    for b in range(n_bands):
        y0, y1 = b * band_rows, min((b + 1) * band_rows, ny)
        h = y1 - y0
        df = pd.DataFrame(
            {
                "clat": np.repeat(cell_lat_all[y0:y1], nx),
                "clon": np.tile(cell_lon, h),
                **{p: layers[p][y0:y1].ravel() for p in PROPERTIES},
            }
        )
        df = df.dropna(subset=["phh2o", "sand"], how="all")
        if len(df):
            # Hücre içi ortalamayı iki adımda al: bant içi toplam+adet, sonda birleştir.
            g = df.groupby(["clat", "clon"]).agg(["sum", "count"])
            g.columns = ["_".join(c) for c in g.columns]
            parts.append(g.reset_index())
        print(f"  bant {b + 1}/{n_bands}  satır {y0}-{y1}", end="\r")
    print()

    if not parts:
        return pd.DataFrame()

    allp = pd.concat(parts, ignore_index=True)
    key = ["clat", "clon"]
    sums = allp.groupby(key, as_index=False).sum(numeric_only=True)

    agg = pd.DataFrame({"clat": sums["clat"], "clon": sums["clon"]})
    for prop in PROPERTIES:
        # Bant sınırında bölünen hücreler burada doğru şekilde birleşiyor:
        # toplamların toplamı / adetlerin toplamı = gerçek ortalama.
        agg[prop] = sums[f"{prop}_sum"] / sums[f"{prop}_count"].replace(0, np.nan)

    out = pd.DataFrame(
        {
            # soil_key: 0,01° ızgara anahtarı. store.ts'teki 0,002°'lik REST
            # anahtarından FARKLIDIR; TS sağlayıcısı bu adımı bilerek kullanır.
            "soil_key": (agg["clat"] * SOIL_GRID_STEP).round(4).map("{:.4f}".format)
            + ":"
            + (agg["clon"] * SOIL_GRID_STEP).round(4).map("{:.4f}".format),
            "lat": (agg["clat"] * SOIL_GRID_STEP).round(4),
            "lng": (agg["clon"] * SOIL_GRID_STEP).round(4),
            "ph": agg["phh2o"].round(1),
            "sand": agg["sand"].round(1),
            "silt": agg["silt"].round(1),
            "clay": agg["clay"].round(1),
            # som_pct = soc(g/kg) * 1,724 / 10   (soil.ts ile birebir)
            "som_pct": (agg["soc"] * 1.724 / 10).round(2),
        }
    )
    # soil.ts kuralı: ph ve sand'in İKİSİ de yoksa profil YOK sayılır.
    return out[~(out["ph"].isna() & out["sand"].isna())].reset_index(drop=True)


def main() -> int:
    ap = argparse.ArgumentParser()
    ap.add_argument("--full", action="store_true")
    args = ap.parse_args()
    bbox = TURKEY_BBOX if args.full else PILOT_BBOX

    print(f"Kapsam: {bbox.name} · ızgara {SOIL_GRID_STEP}°")
    df = build(bbox)

    OUT_DIR.mkdir(parents=True, exist_ok=True)
    out = OUT_DIR / f"soil_grid_{bbox.name}.csv"
    df.to_csv(out, index=False)

    # Kapsama: bbox'taki teorik hücre sayısına göre doluluk.
    theo = int(round((bbox.north - bbox.south) / SOIL_GRID_STEP)) * int(
        round((bbox.east - bbox.west) / SOIL_GRID_STEP)
    )
    print(f"\nSatır      : {len(df):,}")
    print(f"Teorik hücre: {theo:,}")
    print(f"Kapsama    : %{100 * len(df) / theo:.1f}  (eksik %{100 * (1 - len(df) / theo):.1f})")
    print(f"pH dolu    : %{100 * df['ph'].notna().mean():.1f}")
    print(f"Doku dolu  : %{100 * df['sand'].notna().mean():.1f}")
    print(f"Dosya      : {out}  ({out.stat().st_size / 1024 / 1024:.1f} MB)")
    return 0


if __name__ == "__main__":
    raise SystemExit(main())