Análisis Raster con Python para Dummies

Teledetección: técnica para obtener información a distancia. Proceso: Sol → Reflectancia (superficie) → Radiancia (sensor). Sensores pasivos (dependen del sol) vs activos (generan su propia energía). Una imagen satelital es una matriz numérica donde cada píxel representa radiancia. Resolución radiométrica: 8 bits = 256 valores. Resolución espacial: Landsat = 30m/píxel.

Librerías Python

import numpy as np
import matplotlib.pyplot as plt
from PIL import Image

Cargar y Analizar Imagen

imagen = Image.open("ruta.tif")
matriz = np.asarray(imagen)

print(np.max(matriz), np.min(matriz))
print(np.mean(matriz), np.std(matriz))

Visualización e Histograma

plt.imshow(matriz, cmap='nipy_spectral')
plt.colorbar()
plt.show()

plt.hist(matriz.ravel(), bins=50)
plt.show()

Usar clim para ajustar contraste y enfocar el rango de datos relevante.

Una imagen satelital es una matriz donde cada píxel representa radiancia. Usando el modelo RGB asignamos tres bandas a los canales Rojo, Verde y Azul para generar imágenes compuestas a color.

Cargar y Apilar Bandas

banda4 = np.asarray(Image.open(r"B4.tif"))
banda3 = np.asarray(Image.open(r"B3.tif"))
banda2 = np.asarray(Image.open(r"B2.tif"))

rgb = np.dstack((banda4, banda3, banda2))

Normalización

Matplotlib espera valores 0-1. Dividir cada banda por su valor máximo:

b4_norm = banda4 / np.max(banda4)
rgb_norm = np.dstack((b4_norm, b3_norm, b2_norm))
plt.imshow(rgb_norm)

Composiciones de Color (Landsat 9 OLI-2)

Color Natural (4-3-2): Vegetación verde, suelos marrones. Falso Color Estándar (5-4-3): NIR al canal rojo → vegetación en rojo vibrante. Agricultura (6-5-2): SWIR → vegetación verde brillante, suelos magenta. Geología (7-6-2): Estructuras geológicas. Urbano (7-6-4): Discriminar infraestructura.

Firmas Espectrales

Cada material interactúa de manera única con la radiación. El agua absorbe en NIR (aparece oscura). La vegetación refleja alto en NIR (estructura celular). Conocer las firmas permite seleccionar bandas estratégicamente para maximizar contraste.

Asignación de categorías temáticas (agua, bosques, ciudades) a cada píxel usando Machine Learning. El modelo aprende firmas espectrales (curvas de reflectancia únicas) a partir de muestras de entrenamiento.

Librerías

import numpy as np
import rasterio
from sklearn.model_selection import train_test_split
from sklearn import preprocessing
from sklearn.neighbors import KNeighborsClassifier

Preparación de Datos

Las sub-imágenes se "aplanan" (reshape) para convertirlas en vectores de características espectrales:

imagen_aplanada = imagen.reshape(imagen.shape[0], -1).T
pixeles.extend(imagen_aplanada)
etiquetas.extend([cobertura] * imagen_aplanada.shape[0])

Entrenamiento y Clasificación

Se divide en entrenamiento (80%) y prueba (20%) con train_test_split. Se estandarizan los datos con StandardScaler.

Algoritmos

• KNN: Clasifica según los "k" vecinos más cercanos.
• Random Forest: Conjunto de árboles de decisión (más robusto).
• SVM: Busca el hiperplano óptimo de separación.

rf = RandomForestClassifier(n_estimators=100)
rf.fit(pixeles_entrenamiento, etiquetas_entrenamiento)
resultado = rf.predict(imagen_a_clasificar).reshape(alto, ancho)

Visualizar resultados con ListedColormap en matplotlib.

Técnica para agrupar píxeles según su similitud espectral sin muestras previas. Se utiliza el algoritmo K-Means para identificar patrones y segmentar la imagen automáticamente en K clusters definidos por el usuario.

Librerías y Preparación

from sklearn.cluster import KMeans
from sklearn.preprocessing import StandardScaler
import rasterio

Es necesario reestructurar la imagen de (Bandas, Alto, Ancho) a un arreglo bidimensional (Píxeles, Bandas) y normalizar los datos con StandardScaler.

Ejecución de K-Means

kmeans = KMeans(n_clusters=5, random_state=42)
etiquetas = kmeans.fit_predict(datos_escalados)

# Redimensionar al tamaño original
clasificacion = etiquetas.reshape(alto, ancho)

Interpretación

El algoritmo agrupa respuestas espectrales similares pero no asigna nombres temáticos (agua, vegetación, etc.). La interpretación semántica final es responsabilidad del analista, quien asocia cada cluster con una cobertura real.

