Python y Teledetección - Análisis de Imágenes Satelitales

Chapter 06 // Geospatial

Python y Teledetección

Procesamiento de imágenes satelitales con GDAL, NumPy y PIL

La teledetección (Remote Sensing) consiste en recopilar información sobre la Tierra sin contacto físico directo. Típicamente involucra imágenes satelitales, imágenes aéreas, LIDAR (Light Detection and Ranging) y radar de apertura sintética. En teledetección, recorremos cada píxel de una imagen y realizamos consultas o procesos matemáticos.

Una imagen puede verse como un gran array numérico. En teledetección, estos arrays pueden ser enormes, desde decenas de megabytes hasta varios gigabytes. Usaremos PIL para procesamiento de imágenes y NumPy para matemáticas con arrays multidimensionales.

MÓDULO_01

Intercambio de Bandas

Nuestros ojos solo pueden ver colores en el espectro visible como combinaciones de rojo, verde y azul (RGB). Los sensores aerotransportados y espaciales pueden capturar longitudes de onda de energía fuera del espectro visible. Para visualizar estos datos, movemos imágenes que representan diferentes longitudes de onda dentro y fuera de los canales RGB.

Espectro Electromagnético en Sensores

AZUL

VERDE

ROJO

NIR

SWIR

TIR

GDAL incluye el módulo gdal_array que carga y guarda imágenes de teledetección hacia y desde arrays NumPy para fácil manipulación. Las imágenes en NumPy son arrays multidimensionales en el orden: banda, altura, ancho.

Cargar y Reorganizar Bandas

swap_bands.py

from gdal import gdal_array

# Nombre de la imagen fuente
src = "FalseColor.tif"

# Cargar la imagen en un array numpy
arr = gdal_array.LoadFile(src)

# Intercambiar bandas 1 y 2 para imagen en color natural
# Usamos "advanced slicing" de numpy para reordenar
output = gdal_array.SaveArray(
    arr[[1, 0, 2], :], 
    "swap.tif",
    format="GTiff", 
    prototype=src
)

# Liberar referencia para evitar archivo corrupto
output = None

El argumento prototype permite especificar otra imagen de la cual GDAL copia información de referencia espacial. Sin esto, tendríamos una imagen sin georreferenciación.

Transformación de Bandas

Banda 0 (NIR)

Banda 1 (R)

Banda 2 (G)

↓ reorder ↓

R → Canal R

NIR → Canal G

G → Canal B

MÓDULO_02

Creación de Histogramas

Un histograma es la distribución estadística de valores en un dataset. El eje horizontal representa valores únicos mientras el eje vertical representa la frecuencia. Es una operación clave en el procesamiento raster, útil para mejorar contraste, clasificación de objetos y comparación de imágenes.

Beneficios del Análisis

DIAGNÓSTICO

Identificar problemas de exposición y distribución de valores

CLASIFICACIÓN

Base para agrupar píxeles por valores similares

COMPARACIÓN

Evaluar diferencias entre múltiples imágenes

Función de Histograma

histogram.py

from gdal import gdal_array

def histogram(a, bins=list(range(0, 256))):
    """
    Función de histograma para array multidimensional.
    a = array
    bins = rango de números a emparejar
    """
    # Aplanar el array para el histograma
    fa = a.flat
    
    # Ordenar y buscar posiciones
    n = gdal_array.numpy.searchsorted(
        gdal_array.numpy.sort(fa), bins
    )
    n = gdal_array.numpy.concatenate([n, [len(fa)]])
    
    # Calcular frecuencias
    hist = n[1:] - n[:-1]
    return hist

# Cargar imagen y calcular histograma por banda
src = "swap.tif"
arr = gdal_array.LoadFile(src)

histograms = []
for band in arr:
    histograms.append(histogram(band))

Visualización de Histograma (Ejemplo)

MÓDULO_03

Estiramiento de Histograma

Una operación de estiramiento de histograma redistribuye los valores de píxeles a través de toda la escala. Al hacerlo, tenemos más valores en niveles de mayor intensidad y la imagen se vuelve más brillante y con mejor contraste.

