RAVE-SCOUT/src/score.py

"""Combina las capas en dos rankings distintos y extrae candidatos.

1) AISLAMIENTO PURO: ¿qué punto de España está más lejos de cualquier rastro de
   civilización? Es la pregunta interesante, pero su respuesta suele ser una
   cara norte a 2.400 m sin forma de llegar.

2) APTITUD PARA RAVE: aislamiento acotado por lo que la fiesta necesita de
   verdad. Acceso rodado, suelo llano, fuera de espacio protegido y militar, y
   a ser posible en una hondonada, que encierra el sonido y tapa las luces.

Las dos restricciones tiran en sentidos opuestos: cuanto más aislado, peor
acceso. El sitio bueno no es el máximo de ninguna de las dos, sino el máximo de
una con la otra acotada, que es justo lo que hace el modelo.
"""
import pickle
import sys

import numpy as np
import rasterio
from shapely.geometry import Point
from shapely.strtree import STRtree

sys.path.insert(0, str(__file__.rsplit("/", 1)[0]))
import config as C

# --- parámetros del modelo (todo en metros) ---
MIN_D_BUILD = 1500     # umbral duro: por debajo, te oyen desde casa
MAX_D_ACCESS = 600     # más lejos de un vial, no metes el equipo
MIN_D_PROT = 500       # margen al borde del espacio protegido (zona periférica)
MAX_SLOPE = 8.0        # grados; por encima no se monta ni se aparca
SAT_BUILD = 5000       # a partir de aquí el aislamiento ya no suma más
SAT_PLACE = 8000
SAT_ROAD = 5000

N_CAND = 60
MIN_SEP = 15000        # separación mínima entre candidatos, para no devolver
                       # 60 celdas del mismo valle



def load(name):
    return np.load(C.INTERIM / f"{name}.npy")


def norm(a, sat):
    return np.clip(a / sat, 0, 1)


