Metadata-Version: 2.4
Name: cfasig
Version: 0.7.4
Summary: Wrapper simple en español sobre geopandas/shapely/rasterio para tareas de SIG en CONSAEFA.
Author: Antruc
License-Expression: MIT
Project-URL: Homepage, https://github.com/CONSAEFA/cfasig
Project-URL: Repository, https://github.com/CONSAEFA/cfasig
Project-URL: Issues, https://github.com/CONSAEFA/cfasig/issues
Keywords: gis,sig,geopandas,shapely,rasterio,forestal
Classifier: Development Status :: 4 - Beta
Classifier: Intended Audience :: Science/Research
Classifier: Topic :: Scientific/Engineering :: GIS
Classifier: Programming Language :: Python :: 3
Classifier: Programming Language :: Python :: 3.10
Classifier: Programming Language :: Python :: 3.14
Classifier: Operating System :: OS Independent
Classifier: Natural Language :: Spanish
Requires-Python: >=3.10
Description-Content-Type: text/markdown
License-File: LICENSE
Requires-Dist: geopandas>=0.13
Requires-Dist: shapely>=2.0
Requires-Dist: pandas>=1.5
Requires-Dist: numpy>=1.24
Requires-Dist: rasterio>=1.3
Requires-Dist: scipy>=1.10
Requires-Dist: pysheds>=0.3
Requires-Dist: openpyxl>=3.1
Provides-Extra: dev
Requires-Dist: pytest>=7; extra == "dev"
Requires-Dist: ruff>=0.1; extra == "dev"
Dynamic: license-file

# cfasig: Librería SIG, CONSAEFA S.C.

