Curso: CC3084 · Ciencia de Datos · Semestre 02, 2026 · UVG
Vamos a analizar una ventana de 61 × 61 km del altiplano central de Guatemala. Tenemos tres insumos: los límites de 59 municipios, una escena multiespectral de seis bandas a 60 metros por píxel, y 3000 observaciones de campo con la cobertura del suelo que se observó en cada punto.
Cargamos los municipios y vemos cómo se guarda una geometría. Medimos el área de uno de ellos en grados y en metros para ver cuánto cuesta equivocarse de sistema de coordenadas, y lo dibujamos en tres proyecciones distintas. Guardamos la misma capa en shapefile, GeoPackage y GeoJSON para ver qué se pierde en cada formato.
Después abrimos la escena. Sacamos la firma espectral de cada cobertura, comparamos color natural con falso color y calculamos NDVI, NDWI y NDBI. Degradamos la imagen en sus cuatro resoluciones para ver qué información desaparece con cada una, y simulamos un levantamiento LiDAR para separar la altura del dosel del terreno. Luego generamos cuatro fechas con nubes, las enmascaramos, las combinamos en un compuesto de mediana y le corregimos la atmósfera.
Con la imagen limpia hacemos operaciones SIG. Recortamos, unimos puntos con polígonos, generamos buffers, disolvemos municipios en regiones y sacamos estadística zonal. Encadenamos varias operaciones para decidir dónde conviene ubicar un centro de salud nuevo. Medimos autocorrelación espacial con el índice de Moran implementado desde cero, localizamos conglomerados con LISA y Getis-Ord, agrupamos puntos con DBSCAN e interpolamos una superficie con IDW.
La última parte es aprendizaje automático. Entrenamos un Random Forest para clasificar la cobertura de cada píxel de la escena, y después lo evaluamos de cuatro maneras distintas: partiendo los datos al azar y partiéndolos por zonas geográficas cada vez más separadas. Las exactitudes que salen no se parecen entre sí, y esa diferencia es el punto central del cuaderno. Cerramos probando k-means sin etiquetas, agregando variables espaciales derivadas, y pasando del píxel al parche de 7×7 para conectar con las redes convolucionales.
Al final revisamos las trampas del análisis espacial —MAUP, falacia ecológica, efecto de borde— y la ética de trabajar con ubicaciones de personas.
Los datos
Insumo
Qué es
Origen
municipios.gpkg
59 municipios de Guatemala
Reales — geoBoundaries ADM2 (CC-BY 4.0)
escena_multiespectral.tif
6 bandas a 60 m, EPSG:32615
Simulada a partir de firmas espectrales publicadas
cobertura_verdad.tif
La cobertura que generó la escena
Derivado de la simulación
observaciones_campo.geojson
3000 puntos etiquetados
Derivado, con muestreo agrupado
La escena es simulada porque bajar Sentinel-2 o Landsat reales exige credenciales de Copernicus o USGS, y eso rompe la reproducibilidad. La geografía y el CRS sí son reales: los municipios son de Guatemala y están en UTM 15N. A cambio de perder autenticidad en la reflectancia, ganamos algo que ninguna imagen real da: conocemos la clase verdadera de cada píxel, así que podemos medir los errores en vez de estimarlos.
1. Imports y datos
Fijamos la semilla para que el cuaderno sea reproducible y construimos el conjunto de datos si aún no existe.
Código
import sysimport jsonfrom pathlib import Pathimport numpy as npimport pandas as pdimport geopandas as gpdimport matplotlib.pyplot as pltfrom matplotlib.colors import ListedColormapimport rasteriofrom rasterio.features import rasterizefrom shapely.geometry import boxfrom matplotlib.patches import RectangleRAIZ = Path.cwd().parent if Path.cwd().name =="notebooks"else Path.cwd()sys.path.append(str(RAIZ /"src"))import build_dataset as bdnp.random.seed(42)plt.rcParams.update({"figure.dpi": 110,"font.size": 9,"axes.titlesize": 10,"axes.grid": True,"grid.alpha": 0.25,})# Paleta consistente para las 5 clases de cobertura en todo el cuadernoCOLORES = ["#1c7ed6", "#2b8a3e", "#a9e34b", "#e03131", "#d9a441"]CMAP_CLASES = ListedColormap(COLORES)rutas = bd.ensure_dataset(size=1024, points=3000, seed=42)meta = json.loads(Path(rutas["metadata"]).read_text())CLASES, BANDAS = meta["clases"], meta["bandas"]print(f"Escena {meta['tamano_px']}x{meta['tamano_px']} px de {meta['resolucion_m']:.0f} m "f"({meta['extension_km']:.0f} km de lado) en {meta['crs']}")print(f"Municipios: {meta['n_municipios']} — {meta['municipios']}")
Los datos ya existen; no se regenera nada (usa --force para rehacer).
Escena 1024x1024 px de 60 m (61 km de lado) en EPSG:32615
Municipios: 59 — geoBoundaries gbOpen GTM ADM2 (CC-BY 4.0)
2. La ubicación como dato
Qué es un dato geoespacial
Casi todo dato lleva pegada una pregunta implícita: ¿dónde ocurrió? Una venta ocurre en una tienda, un caso de dengue en una aldea, una deforestación en una parcela concreta. Cuando un dato incluye su ubicación en la superficie de la Tierra, decimos que es geoespacial.
Un registro geoespacial tiene tres componentes:
Componente
Pregunta
Ejemplo
Espacial
¿dónde está?
una geometría (punto, línea, polígono) o una celda de una cuadrícula
Temporal
¿cuándo?
la misma parcela no es igual en enero que en junio
Atributos
¿qué hay ahí?
población, temperatura, tipo de cultivo, número de casos
Un dato con las tres dimensiones se llama espacio-temporal, y es el caso más común en la práctica: series de tiempo por ubicación.
La primera ley de la geografía
“Todo está relacionado con todo lo demás, pero las cosas cercanas están más relacionadas que las lejanas.”
Es la idea que sostiene casi todo el análisis espacial:
El precio de una casa se parece al de la casa de al lado.
La temperatura de hoy en un punto se parece a la de un kilómetro más allá.
Un brote de enfermedad tiende a aparecer en vecindarios contiguos.
Esa dependencia entre observaciones cercanas se llama autocorrelación espacial, y rompe un supuesto básico de la estadística clásica: que las observaciones son independientes.
El modelo vectorial
Hay dos maneras de representar el mundo. La primera es el modelo vectorial: geometrías definidas por coordenadas.
Geometría
Qué es
Ejemplo
Punto
una ubicación sin dimensión
un hospital, un árbol, un evento de GPS
Línea
una secuencia de puntos conectados
una carretera, un río, una ruta de bus
Polígono
una línea cerrada que encierra un área
un municipio, una parcela, un lago
Multi-geometría
un registro con varias partes
un archipiélago es un multipolígono
Cada geometría lleva una tabla de atributos asociada. Es el modelo natural para objetos discretos y para datos administrativos.
En Python eso vive en un GeoDataFrame: un DataFrame de pandas con una columna extra, la geometría. Todo lo que ya sabes de pandas sigue funcionando; encima aparecen las operaciones espaciales.
# La columna de geometria no guarda texto: guarda objetos geometricosg = munis.geometry.iloc[0]print("tipo de geometria :", g.geom_type)print("vertices :", len(g.exterior.coords) if g.geom_type =="Polygon"elsesum(len(p.exterior.coords) for p in g.geoms))print("area (km2) :", round(g.area /1e6, 1))print("centroide (UTM) :", (round(g.centroid.x), round(g.centroid.y)))
tipo de geometria : Polygon
vertices : 233
area (km2) : 131.0
centroide (UTM) : (720209, 1609302)
3. CRS, datums y proyecciones
La forma de la Tierra
La Tierra no es una esfera perfecta ni una superficie plana, y esa incomodidad se cuela en todo dato geoespacial. Tres conceptos, en orden de abstracción:
Geoide: la forma real de la Tierra, definida por su campo gravitatorio. Es irregular, con bultos y hundimientos. Nadie calcula sobre ella directamente.
Elipsoide: una aproximación matemática al geoide —una esfera achatada en los polos—, mucho más simple de calcular.
Datum: el elipsoide más su anclaje al planeta. Define el origen desde el cual se miden las coordenadas. El más usado es WGS84, el del GPS.
Por qué importa el datum: dos coordenadas idénticas en datums distintos pueden apuntar a lugares separados por cientos de metros. Un dato de latitud y longitud sin su datum está incompleto.
Coordenadas geográficas
Latitud: ángulo hacia el norte o el sur desde el ecuador, de −90° a 90°.
Longitud: ángulo hacia el este o el oeste desde Greenwich, de −180° a 180°.
Se expresan en grados, no en metros, y ahí está la trampa: un grado de latitud mide unos 111 km en cualquier parte del planeta, pero un grado de longitud mide 111 km en el ecuador y 0 km en los polos, porque los meridianos convergen.
Consecuencia práctica: calcular distancias o áreas restando grados es incorrecto. No es una aproximación tosca — es una operación sin sentido físico. Vamos a comprobarlo con números, y a ver de qué tamaño es el problema.
Tomamos el mismo municipio en dos sistemas de referencia:
EPSG:4326 — WGS84, coordenadas en grados.
EPSG:32615 — UTM zona 15N, coordenadas en metros.
Código
muni = munis.sort_values("area_km2", ascending=False).iloc[[0]]nombre = muni["municipio"].iloc[0]area_metrica = muni.to_crs(32615).geometry.area.iloc[0] /1e6# km2 realesarea_grados = muni.to_crs(4326).geometry.area.iloc[0] # grados^2 (!)# La trampa clasica: convertir grados^2 a km2 como si un grado fuera constanteKM_POR_GRADO =111.32area_ingenua = area_grados * KM_POR_GRADO **2print(f"Municipio: {nombre}\n")print(f" Area en EPSG:32615 (metros) : {area_metrica:10.2f} km2 <- correcta")print(f" Area en EPSG:4326 (grados^2) : {area_grados:10.6f} grados^2 <- sin sentido fisico")print(f" 'Convirtiendo' con 111.32 km/grado: {area_ingenua:6.2f} km2")print(f"\n Error de la conversion ingenua: {(area_ingenua/area_metrica -1) *100:+.1f} %")
Municipio: Escuintla
Area en EPSG:32615 (metros) : 546.71 km2 <- correcta
Area en EPSG:4326 (grados^2) : 0.045770 grados^2 <- sin sentido fisico
'Convirtiendo' con 111.32 km/grado: 567.19 km2
Error de la conversion ingenua: +3.7 %
/var/folders/97/04drgvm13z7bxfszg04r31f40000gn/T/ipykernel_90688/3718448756.py:5: UserWarning: Geometry is in a geographic CRS. Results from 'area' are likely incorrect. Use 'GeoSeries.to_crs()' to re-project geometries to a projected CRS before this operation.
area_grados = muni.to_crs(4326).geometry.area.iloc[0] # grados^2 (!)
Dos cosas que vale la pena separar, porque se confunden todo el tiempo:
El número en grados cuadrados no es un área. Es un producto de dos ángulos. No tiene unidades de superficie, y no se vuelve un área por multiplicarlo por algo.
El atajo de convertir con 111.32 km/grado sí es una aproximación, y aquí resulta bastante buena: se queda a menos del 4 %. Guatemala está cerca del ecuador, donde un grado de longitud todavía mide casi lo que mide allí.
O sea que el atajo no explota en esta latitud. Pero su error depende directamente de dónde estés, y crece rápido.
Fíjate en la advertencia que imprimió GeoPandas arriba: Geometry is in a geographic CRS. Results from 'area' are likely incorrect. La librería está avisando del error exacto que acabamos de cometer a propósito. Cuando esa advertencia aparezca en tu código, no la silencies: reproyecta.
Código
lat = munis.to_crs(4326).geometry.centroid.y.mean()print(f"Latitud media de la ventana: {lat:.2f}° N\n")print(f" 1° de latitud ≈ {110.57:.2f} km (casi constante en todo el planeta)")print(f" 1° de longitud ≈ {111.32* np.cos(np.radians(lat)):.2f} km en esta latitud")print(f" 1° de longitud ≈ {111.32:.2f} km en el ecuador")print(f" 1° de longitud ≈ {0.0:.2f} km en el polo")print("\nError del atajo 'grados^2 x 111.32^2' según la latitud:")for lat_prueba in [0, 15, 30, 45, 60, 75]: error = (1/ np.cos(np.radians(lat_prueba)) -1) *100 marca =" <- Guatemala"if lat_prueba ==15else""print(f" {lat_prueba:2d}° : {error:+7.1f} %{marca}")print(f"\nA {lat:.1f}° N, un grado de longitud ya perdió "f"{(1- np.cos(np.radians(lat))) *100:.1f}% de su longitud ecuatorial.")
Latitud media de la ventana: 14.59° N
1° de latitud ≈ 110.57 km (casi constante en todo el planeta)
1° de longitud ≈ 107.73 km en esta latitud
1° de longitud ≈ 111.32 km en el ecuador
1° de longitud ≈ 0.00 km en el polo
Error del atajo 'grados^2 x 111.32^2' según la latitud:
0° : +0.0 %
15° : +3.5 % <- Guatemala
30° : +15.5 %
45° : +41.4 %
60° : +100.0 %
75° : +286.4 %
A 14.6° N, un grado de longitud ya perdió 3.2% de su longitud ecuatorial.
/var/folders/97/04drgvm13z7bxfszg04r31f40000gn/T/ipykernel_90688/1004881034.py:1: UserWarning: Geometry is in a geographic CRS. Results from 'centroid' are likely incorrect. Use 'GeoSeries.to_crs()' to re-project geometries to a projected CRS before this operation.
lat = munis.to_crs(4326).geometry.centroid.y.mean()
Qué es una proyección
Una proyección es la receta para pasar de la superficie curva de la Tierra a un plano. Toda proyección distorsiona algo; no existe una perfecta. Lo único que puedes elegir es qué sacrificas:
Familia
Preserva
Deforma
Ejemplo
Conformes
los ángulos y las formas locales
las áreas
Mercator: Groenlandia parece del tamaño de África, y es 14 veces más pequeña
Equivalentes
las áreas
las formas
mapas temáticos donde importa comparar superficies
Equidistantes
las distancias
todo lo demás
y solo desde ciertos puntos o líneas
Tres reglas prácticas que resuelven el 90 % de los casos:
Para medir distancias y áreas, usa una proyección métrica local, nunca grados.
UTM divide el mundo en 60 husos de 6° de ancho y da coordenadas en metros dentro de cada uno. Guatemala cae en los husos 15 y 16 norte.
Para mapas web, la proyección de facto es Web Mercator (EPSG:3857): buena para navegar, mala para comparar áreas.
Códigos EPSG
Cada sistema de referencia de coordenadas (CRS) tiene un código EPSG que lo identifica sin ambigüedad:
Código
Qué es
Unidades
EPSG:4326
WGS84, latitud y longitud
grados
EPSG:3857
Web Mercator, la de los mapas web
metros (distorsionados)
EPSG:32615 / 32616
UTM zonas 15N y 16N sobre WGS84
metros
La regla número uno al combinar capas: verificar que todas estén en el mismo CRS. Superponer datos en CRS distintos produce mapas que se ven bien y análisis que están mal. Es la primera fuente de errores en este campo.
La misma verdad, tres proyecciones
Ninguna proyección es neutral. Aquí está exactamente el mismo conjunto de municipios en tres CRS: cambia la forma, cambia el área, cambia la impresión.
Código
proyecciones = [ (4326, "EPSG:4326 — grados\n(no es una proyección: son ángulos)"), (3857, "EPSG:3857 — Web Mercator\n(conforme: infla las áreas)"), (32615, "EPSG:32615 — UTM 15N\n(métrica local: la correcta aquí)"),]fig, axes = plt.subplots(1, 3, figsize=(13, 4.4))for ax, (epsg, titulo) inzip(axes, proyecciones): capa = munis.to_crs(epsg) capa.plot(ax=ax, facecolor="#dbe4ff", edgecolor="#364fc7", linewidth=0.4) capa[capa["municipio"] == nombre].plot(ax=ax, facecolor="#ff8787", edgecolor="#c92a2a", linewidth=0.8) ax.set_title(titulo, fontsize=9) ax.set_xticks([]); ax.set_yticks([]); ax.grid(False)fig.suptitle("El mismo territorio en tres sistemas de referencia", y=1.0)plt.tight_layout(); plt.show()print("Área del municipio resaltado, medida en cada CRS:")for epsg, _ in proyecciones: a = munis[munis["municipio"] == nombre].to_crs(epsg).geometry.area.iloc[0] unidad ="grados^2"if epsg ==4326else"km2" valor = a if epsg ==4326else a /1e6print(f" EPSG:{epsg:<6}{valor:12.4f}{unidad}")
Área del municipio resaltado, medida en cada CRS:
EPSG:4326 0.0458 grados^2
EPSG:3857 585.3212 km2
EPSG:32615 546.7081 km2
/var/folders/97/04drgvm13z7bxfszg04r31f40000gn/T/ipykernel_90688/2703781437.py:20: UserWarning: Geometry is in a geographic CRS. Results from 'area' are likely incorrect. Use 'GeoSeries.to_crs()' to re-project geometries to a projected CRS before this operation.
a = munis[munis["municipio"] == nombre].to_crs(epsg).geometry.area.iloc[0]
El caso de Web Mercator merece atención: da un área en metros, así que parece una respuesta válida. Pero infla sistemáticamente con la latitud, y ese error crece hacia los polos. Es la proyección de los mapas web, y por eso tantos análisis publicados comparan áreas mal.
4. Formatos
Los formatos no son un detalle administrativo: cada uno decide qué se puede guardar, qué se pierde y qué tan rápido se lee. Escribimos la misma capa de municipios en tres formatos vectoriales y comparamos qué sobrevive al viaje.
Código
import tempfile, osTMP = Path(tempfile.mkdtemp(prefix="formatos_"))prueba = munis[["municipio", "codigo", "geometry"]].copy()prueba["indice_de_vegetacion_promedio"] = np.random.rand(len(prueba)) # nombre largo a propositosalidas = {"GeoJSON (.geojson)": (TMP /"munis.geojson", "GeoJSON"),"GeoPackage (.gpkg)": (TMP /"munis.gpkg", "GPKG"),"Shapefile (.shp)": (TMP /"munis.shp", "ESRI Shapefile"),}filas = []for etiqueta, (ruta, driver) in salidas.items(): prueba.to_file(ruta, driver=driver)# Un shapefile son VARIOS archivos: hay que sumarlos todos piezas =list(ruta.parent.glob(ruta.stem +".*")) peso =sum(p.stat().st_size for p in piezas) /1024 leido = gpd.read_file(ruta) filas.append({"formato": etiqueta,"archivos": len(piezas),"KB": round(peso, 1),"CRS conservado": leido.crs == prueba.crs,"nombre de campo": [c for c in leido.columns if c.startswith("indice")][0], })pd.DataFrame(filas)
/var/folders/97/04drgvm13z7bxfszg04r31f40000gn/T/ipykernel_90688/2441975763.py:16: UserWarning: Column names longer than 10 characters will be truncated when saved to ESRI Shapefile.
prueba.to_file(ruta, driver=driver)
/Users/menene/code/data/05_datos_geoespaciales/.venv/lib/python3.13/site-packages/pyogrio/raw.py:733: RuntimeWarning: Normalized/laundered field name: 'indice_de_vegetacion_promedio' to 'indice_de_'
ogr_write(
formato
archivos
KB
CRS conservado
nombre de campo
0
GeoJSON (.geojson)
1
680.5
True
indice_de_vegetacion_promedio
1
GeoPackage (.gpkg)
2
1016.5
True
indice_de_vegetacion_promedio
2
Shapefile (.shp)
7
1234.0
True
indice_de_
Tres cosas que se ven en la tabla y que muerden en la práctica:
El shapefile no es un archivo, son varios (.shp, .shx, .dbf, .prj) que deben viajar juntos. Si alguien te manda solo el .shp, no tienes nada.
El shapefile trunca los nombres de campo a 10 caracteres. Ese indice_de_vegetacion_promedio llegó mutilado al otro lado.
GeoPackage guarda todo en un solo archivo SQLite, sin límite de nombres, y admite varias capas —vectoriales y ráster— dentro del mismo .gpkg. Es el reemplazo moderno recomendado.
GeoJSON es cómodo para la web y para APIs, pero pesa mucho más (es texto) y por especificación se asume en EPSG:4326: no lo uses como formato de trabajo en un CRS métrico.
Falta uno que te vas a encontrar seguro: KML / KMZ, el formato de Google Earth. Está orientado a visualización, no a análisis —lleva estilos, iconos y globos de información— y el KMZ es simplemente un KML comprimido. Sirve para compartir un resultado con alguien que no usa un SIG; no para trabajar.
Código
# WKT: la geometria como texto. Es asi como viaja dentro de una base de datos# (PostGIS, SpatiaLite) y como se lee en una consulta SQL.from shapely import wktg = munis.geometry.iloc[0]texto = g.wktprint("WKT (primeros 120 caracteres):")print(" ", texto[:120], "...")print(f"\nlongitud del WKT : {len(texto):,} caracteres")print(f"longitud del WKB : {len(g.wkb):,} bytes <- la version binaria, ~mitad de tamano")print("ida y vuelta exacta:", wkt.loads(texto).equals(g))# Y el CRS tambien es texto estandarizadoprint("\nCRS en WKT (recortado):")print(" ", munis.crs.to_wkt()[:150], "...")
WKT (primeros 120 caracteres):
POLYGON ((719733.1249511703 1614665.750049108, 719635.1249280736 1614724.2501063899, 719540.0625411432 1614803.125137798 ...
longitud del WKT : 8,784 caracteres
longitud del WKB : 3,741 bytes <- la version binaria, ~mitad de tamano
ida y vuelta exacta: True
CRS en WKT (recortado):
PROJCRS["WGS 84 / UTM zone 15N",BASEGEOGCRS["WGS 84",ENSEMBLE["World Geodetic System 1984 ensemble",MEMBER["World Geodetic System 1984 (Transit)"],MEM ...
Código
# El lado raster: un GeoTIFF lleva su georreferenciacion adentro.with rasterio.open(rutas["escena"]) as src: perfil = src.profileprint("Perfil del GeoTIFF:")for clave in ["driver", "dtype", "count", "width", "height", "crs", "compress"]:print(f" {clave:10}: {perfil.get(clave)}")print(f" {'transform':10}: {tuple(round(v, 2) for v in src.transform[:6])}")print(f"\n ventanas internas (tiles): {src.is_tiled}")print(f" vistas piramidales : {src.overviews(1) or'ninguna'}")
Un COG (Cloud Optimized GeoTIFF) es este mismo archivo, pero organizado en ventanas internas y con vistas piramidales, de modo que un cliente pueda pedir por HTTP solo el pedazo que necesita sin bajar el archivo completo. Es lo que hace posible que un catálogo STAC sirva petabytes de imágenes sin que nadie los descargue: el estándar de búsqueda (STAC) y el de lectura parcial (COG) trabajando juntos.
Para cubos de datos con más de dos dimensiones —x, y, tiempo, variable— el GeoTIFF se queda corto y se usan NetCDF, HDF o Zarr. Y para que dos sistemas distintos se entiendan sin compartir archivos están los servicios del OGC: WMS entrega el mapa ya dibujado, WFS entrega las geometrías, WCS entrega el ráster crudo.
5. El modelo ráster
Discreto contra continuo
La decisión no es de gusto, es del fenómeno:
Fenómenos discretos: objetos con bordes definidos —un municipio, una carretera, un poste—. Se representan bien como vectores.
Fenómenos continuos: magnitudes que existen en todo punto del espacio —elevación, temperatura, precipitación, reflectancia—. Se representan bien como rásteres.
Elegir mal es una fuente frecuente de errores: forzar un fenómeno continuo a polígonos inventa fronteras que no existen; forzar uno discreto a celdas pierde el borde exacto.
La cuadrícula
Un ráster representa el espacio como una cuadrícula regular de celdas, cada una con un valor:
Es esencialmente una matriz, igual que una imagen digital.
La resolución espacial es el tamaño de la celda en el terreno: 10 m, 30 m, 250 m.
Cada banda es una matriz distinta sobre la misma cuadrícula.
Vectorial
Ráster
Bueno para
fronteras, redes, ubicaciones puntuales, datos censales
En la práctica se combinan constantemente: recortar un ráster de NDVI con el polígono de un municipio y sacar el promedio es una operación cotidiana — y es exactamente lo que haremos en la sección 10.
Los sensores remotos
Un sensor remoto es cualquier instrumento que mide algo sin estar en contacto con el objeto medido. En la práctica: observar la Tierra desde el aire o el espacio.
La cadena es siempre la misma:
una fuente de energía ilumina la superficie → la superficie refleja o emite radiación → el sensor la capta → el resultado se convierte en una matriz de números georreferenciada.
Ese último paso es el que nos interesa: lo que llega a tus manos no es una foto, es una matriz con coordenadas.
Ahora el otro modelo del mundo. Un ráster es una matriz georreferenciada: números en una cuadrícula, más los metadatos que la anclan al planeta.
Lo que convierte una matriz cualquiera en un dato geoespacial es la transformación afín (dónde está la esquina y cuánto mide cada celda) y el CRS. Sin eso, es solo una imagen.
Código
with rasterio.open(rutas["escena"]) as src: escena = src.read().astype("float32") # (bandas, alto, ancho) transform, crs_escena = src.transform, src.crs limites = src.bounds descripciones = src.descriptionswith rasterio.open(rutas["verdad"]) as src: verdad = src.read(1)EXT = (limites.left, limites.right, limites.bottom, limites.top)print(f"forma : {escena.shape} (bandas, alto, ancho)")print(f"tipo : {escena.dtype}")print(f"CRS : {crs_escena.to_string()}")print(f"tamaño de celda: {transform.a:.0f} m x {abs(transform.e):.0f} m")print(f"esquina sup-izq: ({transform.c:,.0f}, {transform.f:,.0f}) m UTM")print(f"bandas : {', '.join(descripciones)}")print(f"\nEs un tensor de {escena.size:,} números "f"({escena.nbytes /1e6:.0f} MB en memoria).")
forma : (6, 1024, 1024) (bandas, alto, ancho)
tipo : float32
CRS : EPSG:32615
tamaño de celda: 60 m x 60 m
esquina sup-izq: (700,000, 1,640,000) m UTM
bandas : B02_azul, B03_verde, B04_rojo, B08_nir, B11_swir1, B12_swir2
Es un tensor de 6,291,456 números (25 MB en memoria).
Color natural y falso color
Un sensor multiespectral ve más allá del ojo humano. Con las mismas tres posiciones de color (rojo, verde, azul de la pantalla) podemos mostrar distintas bandas y revelar cosas distintas:
Color natural (R=rojo, G=verde, B=azul): el paisaje como lo veríamos.
Falso color (R=NIR, G=rojo, B=verde): la vegetación sana refleja muchísimo infrarrojo cercano, así que aparece en rojo intenso. Es la composición clásica de agricultura y silvicultura.
Código
def estirar(banda, p=2):'''Estiramiento por percentiles: lo que hace cualquier visor de imagen satelital.''' lo, hi = np.percentile(banda, (p, 100- p))return np.clip((banda - lo) / (hi - lo +1e-12), 0, 1)AZUL, VERDE, ROJO, NIR, SWIR1, SWIR2 =range(6)rgb = np.dstack([estirar(escena[ROJO]), estirar(escena[VERDE]), estirar(escena[AZUL])])falso = np.dstack([estirar(escena[NIR]), estirar(escena[ROJO]), estirar(escena[VERDE])])fig, axes = plt.subplots(1, 3, figsize=(14.5, 5))axes[0].imshow(rgb, extent=EXT); axes[0].set_title("Color natural (R,G,B)")axes[1].imshow(falso, extent=EXT); axes[1].set_title("Falso color (NIR,R,G)\nla vegetación se vuelve roja")im = axes[2].imshow(verdad, cmap=CMAP_CLASES, extent=EXT, vmin=0, vmax=4)axes[2].set_title("Cobertura real (la verdad que simulamos)")for ax in axes: munis.boundary.plot(ax=ax, color="white", linewidth=0.45, alpha=0.75) ax.set_xlim(EXT[0], EXT[1]); ax.set_ylim(EXT[2], EXT[3]) ax.set_xticks([]); ax.set_yticks([]); ax.grid(False)cbar = fig.colorbar(im, ax=axes[2], fraction=0.046, ticks=range(5))cbar.ax.set_yticklabels(CLASES, fontsize=7)plt.tight_layout(); plt.show()
Se distinguen los elementos del paisaje: un lago, un río que baja hacia él, los núcleos urbanos y el contraste entre bosque en las zonas altas y cultivo en el valle. Las líneas blancas son los municipios reales.
Nota cómo en falso color el agua es casi negra (absorbe todo el NIR) mientras la vegetación se enciende. Eso no es un efecto estético: es el comportamiento físico de cada material, y es exactamente lo que vamos a medir en la sección siguiente.
Las seis bandas por separado
Cada banda es una matriz distinta sobre la misma cuadrícula. Mira cómo el agua cambia de brillante a negra al pasar del visible al infrarrojo.
Código
fig, axes = plt.subplots(2, 3, figsize=(13, 8))for i, ax inenumerate(axes.ravel()): ax.imshow(escena[i], cmap="gray", extent=EXT) ax.set_title(f"{BANDAS[i]} [{escena[i].min():.2f}, {escena[i].max():.2f}]") ax.set_xticks([]); ax.set_yticks([]); ax.grid(False)fig.suptitle("Las 6 bandas de la escena — misma cuadrícula, información distinta", y=1.0)plt.tight_layout(); plt.show()
6. Firmas espectrales
El espectro electromagnético
Los sensores no se limitan a lo que ve el ojo humano. Cada región del espectro revela información distinta, y esa es la razón de ser de la teledetección:
Región
Qué revela
Visible (azul, verde, rojo)
lo que vemos; el rojo sirve para analizar vegetación
Infrarrojo cercano (NIR)
la vegetación sana lo refleja intensamente — la banda estrella para agricultura
Infrarrojo de onda corta (SWIR)
humedad del suelo y la vegetación; incendios y minerales
Infrarrojo térmico
temperatura de superficie; islas de calor urbanas, estrés hídrico
Microondas
las usa el radar; atraviesan nubes y lluvia
Nuestra escena tiene 6 bandas: azul, verde, rojo, NIR, SWIR1 y SWIR2 — la configuración típica de Sentinel-2 y Landsat.
La huella digital de cada material
Cada material refleja una proporción distinta de energía en cada longitud de onda. Ese patrón es su firma espectral, y funciona como una huella digital:
Vegetación sana: refleja poco en el rojo (la clorofila lo absorbe) y muchísimo en el NIR.
Vegetación estresada o seca: sube la reflectancia en el rojo y baja en el NIR. Por eso el estrés se detecta antes de que la planta se vea amarilla.
Agua: absorbe casi todo el NIR, por eso se ve muy oscura en esa banda.
Suelo desnudo: reflectancia moderada y creciente hacia el infrarrojo.
Clasificar cobertura terrestre consiste, en el fondo, en distinguir firmas espectrales. Todo lo que viene después —índices, Random Forest, deep learning— son formas cada vez más sofisticadas de hacer exactamente eso.
Aquí la extraemos de la escena: para cada clase de cobertura, la reflectancia media en cada banda, con su rango intercuartílico.
Código
LONGITUDES = [490, 560, 665, 842, 1610, 2190] # nm, centros tipo Sentinel-2fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(13, 4.6))for k, (clase, color) inenumerate(zip(CLASES, COLORES)): mascara = verdad == k muestras = escena[:, mascara][:, ::37] # submuestreo para ir rápido media = muestras.mean(axis=1) q1, q3 = np.percentile(muestras, (25, 75), axis=1) ax1.plot(LONGITUDES, media, "o-", color=color, label=clase, linewidth=2) ax1.fill_between(LONGITUDES, q1, q3, color=color, alpha=0.15)ax1.set_xlabel("Longitud de onda (nm)"); ax1.set_ylabel("Reflectancia")ax1.set_title("Firmas espectrales por cobertura\n(línea = media, banda = rango intercuartílico)")ax1.legend(fontsize=8)for x, etiqueta inzip(LONGITUDES, ["azul", "verde", "rojo", "NIR", "SWIR1", "SWIR2"]): ax1.axvline(x, color="gray", linewidth=0.4, alpha=0.4) ax1.text(x, ax1.get_ylim()[1] *0.97, etiqueta, rotation=90, fontsize=6, ha="right", va="top", color="gray")# El espacio rojo-NIR: donde vive el NDVIfor k, (clase, color) inenumerate(zip(CLASES, COLORES)): mascara = verdad == k ax2.scatter(escena[ROJO][mascara][::400], escena[NIR][mascara][::400], s=3, alpha=0.35, color=color, label=clase)lim = np.linspace(0, 0.35, 10)ax2.plot(lim, lim, "k--", linewidth=0.8, alpha=0.6, label="NDVI = 0")ax2.set_xlabel("Rojo"); ax2.set_ylabel("NIR")ax2.set_title("El espacio rojo–NIR\ncuanto más arriba de la diagonal, mayor NDVI")ax2.legend(fontsize=7, markerscale=3)plt.tight_layout(); plt.show()
Léelo con cuidado, porque explica todo lo que viene después:
El agua cae en picada después del verde: absorbe casi todo el infrarrojo.
El bosque y el cultivo tienen el salto característico de la vegetación: mínimo en el rojo (la clorofila lo absorbe) y máximo en el NIR.
Lo urbano y el suelo desnudo son planos y crecientes hacia el SWIR.
Fíjate también en dónde se solapan las bandas de dispersión: bosque con cultivo en el NIR, urbano con suelo desnudo en el SWIR. Esos dos pares son los que el clasificador va a confundir más adelante, y no es casualidad: distinguir cultivo de bosque, o asfalto de suelo seco, también es difícil en teledetección real.
7. Las cuatro resoluciones
Toda imagen satelital se describe con cuatro resoluciones, y siempre hay que negociar entre ellas. Aquí las degradamos una por una sobre la misma escena para ver qué se pierde exactamente en cada caso.
Resolución
Qué mide
Landsat
Sentinel-2
MODIS
Espacial
tamaño del píxel en el terreno
30 m
10–20 m
250 m
Espectral
cuántas bandas y qué tan angostas
11
13
36
Temporal
cada cuánto vuelve a pasar
16 días
5 días
diaria
Radiométrica
niveles de intensidad que distingue
12 bits
12 bits
12 bits
Código
def degradar_espacial(arr, factor):'''Promedia bloques de factor x factor: es lo que hace un pixel mas grande.''' n = arr.shape[-1] // factor * factorif arr.ndim ==2:return arr[:n, :n].reshape(n // factor, factor, -1, factor).mean(axis=(1, 3))return arr[:, :n, :n].reshape(arr.shape[0], n // factor, factor, -1, factor).mean(axis=(2, 4))factores = [(1, "60 m — la escena original"), (2, "120 m"), (4, "240 m — como MODIS"), (8, "480 m")]fig, axes = plt.subplots(1, 4, figsize=(16, 4.3))for ax, (f, etiqueta) inzip(axes, factores): chico = degradar_espacial(escena, f) ax.imshow(np.dstack([estirar(chico[NIR]), estirar(chico[ROJO]), estirar(chico[VERDE])])) ax.set_title(f"{etiqueta}\n{chico.shape[1]}x{chico.shape[2]} px", fontsize=9) ax.set_xticks([]); ax.set_yticks([]); ax.grid(False)fig.suptitle("Resolución espacial: el mismo territorio, píxeles cada vez más grandes", y=1.02)plt.tight_layout(); plt.show()print("Qué le pasa a la clase URBANA al agrandar el píxel:")for f, etiqueta in factores: urb = degradar_espacial((verdad ==3).astype(float), f)print(f" {etiqueta:26} celdas puras urbanas: {(urb >0.9).mean():6.2%}"f" celdas mezcladas: {((urb >0) & (urb <0.9)).mean():6.2%}")
Qué le pasa a la clase URBANA al agrandar el píxel:
60 m — la escena original celdas puras urbanas: 8.04% celdas mezcladas: 0.00%
120 m celdas puras urbanas: 7.52% celdas mezcladas: 1.06%
240 m — como MODIS celdas puras urbanas: 6.86% celdas mezcladas: 2.87%
480 m celdas puras urbanas: 5.76% celdas mezcladas: 6.51%
Ese es el píxel mixto, y es el problema central de la resolución gruesa: a 480 m casi ninguna celda es urbana pura, todas son una mezcla de techo, calle, árbol y patio. El sensor no te da la clase, te da un promedio ponderado de varias. Por eso a resolución gruesa se trabaja con análisis de mezclas espectrales (qué fracción de cada celda es cada cobertura) en vez de clasificación dura.
A 12 y 8 bits la imagen se ve igual, pero el error que introduce la cuantización en el NDVI crece rápido al bajar. Con 2 bits el índice ya no distingue vegetación vigorosa de vegetación mediana: la información no está perdida en la pantalla, está perdida en el número.
Esto importa sobre todo en superficies oscuras —agua, sombra, bosque denso— donde toda la señal vive en los primeros niveles de la escala.
Código
# Resolucion espectral: cuanta informacion aporta cada banda adicional.# Lo medimos donde importa - en la capacidad de separar coberturas.from sklearn.ensemble import RandomForestClassifierfrom sklearn.model_selection import cross_val_scoreobs_previa = gpd.read_file(rutas["observaciones"])f_o, c_o = obs_previa["fila"].values, obs_previa["col"].valuesespectro = escena[:, f_o, c_o].Ty_o = obs_previa["clase_id"].valuesconjuntos = [ ("solo RGB (3 bandas)", [AZUL, VERDE, ROJO]), ("RGB + NIR (4)", [AZUL, VERDE, ROJO, NIR]), ("RGB + NIR + SWIR1 (5)", [AZUL, VERDE, ROJO, NIR, SWIR1]), ("las 6 bandas", [AZUL, VERDE, ROJO, NIR, SWIR1, SWIR2]),]print("Exactitud de clasificación según cuántas bandas ve el sensor")print("(validación cruzada de 5 pliegues, Random Forest)\n")for etiqueta, cols_b in conjuntos: rf = RandomForestClassifier(n_estimators=120, random_state=42, n_jobs=-1) s = cross_val_score(rf, espectro[:, cols_b], y_o, cv=5).mean()print(f" {etiqueta:28}{s:.3f}")
Exactitud de clasificación según cuántas bandas ve el sensor
(validación cruzada de 5 pliegues, Random Forest)
solo RGB (3 bandas) 0.717
RGB + NIR (4) 0.815
RGB + NIR + SWIR1 (5) 0.868
las 6 bandas 0.882
El salto grande está en añadir el NIR: es la banda que separa vegetación de todo lo demás, y por eso es la banda estrella de la teledetección. El SWIR aporta después sobre todo para distinguir urbano de suelo desnudo, que es exactamente el par que se solapaba en las firmas.
Un sensor hiperespectral lleva esto al extremo: cientos de bandas angostas, con las que se identifican minerales o especies vegetales concretas. El precio es el volumen de datos y la fuerte correlación entre bandas vecinas.
Código
# Resolucion temporal: por que la revisita frecuente vale tanto en el tropico.# Si cada pasada tiene probabilidad p de estar nublada, la probabilidad de# conseguir AL MENOS una imagen limpia en k pasadas es 1 - p^k.p_nube =0.70# tipico de la temporada lluviosa en Guatemaladias =90# una temporada de cultivofig, ax = plt.subplots(figsize=(6.5, 3.8))for revisita, nombre in [(16, "Landsat (16 d)"), (5, "Sentinel-2 (5 d)"), (1, "MODIS (diaria)")]: k = np.arange(0, dias // revisita +1) ax.plot(k * revisita, 1- p_nube ** k, "o-", label=nombre, linewidth=2)ax.axhline(0.95, color="gray", linestyle="--", linewidth=0.8)ax.set_xlabel("días transcurridos"); ax.set_ylabel("P(al menos una imagen sin nubes)")ax.set_title(f"Resolución temporal con {p_nube:.0%} de probabilidad de nube por pasada")ax.legend(fontsize=8); plt.tight_layout(); plt.show()for revisita, nombre in [(16, "Landsat"), (5, "Sentinel-2"), (1, "MODIS")]: k =int(np.ceil(np.log(0.05) / np.log(p_nube)))print(f" {nombre:12} necesita {k} pasadas = {k * revisita:3d} días para un 95% de éxito")
Landsat necesita 9 pasadas = 144 días para un 95% de éxito
Sentinel-2 necesita 9 pasadas = 45 días para un 95% de éxito
MODIS necesita 9 pasadas = 9 días para un 95% de éxito
El compromiso es físico, no de presupuesto: más resolución espacial implica menos área por escena y, en general, menor frecuencia de revisita. Un satélite que ve 30 cm no puede además pasar todos los días por todas partes.
Y en el trópico la resolución temporal deja de ser un lujo: con nubosidad persistente, un sensor de 16 días puede pasarse meses sin una imagen limpia de una parcela. Ahí es donde el radar —que atraviesa nubes— deja de ser una alternativa exótica y pasa a ser la única opción.
Las misiones que vas a usar de verdad
Misión
Resolución
Revisita
Costo
Para qué
Landsat (NASA/USGS)
30 m
16 días
gratis
serie continua desde 1972 — insustituible para cambio a largo plazo
Sentinel-2 (ESA)
10–20 m, 13 bandas
5 días
gratis
el caballo de batalla actual
Sentinel-1 (ESA)
radar SAR
6–12 días
gratis
ve a través de las nubes: inundaciones y trópico
MODIS / VIIRS
250–1000 m
diaria
gratis
clima, incendios, monitoreo global
PlanetScope, WorldView
metros o submétrica
casi diaria
comercial
detalle fino, cuando hay presupuesto
Para un proyecto de curso: Sentinel-2 casi siempre. Es gratis, tiene buena resolución espacial y espectral, y 5 días de revisita.
Plataformas y tipos de sensor
Todo lo anterior fue un sensor pasivo en una plataforma satelital. Es un caso entre varios, y el resto de las combinaciones resuelven problemas que este no puede:
Plataforma
Resolución
Cobertura
Cuándo vuelas
Satélite
metros a cientos de metros
global
cuando pasa
Avión
decenas de centímetros
regional
bajo demanda, caro
Dron
centímetros
pocas hectáreas
cuando quieras
Para un proyecto local —una finca, un deslizamiento, un sitio arqueológico— el dron gana por goleada: resolución de centímetros y control total de la fecha. El satélite gana en cualquier cosa que requiera serie histórica o cobertura amplia.
La otra división es más profunda:
Pasivos: miden energía que no producen ellos, así que dependen del Sol. No ven de noche —salvo los térmicos, que miden la emisión de la superficie— y las nubes los bloquean. Todo lo que hemos hecho hasta aquí es pasivo.
Activos: emiten su propia energía y miden el retorno. Funcionan de noche y, según la longitud de onda, atraviesan nubes. LiDAR y radar son activos.
LiDAR
LiDAR emite pulsos láser y mide cuánto tardan en volver. El resultado no es una imagen sino una nube de puntos en tres dimensiones, y de ahí salen dos superficies distintas que conviene no confundir:
DSM (modelo digital de superficie): lo primero que devuelve el pulso —copas de árboles, techos—.
DTM (modelo digital del terreno): el suelo desnudo, tras filtrar todo lo que está encima.
Restar una de otra da la altura de la vegetación o de las construcciones. Es una operación local de álgebra de mapas, y es la base de media teledetección forestal.
Código
from scipy.ndimage import gaussian_filter# Simulamos un levantamiento LiDAR sobre un recorte de nuestra escena.# La geografia es la real de la escena: el bosque esta donde esta el bosque.rng_l = np.random.default_rng(99)REC =slice(300, 556), slice(300, 556) # 256 x 256 px = 15 x 15 kmcobertura_rec = verdad[REC]# Terreno: relieve suave (esto es el DTM)dtm =1400+260* gaussian_filter(rng_l.standard_normal(cobertura_rec.shape), sigma=22)# Altura de lo que hay ENCIMA del terreno, segun la coberturaaltura = np.zeros_like(dtm)altura[cobertura_rec ==1] = rng_l.normal(22, 5, (cobertura_rec ==1).sum()) # bosquealtura[cobertura_rec ==2] = rng_l.normal(1.5, 0.6, (cobertura_rec ==2).sum()) # cultivoaltura[cobertura_rec ==3] = rng_l.normal(7, 4, (cobertura_rec ==3).sum()) # edificiosaltura = np.clip(gaussian_filter(altura, sigma=0.8), 0, None)dsm = dtm + altura # lo primero que toca el pulsochm = dsm - dtm # modelo de altura del doselfig, axes = plt.subplots(1, 4, figsize=(17, 4.3))im = axes[0].imshow(dtm, cmap="terrain"); axes[0].set_title("DTM — el suelo desnudo")fig.colorbar(im, ax=axes[0], fraction=0.046, label="m s.n.m.")im = axes[1].imshow(dsm, cmap="terrain"); axes[1].set_title("DSM — la primera superficie")fig.colorbar(im, ax=axes[1], fraction=0.046, label="m s.n.m.")im = axes[2].imshow(chm, cmap="YlGn", vmin=0, vmax=35)axes[2].set_title("DSM − DTM = altura del dosel")fig.colorbar(im, ax=axes[2], fraction=0.046, label="m")for k, nombre in [(1, "bosque"), (3, "urbano"), (2, "cultivo")]: axes[3].hist(chm[cobertura_rec == k], bins=45, alpha=0.6, label=nombre, density=True)axes[3].set_xlabel("altura sobre el terreno (m)"); axes[3].set_yticks([])axes[3].set_title("La resta separa coberturas\nque el color no distingue")axes[3].legend(fontsize=8)for ax in axes[:3]: ax.set_xticks([]); ax.set_yticks([]); ax.grid(False)plt.tight_layout(); plt.show()print(f"Altura media del bosque : {chm[cobertura_rec ==1].mean():5.1f} m")print(f"Altura media de lo urbano: {chm[cobertura_rec ==3].mean():5.1f} m")print("\nDos coberturas pueden tener la misma firma espectral y una altura")print("completamente distinta. Por eso LiDAR y óptico se combinan tan bien:")print("miden propiedades independientes del mismo lugar.")
Altura media del bosque : 21.8 m
Altura media de lo urbano: 7.0 m
Dos coberturas pueden tener la misma firma espectral y una altura
completamente distinta. Por eso LiDAR y óptico se combinan tan bien:
miden propiedades independientes del mismo lugar.
Ese último histograma es el argumento entero a favor del LiDAR: la estructura vertical es información que ninguna banda espectral contiene. Es lo que permite estimar biomasa, distinguir bosque primario de secundario, modelar inundaciones con el terreno real —y, quitando digitalmente el dosel, encontrar lo que hay debajo: así aparecieron miles de estructuras mayas bajo la selva del Petén.
Radar y SAR
El radar de apertura sintética emite microondas y mide su retorno:
Atraviesa nubes, humo y lluvia, y funciona de noche. En una región con nubosidad permanente es la única fuente confiable.
Mide rugosidad y humedad, no color. Una imagen SAR no se interpreta como una foto: el agua en calma se ve negra (refleja el pulso lejos del sensor) y una ciudad se ve brillantísima.
Por eso es el sensor de referencia para mapear inundaciones: la lámina de agua aparece como una mancha oscura inconfundible, justo cuando el cielo está cubierto y ningún sensor óptico ve nada.
La interferometría (InSAR) compara la fase entre dos pasadas y detecta deformaciones del terreno de milímetros: hundimientos por extracción de agua, deslizamientos lentos, inflado de un edificio volcánico antes de una erupción. En Guatemala se usa para vigilar Fuego y Pacaya.
8. Preprocesamiento
Una imagen recién bajada no está lista para analizarse. Igual que en cualquier proyecto de ciencia de datos, esta etapa se come la mayor parte del tiempo. La cadena completa:
Paso
Qué hace
¿Lo demostramos?
Corrección radiométrica
convierte los números crudos del sensor a radiancia y luego a reflectancia física, para que dos fechas sean comparables
viene hecha: la escena ya está en reflectancia [0, 1]
Corrección atmosférica
descuenta el efecto de la atmósfera — es la diferencia entre un producto nivel 1 (tope de la atmósfera) y nivel 2 (superficie)
sí, abajo
Corrección geométrica y ortorrectificación
alinea la imagen con el terreno real, corrigiendo la distorsión del relieve
no: los productos modernos ya vienen ortorrectificados
Enmascarado de nubes
descarta píxeles con nube y sombra usando la banda de calidad
sí, abajo
Composición temporal
combina varias fechas (p. ej. la mediana) para obtener una imagen sin nubes
sí, abajo
Mosaico y recorte
une escenas vecinas y recorta al área de interés
recorte en la sección 10
Remuestreo
lleva capas de distinta resolución a una cuadrícula común
sí, abajo
Un detalle sobre la corrección radiométrica que conviene entender aunque no la hagamos: el sensor no mide reflectancia, mide números digitales (DN) — cuentas de un conversor analógico-digital. Sin convertirlos a una magnitud física, dos imágenes del mismo lugar en fechas distintas no son comparables, porque la iluminación solar cambió. Es el equivalente a comparar precios sin ajustar por inflación.
Para verlo con datos, simulamos cuatro fechas de la misma zona, cada una con su propia nubosidad, su bruma atmosférica y su momento fenológico —los cultivos verdean y luego se cosechan—. Es exactamente lo que encontrarías al bajar una serie de Sentinel-2.
Código
from scipy.ndimage import gaussian_filter, shift as desplazar# La bruma es aditiva y se dispersa mas en las longitudes de onda cortas# (dispersion de Rayleigh): pega fuerte en el azul y casi nada en el SWIR.PESOS_BRUMA = np.array([1.00, 0.78, 0.58, 0.32, 0.14, 0.09])def simular_fecha(semilla, fraccion_nube, vigor, bruma):'''Una adquisicion: fenologia + bruma atmosferica + nubes con su sombra.''' rng = np.random.default_rng(semilla) img = escena.copy()# Fenologia: la vegetacion sube o baja su NIR y su rojo segun el momento veg = np.isin(verdad, [1, 2]) img[NIR] = np.where(veg, img[NIR] * vigor, img[NIR]) img[ROJO] = np.where(veg, img[ROJO] * (2- vigor), img[ROJO])# Bruma: aditiva, mas fuerte en longitudes de onda cortas (Rayleigh) img += bruma * PESOS_BRUMA[:, None, None]# Nubes: campo suave umbralizado; la sombra es la nube desplazada campo = gaussian_filter(rng.standard_normal(verdad.shape), sigma=28) campo = (campo - campo.mean()) / campo.std() nube = campo > np.quantile(campo, 1- fraccion_nube) sombra = desplazar(nube.astype(float), (40, 28), order=0) >0.5 sombra &=~nube img = np.where(nube[None], 0.72+0.05* rng.random(img.shape), img) # nube: blanca y brillante img = np.where(sombra[None], img *0.42, img) # sombra: todo mas oscuroreturn np.clip(img, 0, 1).astype("float32"), nube, sombrafechas = [ ("2026-01-12", 12345, 0.10, 0.85, 0.010), # seca, poca nube, vegetacion baja ("2026-03-08", 23456, 0.35, 1.05, 0.055), # bruma fuerte ("2026-06-20", 34567, 0.55, 1.18, 0.030), # lluviosa: mucha nube, vegetacion en su pico ("2026-09-15", 45678, 0.30, 1.10, 0.020),]serie, mascaras = [], []for etiqueta, semilla, fnube, vigor, bruma in fechas: img, nube, sombra = simular_fecha(semilla, fnube, vigor, bruma) serie.append(img); mascaras.append(nube | sombra)serie = np.stack(serie) # (fecha, banda, alto, ancho)mascaras = np.stack(mascaras) # (fecha, alto, ancho)print(f"Serie temporal: {serie.shape} (fechas, bandas, alto, ancho)")for (etiqueta, *_), m inzip(fechas, mascaras):print(f" {etiqueta}: {m.mean():5.1%} de píxeles inservibles (nube o sombra)")print(f"\nPíxeles con las 4 fechas nubladas: {(mascaras.all(axis=0)).mean():.2%}")
Serie temporal: (4, 6, 1024, 1024) (fechas, bandas, alto, ancho)
2026-01-12: 16.8% de píxeles inservibles (nube o sombra)
2026-03-08: 49.8% de píxeles inservibles (nube o sombra)
2026-06-20: 70.9% de píxeles inservibles (nube o sombra)
2026-09-15: 43.1% de píxeles inservibles (nube o sombra)
Píxeles con las 4 fechas nubladas: 3.38%
Código
fig, axes = plt.subplots(2, 4, figsize=(16, 8))for j, ((etiqueta, *_), img, m) inenumerate(zip(fechas, serie, mascaras)): axes[0, j].imshow(np.dstack([estirar(img[ROJO]), estirar(img[VERDE]), estirar(img[AZUL])])) axes[0, j].set_title(f"{etiqueta}\ncolor natural", fontsize=9) axes[1, j].imshow(m, cmap="gray_r") axes[1, j].set_title(f"banda de calidad (QA)\n{m.mean():.0%} enmascarado", fontsize=9)for ax in axes.ravel(): ax.set_xticks([]); ax.set_yticks([]); ax.grid(False)fig.suptitle("Cuatro adquisiciones de la misma zona y su máscara de nube + sombra", y=1.0)plt.tight_layout(); plt.show()
Composición temporal
Ninguna fecha sirve por sí sola. La receta estándar es enmascarar y componer: se descartan los píxeles con nube y sombra usando la banda de calidad, y se toma la mediana de lo que queda en cada píxel. La mediana es la clave: es robusta, así que una nube que se coló en la máscara no arrastra el resultado como sí lo haría un promedio.
Código
import warningsserie_enmascarada = np.where(mascaras[:, None, :, :], np.nan, serie)with warnings.catch_warnings(): # los pixeles nublados en las 4 fechas dan NaN warnings.simplefilter("ignore", RuntimeWarning) compuesto = np.nanmedian(serie_enmascarada, axis=0)# Esos pixeles quedan SIN DATO: hay que declararlos, no rellenarlos en silenciosin_dato = np.isnan(compuesto[0])compuesto = np.where(np.isnan(compuesto), np.nanmedian(compuesto), compuesto)def ndvi_de(img):return (img[NIR] - img[ROJO]) / (img[NIR] + img[ROJO] +1e-10)fig, axes = plt.subplots(1, 4, figsize=(16, 4.4))axes[0].imshow(np.dstack([estirar(serie[2][ROJO]), estirar(serie[2][VERDE]), estirar(serie[2][AZUL])]))axes[0].set_title("Una sola fecha (la más nublada)")axes[1].imshow(np.dstack([estirar(compuesto[ROJO]), estirar(compuesto[VERDE]), estirar(compuesto[AZUL])]))axes[1].set_title("Compuesto de mediana\nlas nubes desaparecieron")im = axes[2].imshow(ndvi_de(serie[2]), cmap="RdYlGn", vmin=-1, vmax=1)axes[2].set_title("NDVI de la fecha nublada"); fig.colorbar(im, ax=axes[2], fraction=0.046)im = axes[3].imshow(ndvi_de(compuesto), cmap="RdYlGn", vmin=-1, vmax=1)axes[3].set_title("NDVI del compuesto"); fig.colorbar(im, ax=axes[3], fraction=0.046)for ax in axes: ax.set_xticks([]); ax.set_yticks([]); ax.grid(False)plt.tight_layout(); plt.show()print(f"Píxeles sin ninguna observación válida: {sin_dato.mean():.2%} "f"(hay que reportarlos, no rellenarlos en silencio)")
Píxeles sin ninguna observación válida: 3.38% (hay que reportarlos, no rellenarlos en silencio)
Corrección atmosférica
La diferencia entre un producto nivel 1 (reflectancia en el tope de la atmósfera) y nivel 2 (reflectancia de superficie) es que al segundo ya le descontaron el efecto de la atmósfera.
El método clásico y más simple es la sustracción del objeto oscuro: se busca en la imagen algo que debería ser casi negro —agua profunda— y se asume que todo lo que ese píxel refleja de más es bruma. Se resta banda por banda.
Como la escena es simulada, tenemos un lujo poco común: sabemos exactamente cuánta bruma le metimos, así que podemos calificar la corrección en vez de confiar en ella.
Código
brumosa = serie[1] # la fecha con mas bruma# El objeto oscuro es el percentil 1 de cada banda... pero SOLO sobre pixeles# validos. Si se incluyen las sombras de nube, lo mas oscuro de la imagen es# una sombra, no un objeto oscuro, y la estimacion se va al suelo.oscuro = np.percentile(brumosa[:, ~mascaras[1]], 1, axis=1)ingenuo = np.percentile(brumosa.reshape(6, -1), 1, axis=1) # sin enmascararcorregida = np.clip(brumosa - oscuro[:, None, None], 0, 1)# Ventaja de simular: sabemos exactamente cuanta bruma le metimos a esta fechabruma_real = fechas[1][4] * PESOS_BRUMAprint(f"{'banda':10}{'real':>8}{'estimada':>10}{'sin enmascarar':>16}")for b, nombre inenumerate(BANDAS):print(f"{nombre:10}{bruma_real[b]:8.4f}{oscuro[b]:10.4f}{ingenuo[b]:16.4f}")print(f"\nError medio de la estimación enmascarada : {np.abs(oscuro - bruma_real).mean():.4f}")print(f"Error medio sin enmascarar : {np.abs(ingenuo - bruma_real).mean():.4f}")print("\nDos cosas: (1) el metodo recupera la bruma casi exactamente, y (2) el")print("patron desciende del azul al SWIR. Eso es dispersion de Rayleigh, y es")print("la firma de que la correccion esta capturando atmosfera y no otra cosa.")fig, axes = plt.subplots(1, 3, figsize=(14, 4.4))axes[0].imshow(np.dstack([estirar(brumosa[ROJO]), estirar(brumosa[VERDE]), estirar(brumosa[AZUL])]))axes[0].set_title("Nivel 1: reflectancia en el tope de la atmósfera")axes[1].imshow(np.dstack([estirar(corregida[ROJO]), estirar(corregida[VERDE]), estirar(corregida[AZUL])]))axes[1].set_title("Nivel 2 aproximado: tras restar el objeto oscuro")agua = verdad ==0axes[2].plot(LONGITUDES, brumosa[:, agua].mean(axis=1), "o-", label="antes", linewidth=2)axes[2].plot(LONGITUDES, corregida[:, agua].mean(axis=1), "o-", label="después", linewidth=2)axes[2].set_xlabel("Longitud de onda (nm)"); axes[2].set_ylabel("Reflectancia del agua")axes[2].set_title("El agua debería ser casi negra\nen todas las bandas")axes[2].legend(fontsize=8); axes[2].grid(True, alpha=0.3)for ax in axes[:2]: ax.set_xticks([]); ax.set_yticks([]); ax.grid(False)plt.tight_layout(); plt.show()
banda real estimada sin enmascarar
B02_azul 0.0550 0.0550 0.0231
B03_verde 0.0429 0.0429 0.0250
B04_rojo 0.0319 0.0319 0.0145
B08_nir 0.0176 0.0179 0.0178
B11_swir1 0.0077 0.0103 0.0136
B12_swir2 0.0049 0.0058 0.0092
Error medio de la estimación enmascarada : 0.0006
Error medio sin enmascarar : 0.0129
Dos cosas: (1) el metodo recupera la bruma casi exactamente, y (2) el
patron desciende del azul al SWIR. Eso es dispersion de Rayleigh, y es
la firma de que la correccion esta capturando atmosfera y no otra cosa.
Compara las dos columnas de la tabla y verás por qué el orden de los pasos importa: si estimas la bruma sin enmascarar antes, lo más oscuro de la imagen no es agua, es una sombra de nube, y el método te devuelve una corrección que no tiene nada que ver con la atmósfera. Enmascarar primero, corregir después.
El método tiene un supuesto fuerte y conviene tenerlo presente: el objeto oscuro tiene que ser realmente negro. Si el agua de tu escena es somera o turbia, la corrección le resta a toda la imagen la reflectancia propia del agua y oscurece de más. Por eso los productos de nivel 2 reales usan modelos de transferencia radiativa (Sen2Cor, LaSRC) y no un percentil — pero la intuición es exactamente esta.
Remuestreo
Cuando se combinan capas de distinta resolución hay que llevarlas a una cuadrícula común, y el método de remuestreo depende del tipo de dato:
Vecino más cercano para datos categóricos. Es el único que no inventa clases.
Bilineal o cúbico para datos continuos, donde promediar sí tiene sentido.
Código
from scipy.ndimage import zoomrecorte = verdad[380:460, 380:460] # un pedazo con varias clasesvecino = zoom(recorte, 3, order=0) # vecino mas cercanobilineal = zoom(recorte.astype(float), 3, order=1) # bilineal (MAL para categorias)fig, axes = plt.subplots(1, 3, figsize=(13.5, 4.6))axes[0].imshow(recorte, cmap=CMAP_CLASES, vmin=0, vmax=4, interpolation="nearest")axes[0].set_title(f"Original\nclases presentes: {sorted(np.unique(recorte))}")axes[1].imshow(vecino, cmap=CMAP_CLASES, vmin=0, vmax=4, interpolation="nearest")axes[1].set_title(f"Vecino más cercano (BIEN)\nclases: {sorted(np.unique(vecino))}")axes[2].imshow(bilineal, cmap=CMAP_CLASES, vmin=0, vmax=4, interpolation="nearest")axes[2].set_title(f"Bilineal (MAL)\nvalores distintos: {len(np.unique(bilineal))}")for ax in axes: ax.set_xticks([]); ax.set_yticks([]); ax.grid(False)plt.tight_layout(); plt.show()print("El remuestreo bilineal promedio la clase 0 (agua) con la 2 (cultivo)")print("y produjo un 1 (bosque) donde no habia bosque. Las etiquetas son")print("numeros, pero NO son cantidades: promediarlas no significa nada.")
El remuestreo bilineal promedio la clase 0 (agua) con la 2 (cultivo)
y produjo un 1 (bosque) donde no habia bosque. Las etiquetas son
numeros, pero NO son cantidades: promediarlas no significa nada.
9. Índices espectrales
Un índice espectral combina bandas para resaltar un fenómeno. Casi siempre es una diferencia normalizada, que acota el resultado a \([-1, 1]\) y lo hace menos sensible a la iluminación.
Esto es álgebra de mapas de tipo local: una operación celda por celda, aplicada a más de un millón de píxeles a la vez, sin escribir un solo ciclo.
Cómo se lee un NDVI
El NDVI es el índice más usado en teledetección, y su lógica es directa: la vegetación sana absorbe rojo y refleja NIR, así que la diferencia se dispara. Va de −1 a 1 y se lee así:
Valor
Qué hay ahí
negativo
agua, nieve, nubes
cerca de 0
suelo desnudo, roca, área construida
0.2 a 0.5
vegetación escasa o pastizales
arriba de 0.6
vegetación densa y sana
Los otros índices que vale la pena conocer
Todos siguen la misma receta de diferencia normalizada, cambiando qué bandas se restan:
Índice
Bandas
Para qué
NDWI
verde y NIR
cuerpos de agua y humedad
NDBI
SWIR y NIR
áreas construidas
EVI
corrige NDVI por suelo y atmósfera
no se satura en vegetación muy densa
NBR
NIR y SWIR
severidad de incendios: la diferencia antes/después cuantifica el área quemada
La razón de que casi todos sean diferencias normalizadas: el cociente acota el resultado a [−1, 1] y lo vuelve menos sensible a la iluminación. Si una ladera está en sombra, ambas bandas bajan proporcionalmente y el índice apenas se mueve — justo lo que quieres.
Código
def diferencia_normalizada(a, b):'''(a - b) / (a + b), protegida contra división entre cero.'''return (a - b) / (a + b +1e-10)ndvi = diferencia_normalizada(escena[NIR], escena[ROJO])ndwi = diferencia_normalizada(escena[VERDE], escena[NIR])ndbi = diferencia_normalizada(escena[SWIR1], escena[NIR])indices = [("NDVI — vegetación", ndvi, "RdYlGn"), ("NDWI — agua", ndwi, "BrBG"), ("NDBI — construido", ndbi, "coolwarm")]fig, axes = plt.subplots(1, 3, figsize=(14.5, 5))for ax, (titulo, arr, cmap) inzip(axes, indices): im = ax.imshow(arr, cmap=cmap, extent=EXT, vmin=-1, vmax=1) ax.set_title(titulo); ax.set_xticks([]); ax.set_yticks([]); ax.grid(False) fig.colorbar(im, ax=ax, fraction=0.046)plt.tight_layout(); plt.show()print(f"NDVI min={ndvi.min():+.2f} media={ndvi.mean():+.2f} max={ndvi.max():+.2f}")
NDVI min=-1.00 media=+0.55 max=+1.00
El NDVI separa lo que el ojo no
Un solo número por píxel, y las coberturas se ordenan solas.
Código
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(13, 4.2))for k, (clase, color) inenumerate(zip(CLASES, COLORES)): ax1.hist(ndvi[verdad == k].ravel()[::17], bins=70, alpha=0.55, color=color, label=clase, density=True)ax1.set_xlabel("NDVI"); ax1.set_ylabel("densidad")ax1.set_title("Distribución del NDVI por cobertura")ax1.legend(fontsize=8)# Umbral simple: clasificar agua solo con NDWIumbral =0.0prediccion_agua = ndwi > umbralreal_agua = verdad ==0vp = (prediccion_agua & real_agua).sum()fp = (prediccion_agua &~real_agua).sum()fn = (~prediccion_agua & real_agua).sum()precision = vp / (vp + fp) if vp + fp else0recall = vp / (vp + fn) if vp + fn else0ax2.imshow(prediccion_agua, cmap="Blues", extent=EXT)ax2.set_title(f"Agua detectada solo con NDWI > {umbral}\n"f"precisión={precision:.1%} · recall={recall:.1%}")ax2.set_xticks([]); ax2.set_yticks([]); ax2.grid(False)plt.tight_layout(); plt.show()
Un umbral sobre un índice, sin ningún modelo entrenado, ya detecta el agua con buena precisión. Para un proyecto de ciencia de datos, estos índices son variables predictoras listas para usar.
10. Operaciones SIG
Qué es un SIG
Un Sistema de Información Geográfica es un sistema para capturar, almacenar, consultar, analizar y visualizar datos con referencia geográfica. La definición importante es la negativa:
Un SIG no es un programa de mapas. Es una base de datos donde la ubicación es un campo de primera clase, con el que puedes consultar igual que consultas por nombre o por fecha.
Un SIG tiene cinco componentes, y solo uno es software:
Componente
Qué incluye
Hardware
computadoras, GPS, servidores
Software
QGIS, ArcGIS, PostGIS, librerías de Python (lo que usamos aquí)
Datos
la parte más cara y más valiosa del sistema
Personas
quienes diseñan el análisis e interpretan los resultados
Métodos
los procedimientos y flujos de trabajo documentados
El modelo de capas
Un SIG organiza el territorio en capas temáticas superpuestas sobre la misma zona: hidrografía, vías, uso del suelo, población, parcelas. Cada capa es independiente y se combina con las demás cuando el análisis lo pide.
Esa separación es lo que permite responder preguntas que ninguna capa contesta por sí sola — y es literalmente lo que haremos en la sección 12.
El catálogo de operaciones
Operación
Qué hace
¿Dónde?
Consulta por atributo
seleccionar por valores, como en SQL
abajo
Consulta espacial
seleccionar por relación geométrica: contiene, intersecta, toca
abajo
Buffer
el área a cierta distancia de una geometría
abajo
Overlay
combinar capas: intersección, unión, diferencia
sección 12
Clip
recortar una capa con la forma de otra
abajo
Disolución
fusionar geometrías vecinas que comparten un atributo
abajo
Unión espacial
pegarle a cada punto los atributos del polígono que lo contiene
abajo — la más común al preparar datos
Estadística zonal
resumir un ráster dentro de cada polígono
abajo
Geocodificación
convertir una dirección de texto en coordenadas (y la inversa)
no aquí
Análisis de redes
rutas más cortas, áreas de servicio
mencionada en la sección 12
Ahora conectamos los dos modelos: vectorial (municipios, puntos) con ráster (la escena). Aquí la ubicación deja de ser una columna y se convierte en un operador de consulta.
Unión espacial
A cada punto de observación le pegamos los atributos del polígono que lo contiene. No hay llave común entre las tablas: la relación es puramente geométrica.
Código
obs = gpd.read_file(rutas["observaciones"])print(f"{len(obs)} observaciones de campo · CRS = {obs.crs.to_string()}")# Recortamos los municipios a la ventana de la escena (operacion 'clip').# Sin esto, un municipio que solo asoma en la esquina distorsionaria las# estadisticas zonales y el mapa.ventana = gpd.GeoDataFrame(geometry=[box(*limites)], crs=munis.crs)munis_v = gpd.clip(munis, ventana).reset_index(drop=True)munis_v = munis_v[munis_v.geometry.geom_type.isin(["Polygon", "MultiPolygon"])]munis_v["area_km2"] = munis_v.geometry.area /1e6munis_v = munis_v[munis_v["area_km2"] >1.0].reset_index(drop=True)print(f"{len(munis_v)} municipios tras recortar a la ventana y descartar astillas")# Union espacial: cada punto hereda el municipio que lo contieneobs_muni = gpd.sjoin(obs, munis_v[["municipio", "geometry"]], predicate="within")print(f"{len(obs_muni)} puntos cayeron dentro de algún municipio "f"({len(obs) -len(obs_muni)} quedaron fuera)")obs_muni[["clase", "municipio", "geometry"]].head()
3000 observaciones de campo · CRS = EPSG:32615
58 municipios tras recortar a la ventana y descartar astillas
3000 puntos cayeron dentro de algún municipio (0 quedaron fuera)
clase
municipio
geometry
0
bosque
Chichicastenango
POINT (708490 1639970)
1
bosque
Siquinalá
POINT (721810 1594250)
2
bosque
Siquinalá
POINT (723790 1587350)
3
suelo desnudo
Santo Domingo Xenacoj
POINT (745570 1624490)
4
bosque
Parramos
POINT (733690 1612070)
Consulta espacial
En SQL seleccionas por el valor de un campo. En un SIG puedes además seleccionar por la relación geométrica entre dos capas: contiene, intersecta, está dentro de, toca, cruza. Es la misma idea de un WHERE, pero el predicado es geometría.
Código
centro = munis_v.geometry.iloc[len(munis_v) //2]vecindario = gpd.GeoDataFrame(geometry=[centro], crs=munis_v.crs)predicados = ["intersects", "touches", "within", "contains"]print("El mismo par de capas, cuatro predicados distintos:\n")for pred in predicados: sel = munis_v[munis_v.geometry.apply(lambda g: getattr(g, pred)(centro))]print(f" {pred:12} -> {len(sel):2d} municipios")# Consulta por atributo (SQL de toda la vida) vs consulta espacialgrandes = munis_v[munis_v["area_km2"] > munis_v["area_km2"].median()]print(f"\n consulta por atributo (area > mediana) -> {len(grandes)} municipios")
El mismo par de capas, cuatro predicados distintos:
intersects -> 8 municipios
touches -> 7 municipios
within -> 1 municipios
contains -> 1 municipios
consulta por atributo (area > mediana) -> 29 municipios
Fíjate en el detalle que confunde a todo el mundo la primera vez: within y contains devuelven 1, el propio municipio. Una geometría está contenida en sí misma. Y intersects devuelve uno más que touches, por la misma razón. En una consulta real hay que excluir explícitamente la unidad de referencia, o el resultado viene contaminado.
Buffer y disolución
Dos operaciones que aparecen en casi todo flujo de trabajo real:
Buffer: generar el área a cierta distancia de una geometría. Es la operación con la que se responden preguntas de proximidad —«qué hay a menos de 500 m de este río»—. Solo tiene sentido en un CRS métrico: en grados, el buffer sale deformado.
Disolución: fusionar geometrías vecinas que comparten un atributo. Es el groupby de la geometría.
Código
from shapely.geometry import LineString# Una capa de vias SINTETICA para la demo: lineas que conectan los tres nucleos# urbanos de la escena. (No hay una capa de carreteras real en este conjunto.)urbano_mask = verdad ==3fu, cu = np.where(urbano_mask)sub = np.random.default_rng(3).choice(len(fu), 3000, replace=False)from sklearn.cluster import KMeanskm_nucleos = KMeans(n_clusters=3, n_init=10, random_state=42).fit( np.column_stack([fu[sub], cu[sub]]))nucleos_rc = km_nucleos.cluster_centers_xs, ys = rasterio.transform.xy(transform, nucleos_rc[:, 0], nucleos_rc[:, 1])nucleos = gpd.GeoDataFrame( {"nucleo": [f"núcleo {i+1}"for i inrange(3)]}, geometry=gpd.points_from_xy(xs, ys), crs=munis_v.crs)vias = gpd.GeoDataFrame( {"via": ["ruta 1", "ruta 2"]}, geometry=[LineString([nucleos.geometry[0], nucleos.geometry[1]]), LineString([nucleos.geometry[1], nucleos.geometry[2]])], crs=munis_v.crs)RADIO_VIA =1200# metrosfranja = gpd.GeoDataFrame(geometry=[vias.union_all().buffer(RADIO_VIA)], crs=munis_v.crs)print(f"Longitud total de vías : {vias.length.sum() /1000:.1f} km")print(f"Área de la franja de {RADIO_VIA} m : {franja.area.iloc[0] /1e6:.1f} km²")# Disolucion: los municipios se fusionan en 'regiones' segun su NDVImunis_v["region"] = np.where(munis_v.geometry.centroid.y > munis_v.geometry.centroid.y.median(),"norte", "sur")regiones = munis_v.dissolve(by="region", aggfunc={"area_km2": "sum"}).reset_index()print(f"\nDisolución: {len(munis_v)} municipios -> {len(regiones)} regiones")print(regiones[["region", "area_km2"]].to_string(index=False))
Longitud total de vías : 67.8 km
Área de la franja de 1200 m : 166.3 km²
Disolución: 58 municipios -> 2 regiones
region area_km2
norte 1558.847548
sur 2215.308369
Código
fig, axes = plt.subplots(1, 3, figsize=(15, 4.8))munis_v.boundary.plot(ax=axes[0], color="#adb5bd", linewidth=0.4)franja.plot(ax=axes[0], color="#ffd8a8", alpha=0.75, edgecolor="#e8590c", linewidth=0.8)vias.plot(ax=axes[0], color="#d9480f", linewidth=1.8)nucleos.plot(ax=axes[0], color="#212529", markersize=45, marker="s")axes[0].set_title(f"Buffer: franja de {RADIO_VIA} m alrededor de las vías")regiones.plot(column="region", ax=axes[1], cmap="Set2", edgecolor="white", linewidth=0.6, legend=True, legend_kwds={"fontsize": 8})munis_v.boundary.plot(ax=axes[1], color="white", linewidth=0.3)axes[1].set_title("Disolución: municipios fusionados en regiones")munis_v.plot(ax=axes[2], facecolor="#e9ecef", edgecolor="white", linewidth=0.4)munis_v[munis_v.geometry.apply(lambda g: g.touches(centro))].plot( ax=axes[2], facecolor="#a5d8ff", edgecolor="#1971c2", linewidth=0.6)gpd.GeoSeries([centro], crs=munis_v.crs).plot(ax=axes[2], facecolor="#e03131", edgecolor="#c92a2a")axes[2].set_title("Consulta espacial: 'touches'\n(el municipio rojo y sus vecinos)")for ax in axes: ax.set_xticks([]); ax.set_yticks([]); ax.grid(False)plt.tight_layout(); plt.show()
Ese último mapa —un polígono y los que lo tocan— es literalmente la matriz de contigüidad que vamos a construir en la sección 14 para medir autocorrelación espacial. La operación SIG y el estadístico espacial son la misma idea vista dos veces.
Rasterización y estadística zonal
Para resumir un ráster dentro de cada polígono, primero rasterizamos los municipios (los quemamos en una cuadrícula idéntica a la escena) y luego agregamos con np.bincount. Es una operación zonal: se resume una capa dentro de las zonas que define otra.
Código
formas = ((geom, i) for i, geom inenumerate(munis_v.geometry, start=1))zonas = rasterize(formas, out_shape=verdad.shape, transform=transform, fill=0, dtype="int32")n_zonas =len(munis_v) +1plano = zonas.ravel()conteo = np.bincount(plano, minlength=n_zonas).astype(float)# NDVI medio por municipiosuma_ndvi = np.bincount(plano, weights=ndvi.ravel(), minlength=n_zonas)munis_v["ndvi_medio"] = (suma_ndvi / np.maximum(conteo, 1))[1:]# Proporcion de cada cobertura por municipiofor k, clase inenumerate(CLASES): suma_k = np.bincount(plano, weights=(verdad == k).ravel().astype(float), minlength=n_zonas) munis_v[f"prop_{clase.replace(' ', '_')}"] = (suma_k / np.maximum(conteo, 1))[1:]munis_v["pixeles"] = conteo[1:].astype(int)munis_v[["municipio", "area_km2", "pixeles", "ndvi_medio","prop_urbano", "prop_bosque"]].sort_values("prop_urbano", ascending=False).head(8)
Con rásteres, el análisis se vuelve aritmética entre matrices. Las cuatro familias de operaciones, y dónde ya usamos cada una:
Familia
Qué mira para calcular una celda
Ejemplo
Dónde en el cuaderno
Local
esa misma celda, en una o varias capas
NDVI
sección 9
Focal
la vecindad de la celda
suavizado, pendiente, textura
aquí
Zonal
todas las celdas de una zona definida por otra capa
NDVI medio por municipio
sección 10
Global
todo el ráster
distancia al elemento más cercano
aquí
Las que faltan son las dos que más rinden en machine learning geoespacial: focal produce textura y contexto, global produce distancias. Volveremos a ellas en la sección 19 para construir variables.
Código
from scipy.ndimage import uniform_filter, distance_transform_edt, sobel# --- FOCAL: el valor de una celda sale de su vecindad -----------------------ventana_px =5# 5x5 pixeles = 300 x 300 mmedia_focal = uniform_filter(ndvi, size=ventana_px)# Desviacion estandar focal = TEXTURA. Alta donde el paisaje es heterogeneo.media_cuad = uniform_filter(ndvi **2, size=ventana_px)textura = np.sqrt(np.maximum(media_cuad - media_focal **2, 0))# Deteccion de bordes con el gradiente de Sobelbordes = np.hypot(sobel(media_focal, axis=0), sobel(media_focal, axis=1))fig, axes = plt.subplots(1, 4, figsize=(16, 4.4))paneles = [ ("NDVI original (local)", ndvi, "RdYlGn", -1, 1), (f"Media focal {ventana_px}x{ventana_px} (suavizado)", media_focal, "RdYlGn", -1, 1), ("Desviación focal = TEXTURA", textura, "magma", 0, 0.25), ("Bordes (gradiente de Sobel)", bordes, "gray_r", 0, 0.4),]for ax, (titulo, arr, cmap, lo, hi) inzip(axes, paneles): im = ax.imshow(arr, cmap=cmap, extent=EXT, vmin=lo, vmax=hi) ax.set_title(titulo, fontsize=9); fig.colorbar(im, ax=ax, fraction=0.046) ax.set_xticks([]); ax.set_yticks([]); ax.grid(False)plt.tight_layout(); plt.show()
La textura es información que ningún píxel tiene por sí solo: un bosque denso y un cultivo pueden tener NDVI parecido, pero el cultivo es un mosaico de parcelas y el bosque es homogéneo. Con un DEM, estas mismas operaciones focales dan pendiente y orientación de ladera, que son de las variables más predictivas en cualquier modelo ambiental.
Código
# --- GLOBAL: distancia al elemento mas cercano ------------------------------# Importante: usamos agua detectada por NDWI, es decir, informacion OBSERVABLE# en la imagen. Calcular la distancia usando el mapa de verdad seria fuga.agua_observada = ndwi >0.0urbano_observado = (ndbi >0.0) & (ndvi <0.2)# distance_transform_edt mide, para cada celda, la distancia al 0 mas cercanodist_agua = distance_transform_edt(~agua_observada) *60.0# metrosdist_urbano = distance_transform_edt(~urbano_observado) *60.0fig, axes = plt.subplots(1, 3, figsize=(15, 4.6))im = axes[0].imshow(dist_agua /1000, cmap="Blues_r", extent=EXT)axes[0].set_title("Distancia al agua más cercana (km)")fig.colorbar(im, ax=axes[0], fraction=0.046)im = axes[1].imshow(dist_urbano /1000, cmap="Oranges_r", extent=EXT)axes[1].set_title("Distancia a lo construido más cercano (km)")fig.colorbar(im, ax=axes[1], fraction=0.046)# El NDVI cae cerca de lo urbano: una relacion que solo existe en el espaciopaso =400axes[2].scatter(dist_urbano.ravel()[::paso] /1000, ndvi.ravel()[::paso], s=2, alpha=0.15, color="#1c7ed6")axes[2].set_xlabel("distancia a lo construido (km)"); axes[2].set_ylabel("NDVI")axes[2].set_title("Variable derivada contra variable objetivo")for ax in axes[:2]: ax.set_xticks([]); ax.set_yticks([]); ax.grid(False)plt.tight_layout(); plt.show()print(f"Distancia media al agua : {dist_agua.mean() /1000:.2f} km "f"(máxima {dist_agua.max() /1000:.1f} km)")print(f"Correlación NDVI ~ distancia a lo construido: "f"{np.corrcoef(dist_urbano.ravel()[::97], ndvi.ravel()[::97])[0, 1]:+.2f}")
Distancia media al agua : 5.00 km (máxima 23.3 km)
Correlación NDVI ~ distancia a lo construido: +0.50
12. Análisis encadenado
Ninguna capa responde sola una pregunta real. El valor de un SIG está en encadenar operaciones hasta que la respuesta salga de la combinación. Este es el ejemplo clásico de análisis de idoneidad, resuelto con nuestros datos:
Partir de la demanda: dónde vive la gente (densidad de superficie urbana).
Restar el área ya cubierta por los centros existentes (buffer).
Exigir accesibilidad: estar dentro de la franja de las vías.
Excluir lo que no es construible: agua y áreas protegidas.
Ordenar los candidatos por población servida.
Código
from rasterio.features import geometry_maskRADIO_COBERTURA =6000# m — area de servicio de un centro existenteVENTANA_DEMANDA =25# px = 1.5 km, para estimar densidad de poblacion# --- 1. Demanda: densidad de superficie urbana (proxy de poblacion) ---------demanda = uniform_filter(urbano_observado.astype(float), size=VENTANA_DEMANDA)# --- 2. Centros existentes y su area de cobertura ---------------------------centros = nucleos.iloc[:2].copy() # dos de los tres nucleoscobertura = centros.geometry.buffer(RADIO_COBERTURA)cubierto =~geometry_mask(cobertura, out_shape=verdad.shape, transform=transform, invert=False)# --- 3. Accesibilidad: la franja de las vias de la seccion 10 ---------------accesible =~geometry_mask(franja.geometry, out_shape=verdad.shape, transform=transform, invert=False)# --- 4. Exclusiones: agua y area protegida (aqui, el bosque) ----------------protegida = verdad ==1excluido = agua_observada | protegida# --- 5. Idoneidad: todas las condiciones a la vez ---------------------------idoneo = (demanda >0.10) & (~cubierto) & accesible & (~excluido)print(f"Población (proxy) total en la ventana : {demanda.sum():,.0f} unidades")print(f" ya cubierta por los centros actuales: {demanda[cubierto].sum() / demanda.sum():.1%}")print(f" sin cubrir : {demanda[~cubierto].sum() / demanda.sum():.1%}")print(f"\nPíxeles idóneos para un centro nuevo : {idoneo.sum():,} "f"({idoneo.mean():.2%} de la escena)")
Población (proxy) total en la ventana : 149,459 unidades
ya cubierta por los centros actuales: 9.2%
sin cubrir : 90.8%
Píxeles idóneos para un centro nuevo : 7,906 (0.75% de la escena)
# Ranking de municipios por poblacion NO cubierta, para priorizar la inversionplano_zonas = zonas.ravel()n_z =len(munis_v) +1pop_total = np.bincount(plano_zonas, weights=demanda.ravel(), minlength=n_z)[1:]pop_sin = np.bincount(plano_zonas, weights=(demanda *~cubierto).ravel(), minlength=n_z)[1:]sitios = np.bincount(plano_zonas, weights=idoneo.ravel().astype(float), minlength=n_z)[1:]ranking = pd.DataFrame({"municipio": munis_v["municipio"],"poblacion_proxy": pop_total.round(0),"no_cubierta": pop_sin.round(0),"% sin cubrir": (pop_sin / np.maximum(pop_total, 1e-9) *100).round(1),"px_idoneos": sitios.astype(int),}).sort_values("no_cubierta", ascending=False)print("Municipios prioritarios: mucha población sin cubrir Y sitios disponibles\n")ranking.head(8)
Municipios prioritarios: mucha población sin cubrir Y sitios disponibles
municipio
poblacion_proxy
no_cubierta
% sin cubrir
px_idoneos
46
San Martín Jilotepeque
34882.0
34882.0
100.0
0
20
Escuintla
16385.0
16385.0
100.0
2581
14
Yepocapa
10840.0
10840.0
100.0
0
36
San Juan Sacatepéquez
9854.0
8108.0
82.3
0
15
Pochuta
6957.0
6699.0
96.3
0
16
Acatenango
9407.0
6357.0
67.6
0
12
Siquinalá
6143.0
6143.0
100.0
0
11
Santa Lucía Cotzumalguapa
5963.0
5963.0
100.0
0
Fíjate en lo que acaba de pasar: la respuesta final —una lista ordenada de municipios— no está en ninguna de las capas de entrada. Sale de encadenar buffer, overlay booleano, rasterización y estadística zonal.
Dos honestidades sobre este análisis:
El área de cobertura real no es un círculo, es un área de servicio de red: los 30 minutos de viaje siguen carreteras, no líneas rectas. Un buffer euclidiano sobrestima la cobertura en terreno montañoso, que es justo el caso de Guatemala.
La «población» aquí es un proxy construido de la imagen. Con datos reales usarías el censo por sector, y entonces aparece el problema de la sección 21: la unidad de agregación cambia el resultado.
La cadena de decisiones es el modelo
Cada umbral que escribimos —demanda > 0.10, 6 km de cobertura, 1200 m de franja— es una decisión de política disfrazada de parámetro. Cambiarlos cambia la lista de municipios. En un análisis de verdad, esos números se justifican y se someten a un análisis de sensibilidad.
13. Cartografía honesta
El mapa es el resultado y también la herramienta de exploración. Antes de las reglas, el catálogo de qué tipo de mapa usar:
Tipo de mapa
Cuándo
Coropleta
una tasa o proporción por unidad administrativa
Símbolos proporcionales
una cantidad absoluta por lugar (el tamaño del símbolo, no el color)
Puntos
eventos individuales, cuando la ubicación exacta importa
Densidad de kernel
muchos puntos, cuando importa la concentración (sección 15)
Isolíneas
una superficie continua: elevación, temperatura
Mapas de flujo
movimiento entre origen y destino
La coropleta es la más usada y la más fácil de arruinar. Dos formas de arruinarla, en orden de frecuencia.
No normalizar
La regla es corta: un mapa de conteos absolutos por municipio es, en el fondo, un mapa de dónde vive la gente —o de qué municipio es más grande—. Usa tasas, densidades o porcentajes. Aquí la versión ráster del mismo error.
Mapeamos dos cosas con los mismos datos: el número absoluto de píxeles urbanos por municipio, y la proporción urbana. Compara los dos mapas.
Código
munis_v["px_urbano"] = (munis_v["prop_urbano"] * munis_v["pixeles"]).round().astype(int)fig, axes = plt.subplots(1, 3, figsize=(15, 4.8))munis_v.plot(column="px_urbano", cmap="OrRd", legend=True, ax=axes[0], edgecolor="white", linewidth=0.4, legend_kwds={"shrink": 0.7})axes[0].set_title("MAL: conteo absoluto de píxeles urbanos")munis_v.plot(column="prop_urbano", cmap="OrRd", legend=True, ax=axes[1], edgecolor="white", linewidth=0.4, legend_kwds={"shrink": 0.7})axes[1].set_title("BIEN: proporción urbana del municipio")munis_v.plot(column="area_km2", cmap="Greys", legend=True, ax=axes[2], edgecolor="white", linewidth=0.4, legend_kwds={"shrink": 0.7})axes[2].set_title("El culpable: el área de cada municipio")for ax in axes: ax.set_xticks([]); ax.set_yticks([]); ax.grid(False)plt.tight_layout(); plt.show()r_area = np.corrcoef(munis_v["px_urbano"], munis_v["area_km2"])[0, 1]r_prop = np.corrcoef(munis_v["prop_urbano"], munis_v["area_km2"])[0, 1]print(f"Correlación con el área del municipio:")print(f" conteo absoluto : r = {r_area:+.2f} <- el mapa de conteos es, en parte, un mapa de áreas")print(f" proporción : r = {r_prop:+.2f}")
Correlación con el área del municipio:
conteo absoluto : r = +0.35 <- el mapa de conteos es, en parte, un mapa de áreas
proporción : r = +0.15
El conteo absoluto correlaciona con el área del municipio más del doble que la proporción. No es una correlación aplastante —los municipios de esta ventana tienen tamaños parecidos, así que el efecto es moderado— pero va en la dirección que delata el problema: parte de lo que muestra el mapa de conteos es simplemente cuánta superficie tiene cada municipio, no cuán urbano es.
Con unidades de tamaños muy dispares (departamentos contra municipios, o países contra ciudades) esa correlación se dispara y el mapa deja de decir nada sobre el fenómeno. Normalizar no es un detalle estético: es lo que separa un resultado de un artefacto de la geometría.
La clasificación
Ya normalizamos. Ahora falta decidir cómo se cortan los intervalos, y esa decisión cambia el mapa tanto como cambiarían los datos. Los tres esquemas habituales, con exactamente los mismos valores:
Intervalos iguales: se parte el rango en trozos del mismo ancho. Intuitivo, pero si la distribución es asimétrica casi todo cae en una clase.
Cuantiles: cada clase tiene el mismo número de unidades. Siempre produce un mapa “bonito” y equilibrado, aunque no haya ninguna diferencia real.
Cortes naturales (Jenks): busca los cortes que minimizan la varianza dentro de cada clase. Es un k-means en una dimensión, y así lo implementamos.
Código
from sklearn.cluster import KMeansK =5valores = munis_v["ndvi_medio"].valuesdef cortes_intervalos_iguales(x, k):return np.linspace(x.min(), x.max(), k +1)def cortes_cuantiles(x, k):return np.quantile(x, np.linspace(0, 1, k +1))def cortes_jenks(x, k):'''Cortes naturales: k-means en una dimension, que es lo que Jenks optimiza.''' km = KMeans(n_clusters=k, n_init=25, random_state=42).fit(x.reshape(-1, 1)) orden = np.argsort(km.cluster_centers_.ravel()) etiquetas = np.argsort(orden)[km.labels_] bordes = [x.min()]for c inrange(k -1): bordes.append((x[etiquetas == c].max() + x[etiquetas == c +1].min()) /2)return np.array(bordes + [x.max()])esquemas = [("Intervalos iguales", cortes_intervalos_iguales), ("Cuantiles", cortes_cuantiles), ("Cortes naturales (Jenks)", cortes_jenks)]fig, axes = plt.subplots(2, 3, figsize=(15, 8), gridspec_kw={"height_ratios": [3, 1]})for j, (nombre, funcion) inenumerate(esquemas): bordes = funcion(valores, K) clase = np.clip(np.digitize(valores, bordes[1:-1]), 0, K -1) munis_v["_clase"] = clase munis_v.plot(column="_clase", cmap="RdYlGn", ax=axes[0, j], vmin=0, vmax=K -1, edgecolor="white", linewidth=0.4) reparto = np.bincount(clase, minlength=K) axes[0, j].set_title(f"{nombre}\nunidades por clase: {list(reparto)}", fontsize=9) axes[0, j].set_xticks([]); axes[0, j].set_yticks([]); axes[0, j].grid(False) axes[1, j].hist(valores, bins=30, color="#adb5bd")for b in bordes[1:-1]: axes[1, j].axvline(b, color="#e03131", linewidth=1.6) axes[1, j].set_xlabel("NDVI medio del municipio"); axes[1, j].set_yticks([])fig.suptitle("El mismo mapa, tres esquemas de clasificación", y=1.0)plt.tight_layout(); plt.show()munis_v = munis_v.drop(columns="_clase")
Los tres mapas son correctos y los tres cuentan una historia distinta. La regla es simple y casi nunca se cumple: declara el esquema que usaste en la leyenda.
Sobre las paletas, tres reglas que evitan la mayoría de los mapas engañosos:
Secuencial (claro → oscuro) para magnitudes que solo crecen.
Divergente (dos colores con un centro neutro) solo cuando existe un punto medio significativo: cero, la media, un umbral de política.
Cualitativa para categorías sin orden.
Nunca arcoíris para datos continuos: el ojo lee saltos donde no los hay, y la paleta se vuelve ilegible en escala de grises y para daltónicos.
Cuando la geometría misma estorba
Queda un problema que ni normalizar ni clasificar resuelven: las unidades grandes dominan visualmente aunque importen poco. Un mapa de Guatemala por departamento le da a Petén más superficie visual que a toda la zona metropolitana, donde vive la gente.
Dos salidas:
Cartograma: se deforma cada unidad para que su tamaño represente la variable (población, votos, casos) en vez de su área geográfica. Se pierde la forma reconocible y se gana honestidad en la proporción.
Mapa de hexágonos: se sustituyen las unidades reales por celdas del mismo tamaño. Cada lugar pesa lo mismo y desaparece el sesgo del área — a cambio de romper las fronteras administrativas.
Y una advertencia final que se aplica a todo lo anterior: cuidado con la proyección. Web Mercator infla las latitudes altas, así que un mapa mundial en esa proyección compara áreas mal por construcción, sin importar cuán bien elegida esté la paleta.
14. Autocorrelación espacial
La primera ley de la geografía dice que las cosas cercanas se parecen más. El índice de Moran lo mide:
donde \(w_{ij}\) es el peso que expresa qué tan vecinas son las unidades \(i\) y \(j\). Se lee como una correlación:
Valor de I
Patrón
Cómo se ve
cerca de +1
agrupado
los valores altos rodeados de altos, los bajos de bajos
cerca de 0
aleatorio
no hay estructura espacial
cerca de −1
disperso
como un tablero de ajedrez: alto, bajo, alto, bajo
Por qué esto no es una curiosidad estadística
Si hay autocorrelación espacial, las observaciones no son independientes — y casi toda la estadística clásica supone que sí lo son. La consecuencia concreta:
Con autocorrelación, los errores estándar salen subestimados y se declaran significativos efectos que no lo son.
Dicho de otro modo: si tienes 58 municipios autocorrelacionados, no tienes 58 observaciones independientes, tienes efectivamente menos. Tu prueba de hipótesis cree tener más información de la que hay.
Esa misma idea, aplicada al machine learning, es lo que hará explotar el modelo en la sección 17.
Vamos a implementarlo sin librerías de estadística espacial: primero la matriz de pesos por contigüidad (dos municipios son vecinos si comparten frontera), luego el índice.
Código
def matriz_contiguidad(gdf):'''Pesos por contigüidad: w[i,j]=1 si los polígonos i y j se tocan.''' n =len(gdf) W = np.zeros((n, n)) arbol = gdf.sindex # índice espacial: evita comparar todos con todosfor i, geom inenumerate(gdf.geometry):for j in arbol.query(geom, predicate="intersects"):if i != j and geom.touches(gdf.geometry.iloc[j]): W[i, j] =1.0return Wdef estandarizar_filas(W):'''Cada fila suma 1: el vecindario de cada unidad pesa lo mismo.''' s = W.sum(axis=1, keepdims=True)return np.divide(W, s, out=np.zeros_like(W), where=s >0)def moran_i(x, W): z = np.asarray(x, float) - np.mean(x) S0 = W.sum()return (len(z) / S0) * (z @ W @ z) / (z @ z)W = matriz_contiguidad(munis_v)Wz = estandarizar_filas(W)vecinos = W.sum(axis=1)print(f"{len(munis_v)} municipios · {int(W.sum() /2)} pares vecinos")print(f"vecinos por municipio: min={vecinos.min():.0f} "f"media={vecinos.mean():.1f} max={vecinos.max():.0f}")
def prueba_permutacion(x, W, n_perm=999, semilla=42):'''¿El patrón observado podría salir de repartir los valores al azar?''' rng = np.random.default_rng(semilla) observado = moran_i(x, W) nulos = np.array([moran_i(rng.permutation(x), W) for _ inrange(n_perm)]) p = (np.sum(np.abs(nulos) >=abs(observado)) +1) / (n_perm +1)return observado, nulos, pvariable = munis_v["ndvi_medio"].valuesI, nulos, p = prueba_permutacion(variable, Wz)print(f"Índice de Moran del NDVI medio: I = {I:+.3f}")print(f"Valor esperado bajo aleatoriedad: E[I] = {-1/ (len(variable) -1):+.3f}")print(f"p-valor (999 permutaciones): {p:.4f}")print("\nInterpretación:", "patrón AGRUPADO — hay autocorrelación espacial"if I >0and p <0.05else"no se distingue de un patrón aleatorio")
Índice de Moran del NDVI medio: I = +0.267
Valor esperado bajo aleatoriedad: E[I] = -0.018
p-valor (999 permutaciones): 0.0020
Interpretación: patrón AGRUPADO — hay autocorrelación espacial
Código
fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(15, 4.4))# 1. Distribucion nula vs observadoax1.hist(nulos, bins=40, color="#adb5bd", edgecolor="white")ax1.axvline(I, color="#e03131", linewidth=2.5, label=f"observado I={I:+.3f}")ax1.axvline(nulos.mean(), color="#495057", linestyle="--", linewidth=1, label=f"media nula={nulos.mean():+.3f}")ax1.set_xlabel("I de Moran"); ax1.set_ylabel("frecuencia")ax1.set_title(f"Prueba de permutaciones\np = {p:.4f}")ax1.legend(fontsize=7)# 2. Diagrama de dispersion de Moranz = (variable - variable.mean()) / variable.std()rezago = Wz @ z # promedio estandarizado del vecindarioax2.scatter(z, rezago, s=28, color="#1c7ed6", alpha=0.75, edgecolor="white")pendiente = np.polyfit(z, rezago, 1)[0]xs = np.linspace(z.min(), z.max(), 10)ax2.plot(xs, pendiente * xs, color="#e03131", linewidth=2, label=f"pendiente = I = {pendiente:+.3f}")ax2.axhline(0, color="gray", linewidth=0.6); ax2.axvline(0, color="gray", linewidth=0.6)ax2.set_xlabel("NDVI estandarizado del municipio")ax2.set_ylabel("NDVI medio de sus vecinos")ax2.set_title("Diagrama de dispersión de Moran")ax2.legend(fontsize=7)# 3. LISA: donde estan los conglomeradosdef moran_local(x, W, n_perm=499, semilla=42): z = (np.asarray(x, float) - np.mean(x)) / np.std(x) Ii = z * (W @ z) rng = np.random.default_rng(semilla) nulos = np.empty((n_perm, len(z)))for k inrange(n_perm): zp = rng.permutation(z) nulos[k] = z * (W @ zp) p = (np.sum(np.abs(nulos) >= np.abs(Ii), axis=0) +1) / (n_perm +1)return Ii, p, z, W @ zIi, p_local, z_std, rezago_std = moran_local(variable, Wz)etiqueta = np.full(len(z_std), "no significativo", dtype=object)sig = p_local <0.05etiqueta[sig & (z_std >0) & (rezago_std >0)] ="alto-alto"etiqueta[sig & (z_std <0) & (rezago_std <0)] ="bajo-bajo"etiqueta[sig & (z_std >0) & (rezago_std <0)] ="alto-bajo"etiqueta[sig & (z_std <0) & (rezago_std >0)] ="bajo-alto"munis_v["lisa"] = etiquetapaleta = {"alto-alto": "#2b8a3e", "bajo-bajo": "#c92a2a","alto-bajo": "#f59f00", "bajo-alto": "#1c7ed6","no significativo": "#e9ecef"}for cat, color in paleta.items(): sub = munis_v[munis_v["lisa"] == cat]iflen(sub): sub.plot(ax=ax3, color=color, edgecolor="white", linewidth=0.4, label=f"{cat} ({len(sub)})")ax3.set_title("LISA — conglomerados locales de NDVI")ax3.legend(fontsize=7, loc="upper left"); ax3.set_xticks([]); ax3.set_yticks([])ax3.grid(False)plt.tight_layout(); plt.show()
El diagrama de dispersión de Moran es la mejor forma de entender el índice: en el eje horizontal va el valor de cada municipio, en el vertical el promedio de sus vecinos, y la pendiente de la recta es el propio índice de Moran.
El mapa LISA muestra dónde están los conglomerados: los verdes son vecindarios consistentemente verdes (alto-alto), los rojos son bolsones de NDVI bajo rodeados de NDVI bajo (bajo-bajo), y los naranjas y azules son valores atípicos espaciales: un municipio que no se parece a sus vecinos.
Este resultado es también la advertencia que viene a continuación: si los municipios vecinos se parecen, las observaciones no son independientes, y casi toda la estadística clásica supone que sí lo son.
Getis-Ord Gi*
Moran responde ¿hay agrupamiento?; LISA responde ¿dónde?. Getis-Ord Gi* responde una pregunta ligeramente distinta y muy usada en la práctica: ¿dónde hay concentraciones de valores altos (puntos calientes) o bajos (puntos fríos) que no se explican por azar?
El resultado es directamente un valor z: por encima de +1.96 es un punto caliente significativo, por debajo de −1.96 uno frío. La estrella del nombre significa que la unidad se incluye a sí misma en su vecindario, que es lo que uno quiere para hablar de “concentración”.
Código
def getis_ord_gi(x, W):'''Gi* como valor z. W binaria; la version estrella se incluye a si misma.''' x = np.asarray(x, float) n =len(x) Ws = W + np.eye(n) # la estrella: cada unidad es su vecina media = x.mean() S = np.sqrt((x **2).mean() - media **2) suma_w = Ws.sum(axis=1) suma_w2 = (Ws **2).sum(axis=1) numerador = Ws @ x - media * suma_w denominador = S * np.sqrt((n * suma_w2 - suma_w **2) / (n -1))return numerador / denominadorgi = getis_ord_gi(munis_v["prop_urbano"].values, W)munis_v["gi_urbano"] = gicalientes = (gi >1.96).sum()frios = (gi <-1.96).sum()print(f"Getis-Ord Gi* sobre la proporción urbana:")print(f" puntos calientes (z > +1.96): {calientes}")print(f" puntos fríos (z < -1.96): {frios}")print(f" no significativos : {len(gi) - calientes - frios}")fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12.5, 4.8))munis_v.plot(column="gi_urbano", cmap="RdBu_r", vmin=-3, vmax=3, legend=True, ax=ax1, edgecolor="white", linewidth=0.4, legend_kwds={"label": "Gi* (valor z)", "shrink": 0.75})ax1.set_title("Getis-Ord Gi*: concentración de superficie urbana")munis_v.plot(column="prop_urbano", cmap="OrRd", legend=True, ax=ax2, edgecolor="white", linewidth=0.4, legend_kwds={"label": "proporción urbana", "shrink": 0.75})ax2.set_title("La variable cruda, para comparar")for ax in (ax1, ax2): ax.set_xticks([]); ax.set_yticks([]); ax.grid(False)plt.tight_layout(); plt.show()
Getis-Ord Gi* sobre la proporción urbana:
puntos calientes (z > +1.96): 6
puntos fríos (z < -1.96): 0
no significativos : 52
La diferencia entre los dos mapas es el punto: el mapa crudo señala municipios con valores altos; Gi* señala vecindarios donde lo alto se acumula. Un municipio muy urbano rodeado de municipios rurales no es un punto caliente —es un valor atípico espacial, que es lo que detecta LISA—.
Es la herramienta estándar para mapear delito, brotes epidémicos o accidentes: lo que interesa no es el punto suelto, sino la zona.
15. Patrones de puntos e interpolación
Hasta aquí analizamos polígonos. Con puntos el juego cambia: no hay unidades predefinidas, y las preguntas son otras. ¿Se agrupan? ¿Dónde está la densidad? ¿Y qué pasa entre los puntos, donde no medimos nada?
DBSCAN
DBSCAN agrupa por densidad: un grupo es una zona donde hay al menos min_samples puntos dentro de un radio eps. No hay que fijar el número de grupos de antemano y —a diferencia de k-means— marca el ruido en vez de forzar cada punto a un grupo.
Código
from sklearn.cluster import DBSCANfrom scipy.ndimage import gaussian_filter as suavizarcoords = np.column_stack([obs.geometry.x.values, obs.geometry.y.values])EPS =1200# metros: dos puntos estan 'densamente conectados' si estan mas cercaMIN_MUESTRAS =8db = DBSCAN(eps=EPS, min_samples=MIN_MUESTRAS).fit(coords)etiquetas_db = db.labels_n_grupos =len(set(etiquetas_db)) - (1if-1in etiquetas_db else0)ruido = (etiquetas_db ==-1).sum()print(f"DBSCAN con eps={EPS} m y min_samples={MIN_MUESTRAS}")print(f" grupos encontrados : {n_grupos}")print(f" puntos marcados como ruido: {ruido} ({ruido /len(coords):.1%})")print(f"\nLas observaciones se generaron alrededor de ~66 cuadrillas de campo con")print(f"una dispersión de ~2 km. DBSCAN recuperó {n_grupos} de esos conglomerados sin")print(f"saber nada del diseño de muestreo: las cuadrillas que se solapan se funden.")print(f"\nY todo depende de eps — no hay un valor 'correcto', hay una escala elegida:")for eps in [600, 1200, 2000, 3000]: l = DBSCAN(eps=eps, min_samples=MIN_MUESTRAS).fit(coords).labels_ n =len(set(l)) - (1if-1in l else0)print(f" eps = {eps:5d} m -> {n:3d} grupos, {(l ==-1).mean():5.1%} de ruido")
DBSCAN con eps=1200 m y min_samples=8
grupos encontrados : 41
puntos marcados como ruido: 911 (30.4%)
Las observaciones se generaron alrededor de ~66 cuadrillas de campo con
una dispersión de ~2 km. DBSCAN recuperó 41 de esos conglomerados sin
saber nada del diseño de muestreo: las cuadrillas que se solapan se funden.
Y todo depende de eps — no hay un valor 'correcto', hay una escala elegida:
eps = 600 m -> 26 grupos, 90.4% de ruido
eps = 1200 m -> 41 grupos, 30.4% de ruido
eps = 2000 m -> 8 grupos, 2.5% de ruido
eps = 3000 m -> 2 grupos, 0.2% de ruido
Código
# Estimacion de densidad por kernel (KDE): de nube de puntos a superficie continuaCELDA =8# agregamos cada 8 px para el conteon_celdas = verdad.shape[0] // CELDAconteo = np.zeros((n_celdas, n_celdas))np.add.at(conteo, (obs["fila"].values // CELDA, obs["col"].values // CELDA), 1)densidad = suavizar(conteo, sigma=2.5)fig, axes = plt.subplots(1, 3, figsize=(15, 4.8))ruido_mask = etiquetas_db ==-1axes[0].scatter(coords[ruido_mask, 0], coords[ruido_mask, 1], s=4, color="#adb5bd", label=f"ruido ({ruido})")axes[0].scatter(coords[~ruido_mask, 0], coords[~ruido_mask, 1], s=5, c=etiquetas_db[~ruido_mask], cmap="tab20")axes[0].set_title(f"DBSCAN: {n_grupos} conglomerados de muestreo")axes[0].legend(fontsize=7, loc="upper right")im = axes[1].imshow(densidad, cmap="inferno", extent=EXT)axes[1].set_title("Densidad de kernel (KDE)\nla nube de puntos como superficie")fig.colorbar(im, ax=axes[1], fraction=0.046)axes[2].imshow(verdad, cmap=CMAP_CLASES, extent=EXT, vmin=0, vmax=4, alpha=0.4)# Coordenadas del centro de cada celda de la rejilla de densidad, en metros UTMgx_m = EXT[0] + (np.arange(n_celdas) * CELDA + CELDA /2) *60gy_m = EXT[3] - (np.arange(n_celdas) * CELDA + CELDA /2) *60axes[2].contour(gx_m, gy_m, densidad, levels=5, colors="black", linewidths=0.7)axes[2].set_title("Isolíneas de densidad sobre la cobertura")for ax in axes: ax.set_xticks([]); ax.set_yticks([]); ax.grid(False) ax.set_xlim(EXT[0], EXT[1]); ax.set_ylim(EXT[2], EXT[3])plt.tight_layout(); plt.show()
Cuidado con el mapa de calor. Ese mapa muestra dónde fuimos a medir, no dónde está el fenómeno. Es el error más común al leer un KDE: un mapa de densidad de casos reportados es, en buena parte, un mapa de dónde hay quien reporte. Sin normalizar por la población en riesgo —o por el esfuerzo de muestreo— no dice lo que parece decir.
Estos conglomerados no son un detalle: son la razón por la que la validación cruzada aleatoria va a fallar en la sección 17. DBSCAN nos acaba de dar, sin saber nada del experimento, los grupos que un modelo honesto debería respetar al partir los datos.
Interpolación
Con mediciones en pocos puntos queremos una superficie completa. IDW (distancia inversa ponderada) es el método más directo: cada punto conocido pesa en proporción inversa a su distancia, elevada a una potencia \(p\).
def idw(xy_conocidos, z, xy_destino, p=2.0, epsilon=1e-9):'''Distancia inversa ponderada, vectorizada.''' d = np.sqrt(((xy_destino[:, None, :] - xy_conocidos[None, :, :]) **2).sum(axis=2)) w =1.0/ (d ** p + epsilon)return (w @ z) / w.sum(axis=1)# Tomamos 120 'estaciones de medicion' del NDVI y reconstruimos toda la superficierng_i = np.random.default_rng(11)muestra = rng_i.choice(len(obs), 120, replace=False)f_m, c_m = obs["fila"].values[muestra], obs["col"].values[muestra]z_m = ndvi[f_m, c_m]xy_m = np.column_stack([c_m.astype(float), f_m.astype(float)])# Rejilla destino gruesa (128x128) y verdad agregada a esa misma rejillaPASO =8gy, gx = np.mgrid[0:verdad.shape[0]:PASO, 0:verdad.shape[1]:PASO]xy_destino = np.column_stack([gx.ravel().astype(float), gy.ravel().astype(float)])verdad_grid = ndvi[::PASO, ::PASO]fig, axes = plt.subplots(1, 4, figsize=(17, 4.4))axes[0].imshow(verdad_grid, cmap="RdYlGn", vmin=-1, vmax=1)axes[0].scatter(c_m / PASO, f_m / PASO, s=8, color="black")axes[0].set_title(f"NDVI real + las {len(muestra)} estaciones")for ax, p inzip(axes[1:], [1.0, 2.0, 6.0]): estimado = idw(xy_m, z_m, xy_destino, p=p).reshape(gy.shape) rmse = np.sqrt(np.mean((estimado - verdad_grid) **2)) ax.imshow(estimado, cmap="RdYlGn", vmin=-1, vmax=1) ax.set_title(f"IDW con p={p:g}\nRMSE = {rmse:.3f}")for ax in axes: ax.set_xticks([]); ax.set_yticks([]); ax.grid(False)fig.suptitle("Interpolación: la misma muestra, distinta potencia", y=1.02)plt.tight_layout(); plt.show()
La potencia \(p\) controla cuánto manda el vecino más cercano:
\(p\) bajo: la superficie es suave y tiende a la media global; se pierden los extremos.
\(p\) alto: cada punto domina su entorno inmediato y aparecen los famosos “ojos de buey” alrededor de cada estación.
Y la limitación de fondo: IDW no dice nada sobre su propia incertidumbre. Devuelve un número en cada celda sin decirte cuánto confiar en él.
Ahí entra el kriging, que modela explícitamente la estructura de la autocorrelación con un variograma —cuánto se parecen dos puntos en función de la distancia que los separa— y entrega, junto a la estimación, su error estándar. Es más trabajo y exige que la estructura espacial sea estable, pero es lo que se usa cuando la incertidumbre importa: contaminación, minería, clima.
Entre ambos hay alternativas que suavizan de otra forma —vecino natural, que pondera por el área que cada punto “cede” en una teselación de Voronoi, y los splines, que ajustan una superficie flexible que pasa por los puntos—. La elección cambia el aspecto del mapa más de lo que la mayoría admite.
Regla de decisión rápida. ¿Necesitas la incertidumbre? Kriging. ¿Solo necesitas una superficie razonable rápido? IDW. ¿La variable tiene saltos bruscos (fallas, fronteras)? Ninguno de los dos: ningún interpolador suave respeta una discontinuidad.
16. Clasificación supervisada de cobertura terrestre
La tarea clásica de la teledetección
Clasificación de cobertura terrestre: asignar a cada píxel una clase de uso del suelo —bosque, agua, cultivo, urbano, suelo desnudo—. Hay dos caminos:
Supervisada (esta sección): se etiquetan a mano regiones de entrenamiento y se ajusta un clasificador. Random Forest y gradient boosting funcionan muy bien con bandas e índices como variables.
No supervisada (sección 18): se agrupan los píxeles por similitud espectral y después se interpreta cada grupo.
Qué usar como variables
Las bandas crudas.
Los índices espectrales — que son bandas combinadas, pero condensan mejor la información (sección 9).
La textura y el contexto del vecindario (sección 19).
Y la más poderosa cuando existe: la serie temporal completa del año. Los cultivos se distinguen por su fenología —cómo cambia su NDVI a lo largo de la temporada— mucho mejor que por una sola fecha. Maíz y caña se ven casi iguales en marzo y muy distintos en el ciclo completo.
Aquí usamos 6 bandas + 3 índices = 9 variables por píxel, y entrenamos un Random Forest con las observaciones de campo para predecir la cobertura de cada píxel de la escena.
def bosque():return RandomForestClassifier(n_estimators=200, random_state=42, n_jobs=-1)modelo = bosque().fit(X, y)# Predecir TODA la escena: mas de un millon de pixelestodos = construir_variables(escena.reshape(6, -1).T)mapa_pred = modelo.predict(todos).reshape(verdad.shape)exactitud_total = (mapa_pred == verdad).mean()print(f"Exactitud sobre la escena completa: {exactitud_total:.1%}")print(f"({mapa_pred.size:,} píxeles clasificados)")
Exactitud sobre la escena completa: 82.5%
(1,048,576 píxeles clasificados)
Mira el mapa de errores: no están repartidos al azar. Se concentran en regiones enteras y en los bordes entre coberturas. Ese agrupamiento de los errores es autocorrelación espacial otra vez, y es la pista de lo que viene.
Código
mc = confusion_matrix(verdad.ravel()[::7], mapa_pred.ravel()[::7], normalize="true")fig, ax = plt.subplots(figsize=(5.6, 4.8))im = ax.imshow(mc, cmap="Blues", vmin=0, vmax=1)ax.set_xticks(range(5)); ax.set_yticks(range(5))ax.set_xticklabels(CLASES, rotation=45, ha="right", fontsize=8)ax.set_yticklabels(CLASES, fontsize=8)ax.set_xlabel("predicho"); ax.set_ylabel("real")ax.set_title("Matriz de confusión (normalizada por fila)")for i inrange(5):for j inrange(5): ax.text(j, i, f"{mc[i, j]:.2f}", ha="center", va="center", fontsize=8, color="white"if mc[i, j] >0.5else"#212529")ax.grid(False); fig.colorbar(im, fraction=0.046)plt.tight_layout(); plt.show()
Código
importancias = pd.Series(modelo.feature_importances_, index=NOMBRES_VAR).sort_values()fig, ax = plt.subplots(figsize=(6.5, 3.6))importancias.plot.barh(ax=ax, color="#1c7ed6")ax.set_title("Importancia de cada variable en el Random Forest")ax.set_xlabel("importancia")plt.tight_layout(); plt.show()
Los índices espectrales suelen pelear los primeros puestos con las bandas crudas: una división bien elegida entre dos bandas condensa más información sobre la cobertura que cualquiera de las dos por separado.
17. Fuga de información espacial
Aquí está el punto más importante del cuaderno, y la razón por la que el aprendizaje automático con datos geoespaciales necesita su propio capítulo.
Vamos a evaluar el mismo modelo, con los mismos datos, cambiando una sola cosa: cómo se parte el conjunto en entrenamiento y prueba.
Validación cruzada aleatoria — la de siempre: se barajan los puntos y se reparten en 5 pliegues.
Validación cruzada espacial por bloques — la escena se divide en una cuadrícula, y bloques enteros van a prueba. El modelo tiene que predecir en territorio que nunca vio.
Recuerda cómo se tomaron las observaciones: agrupadas, como en campo real. Con una partición aleatoria, un punto de prueba casi siempre tiene un vecino suyo en el entrenamiento.
Código
LADO_BLOQUE =4# cuadrícula de 4x4 = 16 bloquesn_px = verdad.shape[0]tam = n_px // LADO_BLOQUEbloque = (filas // tam) * LADO_BLOQUE + (cols // tam)fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12.5, 5.2))ax1.imshow(verdad, cmap=CMAP_CLASES, extent=EXT, vmin=0, vmax=4, alpha=0.35)ax1.scatter(obs.geometry.x, obs.geometry.y, s=2.5, c="#212529", alpha=0.55)for i inrange(1, LADO_BLOQUE): ax1.axvline(EXT[0] + i * tam *60, color="#e03131", linewidth=1.2) ax1.axhline(EXT[2] + i * tam *60, color="#e03131", linewidth=1.2)ax1.set_title("Observaciones agrupadas y los 16 bloques espaciales")ax1.set_xticks([]); ax1.set_yticks([]); ax1.grid(False)ax2.scatter(obs.geometry.x, obs.geometry.y, s=4, c=bloque, cmap="tab20", alpha=0.85)ax2.set_title("Cada color es un bloque: en la CV espacial\nun bloque entero va a prueba")ax2.set_xticks([]); ax2.set_yticks([]); ax2.grid(False)plt.tight_layout(); plt.show()print(f"{len(np.unique(bloque))} bloques con observaciones")
16 bloques con observaciones
Código
def evaluar(particiones, etiqueta): exactitudes = []for entrena, prueba in particiones: m = bosque().fit(X[entrena], y[entrena]) exactitudes.append(accuracy_score(y[prueba], m.predict(X[prueba]))) exactitudes = np.array(exactitudes)print(f"{etiqueta:<34}{exactitudes.mean():.3f} ± {exactitudes.std():.3f} "f"pliegues: {np.round(exactitudes, 3)}")return exactitudesprint("Exactitud por estrategia de validación\n"+"-"*78)acc_aleatoria = evaluar( KFold(n_splits=5, shuffle=True, random_state=42).split(X),"CV aleatoria (5 pliegues)")acc_espacial = evaluar( GroupKFold(n_splits=5).split(X, y, groups=bloque),"CV espacial por bloques (5)")brecha = (acc_aleatoria.mean() - acc_espacial.mean()) *100print("-"*78)print(f"BRECHA: {brecha:.1f} puntos porcentuales de exactitud imaginaria.")
Exactitud por estrategia de validación
------------------------------------------------------------------------------
CV aleatoria (5 pliegues) 0.884 ± 0.017 pliegues: [0.882 0.887 0.853 0.903 0.897]
CV espacial por bloques (5) 0.731 ± 0.080 pliegues: [0.734 0.805 0.669 0.832 0.618]
------------------------------------------------------------------------------
BRECHA: 15.3 puntos porcentuales de exactitud imaginaria.
Código
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12.5, 4.2), gridspec_kw={"width_ratios": [1, 1.3]})medias = [acc_aleatoria.mean(), acc_espacial.mean()]errores = [acc_aleatoria.std(), acc_espacial.std()]barras = ax1.bar(["CV aleatoria\n(optimista)", "CV espacial\n(honesta)"], medias, yerr=errores, capsize=6, color=["#e03131", "#2b8a3e"], width=0.55)for barra, valor inzip(barras, medias): ax1.text(barra.get_x() + barra.get_width() /2, valor +0.025, f"{valor:.1%}", ha="center", fontsize=11, fontweight="bold")ax1.set_ylim(0, 1.08); ax1.set_ylabel("exactitud")ax1.set_title(f"La misma exactitud, medida de dos formas\nbrecha = {brecha:.1f} puntos")ax1.annotate("", xy=(0, medias[0]), xytext=(1, medias[1]), arrowprops=dict(arrowstyle="<->", color="#495057", linewidth=1.2))ax2.boxplot([acc_aleatoria, acc_espacial], tick_labels=["aleatoria", "espacial"], widths=0.5, patch_artist=True, boxprops=dict(facecolor="#dbe4ff"), medianprops=dict(color="#c92a2a"))ax2.scatter(np.ones(len(acc_aleatoria)), acc_aleatoria, color="#e03131", zorder=3, s=28)ax2.scatter(np.full(len(acc_espacial), 2), acc_espacial, color="#2b8a3e", zorder=3, s=28)ax2.set_ylabel("exactitud por pliegue")ax2.set_title("Dispersión entre pliegues\nla CV espacial es más baja Y más variable")plt.tight_layout(); plt.show()
Qué acaba de pasar
El modelo no cambió. Los datos no cambiaron. Lo único que cambió fue la pregunta que le hicimos a la validación:
La CV aleatoria pregunta: ¿sabes predecir un punto que está a 200 metros de uno que ya viste? Y la respuesta, casi siempre, es que sí — pero eso no es generalizar, es interpolar entre vecinos.
La CV espacial pregunta: ¿sabes predecir en una región donde nunca estuviste? Que es exactamente lo que le vas a pedir al modelo cuando lo pongas en producción sobre territorio nuevo.
La brecha tiene dos causas, y las dos están en los datos:
Las observaciones están agrupadas. Puntos vecinos son casi duplicados.
La misma etiqueta no significa lo mismo en todo el territorio. En la simulación, cada clase tiene una deriva regional en su firma espectral —el equivalente a que el “cultivo” del norte sea maíz y el del sur caña. El modelo aprende la versión local y falla en la versión de al lado.
Honestidad sobre el número. La magnitud exacta de la brecha en este cuaderno depende de cuánta deriva regional metimos en la simulación, que está calibrada para que el efecto se vea en clase. Lo que no es un artefacto es el fenómeno: en teledetección real, publicaciones que reevalúan modelos de cobertura con validación espacial reportan caídas del mismo orden. La lección es la dirección y el mecanismo, no el decimal.
Y una advertencia que se aplica a cualquier dato con dimensión espacial: si tus datos tienen coordenadas y estás usando train_test_split con shuffle=True, tu métrica está inflada. No un poco: potencialmente así.
Código
# ¿Y si medimos contra la verdad de TODA la escena, incluida la que nunca vimos?# Es la evaluacion mas honesta posible, y solo es posible porque simulamos.bloques_prueba = [0, 5, 10, 15] # diagonal de la cuadriculaentrena =~np.isin(bloque, bloques_prueba)m_local = bosque().fit(X[entrena], y[entrena])pred_global = m_local.predict(todos).reshape(verdad.shape)mascara_prueba = np.zeros_like(verdad, dtype=bool)for b in bloques_prueba: r0, c0 = (b // LADO_BLOQUE) * tam, (b % LADO_BLOQUE) * tam mascara_prueba[r0:r0 + tam, c0:c0 + tam] =Trueacc_vistas = (pred_global[~mascara_prueba] == verdad[~mascara_prueba]).mean()acc_nuevas = (pred_global[mascara_prueba] == verdad[mascara_prueba]).mean()print(f"Entrenado solo con 12 de los 16 bloques:\n")print(f" exactitud en regiones YA VISTAS : {acc_vistas:.1%}")print(f" exactitud en regiones NUEVAS : {acc_nuevas:.1%}")print(f" caída : {(acc_vistas - acc_nuevas) *100:.1f} puntos")def marcar_bloques(ax):for b in bloques_prueba: r0, c0 = (b // LADO_BLOQUE) * tam, (b % LADO_BLOQUE) * tam ax.add_patch(Rectangle((EXT[0] + c0 *60, EXT[3] - (r0 + tam) *60), tam *60, tam *60, fill=False, edgecolor="black", linewidth=1.8))fig, axes = plt.subplots(1, 2, figsize=(11.5, 5))axes[0].imshow(pred_global, cmap=CMAP_CLASES, extent=EXT, vmin=0, vmax=4)marcar_bloques(axes[0])axes[0].set_title("Predicción entrenando sin los bloques marcados")err = (pred_global != verdad).astype(float)axes[1].imshow(err, cmap="Reds", extent=EXT)marcar_bloques(axes[1])axes[1].set_title("Errores: se acumulan dentro de los bloques nunca vistos")for ax in axes: ax.set_xticks([]); ax.set_yticks([]); ax.grid(False)plt.tight_layout(); plt.show()
Entrenado solo con 12 de los 16 bloques:
exactitud en regiones YA VISTAS : 85.1%
exactitud en regiones NUEVAS : 66.3%
caída : 18.8 puntos
La escalera de exactitudes
Podría objetarse que los bloques de 4×4 son arbitrarios. Lo son. Así que midamos con cuatro particiones distintas —aleatoria, por conglomerado de muestreo (los de DBSCAN de la sección 15), por bloques de 4×4 y por bloques de 2×2— y grafiquemos la exactitud contra algo medible: qué tan lejos queda el punto de prueba del dato de entrenamiento más cercano.
Código
from sklearn.model_selection import cross_val_scorefrom scipy.spatial import cKDTreegrupos_db = etiquetas_db.copy()es_ruido = grupos_db ==-1# Cada punto marcado como ruido es un grupo de una sola observaciongrupos_db[es_ruido] = grupos_db.max() +1+ np.arange(es_ruido.sum())bloque_grande = (filas // (n_px //2)) *2+ (cols // (n_px //2)) # cuadricula 2x2def evaluar_con_distancia(particiones):'''Exactitud por pliegue Y distancia al vecino de entrenamiento mas cercano.''' exactitudes, distancias = [], []for entrena, prueba in particiones: m = bosque().fit(X[entrena], y[entrena]) exactitudes.append(accuracy_score(y[prueba], m.predict(X[prueba]))) d, _ = cKDTree(coords[entrena]).query(coords[prueba]) distancias.append(np.median(d))return np.mean(exactitudes), np.mean(distancias)estrategias = [ ("aleatoria", KFold(5, shuffle=True, random_state=42).split(X), "#e03131"), ("por conglomerado", GroupKFold(5).split(X, y, groups=grupos_db), "#f59f00"), ("bloques 4x4", GroupKFold(5).split(X, y, groups=bloque), "#1c7ed6"), ("bloques 2x2", GroupKFold(4).split(X, y, groups=bloque_grande), "#2b8a3e"),]print(f"{'partición':20}{'exactitud':>10}{'distancia mediana al entrenamiento':>36}")print("-"*70)puntos = []for nombre, particiones, color in estrategias: acc, dist = evaluar_con_distancia(particiones) puntos.append((nombre, acc, dist, color))print(f"{nombre:20}{acc:10.3f}{dist:32,.0f} m")
partición exactitud distancia mediana al entrenamiento
----------------------------------------------------------------------
aleatoria 0.884 432 m
por conglomerado 0.820 1,752 m
bloques 4x4 0.731 4,588 m
bloques 2x2 0.614 9,994 m
Código
fig, ax = plt.subplots(figsize=(7.5, 4.6))for nombre, acc, dist, color in puntos: ax.scatter(dist /1000, acc, s=140, color=color, zorder=3, edgecolor="white") ax.annotate(f"{nombre}\n{acc:.1%}", (dist /1000, acc), textcoords="offset points", xytext=(0, 14), ha="center", fontsize=8)xs = np.array([p[2] for p in puntos]) /1000ys = np.array([p[1] for p in puntos])orden = np.argsort(xs)ax.plot(xs[orden], ys[orden], color="#adb5bd", linestyle="--", zorder=1)ax.set_xlabel("distancia mediana del punto de prueba al entrenamiento (km)")ax.set_ylabel("exactitud reportada")ax.set_ylim(min(ys) -0.08, max(ys) +0.09)ax.set_title("La exactitud que reportas es una función de\ncuán lejos pusiste el conjunto de prueba")plt.tight_layout(); plt.show()
Este es el resultado que conviene llevarse de todo el cuaderno: no existe “la” exactitud del modelo. Existe una exactitud por cada distancia de predicción, y la curva baja monótonamente.
Con partición aleatoria, la mediana de la distancia al vecino de entrenamiento es de unos pocos cientos de metros: se le está preguntando al modelo por un sitio que prácticamente ya vio. A medida que alejamos la prueba, la métrica cae.
Entonces la pregunta correcta no es ¿qué exactitud tiene mi modelo? sino:
¿A qué distancia de mis datos de entrenamiento lo voy a usar?
Si vas a predecir en las mismas parcelas donde muestreaste, la CV aleatoria no está mintiendo. Si vas a extrapolar a otro departamento, la cifra honesta es la del extremo derecho de la gráfica. Elegir la partición es declarar para qué sirve el modelo.
La regla práctica: si tus datos tienen coordenadas, la unidad de partición no es la observación, es el lugar.
18. Aprendizaje no supervisado
No siempre hay etiquetas. La alternativa clásica en teledetección es la clasificación no supervisada: agrupar los píxeles por similitud espectral con k-means y después interpretar qué es cada grupo mirando su firma.
Es el flujo que se usa cuando no hay presupuesto para trabajo de campo, y sigue siendo el primer paso razonable ante una imagen de la que no se sabe nada.
Código
from sklearn.cluster import MiniBatchKMeansfrom sklearn.preprocessing import StandardScalerK_GRUPOS =5espectral = escena.reshape(6, -1).T # (pixeles, bandas)escalador = StandardScaler().fit(espectral[::50])km = MiniBatchKMeans(n_clusters=K_GRUPOS, random_state=42, n_init=10, batch_size=4096).fit(escalador.transform(espectral[::50]))grupos = km.predict(escalador.transform(espectral)).reshape(verdad.shape)# El paso que NADIE puede saltarse: interpretar cada grupo.# Aqui hacemos trampa a proposito (usamos la verdad) para poder medir; en la# vida real se interpreta mirando la firma espectral de cada centroide.traduccion = {}for g inrange(K_GRUPOS): clases_en_grupo = verdad[grupos == g] traduccion[g] = np.bincount(clases_en_grupo, minlength=5).argmax()grupos_traducidos = np.vectorize(traduccion.get)(grupos)acc_kmeans = (grupos_traducidos == verdad).mean()print("Interpretación de cada grupo (por mayoría contra la verdad):")for g, k in traduccion.items(): tam = (grupos == g).mean() pureza = (verdad[grupos == g] == k).mean()print(f" grupo {g} -> {CLASES[k]:14}{tam:6.1%} de la escena, pureza {pureza:.1%}")print(f"\nExactitud del k-means traducido : {acc_kmeans:.1%}")print(f"Exactitud del Random Forest : {exactitud_total:.1%} <- con etiquetas")
Interpretación de cada grupo (por mayoría contra la verdad):
grupo 0 -> cultivo 44.8% de la escena, pureza 61.7%
grupo 1 -> cultivo 20.0% de la escena, pureza 57.3%
grupo 2 -> suelo desnudo 7.0% de la escena, pureza 82.3%
grupo 3 -> bosque 20.0% de la escena, pureza 57.7%
grupo 4 -> agua 8.1% de la escena, pureza 45.6%
Exactitud del k-means traducido : 60.1%
Exactitud del Random Forest : 82.5% <- con etiquetas
Código
fig, axes = plt.subplots(1, 3, figsize=(15.5, 5))axes[0].imshow(verdad, cmap=CMAP_CLASES, extent=EXT, vmin=0, vmax=4)axes[0].set_title("Cobertura real")axes[1].imshow(grupos, cmap="tab10", extent=EXT, vmin=0, vmax=9)axes[1].set_title(f"k-means con k={K_GRUPOS}\n(sin usar ninguna etiqueta)")axes[2].imshow(grupos_traducidos, cmap=CMAP_CLASES, extent=EXT, vmin=0, vmax=4)axes[2].set_title(f"Grupos ya interpretados ({acc_kmeans:.1%})")for ax in axes: ax.set_xticks([]); ax.set_yticks([]); ax.grid(False)plt.tight_layout(); plt.show()# Las firmas de los centroides: esto es lo que un analista mira para interpretarcentros = escalador.inverse_transform(km.cluster_centers_)fig, ax = plt.subplots(figsize=(7, 4))for g inrange(K_GRUPOS): ax.plot(LONGITUDES, centros[g], "o-", linewidth=2, label=f"grupo {g} -> parece {CLASES[traduccion[g]]}")ax.set_xlabel("Longitud de onda (nm)"); ax.set_ylabel("Reflectancia")ax.set_title("Firma espectral de cada centroide\nasí se interpreta una clasificación no supervisada")ax.legend(fontsize=7); plt.tight_layout(); plt.show()
Dos lecciones que se ven mejor aquí que en cualquier explicación:
k-means encuentra estructura, no significado. Mira la traducción de los grupos: aunque pedimos exactamente 5 —el número real de clases— los grupos no se reparten uno por clase. Alguna clase se queda sin grupo propio y otra se parte en dos, porque k-means sigue la varianza espectral, no nuestras categorías. Lo urbano, que es minoritario y espectralmente parecido al suelo desnudo, es siempre el primero en desaparecer.
El número de grupos es una decisión tuya. Con \(k\) mayor que el número de clases reales, los grupos se subdividen y luego se re-agrupan a mano; es una práctica común y perfectamente válida.
Y la comparación con el Random Forest pone precio a las etiquetas: lo que cuestan las visitas de campo se traduce, aquí, en puntos de exactitud.
19. Ingeniería de variables espaciales
Hay una frase que se repite mucho en este campo: el feature engineering espacial suele mejorar más un modelo que cambiar de algoritmo. Vamos a comprobarla — y sobre todo a comprobar bajo qué validación mejora.
Añadimos cuatro variables construidas con el álgebra de mapas de la sección 11:
Variable
Familia
Qué aporta
NDVI focal 5×5
focal
el contexto: qué hay alrededor, no solo aquí
Textura del NDVI 5×5
focal
heterogeneidad: cultivo en mosaico vs. bosque uniforme
Distancia al agua
global
proximidad a ríos y lagos
Distancia a lo construido
global
gradiente urbano–rural
La trampa que hay que evitar. Las cuatro se calculan a partir de la imagen, no del mapa de verdad. Una variable como “distancia al bosque más cercano” calculada sobre las etiquetas sería fuga descarada: en producción no tendrías ese mapa —es justo lo que estás tratando de predecir—.
Código
# Las cuatro capas derivadas, ya calculadas en la seccion 11capas_extra = np.stack([media_focal, textura, dist_agua /1000.0, dist_urbano /1000.0])NOMBRES_EXTRA = ["NDVI_focal", "textura", "dist_agua_km", "dist_urbano_km"]X_extra = np.column_stack([X, capas_extra[:, filas, cols].T])print(f"X original : {X.shape[1]} variables")print(f"X ampliada : {X_extra.shape[1]} variables ({', '.join(NOMBRES_EXTRA)})")
X original : 9 variables
X ampliada : 13 variables (NDVI_focal, textura, dist_agua_km, dist_urbano_km)
Código
resultados = []for etiqueta, matriz in [("solo espectrales (9)", X), ("+ espaciales (13)", X_extra)]: ale = cross_val_score(bosque(), matriz, y, cv=KFold(5, shuffle=True, random_state=42)) esp = cross_val_score(bosque(), matriz, y, cv=GroupKFold(5), groups=bloque) resultados.append({"variables": etiqueta,"CV aleatoria": round(ale.mean(), 3),"CV espacial": round(esp.mean(), 3),"brecha (pp)": round((ale.mean() - esp.mean()) *100, 1)})tabla = pd.DataFrame(resultados)print(tabla.to_string(index=False))print()d_ale = tabla["CV aleatoria"][1] - tabla["CV aleatoria"][0]d_esp = tabla["CV espacial"][1] - tabla["CV espacial"][0]print(f"Las variables espaciales cambian la CV aleatoria en {d_ale *100:+.1f} pp")print(f"y la CV espacial en {d_esp *100:+.1f} pp.")print("\nLa que importa es la segunda: es la que estima el desempeño en territorio nuevo.")
variables CV aleatoria CV espacial brecha (pp)
solo espectrales (9) 0.884 0.731 15.3
+ espaciales (13) 0.921 0.726 19.4
Las variables espaciales cambian la CV aleatoria en +3.7 pp
y la CV espacial en -0.5 pp.
La que importa es la segunda: es la que estima el desempeño en territorio nuevo.
Y aquí está la sorpresa —o la lección, según cómo se mire.
Las cuatro variables espaciales suben la exactitud aleatoria varios puntos y no mejoran la espacial. Traducido: no aprendieron nada nuevo sobre las coberturas; aprendieron dónde están las cosas en este mapa concreto.
Tiene todo el sentido. La media focal y la textura del NDVI resumen «cómo es el vecindario», y en la CV aleatoria el vecindario del punto de prueba está en el entrenamiento. Es fuga espacial otra vez, ahora escondida dentro de una variable en vez de dentro de la partición. Las variables de distancia son aún más directas: son, en el fondo, una forma suave de codificar la coordenada.
Ojo con la conclusión, que no es «no uses variables espaciales»:
Si vas a predecir dentro de la misma zona —rellenar huecos de nube, mapear las parcelas no visitadas de la misma finca— estas variables sí ayudan de verdad, y la CV aleatoria es la métrica adecuada.
Si vas a extrapolar a territorio nuevo, la ganancia se evapora, y solo la CV espacial te lo dice.
Lo fácil habría sido publicar la cifra aleatoria —la más alta de todo el cuaderno— y sentirse bien.
Código
modelo_extra = bosque().fit(X_extra, y)importancias_extra = pd.Series(modelo_extra.feature_importances_, index=NOMBRES_VAR + NOMBRES_EXTRA).sort_values()fig, ax = plt.subplots(figsize=(7, 4.4))colores_barras = ["#e8590c"if n in NOMBRES_EXTRA else"#1c7ed6"for n in importancias_extra.index]importancias_extra.plot.barh(ax=ax, color=colores_barras)ax.set_title("Importancia de variables — naranja = derivadas del espacio")ax.set_xlabel("importancia")plt.tight_layout(); plt.show()
Cuidado al leer esa gráfica: una variable espacial con mucha importancia no es automáticamente una buena noticia. Si la variable codifica dónde está el píxel más que qué es, el modelo aprende a memorizar el mapa —y eso funciona espectacularmente bien en la CV aleatoria y se derrumba en territorio nuevo—.
La coordenada cruda (x, y) es el caso extremo: como variable predictora es casi siempre una trampa. El modelo la usa para ubicarse en el mapa de entrenamiento, no para entender el fenómeno.
La prueba de fuego es siempre la misma: si una variable sube la métrica aleatoria y no la espacial, no aprendió nada — memorizó el mapa.
20. Deep learning sobre imágenes satelitales
Todo lo que ya sabes de redes convolucionales para imágenes aplica aquí, con tres matices que cambian el código:
Hay más de tres canales. Una CNN preentrenada en fotos espera 3; aquí tenemos 6, y en un hiperespectral, cientos. Hay que adaptar la primera capa.
Los valores no están en [0, 255]. Son reflectancias en [0, 1], y la normalización estándar de ImageNet no aplica.
El contexto geográfico y temporal importa. Un parche no es una foto suelta: tiene vecinos y tiene historia.
El cambio conceptual es de unidad de análisis: del píxel al parche. Hasta ahora cada observación era un vector de 9 o 13 números. Ahora es un tensor.
Código
RADIO =3# parche de 7x7 pixeles = 420 x 420 mLADO =2* RADIO +1acolchada = np.pad(escena, ((0, 0), (RADIO, RADIO), (RADIO, RADIO)), mode="edge")parches = np.stack([acolchada[:, f:f + LADO, c:c + LADO]for f, c inzip(filas, cols)]) # (n, 6, 7, 7)print(f"Píxel suelto : {X.shape} <- lo que usó el Random Forest")print(f"Parches : {parches.shape} <- (observaciones, canales, alto, ancho)")print(f"\nEs exactamente la forma que espera una CNN en PyTorch: (N, C, H, W).")print(f"La diferencia con una foto: C=6 en vez de 3, y los valores van de "f"{parches.min():.2f} a {parches.max():.2f}, no de 0 a 255.")
Píxel suelto : (3000, 9) <- lo que usó el Random Forest
Parches : (3000, 6, 7, 7) <- (observaciones, canales, alto, ancho)
Es exactamente la forma que espera una CNN en PyTorch: (N, C, H, W).
La diferencia con una foto: C=6 en vez de 3, y los valores van de 0.00 a 0.74, no de 0 a 255.
Código
# Un vistazo a los parches, en falso colorfig, axes = plt.subplots(2, 6, figsize=(13, 4.6))for k inrange(5): idx = np.where(y == k)[0][:2]for j, i inenumerate(idx): ax = axes[j, k] p = parches[i] ax.imshow(np.dstack([estirar(p[NIR]), estirar(p[ROJO]), estirar(p[VERDE])]), interpolation="nearest") ax.set_title(CLASES[k], fontsize=8) ax.set_xticks([]); ax.set_yticks([]); ax.grid(False)for j inrange(2): axes[j, 5].axis("off")fig.suptitle("Parches de 7×7 en falso color: la unidad de análisis de una CNN", y=1.02)plt.tight_layout(); plt.show()
¿Vale la pena el contexto?
Antes de montar una CNN conviene contestar la pregunta barata: ¿aporta algo mirar la vecindad? Lo medimos con un perceptrón multicapa sobre el parche aplanado. No es una CNN —no comparte pesos ni es invariante a traslaciones— pero sí responde si la información de contexto existe.
Y lo medimos, por supuesto, con validación espacial.
Código
from sklearn.neural_network import MLPClassifierfrom sklearn.pipeline import make_pipelineX_parches = parches.reshape(len(parches), -1) # 6*7*7 = 294 variablesdef red():return make_pipeline( StandardScaler(), MLPClassifier(hidden_layer_sizes=(64,), max_iter=400, random_state=42, early_stopping=True, n_iter_no_change=15))print("Exactitud con validación cruzada ESPACIAL por bloques\n"+"-"*58)for etiqueta, matriz, estimador in [ ("Random Forest, píxel suelto (9 var.)", X, bosque), ("Random Forest, + espaciales (13 var.)", X_extra, bosque), ("Red neuronal sobre parches 7x7 (294)", X_parches, red),]:with warnings.catch_warnings(): warnings.simplefilter("ignore") s = cross_val_score(estimador(), matriz, y, cv=GroupKFold(5), groups=bloque)print(f" {etiqueta:40}{s.mean():.3f} ± {s.std():.3f}")
Exactitud con validación cruzada ESPACIAL por bloques
----------------------------------------------------------
Random Forest, píxel suelto (9 var.) 0.731 ± 0.080
Random Forest, + espaciales (13 var.) 0.726 ± 0.107
Red neuronal sobre parches 7x7 (294) 0.826 ± 0.094
Esta vez el contexto sí aportó. La red sobre parches gana varios puntos en validación espacial —la honesta—, justo donde las variables espaciales de la sección anterior no habían ganado nada.
La diferencia entre los dos casos vale la pena entenderla, porque es sutil:
Las variables de la sección 19 resumían el vecindario en cuatro números, y dos de ellos (las distancias) eran básicamente coordenadas disfrazadas.
El parche entrega la estructura espacial completa: la disposición de los valores, no solo su promedio. Un cultivo en mosaico y un bosque uniforme pueden tener la misma media y la misma textura, y verse completamente distintos en 7×7.
Eso es exactamente lo que una CNN explota, y por qué la segmentación semántica domina el mapeo de cobertura desde hace una década.
Lo que sigue después de este experimento, en un proyecto real:
Clasificación de parches con una CNN: la misma idea, pero con convoluciones que comparten pesos y capturan patrones espaciales a varias escalas. Detecta edificios, embarcaciones, piscinas, parcelas deforestadas.
Segmentación semántica con U-Net: en vez de una etiqueta por parche, una etiqueta por píxel. Es el estándar para mapear cobertura y extraer huellas de edificios; la arquitectura encoder–decoder con conexiones de salto conserva el detalle fino que el encoder comprime.
Series temporales: un cubo (tiempo, banda, alto, ancho) con LSTM o transformers. Los cultivos se distinguen por su fenología —cómo cambia su NDVI a lo largo del año— mucho mejor que por una sola fecha.
Y una advertencia que vale doble con deep learning: todo lo de la sección 17 sigue aplicando. Partir parches al azar es peor que partir píxeles al azar, porque dos parches vecinos literalmente comparten píxeles. La partición debe ser por regiones o por escenas completas.
21. Trampas del análisis espacial
MAUP
El problema de la unidad de área modificable dice que los resultados cambian según el tamaño y el trazo de las unidades espaciales. No es una sutileza teórica: con exactamente los mismos píxeles, la relación entre dos variables puede cambiar de magnitud e incluso de signo.
Medimos la correlación entre la proporción urbana y el NDVI medio, agregando a distintas escalas.
Código
def agregar_en_cuadricula(lado_celda_px):'''Agrega la escena en celdas cuadradas y devuelve (prop_urbana, ndvi_medio).''' n = verdad.shape[0] // lado_celda_px * lado_celda_px urb = (verdad[:n, :n] ==3).reshape( n // lado_celda_px, lado_celda_px, -1, lado_celda_px).mean(axis=(1, 3)) veg = ndvi[:n, :n].reshape( n // lado_celda_px, lado_celda_px, -1, lado_celda_px).mean(axis=(1, 3))return urb.ravel(), veg.ravel()escalas = [(8, "0.5 km"), (32, "2 km"), (128, "7.7 km"), (256, "15 km")]filas_tabla = []fig, axes = plt.subplots(1, 5, figsize=(17.5, 3.5))for ax, (lado, etiqueta) inzip(axes, escalas): u, v = agregar_en_cuadricula(lado) r = np.corrcoef(u, v)[0, 1] filas_tabla.append({"unidad": f"cuadrícula {etiqueta}", "n": len(u), "r": r}) ax.scatter(u, v, s=8, alpha=0.4, color="#1c7ed6") ax.set_title(f"cuadrícula {etiqueta}\nn={len(u)} r={r:+.2f}", fontsize=9) ax.set_xlabel("proporción urbana")axes[0].set_ylabel("NDVI medio")# Y ahora con las unidades administrativas realesr_muni = np.corrcoef(munis_v["prop_urbano"], munis_v["ndvi_medio"])[0, 1]filas_tabla.append({"unidad": "municipios reales", "n": len(munis_v), "r": r_muni})axes[4].scatter(munis_v["prop_urbano"], munis_v["ndvi_medio"], s=26, alpha=0.75, color="#e03131")axes[4].set_title(f"municipios reales\nn={len(munis_v)} r={r_muni:+.2f}", fontsize=9)axes[4].set_xlabel("proporción urbana")fig.suptitle("MAUP — la misma relación, cinco unidades de agregación", y=1.04)plt.tight_layout(); plt.show()pd.DataFrame(filas_tabla)
unidad
n
r
0
cuadrícula 0.5 km
16384
-0.106369
1
cuadrícula 2 km
1024
-0.066803
2
cuadrícula 7.7 km
64
-0.012309
3
cuadrícula 15 km
16
0.025299
4
municipios reales
58
-0.038466
La correlación se fortalece al agregar en celdas más grandes: promediar elimina el ruido de píxel y deja solo la señal de fondo. Es un efecto conocido —la agregación infla las correlaciones— y explica por qué un estudio a nivel departamento puede reportar una relación mucho más fuerte que el mismo estudio a nivel de manzana.
Los municipios reales no caen en la tendencia de las cuadrículas, y esa es la otra mitad del MAUP: no solo importa el tamaño de la unidad, sino dónde se trazan las fronteras. Las fronteras administrativas no se dibujaron pensando en tu análisis.
La consecuencia práctica: la unidad de agregación es una decisión metodológica, tan importante como la elección del modelo, y hay que declararla.
Falacia ecológica
El MAUP tiene un pariente cercano y más peligroso: concluir algo sobre individuos a partir de datos agregados por zona. Aquí lo vemos con píxeles, que son los “individuos” de esta escena.
Código
# A nivel de MUNICIPIO: relacion entre proporcion urbana y NDVI medior_agregado = np.corrcoef(munis_v["prop_urbano"], munis_v["ndvi_medio"])[0, 1]# A nivel de PIXEL, dentro de cada municipio por separadocorrelaciones_internas = []for i inrange(len(munis_v)): dentro = zonas == (i +1)if dentro.sum() <500:continue es_urbano = (verdad[dentro] ==3).astype(float) ndvi_px = ndvi[dentro]if es_urbano.std() >0: correlaciones_internas.append(np.corrcoef(es_urbano, ndvi_px)[0, 1])print(f"Correlación urbano–NDVI entre MUNICIPIOS (agregado): r = {r_agregado:+.2f}")print(f"Correlación urbano–NDVI entre PÍXELES (individual) : r = "f"{np.mean(correlaciones_internas):+.2f} (promedio de {len(correlaciones_internas)} municipios)")fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 4.2))ax1.scatter(munis_v["prop_urbano"], munis_v["ndvi_medio"], s=30, color="#e03131")ax1.set_xlabel("proporción urbana del municipio"); ax1.set_ylabel("NDVI medio")ax1.set_title(f"Nivel agregado: {len(munis_v)} municipios\nr = {r_agregado:+.2f}")ax2.hist(correlaciones_internas, bins=20, color="#1c7ed6", edgecolor="white")ax2.axvline(r_agregado, color="#e03131", linewidth=2, label="r agregado")ax2.axvline(np.mean(correlaciones_internas), color="#1c7ed6", linestyle="--", linewidth=2, label="r individual medio")ax2.set_xlabel("correlación a nivel de píxel, dentro de cada municipio")ax2.set_title("Nivel individual: una correlación por municipio")ax2.legend(fontsize=8)plt.tight_layout(); plt.show()
Correlación urbano–NDVI entre MUNICIPIOS (agregado): r = -0.04
Correlación urbano–NDVI entre PÍXELES (individual) : r = -0.34 (promedio de 41 municipios)
Las dos correlaciones no tienen por qué coincidir, y en general no coinciden: son preguntas distintas medidas sobre unidades distintas. Que un municipio con más superficie urbana tenga menos vegetación no permite afirmar nada sobre qué pasa dentro de un municipio concreto.
El ejemplo canónico: que un municipio tenga alto ingreso y alta mortalidad no implica que los ricos mueran más. Puede ser exactamente al revés — que los pobres de los municipios ricos mueran mucho.
Efecto de borde
Las unidades del borde del área de estudio tienen menos vecinos, simplemente porque el vecino de al lado quedó fuera del recorte. Todo estadístico de vecindad —Moran, LISA, Gi*, cualquier operación focal— queda sesgado ahí.
Código
# Marcamos los municipios que tocan el borde de la ventana de estudiofrontera = gpd.GeoSeries([box(*limites).boundary], crs=munis_v.crs).iloc[0]en_borde = munis_v.geometry.apply(lambda g: g.intersects(frontera)).valuesn_vecinos = W.sum(axis=1)print(f"Municipios en el borde de la ventana : {en_borde.sum()} de {len(munis_v)}")print(f" vecinos promedio en el borde : {n_vecinos[en_borde].mean():.1f}")print(f" vecinos promedio en el interior: {n_vecinos[~en_borde].mean():.1f}")# Cuanto cambia el indice de Moran si excluimos el bordeI_todos = moran_i(munis_v["ndvi_medio"].values, Wz)interior =~en_bordeW_int = estandarizar_filas(W[np.ix_(interior, interior)])I_interior = moran_i(munis_v["ndvi_medio"].values[interior], W_int)print(f"\nMoran con todos los municipios : I = {I_todos:+.3f}")print(f"Moran solo con el interior : I = {I_interior:+.3f}")fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 4.6))munis_v.assign(borde=en_borde).plot(column="borde", ax=ax1, cmap="Set1", edgecolor="white", linewidth=0.4, legend=True, legend_kwds={"fontsize": 8})ax1.set_title("Municipios que tocan el borde de la ventana")ax1.set_xticks([]); ax1.set_yticks([]); ax1.grid(False)ax2.boxplot([n_vecinos[en_borde], n_vecinos[~en_borde]], tick_labels=["borde", "interior"], widths=0.5, patch_artist=True, boxprops=dict(facecolor="#dbe4ff"), medianprops=dict(color="#c92a2a"))ax2.set_ylabel("número de vecinos")ax2.set_title("El borde tiene menos vecinos por construcción")plt.tight_layout(); plt.show()
Municipios en el borde de la ventana : 27 de 58
vecinos promedio en el borde : 4.3
vecinos promedio en el interior: 6.2
Moran con todos los municipios : I = +0.267
Moran solo con el interior : I = +0.334
Las mitigaciones habituales: usar una zona de amortiguamiento (traer datos más allá del área de estudio aunque no se reporten), corregir los pesos de las unidades de borde, o simplemente declarar que los resultados del borde son menos confiables.
Y la cuarta trampa, la que no se arregla con estadística: el sesgo de cobertura. Las zonas rurales y pobres están sistemáticamente peor mapeadas, y esa ausencia de datos se confunde con ausencia del fenómeno. Un modelo entrenado sobre cobertura desigual perpetúa la desigualdad.
22. Ética
Cerramos con lo que más importa fuera del aula. La ubicación es un dato personal particularmente sensible, y por seis razones distintas:
Reidentificación. Bastan cuatro puntos espacio-temporales para identificar de forma única a la mayoría de las personas en un conjunto de datos de movilidad. Casa, trabajo y dos paradas más ya son una huella digital. Anonimizar quitando el nombre no es suficiente.
Agregación y desplazamiento. Publicar datos sensibles exige agregar a unidades suficientemente grandes o desplazar las coordenadas de forma controlada. Publicar la ubicación exacta de un caso de enfermedad o de una vivienda es revelar a una persona.
Vigilancia. El rastreo masivo de ubicación plantea problemas de derechos que la técnica no resuelve. Que se pueda hacer no lo hace legítimo.
Sesgo de datos. Los mapas colaborativos están mucho más completos en zonas ricas y urbanas. Un modelo entrenado sobre esa cobertura desigual perpetúa la desigualdad: lo que no está mapeado no existe para el algoritmo.
El poder del mapa. Todo mapa es una decisión sobre qué se muestra y qué se omite. Las fronteras, los nombres y las categorías de un mapa nunca son neutrales — y quien dibuja el mapa ejerce ese poder, sepa o no que lo tiene.
Consentimiento. La ubicación se recolecta con frecuencia de forma pasiva desde teléfonos y aplicaciones, sin que la persona entienda el alcance de lo que está cediendo.
Vamos a ponerle números al primero. Supongamos que las observaciones no son parcelas sino casos de una enfermedad. Publicar el mapa de puntos exactos es publicar dónde vive cada persona.
Código
rng = np.random.default_rng(7)casos = obs.sample(320, random_state=7).copy()# Estrategia 1: desplazamiento aleatorio controlado (jitter)RADIO =1500# metrosangulo = rng.uniform(0, 2* np.pi, len(casos))radio = RADIO * np.sqrt(rng.uniform(0, 1, len(casos)))casos_desplazados = casos.copy()casos_desplazados["geometry"] = gpd.points_from_xy( casos.geometry.x + radio * np.cos(angulo), casos.geometry.y + radio * np.sin(angulo), crs=casos.crs)# Estrategia 2: agregacion a la unidad administrativapor_muni = gpd.sjoin(casos, munis_v[["municipio", "geometry"]], predicate="within")tasas = por_muni.groupby("municipio").size().rename("casos").reset_index()mapa_tasas = munis_v.merge(tasas, on="municipio", how="left")mapa_tasas["casos"] = mapa_tasas["casos"].fillna(0)mapa_tasas["tasa"] = mapa_tasas["casos"] / mapa_tasas["area_km2"]fig, axes = plt.subplots(1, 3, figsize=(15, 4.8))munis_v.boundary.plot(ax=axes[0], color="#adb5bd", linewidth=0.4)casos.plot(ax=axes[0], color="#e03131", markersize=6)axes[0].set_title("NO PUBLICAR: ubicación exacta\ncada punto es una persona")munis_v.boundary.plot(ax=axes[1], color="#adb5bd", linewidth=0.4)casos_desplazados.plot(ax=axes[1], color="#f59f00", markersize=6)axes[1].set_title(f"Mitigación 1: desplazamiento aleatorio\nradio de {RADIO} m")mapa_tasas.plot(column="tasa", cmap="OrRd", legend=True, ax=axes[2], edgecolor="white", linewidth=0.4, legend_kwds={"label": "casos por km²", "shrink": 0.7})axes[2].set_title("Mitigación 2: agregación\ncasos por km² (normalizado)")for ax in axes: ax.set_xticks([]); ax.set_yticks([]); ax.grid(False)plt.tight_layout(); plt.show()d = np.hypot(casos_desplazados.geometry.x - casos.geometry.x, casos_desplazados.geometry.y - casos.geometry.y)print(f"Desplazamiento aplicado: mediana {d.median():.0f} m, máximo {d.max():.0f} m")print(f"Agregación: {len(casos)} ubicaciones individuales -> "f"{(mapa_tasas['casos'] >0).sum()} valores municipales")
Desplazamiento aplicado: mediana 1061 m, máximo 1499 m
Agregación: 320 ubicaciones individuales -> 39 valores municipales
Ninguna de las dos mitigaciones es gratis:
El desplazamiento conserva el patrón general pero destruye el análisis a escala fina, y si el radio es pequeño no protege gran cosa.
La agregación protege bien, pero introduce de lleno el MAUP y la falacia ecológica de la sección anterior.
Ese es el intercambio real de la privacidad geoespacial: utilidad contra protección, y no existe un punto óptimo universal. Lo que sí existe es la obligación de elegirlo explícitamente y declararlo, en vez de publicar los puntos crudos porque era lo más fácil.
Y una advertencia sobre el desplazamiento aleatorio: no protege si el atacante tiene el dato repetido. Si publicas la misma persona desplazada varias veces, el promedio de los desplazamientos converge a su ubicación real. El ruido hay que aplicarlo una vez y de forma consistente, no en cada publicación.
23. Aplicaciones y cierre
Dónde se usa todo esto
Dominio
Qué se hace con lo de este cuaderno
Agricultura de precisión
monitoreo con NDVI, estimación de rendimiento, detección temprana de estrés hídrico, aplicación variable de insumos
Gestión de desastres
mapeo de inundaciones con SAR (ve a través de la nube), evaluación de daños tras un sismo, zonas de riesgo por lahares
Ambiente y clima
deforestación, cambio de cobertura, islas de calor urbanas, monitoreo de glaciares
Salud pública
mapeo de brotes (Gi*, KDE), accesibilidad a servicios, identificación de zonas desatendidas
Planificación urbana
crecimiento de la mancha urbana, catastro, transporte, ubicación óptima de servicios
Logística y retail
optimización de rutas, áreas de cobertura, selección de sitios
Seguridad alimentaria
alerta temprana combinando lluvia, NDVI y precios
En Guatemala son casos vivos el monitoreo de la deforestación en el Petén, la vigilancia del complejo volcánico de Fuego y Pacaya con InSAR, el mapeo de riesgo por deslizamientos en la cuenca del Amatitlán y el LiDAR arqueológico que reveló miles de estructuras mayas bajo el dosel de la selva.
El recorrido completo
Sección
Concepto
2
La ubicación como dato: geometría + atributos + tiempo
3
CRS, datums y proyecciones; medir en grados es un error
Consulta espacial, buffer, disolución, unión espacial, estadística zonal
11
Álgebra de mapas focal y global
12
Análisis encadenado de idoneidad
13
Cartografía honesta: normalizar y clasificar
14
Autocorrelación espacial: Moran, LISA, Getis-Ord
15
Patrones de puntos: DBSCAN, KDE e interpolación IDW
16
Clasificación supervisada de cobertura terrestre
17
Fuga de información espacial
18
Aprendizaje no supervisado: agrupamiento espectral
19
Ingeniería de variables espaciales
20
Deep learning sobre imágenes satelitales
21
MAUP, falacia ecológica y efecto de borde
22
Privacidad de la ubicación
Las cuatro ideas que valen más que el código
Un CRS no es metadato decorativo. Es la diferencia entre un número y un sinsentido, y combinar capas sin verificarlo es la primera fuente de errores.
La autocorrelación espacial rompe el supuesto de independencia. Aparece tres veces en este cuaderno: cuando la medimos (Moran, Gi*), cuando la explotamos (interpolación, variables de contexto) y cuando nos habría engañado (validación cruzada).
La validación espacial no es opcional. Si tus datos tienen coordenadas, la unidad de partición no es la observación, es el lugar. Todo lo demás es una métrica inflada.
La ubicación identifica personas. El cuidado ético no es un apéndice al final del análisis, es una decisión de diseño desde la primera línea.
Las herramientas que usamos, y las que faltan
Todo este cuaderno corre con GeoPandas (vectorial), Rasterio (ráster), Shapely (geometría), pyproj (CRS), SciPy (álgebra de mapas focal y global) y scikit-learn (modelos). Deliberadamente implementamos Moran y Getis-Ord a mano, pero en producción se usa PySAL, que trae eso y mucho más.
Lo que queda fuera y conviene conocer:
rioxarray / xarray: cubos de datos multidimensionales (x, y, tiempo, banda). Es la forma cómoda de trabajar la serie temporal de la sección 8.
Google Earth Engine: catálogo planetario con cómputo en la nube. Evita descargar terabytes; es como se hace teledetección a escala nacional.
QGIS: SIG de escritorio libre, con consola de Python. Imprescindible para inspeccionar datos antes de programar nada.
PostGIS: la extensión espacial de PostgreSQL. Consultas SQL con predicados geométricos e índices espaciales; es la forma seria de guardar datos geoespaciales, no en carpetas de shapefiles.
Folium / Kepler.gl / plotly: mapas interactivos para comunicar resultados.
Retos para el estudiante
Cambia la semilla del generador y vuelve a correr todo. ¿Se mantiene la brecha entre validación aleatoria y espacial? ¿Cuánto varía?
Sube y baja la deriva regional (AMPLITUD_DERIVA en src/build_dataset.py) y grafica la brecha en función de ella. ¿Dónde deja de importar la validación espacial?
Cambia el tamaño de bloque de la CV espacial (2×2, 4×4, 8×8). ¿Qué pasa con la brecha cuando los bloques son muy pequeños? ¿Por qué?
Mete las coordenadas como variables (x, y) en el Random Forest de la sección 19. Mide la exactitud con CV aleatoria y con CV espacial. La diferencia entre ambas es la definición operativa de “memorizar el mapa”.
Clasifica con la serie temporal de la sección 8 en vez de una sola fecha: 4 fechas × 6 bandas = 24 variables. ¿Cuánto ayuda la fenología?
Sustituye el índice de Moran por vecinos de k-más-cercanos en vez de contigüidad. ¿Cambia la conclusión sobre el patrón del NDVI?
Repite el análisis de idoneidad de la sección 12 cambiando los umbrales (radio de cobertura, franja de vías, umbral de demanda). ¿Cuánto cambia la lista de municipios prioritarios? Eso es un análisis de sensibilidad.
Trae datos reales. Descarga una escena Sentinel-2 de Copernicus para la misma zona y repite las secciones 5 a 9 con reflectancia auténtica.