Aritmética de bandas para resaltar características terrestres (vegetación, humedad, minerales) mediante índices espectrales. Usaremos Rasterio para leer imágenes y normalizar los valores (0-1).

Normalización de Bandas (Landsat 9)

banda_rojo = imagen.read(4) / imagen.read(4).max()
banda_nir = imagen.read(5) / imagen.read(5).max()

Índice NDVI (Vegetación)

Relaciona el Infrarrojo Cercano y el Rojo. Valores cercanos a 1 indican vegetación vigorosa.

ndvi = (banda_nir - banda_rojo) / (banda_nir + banda_rojo)
plt.imshow(ndvi, cmap='nipy_spectral')

Índices Hídricos (MNDWI)

El agua refleja en el verde y absorbe en el infrarrojo. Valores altos resaltan agua o humedad.

mndwi = (banda_verde - banda_nir) / (banda_verde + banda_nir)
plt.imshow(mndwi, cmap='ocean')

Aplicaciones Geológicas (Arcilla)

Identificación de minerales arcillosos mediante la razón entre bandas SWIR.

indice_arcilla = banda_swir1 / banda_swir2
plt.imshow(indice_arcilla, cmap='terrain')

Análisis de cambios multitemporal aplicado a la tragedia de Armero (1985) usando imágenes Landsat 5 (escenarios previo y posterior). Flujo completo con Rasterio, NumPy y Matplotlib.

Exploración y Normalización

Lectura de las 6 bandas espectrales y normalización (escala 0-1) para facilitar la visualización y operaciones matemáticas.

banda_norm = banda / banda.max()

Composiciones RGB y Comparación

Uso de np.dstack para crear composiciones de Falso Color Estándar y composiciones hídricas que resaltan flujos de lodo y escombros. Comparación visual directa usando subplots.

Detección de Cambios con MNDWI

Cálculo del Índice Diferencial de Agua Normalizado Modificado:

mndwi = (verde - infrarrojo) / (verde + infrarrojo)
diferencia = mndwi_despues - mndwi_antes

Umbral Estadístico y Máscara

Aislamiento de zonas afectadas mediante un umbral calculado como promedio + desviación estándar de la diferencia. Manejo de nulos (NaN) reemplazándolos con ceros.

Clasificación y Cuantificación

Generación de máscara binaria (1 para zona de desastre, 0 para el resto). Estimación del área total multiplicando píxeles por resolución espacial (30m x 30m) y conversión a hectáreas.

Avanzamos de los algoritmos tradicionales hacia el Deep Learning (aprendizaje profundo). Una red neuronal se inspira en el cerebro biológico, usando capas de entrada, capas ocultas densas y una capa de salida que entrega un vector de probabilidades por cada píxel.

Preparación de Datos

from sklearn.preprocessing import LabelBinarizer
import tensorflow as tf
from tensorflow.keras import layers, models

binarizer = LabelBinarizer()
etiquetas_bin = binarizer.fit_transform(etiquetas)

Transformamos las categorías en vectores binarios (ej. [0, 1, 0]) para que la red pueda interpretarlos.

Arquitectura del Modelo

Usamos un modelo Sequential. Las capas ocultas usan activación ReLU para introducir no linealidades, y la capa de salida usa Softmax para generar la distribución de probabilidad entre las clases.

model = models.Sequential([
    layers.Input(shape=(n_bandas,)),
    layers.Dense(64, activation='relu'),
    layers.Dense(n_clases, activation='softmax')
])

Compilación y Entrenamiento

Optimizador Adam y función de pérdida categorical_crossentropy. Se definen epochs (recorridos totales) y batch_size (muestras por grupo) para evitar el sobreajuste (overfitting).

Predicción y Reconstrucción

Aplicamos predict() y usamos inverse_transform() para volver a las etiquetas originales. Finalmente, redimensionamos con reshape() al tamaño original de la imagen y exportamos a GeoTIFF con metadatos espaciales.

Implementación de Inteligencia Artificial para segmentación geoespacial con SAM Geo (Segment Anything Model de Meta). A diferencia de la clasificación de píxeles, la segmentación reconoce objetos completos como edificios o cuerpos de agua.

Entorno y Configuración

pip install segment-geospatial leafmap geopandas

Descarga de Datos

Usamos leafmap para generar un mapa interactivo, centrar la vista y descargar las teselas satelitales en formato GeoTIFF local.

Segmentación Automática

from samgeo import SamGeo

sam = SamGeo(model_type="vit_h", checkpoint="sam_vit_h_4b8939.pth")
sam.generate(input_raster, output_mask)

Vectorización y Exportación

Convertimos la máscara ráster resultante a polígonos georreferenciados en formatos estándar como GeoPackage (.gpkg) o Shapefile (.shp) para su uso en SIG.

sam.tiff_to_vector(output_mask, output_vector)

import geopandas as gpd
gdf = gpd.read_file(output_vector)
gdf.plot()