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

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.6.0 | 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
• 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.

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).
• 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="tin", 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: "tin" (triangulación de Delaunay, respeta quiebres de pendiente), "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) o "idw" (distancia inversa, respaldo si el TIN falla por geometría degenerada). Fuera de la envolvente convexa devuelve NaN (no extrapola; rellena luego con rellenar_nodata). Descarta curvas con cota NaN. resolucion por defecto = la equidistancia detectada (los metros de la carta = el intervalo entre curvas); equidistancia se autodetecta (moda de las diferencias entre cotas). 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="tin", 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) y compara la elevación real contra la del MDE en los puntos reservados. semilla fija el reparto para reproducibilidad. 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. 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): vector → máscara raster.
• 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, suavizar=3, nodata=-9999): orientación cardinal por píxel (0=plano/sin orientación, 1=N, 2=E, 3=S, 4=O). El suavizado promedia solo vecinos válidos (el NoData no crea acantilados falsos en los bordes) y los píxeles NoData quedan como 0.
• 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.

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.

RESOLUCION = 20          # m/px (None = equidistancia de la carta)
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).

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.

"""
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

# 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).

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.

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, que ya se importa en raster.py, + 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.