[![Code style: black](https://img.shields.io/badge/code%20style-black-000000.svg)](https://github.com/psf/black)

Wrapper en español sobre geopandas/shapely (vector) y rasterio/pysheds/scipy (raster) que simplifica procesos SIG (cargar de shapefile o de puntos CSV/Excel, recortar, disolver, buffer, combinar, convertir entre puntos/líneas/polígonos, guardar; remuestrear, hidrología, aspecto, pendiente, reclasificación por rangos). La idea es tener una sola función clara por operación en lugar de repetir la sintaxis de cada librería en cada script.

**Versión:** 0.7.4 | **Fecha:** Julio 2026
**Paquete:** `cfasig` (Python ≥ 3.10; probado en 3.14)
**Ruta local:** `%USERPROFILE%\Downloads\cfasig\`

## Stack

- **Base vector:** geopandas ≥ 0.13 (GeoDataFrame como estructura central).
- **Geometría:** shapely ≥ 2.0 (`unary_union`, `box`, `buffer`, `difference`, `split`).
- **Tablas:** pandas ≥ 1.5 (concatenación de capas), openpyxl ≥ 3.1 (lectura de .xlsx en `cargar_puntos`, escritura en `exportar_tabla`).
- **Raster:** rasterio ≥ 1.3 (recorte, remuestreo, rasterización, poligonización), numpy ≥ 1.24 (arrays), scipy ≥ 1.10 (`ndimage`: etiquetado, relleno, suavizado; `interpolate`/`spatial`: interpolación de MDE por TIN, spline de placa delgada e IDW), pysheds ≥ 0.3 (dirección/acumulación de flujo).
- **Empaquetado:** setuptools + pyproject.toml. Instalable en modo editable con `pip install -e .`. Licencia MIT. Expone el comando de consola `cfasig` (ver sección CLI). Dev: pytest ≥ 7, ruff ≥ 0.1.
- **Idioma:** API, funciones y docstrings en español; distancias en las unidades del CRS (metros si es UTM).
- **Filosofía:** wrapper delgado. Cada función es una operación conocida de las librerías base con nombre simple, valores por defecto sensatos y mensajes de aviso en consola. No reinventa; ordena. No está limitado a geopandas/shapely: cubre SIG en general (vector y raster).

## Estructura de archivos

Layout `src/`: el paquete vive bajo `src/cfasig/`, así que hay que instalarlo (`pip install -e .`) para importarlo; evita que los tests importen el código desde la carpeta en vez del instalado.

- **Raíz:** pyproject.toml / README.md / LICENSE / CLAUDE.md. Además, dos scripts de trabajo real (no son parte del paquete, viven aquí porque son del predio en curso y traen rutas absolutas): `predio_completo.py`, las cinco etapas encadenadas (caminos → hidrología → MDE → pendientes → clasificación), y `prueba_mde.py`, el barrido que compara métodos e intervalos de muestreo del MDE por superficie en `0-5`, `>100` y sobrepaso de cotas.
- **ejemplos/**: scripts de referencia: `caminos.py` (buffer jerárquico), `hidrologia.py` (buffer por condición), `mde.py` (MDE desde curvas, quemado y acondicionado). Rutas de prueba relativas (`datos/` de entrada, `salidas/` de salida); `guardar`/`guardar_raster` crean la carpeta de salida solos.
- **src/cfasig/**
  - `__init__.py`: expone toda la API pública (`import cfasig as sig`)
  - `archivo.py`: entrada/salida de capas
  - `proyeccion.py`: sistemas de coordenadas (CRS)
  - `geometria.py`: operaciones geométricas (vector) y utilidades de columnas
  - `raster.py`: operaciones raster (recorte, remuestreo, hidrología, aspecto, pendiente, reclasificación)
  - `cli.py`: comando de consola `cfasig` (convierte entre formatos)
  - `py.typed`: marca el paquete como tipado (PEP 561); las firmas llevan anotaciones y el editor las aprovecha sin stubs aparte
- **tests/**: suite pytest (`conftest.py` con fixtures + un `test_*.py` por módulo)

## CRS: quién alinea solo

`recortar`, `intersectar`, `quitar_solape` y `unir_atributos` reproyectan el **segundo** argumento al CRS del primero por su cuenta. La capa que procesas nunca se mueve; se mueve la máscara. **No escribas `reproyectar` antes de llamarlas**: la línea sobra y es donde se cuelan los errores.

De ahí el orden barato: recorta primero en el CRS de origen y reproyecta el recorte después, no al revés. Mueves menos geometrías.

`combinar` no alinea a propósito (concatenar dos sistemas da una capa mal ubicada en silencio): ahí sí, `asegurar_crs` antes.

Nunca pases `capa.crs.to_epsg()` como destino. Devuelve `None` cuando el CRS no tiene código EPSG registrado, que es el caso de la cartografía INEGI (Cónica Conforme de Lambert ITRF92, viene como WKT suelto). Pasa `capa.crs` directo.

```python
cauces = sig.recortar(cauces, cuadro)                    # sí
cauces = sig.recortar(cauces, sig.reproyectar(cuadro, cauces.crs.to_epsg()))  # no
```

## API pública (`import cfasig as sig`)

**archivo.py**

- `cargar(ruta, capa=None, mostrar=True)`: lee una capa vectorial (.shp/.gpkg/.geojson/.gpx) como GeoDataFrame; `capa` elige la capa en formatos multicapa (GPX: waypoints/routes/tracks; GeoPackage). Opcionalmente imprime nº de entidades y CRS.
- `cargar_puntos(ruta, x="x", y="y", epsg=None, orden=None, mostrar=True)`: arma una capa de puntos desde un CSV o Excel (.csv/.xlsx/.xls) con columnas de coordenadas. `x`/`y` nombran las columnas de coordenada; `epsg` fija el CRS (el archivo de texto no lo trae, sin él no se puede reproyectar ni medir áreas; avisa si falta); `orden` ordena las filas por una columna de secuencia antes de armar la capa (importa si luego conviertes a línea/polígono). Todas las columnas del archivo quedan como atributos.
- `convertir(entrada, salida, capa=None, gpx_como="track", mostrar=True)`: convierte de un formato a otro (`cargar` + `guardar`) en una línea; útil para bucles de lote. Si la entrada es GPX y no se da `capa`, autodetecta la capa con datos (un GPX de solo tracks leído a secas sale vacío porque la capa por defecto es waypoints).
- `guardar(gdf, ruta, gpx_como="track", mostrar=True)`: escribe la capa; el formato se deduce de la extensión. Para .gpx/.kml/.kmz reproyecta automáticamente a EPSG:4326 (esos formatos solo aceptan lon/lat) y usa el driver adecuado (GPX / LIBKML). Convertir = `cargar` + `guardar`: p. ej. `guardar(cargar("predio.shp"), "predio.kmz")`. GPX no admite polígonos: se exporta su contorno como línea. `gpx_como` decide cómo se escriben las líneas en GPX: `"track"` (por defecto, como trabaja el GPS aquí) o `"route"`. Antes de escribir avisa si una columna de área (`SUP`, `superficie_ha`...) ya no cuadra con la geometría (pista de que se calculó antes de recortar/reproyectar y quedó obsoleta).
- `exportar_tabla(gdf, ruta, mostrar=True)`: escribe solo la tabla de atributos (sin geometría) a `.csv` o `.xlsx`. Para llevar los datos a Excel sin pasar por un SIG. Valida la extensión antes de escribir y avisa de columnas de área obsoletas igual que `guardar`. `# ponytail:` función aparte en vez de aceptar `.csv` en `guardar` (que escribe capas, no tablas); sin opciones de separador, encoding ni selección de columnas: filtra con `gdf[cols]` antes de llamar.
- `cargar_waypoints(ruta, mostrar=True)`: lee los waypoints de un GPX con la tabla de atributos limpia (columnas NAME, LAYER, ELEVATION, time en ISO UTC, sym), descartando las ~18 columnas vacías del esquema fijo de GPX. Estilo Global Mapper.
- `cargar_tracks(ruta, utc_offset=-6, mostrar=True)`: lee los tracks de un GPX, una línea por tramo (`<trkseg>`), con columnas limpias (NAME, LAYER, gpxx_DisplayColor, START_TIME, END_TIME). Los tiempos salen del primer/último punto del tramo convertidos a hora local (`utc_offset`, por defecto -6 = Jalisco) con formato español. `convertir` usa estos dos lectores automáticamente para GPX (waypoints/tracks); las routes u otras capas caen al lector genérico.

**proyeccion.py**

- `reproyectar(gdf, destino)`: reproyecta siempre al CRS indicado (`to_crs`). `destino` acepta EPSG (int), cadena (`'EPSG:32613'`, WKT, proj4) u objeto CRS. Para alinear con una capa cuyo CRS no tiene código EPSG registrado (cartografía INEGI en Cónica Conforme de Lambert ITRF92, que viene como WKT suelto), pasa `otra.crs` en vez de `otra.crs.to_epsg()`: eso último da `None` y truena dentro de geopandas. Falla con `ValueError` claro si el destino es `None`.
- `asegurar_crs(gdf, destino, nombre="capa")`: reproyecta solo si la capa no está ya en ese CRS; avisa cuando lo hace. Evita reproyecciones innecesarias. Mismo `destino` que `reproyectar`.

**geometria.py**

- `recortar(gdf, mascara)`: recorta (`clip`) por otra capa (o por una geometría shapely suelta); reproyecta la máscara al CRS de `gdf` si difieren. Conserva solo el tipo de geometría de entrada (`keep_geom_type`): descarta las esquirlas punto/línea que el clip genera al tocar un vértice, así el resultado guarda a shapefile sin error.
- `disolver(gdf, campo=None, valor=None)`: une todas las geometrías en una sola (`unary_union`); opcionalmente etiqueta el resultado con `campo=valor`. La capa que devuelve trae solo `geometry` (más `campo` si se pasa): los demás atributos se pierden, porque de N entidades sale una y no hay a quién asignárselos. Si los necesitas, usa `disolver_por_grupo`, que sí los conserva (primer valor de cada grupo). De una capa vacía devuelve igual una fila, con la geometría vacía, y avisa por pantalla: esa fila parece válida, viaja aguas abajo y no revienta hasta el resumen final, o nunca.
- `disolver_por_grupo(gdf, campo)`: un polígono por cada valor único de `campo` (dissolve por grupo, equivalente al "Disolver por gridcode" de ArcMap); a diferencia de `disolver`, `campo` agrupa, no solo etiqueta. Conserva el CRS.
- `buffer(gdf, distancia)`: área de influencia por geometría (unidades del CRS); devuelve copia, no muta el original.
- `quitar_solape(gdf, otro, limpiar=True)`: resta de `gdf` lo que pise `otro` (capa o geometría shapely suelta) con `difference`, dando prioridad a `otro`; limpia vacías por defecto. Reproyecta `otro` al CRS de `gdf` si difieren, igual que `recortar`. Si el recorte deja una mezcla de área y línea (bordes que se rozan), conserva solo la de mayor dimensión: esas esquirlas no son solape y revientan al guardar a shapefile.
- `resolver_por_prioridad(capas)`: resuelve los solapes de una lista ordenada de mayor a menor prioridad y devuelve una sola capa sin encimados. Cada capa se recorta contra la unión de **todas** las anteriores, no solo contra la inmediata superior: con 3+ capas, recortar solo contra la vecina deja pasar el solape de la tercera con la primera allí donde la segunda no cubre.
- `resolver_por_cercania(gdf_a, gdf_b, resolucion=1.0, nombres=("A","B"), max_pixeles=20_000_000)`: para dos capas que no deberían solaparse y no tienen jerarquía entre sí, parte la zona en conflicto y le da cada pedazo a la capa cuyo núcleo (su parte no conflictiva) esté más cerca. Es **Asignación Euclidiana** (*Euclidean Allocation*, Spatial Analyst de ArcGIS) con dos matices: las fuentes son solo dos (los dos núcleos) y la asignación se aplica solo dentro del solape, no a todo el ráster; el resultado equivale a partir la zona en conflicto con polígonos de Thiessen respecto a los dos núcleos, y `resolucion` es la malla temporal que aproxima ese borde. Ambas capas deben traer una sola entidad. Devuelve `(gdf_a, gdf_b, conflicto)`; `conflicto` trae la columna `gana` para revisar el reparto a ojo. En empate gana `gdf_a`. Si la malla pasa de `max_pixeles` falla pidiendo subir `resolucion`: ojo, el tope mide el **rectángulo envolvente** del conflicto, no su área, así que unas pocas esquirlas en esquinas opuestas del predio dan un rectángulo enorme con casi nada dentro. Ante ese error, sube `resolucion` a 2 o 5 (el borde entre capas queda definido a esos metros) antes de sospechar del solape.
- `combinar(capas)`: concatena una lista de capas en una sola (toma el CRS de la primera). Falla con `ValueError` si alguna capa viene en otro CRS: no reproyecta sola, porque concatenar coordenadas de dos sistemas da una capa silenciosamente mal ubicada. Reproyecta antes con `asegurar_crs`.
- `puntos_a_linea(gdf, campo=None)`: une los puntos en una línea siguiendo el orden de las filas; con `campo` genera una línea por cada valor distinto (agrupación). Requiere ≥2 puntos por línea.
- `puntos_a_poligono(gdf, campo=None, envolvente=False)`: une los puntos en un polígono usándolos como vértices en orden de filas; con `envolvente=True` usa la envolvente convexa (útil si no vienen ordenados por el contorno); con `campo`, un polígono por grupo. Requiere ≥3 puntos.
- `linea_a_poligono(gdf)`: cierra cada línea en un polígono (sus vértices como contorno); si la línea está abierta, une el último punto con el primero. Conserva atributos.
- `crear_cuadro(gdf, margen=0)`: rectángulo (bounding box) alrededor de la capa, ampliado `margen`; útil como máscara de recorte previo rápido.
- `cuadros_arcmap(predio, resolucion, margen_borde_celdas=10, margen_salida=0.0)`: genera `(cuadro_recorte, cuadro_salida)` para el patrón de extent de Topo to Raster (ArcMap). `cuadro_salida` = predio (ampliado `margen_salida`), el extent que conservas; `cuadro_recorte` = `cuadro_salida` ampliado `margen_borde_celdas × resolucion`, al que se recortan curvas y cauces. Así la entrada sobresale N celdas respecto a la salida y las celdas de borde se interpolan con datos a ambos lados (el extent de salida queda metido dentro de los datos, como recomienda ArcMap). El borde escala con la resolución: 200 m a 20 m/px con 10 celdas.
- `limpiar_vacias(gdf)`: elimina entidades con geometría vacía o nula.
- `reparar_geometrias(gdf)`: `buffer(0)` + descarta inválidas/vacías; el patrón de reparación que se repite tras cada operación pesada.
- `calcular_superficie(gdf, campo="superficie_ha", en_hectareas=True, decimales=2)`: añade una columna de área (m² del CRS, o ha si `en_hectareas`). Redondea a **2 decimales** por defecto (0.01 ha = 100 m², suficiente para uso forestal/catastral y sin el falso detalle de 4); súbelo con `decimales=` si necesitas más. El valor es una instantánea de la geometría actual: calcúlalo como último paso, después de reproyectar y recortar, o quedará obsoleto (`guardar` avisa si no cuadra, con holgura de 0.1% + 0.005 para no dar falsas alarmas por ese redondeo).
- `intersectar(gdf, otro)`: intersección con `overlay`; corta `gdf` al solape con `otro` generando una parte por cada par de entidades que se pisan (distinto de `recortar`, que solo recorta). Conserva los atributos de `gdf`; de `otro` solo usa la forma, no trae sus columnas (para eso, `unir_atributos`). Conserva solo el tipo de geometría de `gdf` (`keep_geom_type=True`): dos polígonos que comparten borde intersecan también en líneas y puntos, y esas esquirlas de contacto no son solape. Se pasa explícito para que geopandas no emita un `UserWarning` por cada capa con bordes compartidos.
- `unir_atributos(gdf, otro, como="left", predicado="intersects")`: spatial join: pega las columnas de `otro` a `gdf` según su posición sin cortar geometrías (a diferencia de `intersectar`). `predicado` = 'intersects'/'within'/'contains'... Reproyecta `otro` si difiere el CRS.
- `calcular_longitud(gdf, campo=None, en_km=False, decimales=2)`: añade una columna con la longitud de cada geometría (unidades del CRS, o km si `en_km`); sin `campo` el nombre es `longitud_m`, o `longitud_km` con `en_km`. Gemelo de `calcular_superficie` para líneas, mismo redondeo a 2 decimales.
- `simplificar(gdf, tolerancia, conservar_topologia=True)`: reduce vértices (Douglas-Peucker); `tolerancia` en unidades del CRS. Con `conservar_topologia` (por defecto) evita auto-cruces y no rompe bordes compartidos entre polígonos vecinos.
- `cerrar_microhuecos(gdf, distancia, estilo_junta=2)`: closing morfológico (expandir/contraer) que cierra huecos menores a `distancia`.
- `subdividir_por_area(gdf, area_max_ha, campo_area="superficie_ha", min_esquirla=100, max_prof=8)`: parte los polígonos que superen `area_max_ha` por bisección recursiva del eje más largo.
- `fusionar_menores(gdf, area_min_ha, area_max_ha, campo_area="superficie_ha", max_iter=30)`: fusiona iterativamente los polígonos bajo el mínimo con su mejor vecino, sin superar el máximo.
- `campo_seguro(gdf, campo)`: devuelve un nombre de columna que no choque con los existentes (si existe, genera una variante única `campo_ab12`).

**raster.py**

- `cargar_raster(ruta, mostrar=True)`: abre un raster de una banda; el nodata declarado en el archivo (p. ej. -9999) entra como NaN, que es como el módulo marca NoData. Devuelve `(array, perfil, transform, crs, resolucion_m)`.
- `guardar_raster(arr, ruta, perfil, mostrar=True)`: escribe un array 2D a raster (float32). `arr` puede ser cualquier cosa convertible a array 2D (p. ej. el `Raster` de pysheds que devuelve `acondicionar_mde`): la conversión va dentro, el script no necesita `numpy`.
- `perfil_raster(transform, crs, nodata=nan, driver="GTiff")`: arma el perfil (metadatos) que pide `guardar_raster`, en vez de escribir el dict a mano en cada script. `height`/`width` no van aquí: `guardar_raster` los toma del array. `nodata` es NaN por defecto (como marca NoData el módulo), así el .tif sale marcado. Para reusar el perfil de un archivo existente no hace falta: `cargar_raster` ya lo devuelve.
- `recortar_raster(ruta, mascara, nodata=nan)`: recorta un raster al contorno de una capa vector (reproyecta la máscara); el nodata declarado del archivo también entra como `nodata`, igual que `cargar_raster`. Devuelve `(array, perfil, transform, crs, resolucion_m)`.
- `rellenar_nodata(arr)`: rellena NaN con el valor del píxel válido más cercano (sin interpolar).
- `remuestrear(arr, transform, crs, res_destino, metodo=bilinear)`: cambia la resolución; devuelve `(array, transform)`.
- `interpolar_mde(curvas, campo_elevacion, resolucion=None, equidistancia=None, intervalo_muestreo=None, cuadro=None, metodo="spline", suavizado=0.0, vecinos=None, mostrar=True)`: genera un MDE continuo interpolando curvas de nivel vectoriales. Muestrea puntos a intervalos regulares sobre cada curva (evita el sesgo de densidad de vértices) e interpola sobre una malla regular: `"spline"`/`"tps"` (spline de placa delgada local, superficie suave que evita el escalonado de los triángulos planos entre curvas; deduplica los puntos XY coincidentes, p. ej. el inicio == fin de una curva cerrada, que dejarían singular el sistema RBF), `"tin"` (triangulación de Delaunay, respeta quiebres de pendiente) o `"idw"` (distancia inversa, respaldo si el TIN falla por geometría degenerada). El TIN produce pendientes falsas: un triángulo con sus tres vértices sobre la misma curva sale plano, y eso ocurre justo en cimas, lomos y fondos de vaguada. Medido sobre un predio de 4800 ha en sierra (equidistancia 20 m, píxel 10 m), el TIN reportaba 736 ha con menos de 5% de pendiente contra 123-164 ha del spline: ~600 ha de plano inventado. Por eso el default es `"spline"`; `"tin"` sigue siendo lo correcto si la cartografía trae líneas de quiebre y quieres que se respeten tal cual. El precio del spline es ondulación entre curvas, que infla el rango >100%; se controla con tres perillas: `intervalo_muestreo` (indirecta: subirlo hasta la separación horizontal media entre curvas bajó el >100% de 367 a 252 ha), y `suavizado`/`vecinos` (directas: `suavizado=0` interpola exacto y >0 relaja el ajuste; no está en metros, barre por órdenes de magnitud. `vecinos=None` = el default de cada método, 48 en spline y 8 en IDW, que no comparten escala; menos vecinos = menos rizo y menos margen de datos necesita el cuadro de recorte. `"tin"` ignora ambos y avisa). No calibres `suavizado` con `validar_mde`: suavizar es dejar de pasar por los puntos que ese hold-out reserva, así que su RMS empeora siempre y optimizarlo te devolvería 0 cada vez. Falta una cuarta que no es un parámetro de esta función: **el tamaño de píxel**. El rizo tiene la longitud de onda de la separación entre curvas (λ) y una diferencia finita de paso `h` lo convierte en pendiente espuria `A·sin(2πh/λ)/h`, máxima cuando `h ≈ λ/4`. En el mismo predio (λ=37.8 m), `resolucion=10` era el peor valor posible (~57% de pendiente inventada, sobre una pendiente real media de 53%) y `resolucion=20` lo canceló casi entero (~5%), bajando el >100% de 252 a 121 ha: la diferencia abarca una onda completa. El default (mitad de la separación) cae del lado bueno; si fijas `resolucion` a mano, haz la cuenta. Fuera de la envolvente convexa devuelve NaN (no extrapola; rellena luego con `rellenar_nodata`). Descarta curvas con cota NaN. `resolucion` e `intervalo_muestreo` por defecto salen de la **separación horizontal media entre curvas** (`area / largo_total`, redondeada a un valor limpio): `resolucion` = la mitad de esa separación, `intervalo_muestreo` = la separación. La equidistancia no sirve para esto: es una distancia vertical, 20 m de salto de cota no dicen cuántos metros hay de una curva a otra en el suelo. Muestrear más fino que la separación deja la nube de puntos más densa a lo largo de la curva que entre curvas, y esa anisotropía es lo que hace ondular al spline. `equidistancia` se autodetecta (moda de las diferencias entre cotas) y sirve para el reporte y para avisar si la que pasas no coincide. La línea de `mostrar` incluye `separacion:`, para poder juzgar los parámetros que fijaste a mano. Devuelve `(array, transform)`; encaja entre el flujo vector (`recortar`/`disolver`) y `acondicionar_mde`/`quemar_cauces`.
- `validar_mde(curvas, campo_elevacion, resolucion=None, equidistancia=None, intervalo_muestreo=None, cuadro=None, metodo="spline", suavizado=0.0, vecinos=None, porcentaje_reserva=0.2, semilla=None, mostrar=True)`: control de calidad de `interpolar_mde` por reserva de puntos (*hold-out*, como el RMS de ANUDEM/Topo to Raster): muestrea las curvas igual que `interpolar_mde`, aparta al azar `porcentaje_reserva` (20% por defecto) de los puntos, interpola con el resto (mismos `metodo`/parámetros; pásale también los mismos `suavizado`/`vecinos` o validarás otra superficie) y compara la elevación real contra la del MDE en los puntos reservados. `semilla` fija el reparto para reproducibilidad. Dos cosas que este número NO puede hacer: tunear `suavizado` (suavizar es dejar de pasar por los puntos que reserva, así que empeora siempre) y comparar resoluciones, porque muestrea el **centro de celda** más cercano al punto y el residual sale con dos componentes sumadas, el error del interpolador y el de discretizar la malla, y la segunda crece con el píxel. Medido en un predio de sierra, pasar de 10 a 20 m/px subió el RMS de 7.16 a 11.59 sin que la superficie cambiara (el spline se ajusta con los mismos puntos, solo cambia dónde se evalúa); despejadas las dos componentes, el error propio del interpolador eran ~4.8 m en ambos casos. Para comparar resoluciones, compara superficie por rango. Devuelve `(metricas, puntos_validacion)`: `metricas` es un dict (`rms`, `error_medio`, `error_abs_medio`, `p50`/`p90`/`p95`/`max` del error absoluto, `n_puntos_prueba`) y `puntos_validacion` un GeoDataFrame de los puntos reservados con columna `residual` (real − interpolada) para mapear dónde se concentra el error y `curva` con el índice de la curva de la que salió cada punto: con esa columna, cazar una curva con la cota mal capturada es un `groupby("curva")["residual"].mean()` y la curva mala salta con la media sesgada, sin simbolizar mapas. Función pura y aparte de la cadena de producción: verifica el MDE sin alterar `interpolar_mde`. Falla (`ValueError`) si falta CRS, el campo no existe, las curvas están vacías, `porcentaje_reserva` sale de `(0,1)` o no hay puntos para reservar al menos 10 de prueba.
- `rasterizar(gdf, forma, transform, valor=1, relleno=0, tipo=uint8, todo_tocado=False)`: vector → máscara raster. `todo_tocado=False` (igual que ArcMap) pinta el píxel solo si su **centro** cae dentro de la geometría; `True` pinta cualquier píxel que la geometría toque. Usa `True` cuando la máscara se vaya a recortar en vector después: sobra por fuera y el `recortar` la deja en el límite exacto. Con `False`, la celda mordida por el lindero cuyo centro quedó afuera no se poligoniza nunca y el recorte posterior ya no puede devolverla, así que la capa sale con menos superficie que el polígono de referencia y nadie se entera: la fuga vale ≈ `perímetro · px/8` y por tanto se duplica al duplicar el tamaño de píxel (medido: 14.25 ha sobre un predio de 4806 ha a 20 m/px). Déjalo en `False` cuando la máscara **es** el resultado, como al quemar cauces.
- `quemar_cauces(mde, cauces, transform, profundidad)`: *stream burning*: baja la elevación del MDE en los cauces; devuelve `(array, nº_píxeles_quemados)`.
- `combinar_categorias(a, b, factor=10)`: empaqueta dos rasters categóricos en un ID único (`a*factor+b`); recuperar con `//factor` y `%factor`. Falla si `factor` no supera al máximo de `b` (los IDs se pisarían) o si alguno trae valores negativos: el `nodata_valor=-1` de `reclasificar_rangos` es el caso típico, reemplázalo antes de combinar.
- `poligonizar(arr, transform, crs, mascara=None, campo="valor")`: raster → polígonos vector (acepta máscara booleana o entera). Ignora el valor 0, que se toma como fondo: si es una de tus clases, recodifícala antes. **Usa `mascara` para no vectorizar lo que está fuera del área de interés**: recortar en raster antes es mucho más barato que poligonizar de más y recortar el vector después. Con la salida vacía devuelve un GeoDataFrame con columna `geometry` y CRS, usable sin comprobar el largo primero.
- `limpiar_moteado(arr, min_pixeles, conectividad=4, mascara=None)`: absorbe las regiones menores a `min_pixeles` en su vecino más grande. Es el **Sieve** de ArcMap/GDAL: quita el moteado de un raster categórico (el píxel suelto, la mancha de tres celdas) que no significa nada en campo y que al vectorizar se vuelve miles de polígonos diminutos. Va entre `reclasificar_rangos` y `poligonizar`. `min_pixeles` es una decisión de campo, no técnica (a 10 m/px, 4 píxeles = 400 m²); `conectividad` 4 (solo lados) u 8 (también diagonales); pásale `mascara` con el NoData excluido (`codigos > 0`) o las regiones del borde se absorben en él. Devuelve array nuevo, no muta la entrada.
- `acondicionar_mde(mde_path)`: pysheds: rellena pits, depresiones y zonas planas; devuelve `(grid, mde_acondicionado)`.
- `direccion_flujo(grid, dem)`: dirección de flujo D8.
- `acumulacion_flujo(grid, fdir)`: acumulación de flujo (array numpy).
- `etiquetar_cuencas(acumulacion, umbral)`: laderas de no-cauce (acumulación ≤ umbral) etiquetadas como cuencas; devuelve `(array, n_cuencas)`.
- `calcular_aspecto(mde, res, nodata=-9999)`: orientación cardinal por píxel (0=plano/sin orientación, 1=N, 2=E, 3=S, 4=O). Deriva con Horn sobre la ventana 3×3, la misma derivada que `calcular_pendiente`, y aplica la misma regla de NoData de Esri (centro NoData o <7 vecinos válidos → 0). No lleva parámetro de suavizado: la ventana de Horn ya promedia los 8 vecinos, y suavizar aparte hace que pendiente y aspecto no coincidan en qué es plano.
- `calcular_pendiente(mde, res, unidad="porcentaje", nodata=-9999)`: pendiente por píxel con el algoritmo Planar/Horn (1981), el método por defecto de ArcMap 10.3.1 (diferencias finitas de 3er orden sobre ventana 3×3, réplica exacta). `res` acepta un número (celda cuadrada) o `(res_x, res_y)`; `unidad` = `"porcentaje"` o `"grados"`. Aplica la regla de NoData de Esri (centro NoData → NoData; <7 vecinos válidos → NoData, marcados como NaN). Devuelve array float32.
- `reclasificar_rangos(arr, rangos, nodata_valor=-1)`: reclasifica un raster continuo en códigos enteros por rangos (genérica, no solo pendientes). `rangos` es una lista de `(etiqueta, mínimo, máximo)`; convención de límites igual que Reclassify de ArcGIS: cada rango es `(mín, máx]`, salvo el primero que es `[mín, máx]`, así el límite compartido cae en el rango de abajo. Los NoData/NaN quedan en `nodata_valor`. Devuelve `(array_códigos int32, dict {código: etiqueta})` para mapear luego a una columna de texto.

## Flujo de pendientes (ArcMap → cfasig)

Reimplementación del flujo de ArcMap (Superficie > Pendiente → Reclasificar → Sieve → Raster a polígono → recorte → Disolver por gridcode). `SUP` se calcula **después** de recortar para reflejar el área real recortada al predio; `codigo` se conserva en la salida junto con `PENDIENTES` y `SUP`.

**Toda la limpieza va en raster, antes de vectorizar.** Es la diferencia entre poligonizar cientos de miles de esquirlas o unos miles de polígonos con significado. Dos pasos, en este orden:

1. **Recorte al predio con `mascara`.** El MDE se interpola sobre el *bounding box* del predio, así que sin máscara se vectorizan las esquinas que después se tiran.
2. **`limpiar_moteado` (Sieve).** El moteado sub-mínimo no significa nada en campo y es el grueso del conteo de polígonos.

El recorte vectorial (`recortar`) **no se quita**: la máscara raster corta por bordes de píxel y el clip da el límite exacto del predio, que es lo que hace cuadrar la superficie.

```python
mde, perfil, transform, crs, res = sig.cargar_raster(RUTA_MDE)
pendiente = sig.calcular_pendiente(mde, res)
codigos, etiquetas = sig.reclasificar_rangos(pendiente, RANGOS_PENDIENTE)

predio = sig.asegurar_crs(sig.cargar(RUTA_PREDIO), EPSG_UTM, "predio")
mascara = (sig.rasterizar(predio, codigos.shape, transform) == 1) & (codigos > 0)
codigos = sig.limpiar_moteado(codigos, MIN_PIXELES, mascara=mascara)

gdf = sig.poligonizar(codigos, transform, crs, mascara=mascara, campo="codigo")
gdf["PENDIENTES"] = gdf["codigo"].map(etiquetas)
gdf = sig.recortar(gdf, predio)                                    # borde exacto
gdf = sig.disolver_por_grupo(gdf, campo="codigo")
gdf = sig.calcular_superficie(gdf, campo="SUP", en_hectareas=True)  # 2 decimales
sig.guardar(gdf, RUTA_SALIDA)
```

**La pendiente se calcula sobre el MDE interpolado, no sobre el quemado.** `quemar_cauces` baja los cauces varios metros a propósito y `acondicionar_mde` (`resolve_flats`) mete micro-pendientes en los planos: los dos son correctos para análisis de flujo y los dos arruinan la pendiente. Un quemado de 5 m a 10 m/px es un escalón de 50 % en cada orilla de cauce, así que el rango `>100` se llena de líneas de cauce que no existen. Guarda el MDE antes de quemar y úsalo aquí.

## Flujo de MDE (curvas → raster, extent fiel a ArcMap)

Construcción de MDE desde curvas de nivel replicando el manejo de extent de Topo to Raster (ArcMap): las curvas y cauces se recortan a un cuadro amplio (predio + borde de N celdas) y el MDE se interpola sobre el extent de salida (el predio), que queda metido dentro de los datos para que las celdas de borde se interpolen con datos a ambos lados. `validar_mde` corre con los mismos parámetros (resolución, cuadro, método) para que el *hold-out* corresponda al raster final.

```python
RESOLUCION = 20          # m/px (None = media separación horizontal entre curvas)
cuadro_recorte, cuadro_salida = sig.cuadros_arcmap(
    predio, resolucion=RESOLUCION, margen_borde_celdas=10  # borde de datos = 200 m a 20 m/px
)
curvas = sig.recortar(curvas, cuadro_recorte)
cauces = sig.recortar(cauces, cuadro_recorte)

metricas, puntos = sig.validar_mde(curvas, campo_elevacion="ELEVACION",
                                   resolucion=RESOLUCION, cuadro=cuadro_salida, metodo="spline")
mde, transform = sig.interpolar_mde(curvas, campo_elevacion="ELEVACION",
                                    resolucion=RESOLUCION, cuadro=cuadro_salida, metodo="spline")
perfil = sig.perfil_raster(transform, EPSG_UTM)      # nodata=NaN declarado

# rellenar ANTES de quemar: el interpolador deja NaN fuera de la envolvente
# convexa, ahí NaN - profundidad sigue siendo NaN y el relleno posterior
# copiaría el vecino SIN quemar, perdiendo el cauce en silencio.
mde = sig.rellenar_nodata(mde)
sig.guardar_raster(mde, "MDE_BASE.tif", perfil)  # el de pendientes: sin quemar
mde_final, n = sig.quemar_cauces(mde, sig.disolver(cauces), transform, profundidad=5)
if n == 0:
    raise ValueError("ningún píxel quemado: revisa el CRS de la hidrología")

# pysheds lee del disco: guardar la base quemada y sobreescribir ya acondicionada
sig.guardar_raster(mde_final, "MDE_HIDRO.tif", perfil)
grid, dem = sig.acondicionar_mde("MDE_HIDRO.tif")
sig.guardar_raster(dem, "MDE_HIDRO.tif", perfil)     # dem es Raster de pysheds, entra tal cual
```

## Flujo de clasificación de superficies (capas que se pisan)

Dos maneras de resolver un solape, y el orden entre ellas importa. Con jerarquía clara manda el orden de la lista (`resolver_por_prioridad`); sin jerarquía gana el núcleo más cercano (`resolver_por_cercania`). El polígono general entra al final como categoría de relleno, así la suma cuadra con el predio por construcción.

**Primero resta, luego reparte.** A las capas sin jerarquía se les quita antes lo que se llevan las capas de mayor prioridad, y hasta entonces se reparte por cercanía. Si se hace al revés, el reparto afina al metro un borde que la cascada borra después: las hectáreas que reporta el reparto no son las que acaban en el shapefile. Dos casos que solo aparecen con este orden y conviene contemplar: una capa que queda vacía porque las superiores la cubren entera (sáltala, una capa vacía no disputa nada), y un solape que se rompe en esquirlas dispersas (ver el tope de malla en `resolver_por_cercania`).

```python
superiores = sig.disolver(sig.combinar([caminos, hidrologia, mayor100]))
produccion = sig.quitar_solape(produccion, superiores)
restauracion = sig.quitar_solape(restauracion, superiores)

produccion, restauracion, conflicto = sig.resolver_por_cercania(
    produccion, restauracion, nombres=("PRODUCCION", "RESTAURACION")
)
sig.guardar(conflicto, RUTA_CONFLICTOS)      # revisar el reparto a ojo

conservacion = sig.disolver(predio, campo="CLASS_SUP", valor="CONSERVACION")
combinado = sig.resolver_por_prioridad(CAPAS + [conservacion])   # orden = prioridad
combinado = sig.recortar(combinado, predio)  # recorte de seguridad
gdf = sig.disolver_por_grupo(combinado, campo="CLASS_SUP")
gdf = sig.calcular_superficie(gdf, campo="SUP", en_hectareas=True)  # área al final
```

## Ejemplos de uso (`ejemplos/`)

`ejemplos/caminos.py` como referencia. Los otros dos (`ejemplos/hidrologia.py`, `ejemplos/mde.py`) siguen el mismo patrón; `mde.py` es la versión ejecutable del "Flujo de MDE" de más arriba.

```python
"""
Ejemplo: buffer jerárquico de caminos (primario > secundario > saca).

Cada nivel recorta al inferior para que no se solapen. Muestra el uso
de cfasig con un patrón repetido resuelto en un bucle.
"""

import cfasig as sig

# ── CONFIG ───────────────────────────────────────────────────
EPSG_UTM = 32613

# nombre, ruta, buffer (en metros). Orden = prioridad (mayor a menor).
NIVELES = [
    ("PRIMARIO", "datos/camino_primario.shp", 5.0),
    ("SECUNDARIO", "datos/camino_secundario.shp", 3.0),
    ("SACA", "datos/camino_saca.shp", 1.75),
]
CAMPO_TIPO = "camino"

SALIDA_COMBINADO = "salidas/CAMINOS_BUFFER_COMBINADO.shp"
SALIDA_DISUELTO = "salidas/CAMINOS_BUFFER_DISUELTO.shp"

# ── 1. CARGAR, ASEGURAR CRS, DISOLVER Y BUFFER CADA NIVEL ────
capas = []
for nombre, ruta, dist in NIVELES:
    gdf = sig.cargar(ruta)
    gdf = sig.asegurar_crs(gdf, EPSG_UTM, nombre)
    gdf = sig.disolver(gdf, campo=CAMPO_TIPO, valor=nombre)
    gdf = sig.buffer(gdf, dist)
    capas.append(gdf)

# ── 2. QUITAR SOLAPE: cada nivel pierde contra TODOS los superiores ──
# (recortar solo contra el nivel de arriba dejaba pasar el solape saca-primario
#  donde no hay secundario en medio, que es justo cada entronque)
combinado = sig.resolver_por_prioridad(capas)

# ── 3. GUARDAR ───────────────────────────────────────────────
sig.guardar(combinado, SALIDA_COMBINADO)

# ── 4. VERSIÓN DISUELTA (una sola geometría) ─────────────────
sig.guardar(sig.disolver(combinado, campo="caminos", valor="caminos"), SALIDA_DISUELTO)
```

**Estilo y estructura de un script** (confirmado contra los tres scripts de `ejemplos/`):

- Docstring corto al inicio: qué hace y por qué (una frase de "qué", una de "por qué" si no es obvio).
- Bloque `CONFIG` arriba de todo, en mayúsculas: rutas, EPSG, distancias, campos. Nada de esto va suelto en medio de la lógica; si hay que tocar un valor para correrlo en otro predio, se toca aquí y solo aquí.
- Cuando hay una lista de "cosas parecidas" y vale la pena (3+ niveles/casos), se modela como lista de tuplas (`NIVELES` en `caminos.py`) y se recorre con un `for`, en vez de repetir el bloque una vez por elemento; el orden puede codificar prioridad. Con solo 2 casos (perenne/intermitente en `hidrologia.py`) no se generaliza: se escriben los dos bloques directo, sin loop ni lista, porque el loop no ahorra nada con dos ramas. Loop si repites de verdad, directo si no (YAGNI).
- Pasos numerados con comentarios `# ── N. VERBO EN MAYÚSCULAS ──`: cada bloque es una etapa clara del flujo (cargar, quitar solape, combinar, guardar). Sirve para ubicarse en scripts de 40-100 líneas sin funciones propias.
- El script encadena funciones de `cfasig` (`sig.cargar`, `sig.buffer`, ...); no reimplementa nada que la librería ya resuelva. La única lógica que vive en el script es la específica del caso (qué niveles hay, qué prioridad tienen, qué profundidad de quemado o margen de extent pide la cartografía).
- `print()` como reporte, no logging: además de los avisos que ya imprimen funciones como `cargar` (nº de entidades, CRS), el script puede imprimir sus propios diagnósticos a media ejecución (valores únicos de un campo, tamaño de un raster, píxeles quemados) y, si el resultado tiene métricas que valen la pena, un resumen al final (conteos, sumas, rangos): ver `mde.py` e `hidrologia.py`. Si el script ya va a imprimir lo suyo, usa `mostrar=False` en `cargar`/`cargar_raster` para no duplicar el aviso.
- Sin funciones ni clases propias salvo que el script se vuelva a llamar con distintos parámetros; un script de un solo uso es lineal de arriba a abajo, reutilizando el mismo nombre de variable al reasignar (`gdf = sig.algo(gdf)`) en vez de encadenar nombres nuevos por paso.
- Nombres de variable en español, cortos y descriptivos (`gdf`, `capas`, `combinado`); constantes de config en mayúsculas.

## Instalación / uso

```bash
# desde la carpeta que contiene pyproject.toml
pip install -e .

# luego, desde cualquier script:
import cfasig as sig
```

Con layout `src/`, la instalación (`pip install -e .`) es obligatoria: sin ella `import cfasig` no encuentra el paquete.

## CLI (`cfasig`)

La instalación registra el comando de consola `cfasig` (`project.scripts` → `cli.main`). Hoy hace una sola cosa: convertir archivos entre formatos, uno o varios en la misma llamada (proceso por lote). La salida se autonombra (misma carpeta y nombre, extensión nueva).

```bash
cfasig predio.shp kmz            # crea predio.kmz
cfasig ruta.gpx shp              # crea ruta.shp
cfasig a.shp b.shp c.shp kmz     # lote: crea a.kmz, b.kmz, c.kmz
cfasig ayuda                     # muestra la ayuda
```

Uso: `cfasig ARCHIVO [ARCHIVO ...] TIPO`. Para el lote se pasan varios archivos y el TIPO al final; pensado para arrastrar los archivos a la terminal (quedan separados por espacios), sin glob ni comodines porque el usuario no técnico los arrastra. Con varios archivos imprime un resumen final (`Listo: N convertidos, M con error`) y sigue con los demás si uno falla. Un archivo que ya está en el formato destino se salta.

Tipos válidos: `shp`, `gpkg`, `geojson`, `gpx`, `kml`, `kmz`. Opciones avanzadas: `--capa NOMBRE` (capa a leer en archivos multicapa), `--gpx-como {track,route}` y `--si` (sobrescribe salidas existentes sin preguntar; útil en lotes o scripts). Sin `--si` pide confirmación antes de sobrescribir cada archivo. Ayuda y errores en español; la ayuda no carga geopandas (import perezoso). `# ponytail:` sin subcomandos porque solo convierte; migra a subparsers si crecen las operaciones.

## Tests

Suite en `tests/`, un archivo por módulo/tema (`test_archivo`, `test_proyeccion`, `test_geometria`, `test_raster`, más `test_conversion` para el ida y vuelta entre formatos, `test_puntos` para carga de puntos CSV/Excel y su conversión a línea/polígono, `test_gpx` para los lectores limpios de waypoints/tracks, y `test_cli` para el comando de consola), con fixtures compartidas en `conftest.py`. Cubre la lógica determinista de vector y raster; los helpers de flujo de pysheds (`acondicionar_mde`, `direccion_flujo`, `acumulacion_flujo`) quedan pendientes hasta tener un MDE de prueba.

```bash
pip install -e ".[dev]"   # instala pytest
pytest
```

## Convenciones y decisiones

- Wrapper delgado: una función = una operación de la librería base (vector o raster), con nombre simple en español.
- La lógica de dominio (claves de rodal, valores calibrados, umbrales del predio) se queda en los scripts; la librería solo aporta operaciones SIG genéricas.
- Raster: las funciones trabajan sobre arrays 2D (numpy) más su `transform`/`crs`, que es como rasterio maneja los datos. Los scripts no necesitan importar numpy para el flujo normal.
- Las funciones que devuelven capa (`buffer`, `quitar_solape`, `limpiar_vacias`) trabajan sobre una copia; no mutan la entrada.
- Distancias y márgenes en unidades del CRS activo (usar UTM para metros).
- `disolver` y `combinar` conservan el CRS de origen; `recortar` alinea CRS automáticamente.
- Prioridad entre capas modelada con `quitar_solape` (la capa `otro` gana el solape). Con tres o más capas, `resolver_por_prioridad` en vez del bucle a mano: recorta contra la unión acumulada, no contra la vecina.
- Cuando dos capas se solapan pero ninguna manda sobre la otra, `resolver_por_cercania` reparte el conflicto en vez de que una gane entera. Es la única función de vector que baja a raster por dentro (`geometria` importa `rasterizar`/`poligonizar` de `raster`); la malla es temporal y no sale de la función.
- Mensajes informativos por consola (`print`); sin logging formal todavía.
- Seguros implementados dentro de cada función: un helper compartido, `_exigir_crs(gdf, nombre)` (proyeccion.py, usado en proyeccion/geometria/raster), más chequeos en línea (campo inexistente, ruta inexistente, lista/capa vacía, factor que colisiona, máscara que no intersecta). Fallan con `raise ValueError` de mensaje claro.
- EPSG de trabajo típico: 32613 (UTM 13N, WGS84) para la zona de operación.

## Pendientes / ideas para crecer

- Completar los tests de flujo hidrológico (pysheds) con un MDE pequeño de prueba: es lo único de la API que aún no tiene cobertura.
- **Detección de áreas sin vegetación desde imagen digital** (hoy se digitaliza a mano en ArcMap 10.3.1). La cadena ya existe casi entera: `recortar_raster` → *[índice]* → `poligonizar` → `cerrar_microhuecos` → `fusionar_menores` → `simplificar` → `calcular_superficie` → `guardar`. Faltarían dos funciones cortas, sin dependencias nuevas (rasterio/numpy/scipy ya están):
  - `cargar_bandas(ruta, bandas)`: `cargar_raster` solo lee la banda 1, aquí hacen falta 2-3.
  - `indice_vegetacion(bandas, tipo)` + umbral → máscara booleana. NDVI si la imagen trae NIR; VARI/ExG si es RGB puro (captura de ArcMap, Google/Bing, dron sin NIR). Umbral fijo o por Otsu.

  A debatir antes de implementar: el código es trivial, el problema es la clasificación. Sombra de arbolado y pasto seco entran como suelo desnudo; el umbral se mueve con fecha, hora y tipo de suelo; sin NIR la separación es notoriamente peor; y Sentinel-2 a 10 m/px queda grueso para un predio chico. Sirve como **primer borrador** que el técnico corrige (ahí está el ahorro de tiempo de digitalización), no como salida final sin revisar. La descarga de la imagen queda fuera de la librería: depende de la fuente.

  **Sombras: probar primero la regla de brillo.** El suelo desnudo es brillante y la sombra es oscura; ambos dan índice de vegetación bajo, y por eso se confunden, pero en brillo se separan casi limpio. Dos condiciones (`indice < umbral_veg` **y** `brillo > umbral_brillo`) en vez de una, cero dependencias nuevas. Debe probarse esto antes de considerar cualquier modelo. Dos matices: NDVI aguanta la sombra mejor que VARI/ExG porque al ser cociente se auto-normaliza parcialmente contra la iluminación (otro argumento para conseguir imágenes con NIR); y si la sombra es de relieve y no de arbolado, ya se puede modelar en casa con hillshade desde el MDE usando azimut y elevación solar de la fecha.

  **Opción Random Forest (clasificación supervisada) — solo si la regla de brillo no alcanza.** El técnico digitaliza 20-50 polígonos de muestra por clase (`vegetación`, `suelo desnudo`, `sombra`, `camino`, `agua`), cada píxel se vuelve un vector de variables (bandas + índice + brillo + textura con `uniform_filter` de `scipy.ndimage`, + pendiente de `calcular_pendiente`), y el clasificador sustituye al umbral. La cadena posterior (`poligonizar` → `cerrar_microhuecos` → `fusionar_menores` → ...) no cambia. A favor: mata el problema del umbral (el técnico enseña con ejemplos en ese predio y esa imagen, en vez de cazar un número mágico); la sombra pasa a ser una clase más en lugar de un error que corregir; mezcla variables de distinta naturaleza sin normalizar; da probabilidad por píxel, lo que permite marcar como "revisar" solo lo que salga con confianza baja; corre en CPU en segundos.

  Pero abre varios problemas nuevos, y por eso queda como opción a debatir, no como plan:

  - Dependencia nueva (`scikit-learn`) más ~50 líneas de wrapper, contra las ~40 líneas y cero dependencias de la versión con umbral.
  - Clasifica píxel a píxel, sin contexto espacial: produce sal y pimienta. Mitigable con las variables de textura y con `fusionar_menores`, pero es trabajo extra.
  - **No transfiere entre imágenes.** Un modelo entrenado en una fecha se degrada en otra. Implica reentrenar por predio, y por lo tanto definir dónde viven los polígonos de entrenamiento, con qué convención de clases y quién los mantiene. Eso es flujo de trabajo, no una función.
  - Trampa de validación: apartar píxeles al azar da exactitudes falsamente altas porque los píxeles vecinos del mismo polígono quedan de los dos lados. Hay que reservar **polígonos completos**, en la línea de lo que hace `validar_mde` con su reserva. Es una función de validación aparte que también habría que escribir.
  - Entrenamiento sucio, salida sucia: la calidad depende de la pureza de los polígonos de muestra, que es responsabilidad del técnico y no del código.

  Orden sugerido si se retoma: índice + umbral → regla de brillo → medir cuánto falta → recién ahí decidir sobre Random Forest.