def main():
    print("cargando capas…", flush=True)
    d_build = load("d_build")
    d_major = load("d_major")
    d_minor = load("d_minor")
    d_track = load("d_track")
    d_place = load("d_place_any")
    dens = load("dens_build_5km")
    slope = load("slope_deg")
    tpi = load("tpi_2km")
    spain = load("mask_spain").astype(bool)
    protected = load("mask_protected").astype(bool)
    military = load("mask_military").astype(bool)
    water = load("mask_water").astype(bool)

    d_access = np.minimum(np.minimum(d_major, d_minor), d_track)

    # ---------- 1) aislamiento puro ----------
    # Media geométrica de las tres distancias que definen "lejos de todo".
    # Geométrica y no aritmética: así un sitio a 8 km de una casa pero a 200 m
    # de una nacional no puntúa; hay que estar lejos de las tres cosas a la vez.
    iso = np.cbrt(np.maximum(d_build, 1) *
                  np.maximum(np.minimum(d_major, d_minor), 1) *
                  np.maximum(d_place, 1)).astype(np.float32)
    iso[~spain] = np.nan
    np.save(C.INTERIM / "iso_pure.npy", iso)

    valid_iso = spain & np.isfinite(iso)
    print(f"\n=== AISLAMIENTO PURO ===")
    print(f"  celdas válidas: {valid_iso.sum():,}")
    for q in (50, 90, 99, 99.9):
        print(f"  p{q:<5} {np.nanpercentile(iso[valid_iso], q):7.0f} m")
    print(f"  máximo  {np.nanmax(iso):7.0f} m")
    print(f"  distancia máxima a un edificio en España: "
          f"{np.nanmax(np.where(spain, d_build, np.nan))/1000:.1f} km")
    print(f"  distancia máxima a una carretera:         "
          f"{np.nanmax(np.where(spain, np.minimum(d_major, d_minor), np.nan))/1000:.1f} km")

    # ---------- 2) aptitud para rave ----------
    d_prot = load("d_protected")
    d_mil = load("d_military")
    ok = (spain & ~protected & ~military & ~water &
          (d_build >= MIN_D_BUILD) &
          (d_access <= MAX_D_ACCESS) &
          (slope <= MAX_SLOPE) &
          (d_prot >= MIN_D_PROT) &
          np.isfinite(slope))
    print(f"\n=== APTITUD PARA RAVE ===")
    print(f"  celdas que pasan todas las restricciones: {ok.sum():,} "
          f"({ok.sum()/spain.sum()*100:.2f}% del país, "
          f"{ok.sum()*0.01:,.0f} km2)")

    score = (0.34 * norm(d_build, SAT_BUILD) +
             0.16 * norm(d_place, SAT_PLACE) +
             0.20 * norm(np.minimum(d_major, d_minor), SAT_ROAD) +
             0.12 * (1 - np.clip(dens / 200.0, 0, 1)) +
             0.10 * np.clip(-tpi / 30.0, 0, 1) +      # hondonada
             0.08 * (1 - np.clip(slope / MAX_SLOPE, 0, 1)))
    score = score.astype(np.float32)
    score[~ok] = np.nan
    np.save(C.INTERIM / "score_rave.npy", score)

    with rasterio.open(
            C.OUT / "score_rave.tif", "w", driver="GTiff",
            height=C.HEIGHT, width=C.WIDTH, count=1, dtype="float32",
            crs=C.CRS_GRID, transform=C.TRANSFORM, nodata=np.nan,
            compress="deflate", tiled=True) as dst:
        dst.write(score, 1)
    print(f"  ráster guardado en out/score_rave.tif")

    # ---------- extracción de candidatos ----------
    print(f"\nextrayendo {N_CAND} candidatos con {MIN_SEP/1000:.0f} km de separación…",
          flush=True)
    with open(C.INTERIM / "area_admin_municipality.pkl", "rb") as fh:
        munis = pickle.load(fh)
    with open(C.INTERIM / "area_admin_region.pkl", "rb") as fh:
        regs = [r for r in pickle.load(fh) if r["geom"].area > 1e6]
    muni_tree = STRtree([m["geom"] for m in munis])
    reg_tree = STRtree([r["geom"] for r in regs])

    work = np.where(np.isfinite(score), score, -1.0).astype(np.float32)
    rad = int(MIN_SEP / C.RES)
    rows = []
    for i in range(N_CAND):
        idx = int(np.argmax(work))
        r, c = divmod(idx, C.WIDTH)
        if work[r, c] <= 0:
            break
        x = C.X_MIN + (c + 0.5) * C.RES
        y = C.Y_MAX - (r + 0.5) * C.RES
        lon, lat = C.to_wgs(x, y)
        pt = Point(x, y)

        def lookup(tree, recs):
            for j in tree.query(pt):
                if recs[j]["geom"].contains(pt):
                    return recs[j]["name"]
            return ""

        rows.append({
            "rank": i + 1,
            "lat": round(float(lat), 5),
            "lon": round(float(lon), 5),
            "municipio": lookup(muni_tree, munis),
            "comunidad": lookup(reg_tree, regs),
            "score": round(float(score[r, c]), 4),
            "d_edificio_m": int(d_build[r, c]),
            "d_carretera_m": int(min(d_major[r, c], d_minor[r, c])),
            "d_acceso_m": int(d_access[r, c]),
            "d_nucleo_m": int(d_place[r, c]),
            "d_protegido_m": int(d_prot[r, c]),
            "d_militar_m": int(d_mil[r, c]),
            "cota_m": int(np.nan_to_num(load_dem_at(r, c))),
            "pendiente_deg": round(float(slope[r, c]), 1),
            "tpi_m": round(float(tpi[r, c]), 1),
            "edif_5km": int(dens[r, c]),
        })

        # suprimimos un disco alrededor para el siguiente candidato
        r0, r1 = max(0, r - rad), min(C.HEIGHT, r + rad + 1)
        c0, c1 = max(0, c - rad), min(C.WIDTH, c + rad + 1)
        yy, xx = np.ogrid[r0:r1, c0:c1]
        work[r0:r1, c0:c1][((yy - r) ** 2 + (xx - c) ** 2) <= rad * rad] = -1.0

    import csv
    cols = list(rows[0].keys())
    with open(C.OUT / "candidatos.csv", "w", newline="") as fh:
        w = csv.DictWriter(fh, fieldnames=cols)
        w.writeheader()
        w.writerows(rows)
    print(f"  {len(rows)} candidatos -> out/candidatos.csv")
    for r_ in rows[:15]:
        print(f"  #{r_['rank']:<3} {r_['lat']:8.4f},{r_['lon']:9.4f}  "
              f"{r_['municipio'][:26]:26s} {r_['comunidad'][:18]:18s} "
              f"casa={r_['d_edificio_m']:>6}m  ctra={r_['d_carretera_m']:>5}m  "
              f"acc={r_['d_acceso_m']:>4}m")


_dem = None


def load_dem_at(r, c):
    global _dem
    if _dem is None:
        _dem = np.load(C.INTERIM / "dem.npy", mmap_mode="r")
    v = _dem[r, c]
    return 0.0 if not np.isfinite(v) else float(v)


if __name__ == "__main__":
    main()