Algoritmo de Estiramiento

stretch.py

import gdal_array
import operator
from functools import reduce

def histogram(a, bins=list(range(0, 256))):
    fa = a.flat
    n = gdal_array.numpy.searchsorted(
        gdal_array.numpy.sort(fa), bins
    )
    n = gdal_array.numpy.concatenate([n, [len(fa)]])
    hist = n[1:] - n[:-1]
    return hist

def stretch(a):
    """
    Realiza estiramiento de histograma en una imagen.
    """
    hist = histogram(a)
    lut = []
    
    for b in range(0, len(hist), 256):
        # Tamaño de paso - crear bins de intervalo igual
        step = reduce(operator.add, hist[b:b+256]) / 255
        
        # Crear tabla de ecualización
        n = 0
        for i in range(256):
            lut.append(n / step)
            n = n + hist[i + b]
    
    gdal_array.numpy.take(lut, a, out=a)
    return a

# Aplicar estiramiento
src = "swap.tif"
arr = gdal_array.LoadFile(src)
stretched = stretch(arr)

output = gdal_array.SaveArray(
    arr, "stretched.tif",
    format="GTiff", prototype=src
)
output = None

Después del estiramiento, todas las bandas se distribuyen uniformemente. Su distribución relativa entre sí es la misma, pero dentro de la imagen ahora están dispersas a través de todo el espectro.

MÓDULO_04

Recorte de Imágenes

Raramente un analista está interesado en una escena satelital completa, que puede cubrir cientos de millas cuadradas. Dado el tamaño de los datos satelitales, estamos motivados a reducir el tamaño de la imagen solo a nuestra área de interés usando shapefiles como límite.

Pipeline de Recorte

Cargar imagen en array usando gdal_array

↓

Crear lector de shapefile con PyShp

↓

Rasterizar shapefile a imagen georreferenciada

↓

Convertir imagen shapefile en máscara binaria

↓

Filtrar imagen satelital a través de la máscara

↓

Guardar imagen recortada como clip.tif

Funciones de Conversión

clip_utils.py

from osgeo import gdal, gdal_array, osr
import shapefile
from PIL import Image, ImageDraw

def world2Pixel(geoMatrix, x, y):
    """
    Convierte coordenadas geográficas a coordenadas de píxel.
    Usa el objeto GeoTransform de GDAL.
    """
    ulX = geoMatrix[0]
    ulY = geoMatrix[3]
    xDist = geoMatrix[1]
    yDist = geoMatrix[5]
    
    pixel = int((x - ulX) / xDist)
    line = int((ulY - y) / abs(yDist))
    
    return (pixel, line)

def imageToArray(img):
    """
    Convierte una imagen PIL en un array NumPy.
    """
    a = gdal_array.numpy.fromstring(
        img.tobytes(), 'b'
    )
    a.shape = img.im.size[1], img.im.size[0]
    return a

Recorte con Máscara

clip_image.py

# Cargar imagen fuente y shapefile
raster = "stretched.tif"
shp = "hancock"
output = "clip"

srcArray = gdal_array.LoadFile(raster)
srcImage = gdal.Open(raster)
geoTrans = srcImage.GetGeoTransform()

r = shapefile.Reader("{}.shp".format(shp))

# Convertir extensión del shapefile a coordenadas de imagen
minX, minY, maxX, maxY = r.bbox
ulX, ulY = world2Pixel(geoTrans, minX, maxY)
lrX, lrY = world2Pixel(geoTrans, maxX, minY)

# Calcular tamaño de la nueva imagen
pxWidth = int(lrX - ulX)
pxHeight = int(lrY - ulY)
clip = srcArray[:, ulY:lrY, ulX:lrX]

# Crear máscara del polígono
pixels = []
for p in r.shape(0).points:
    pixels.append(world2Pixel(geoTrans, p[0], p[1]))

rasterPoly = Image.new("L", (pxWidth, pxHeight), 1)
rasterize = ImageDraw.Draw(rasterPoly)
rasterize.polygon(pixels, 0)

# Convertir máscara PIL a NumPy y aplicar
mask = imageToArray(rasterPoly)
clip = gdal_array.numpy.choose(
    mask, (clip, 0)
).astype(gdal_array.numpy.uint8)

# Guardar resultado
gdal_array.SaveArray(
    clip, "{}.tif".format(output),
    format="GTiff", prototype=raster
)

Las áreas fuera del límite del condado que aparecen en negro son valores NoData, lo que significa que no hay información en esa ubicación y son ignorados por la mayoría del software geoespacial.

MÓDULO_05

Clasificación de Imágenes

La teledetección automatizada (ARS) raramente se hace en el espectro visible. Cada píxel tiene un valor de reflectancia que es información, no solo color. Python puede ver esos valores y distinguir características agrupando valores de píxeles relacionados. Esta técnica se llama clasificación.

Tipos de Clasificación

NO SUPERVISADA

Sin entrada adicional, solo la imagen. Agrupación por algoritmo de distribución

SUPERVISADA

Usa datos de entrenamiento para guiar al computador

MACHINE LEARNING

Métodos avanzados con IA y aprendizaje automático

Clasificación No Supervisada

classify.py

import gdal
from gdal import gdal_array, osr

# Archivos de entrada y salida
src = "thermal.tif"
tgt = "classified.jpg"

# Cargar imagen en array numpy
srcArr = gdal_array.LoadFile(src)

# Obtener histograma y rangos de clases
classes = 20
hist, bins = gdal_array.numpy.histogram(srcArr, bins=classes)

# Tabla de colores arbitraria para 20 clases
lut = [
    [31, 120, 180], [178, 223, 138],
    [51, 160, 44], [251, 154, 153],
    [227, 26, 28], [253, 191, 111],
    [255, 127, 0], [202, 178, 214],
    [106, 61, 154], [255, 255, 153],
    # ... más colores ...
]

# Crear imagen RGB de salida
rgb = gdal_array.numpy.zeros(
    (3, srcArr.shape[0], srcArr.shape[1]),
    gdal_array.numpy.uint8
)

# Asignar colores a cada clase
for i in range(len(bins) - 1):
    mask = gdal_array.numpy.logical_and(
        bins[i] <= srcArr, srcArr < bins[i+1]
    )
    for j in range(3):
        rgb[j] = gdal_array.numpy.choose(
            mask, (rgb[j], lut[i % len(lut)][j])
        )

# Guardar imagen clasificada
output = gdal_array.SaveArray(
    rgb, tgt, format="JPEG"
)

Paleta de Clasificación (Ejemplo 7 clases)

C1
Suelo

C7
Sano

MÓDULO_06

Extracción de Características

La extracción de características es la digitalización automática de elementos en una imagen hacia puntos, líneas o polígonos. GDAL incluye una función llamada Polygonize() que automáticamente crea un shapefile de polígonos vectoriales a partir de celdas raster contiguas con el mismo valor.

Vectorización Automática

extract_features.py

from osgeo import gdal, ogr, osr

# Archivos de entrada y salida
src = "classified.tif"
tgt = "extract.shp"
tgtLayer = "extract"

# Abrir raster de entrada
srcDS = gdal.Open(src)
band = srcDS.GetRasterBand(1)
mask = band  # Usar banda como máscara

# Configurar shapefile de salida
driver = ogr.GetDriverByName("ESRI Shapefile")
shp = driver.CreateDataSource(tgt)

# Copiar referencia espacial
srs = osr.SpatialReference()
srs.ImportFromWkt(srcDS.GetProjectionRef())
layer = shp.CreateLayer(tgtLayer, srs=srs)

# Configurar atributos
fd = ogr.FieldDefn("DN", ogr.OFTInteger)
layer.CreateField(fd)
dst_field = 0

# ¡Extraer polígonos automáticamente!
extract = gdal.Polygonize(
    band, mask, layer, dst_field, [], None
)

Proceso de Vectorización

Raster Clasificado

→

gdal.Polygonize()

→

Shapefile Vector

La extracción de características funciona bien con cuerpos de agua, islas, carreteras, campos agrícolas, edificios y otras características que tienden a tener valores de píxeles de alto contraste con su fondo.

MÓDULO_07

Detección de Cambios

La detección de cambios es el proceso de tomar dos imágenes georreferenciadas de la misma área exacta en dos fechas diferentes e identificar automáticamente las diferencias. Es otra forma de clasificación de imágenes. La técnica básica es simplemente restar la primera imagen de la segunda.

Algoritmo de Diferencia Simple

change_detection.py

from osgeo import gdal
from gdal import gdal_array
import numpy as np

# Imágenes antes y después
img1 = "before.tif"
img2 = "after.tif"

# Cargar imágenes en arrays
arr1 = gdal_array.LoadFile(img1).astype(np.int16)
arr2 = gdal_array.LoadFile(img2).astype(np.int16)

# Calcular diferencia
diff = arr2 - arr1

# Crear imagen de cambio con colores
# Verde = algo añadido (valor positivo)
# Rojo = algo removido (valor negativo)
classes = 5
lut = [
    [255, 0, 0],    # Rojo: muy negativo
    [255, 128, 128], # Rosa: negativo
    [220, 220, 220], # Gris: sin cambio
    [128, 255, 128], # Verde claro: positivo
    [0, 255, 0]      # Verde: muy positivo
]

# Histograma para definir clases
hist, bins = np.histogram(diff, bins=classes)

# Crear imagen RGB de salida
rgb = np.zeros((3, diff.shape[0], diff.shape[1]), np.uint8)

# Asignar colores basados en rangos
for i in range(len(bins) - 1):
    mask = np.logical_and(
        bins[i] <= diff, diff < bins[i+1]
    )
    for j in range(3):
        rgb[j] = np.choose(mask, (rgb[j], lut[i][j]))

# Guardar resultado
im2 = gdal.Open(img2)
output = gdal_array.SaveArray(
    rgb, "change.tif",
    format="GTiff", prototype=im2
)

Fórmula de Detección de Cambios

Δ = Imagen_después - Imagen_antes

Interpretación de Resultados

Color	Significado	Ejemplo
Verde	Algo fue añadido (valores más brillantes)	Nueva construcción, vehículos
Rojo	Algo fue removido (valores más oscuros)	Demolición, erosión
Gris	Sin cambio significativo	Terreno estable

El concreto es un indicador importante en detección de cambios porque es muy brillante bajo la luz solar y usualmente indica nuevo desarrollo. Inversamente, si un edificio es demolido, la diferencia también es fácil de identificar.

Resumen del Capítulo

01 // Intercambio de Bandas

Reorganizar canales RGB para visualizar datos multiespectrales de forma más intuitiva.

02 // Histogramas

Analizar la distribución estadística de valores para diagnóstico y clasificación.

03 // Estiramiento

Redistribuir valores de píxeles para mejorar contraste y visualización.

04 // Recorte

Reducir imágenes al área de interés usando shapefiles como máscaras.

05 // Clasificación

Agrupar píxeles por valores similares para identificar características.

06 // Extracción

Convertir datos raster clasificados en geometrías vectoriales.

Librerías Principales

GDAL / gdal_array

Biblioteca de abstracción de datos geoespaciales. Carga y guarda imágenes en arrays NumPy.

NumPy

Matemáticas con arrays multidimensionales. Procesamiento de datos raster a alta velocidad.

PIL / Pillow

Python Imaging Library. Manipulación de imágenes y creación de máscaras.

PyShp (shapefile)

Lectura y escritura de shapefiles ESRI para datos vectoriales.

Alejandro Gutiérrez Marín

Archivo

Categorias

Visitas al Blog

Acerca de Mi

jueves, 22 de enero de 2026

Análisis geoespacial con Python - Parte 6 - Teledetección