Python // Geospatial Analysis
Modelado Geoespacial Avanzado
Capítulo 8 — Algoritmos y productos de información del mundo real
El modelado geoespacial avanzado nos permite transformar datos crudos en productos de información que responden preguntas específicas del mundo real. En este capítulo construiremos sobre los conceptos de procesamiento de datos aprendidos para crear sistemas completos utilizados en agricultura, gestión de emergencias, logística y otras industrias.
Un modelo geoespacial es una representación simplificada de algún aspecto del mundo real que nos ayuda a responder una o más preguntas sobre un proyecto o problema. Aprenderemos a combinar múltiples conjuntos de datos procesados para construir estos modelos.
NumPy
GDAL
PIL
PyShp
NetworkX
Rasterio
OVERVIEW_00
Productos de Información
Los ejemplos en este capítulo son procesos reales que utilizarías como analista geoespacial profesional. Cada algoritmo, aunque diseñado para una tarea específica, tiene aplicaciones amplias en el análisis geoespacial.
Mapa de Salud de Cultivos
NDVI para agricultura de precisión
Modelo de Inundación
Predicción de impacto hidrológico
Hillshade Colorizado
Visualización avanzada de terreno
Ruta de Mínimo Costo
Navegación A* sobre terreno
Routing de Calles
Camino más corto en grafos
Geolocalización de Fotos
Extracción de coordenadas EXIF
MODULE_01
Índice de Vegetación NDVI
El Normalized Difference Vegetative Index (NDVI) muestra la salud relativa de las plantas en un área de interés. Utiliza imágenes satelitales o aéreas para destacar la densidad de clorofila, empleando únicamente las bandas roja e infrarroja cercana.
Fórmula NDVI
NDVI =
( Infrarrojo - Rojo )
/
( Infrarrojo + Rojo )
Flujo de Procesamiento NDVI
Imagen Multiespectral
→
Bandas R + IR
→
Fórmula NDVI
→
Clasificación
Configuración del Framework
ndvi_setup.py
from osgeo import gdal, gdal_array, ogr from PIL import Image, ImageDraw import numpy as np def imageToArray(i): """Convierte imagen PIL a array NumPy""" a = gdal_array.numpy.fromstring(i.tobytes(), 'b') a.shape = i.im.size[1], i.im.size[0] return a def world2Pixel(geoMatrix, x, y): """Convierte coordenadas geoespaciales a píxeles""" 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) # Cargar imagen multiespectral source = "farm.tif" target = "ndvi.tif" srcArray = gdal_array.LoadFile(source) srcImage = gdal.Open(source) geoTrans = srcImage.GetGeoTransform() # Extraer bandas roja e infrarroja r = srcArray[1] # Banda roja ir = srcArray[2] # Banda infrarroja
Rasterizar el Shapefile de Recorte
ndvi_rasterize.py
# Abrir shapefile del campo field = ogr.Open("field.shp") lyr = field.GetLayer("field") poly = lyr.GetNextFeature() # Convertir extensión a coordenadas de píxel minX, maxX, minY, maxY = lyr.GetExtent() 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) # Recortar bandas al bounding box rClip = r[ulY:lrY, ulX:lrX] irClip = ir[ulY:lrY, ulX:lrX] # Crear máscara como imagen PIL rasterPoly = Image.new("L", (pxWidth, pxHeight), 1) rasterize = ImageDraw.Draw(rasterPoly) rasterize.polygon(pixels, 0) # Convertir máscara a array NumPy mask = imageToArray(rasterPoly)
Aplicar Fórmula NDVI
ndvi_calculate.py
# Recortar bandas usando la máscara rClip = gdal_array.numpy.choose(mask, (rClip, 0)).astype(gdal_array.numpy.uint8) irClip = gdal_array.numpy.choose(mask, (irClip, 0)).astype(gdal_array.numpy.uint8) # Ignorar warnings de NaN gdal_array.numpy.seterr(all="ignore") # Fórmula NDVI: (infrarrojo - rojo) / (infrarrojo + rojo) ndvi = 1.0 * ((irClip - rClip) / (irClip + rClip + 1.0)) # Convertir NaN a 0 ndvi = gdal_array.numpy.nan_to_num(ndvi) # Guardar como GeoTIFF gtiff = gdal.GetDriverByName('GTiff') gtiff.CreateCopy(target, copy_geo(ndvi, prototype=source))
Clasificar y Colorizar NDVI
ndvi_classify.py
import gdal_array as gd import operator from functools import reduce # Cargar NDVI y crear imagen RGB de salida ndvi = gd.LoadFile("ndvi.tif").astype(gd.numpy.uint8) rgb = gd.numpy.zeros((3, len(ndvi), len(ndvi[0])), gd.numpy.uint8) # Definir 7 clases con rangos de valores classes = [58, 73, 110, 147, 184, 220, 255] # Look-up table: marrón oscuro → verde oscuro lut = [ [120, 69, 25], # Suelo desnudo [255, 178, 74], # Vegetación baja [255, 237, 166], # Vegetación media-baja [173, 232, 94], # Vegetación media [135, 181, 64], # Vegetación media-alta [3, 156, 0], # Vegetación alta [1, 100, 0] # Vegetación muy saludable ] # Clasificar cada rango start = 1 for i in range(len(classes)): mask = gd.numpy.logical_and(start <= ndvi, ndvi <= classes[i]) for j in range(len(lut[i])): rgb[j] = gd.numpy.choose(mask, (rgb[j], lut[i][j])) start = classes[i] + 1 # Guardar imagen clasificada gd.SaveArray(rgb, "ndvi_color.tif", format="GTiff")
Los agricultores usan estos datos para determinar cómo irrigar y aplicar fertilizantes de manera
eficiente y localizada. Las clases pueden convertirse en shapefiles vectoriales para sistemas GPS de aplicación variable.
MODULE_02
Modelo de Inundación
Un modelo de inundación nos ayuda a entender el impacto de una crecida. Utiliza el algoritmo Flood Fill (conocido del juego Buscaminas y herramientas de relleno en editores gráficos) para simular la extensión máxima de agua que puede contener una cuenca a cierta elevación.
Algoritmo Flood Fill
Punto Inicial
→
Filtrar Elevación
→
Cola de Celdas
→
Máscara de Inundación
Función Flood Fill
flood_fill.py
import numpy as np from linecache import getline def floodFill(x, y, arr): """ Flood fill basado en cola - no recursivo. Procesa celdas únicas usando sets de Python. """ # Crear array de salida lleno de ceros fld = np.zeros_like(arr) # Sets para tracking: fill (por procesar), filled (procesadas) fill = set() filled = set() # Agregar punto inicial fill.add((x, y)) # Dimensiones del array height, width = arr.shape while fill: # Obtener siguiente celda x, y = fill.pop() # Saltar si ya procesada o fuera de límites if (x, y) in filled: continue if x < 0 or x >= width: continue if y < 0 or y >= height: continue # Verificar si la celda está bajo el umbral if arr[y][x] == 1: # Marcar como inundada fld[y][x] = 1 filled.add((x, y)) # Agregar vecinos a la cola (4 direcciones) fill.add((x-1, y)) fill.add((x+1, y)) fill.add((x, y-1)) fill.add((x, y+1)) else: filled.add((x, y)) return fld
Predicción de Inundación
predict_flood.py
# Archivos de entrada/salida source = "dem.asc" target = "flood.asc" # Cargar DEM print("Loading DEM...") img = np.loadtxt(source, skiprows=6) # Filtrar: elevaciones < 70 metros = inundables wet = np.where(img < 70, 1, 0) print("Image masked") # Parsear header del ASCII Grid hdr = [getline(source, i) for i in range(1, 7)] values = [float(h.split(" ")[-1].strip()) for h in hdr] cols, rows, lx, ly, cell, nd = values # Punto inicial en el lecho del río (coordenadas de píxel) sx = 2582 sy = 2057 # Ejecutar flood fill print("Beginning flood fill...") fld = floodFill(sx, sy, wet) print("Finished flood fill") # Guardar resultado como ASCII Grid header = "" for i in range(6): header += hdr[i] with open(target, "wb") as f: f.write(bytes(header, 'UTF-8')) np.savetxt(f, fld, fmt="%1i") print("Done!")
Cambiando el valor de elevación (70, 90, etc.) puedes crear múltiples capas de inundación mostrando diferentes escenarios. El modelo puede extenderse con
GDAL.polygonize() para crear polígonos vectoriales.
MODULE_03
Hillshade Colorizado
Combinamos la técnica de shaded relief con clasificación por colores para crear visualizaciones de terreno que sirven como fondos para mapas GIS. Este proceso mezcla un DEM colorizado con un hillshade utilizando PIL.
Pipeline de Visualización
DEM
Relief
Color Classification
PIL Blend (40%)
Hillshade Final
Crear Hillshade Colorizado
color_hillshade.py
import gdal_array as gd from PIL import Image # Archivos de entrada relief = "relief.asc" # Shaded relief pre-calculado dem = "dem.asc" # Digital Elevation Model target = "hillshade.tif" # Cargar relief como fondo bg = gd.numpy.loadtxt(relief, skiprows=6) # Cargar DEM como foreground fg = gd.numpy.loadtxt(dem, skiprows=6)[:-2, :-2] # Crear imagen RGB vacía rgb = gd.numpy.zeros((3, len(fg), len(fg[0])), gd.numpy.uint8) # Clases por rangos de elevación (metros) classes = [356, 649, 942, 1235, 1528, 1821, 2114, 2300, 2700] # Paleta de colores: verde → marrón → gris (nieve) lut = [ [63, 159, 152], # Verde agua (bajo) [96, 235, 155], # Verde claro [100, 246, 174], # Verde brillante [248, 251, 155], # Amarillo [246, 190, 39], # Naranja [242, 155, 39], # Naranja oscuro [165, 84, 26], # Marrón [236, 119, 83], # Terracota [203, 203, 203] # Gris (cumbres) ] # Clasificar start = 1 for i in range(len(classes)): mask = gd.numpy.logical_and(start <= fg, fg <= classes[i]) for j in range(len(lut[i])): rgb[j] = gd.numpy.choose(mask, (rgb[j], lut[i][j])) start = classes[i] + 1 # Convertir arrays a imágenes PIL im1 = Image.fromarray(bg).convert('RGB') im2 = Image.fromarray(rgb.transpose(1, 2, 0)).convert('RGB') # Mezclar con 40% alpha hillshade = Image.blend(im1, im2, .4) # Guardar resultado hillshade.save(target)
MODULE_04
Análisis de Ruta de Mínimo Costo
El algoritmo A* (A-star) busca la ruta óptima considerando tanto la distancia como el costo del terreno. A diferencia del algoritmo Dijkstra tradicional, A* es ideal para navegar sobre grids como un DEM, balanceando múltiples factores.
TERRENO
Peso principal: evita elevaciones altas
DISTANCIA
Progreso hacia el destino
EFICIENCIA
8 direcciones de movimiento
Función de Puntuación
scoring.py
import math def e_dist(p1, p2): """Distancia Euclidiana entre dos puntos""" y1, x1 = p1 y2, x2 = p2 return math.sqrt((x2 - x1)**2 + (y2 - y1)**2) def score(cur, node, h, start, end): """ Calcula puntuación para elegir el siguiente nodo. Mayor peso al terreno, luego progreso hacia el objetivo. """ score = 0 # Elevación y distancias del nodo actual cur_h = h[cur] cur_g = e_dist(cur, end) # Distancia al destino cur_d = e_dist(cur, start) # Distancia desde inicio # Elevación y distancias del nodo vecino node_h = h[node] node_g = e_dist(node, end) node_d = e_dist(node, start) # Comparar valores con pesos if node_h < cur_h: score += cur_h - node_h # Bonus por bajar elevación if node_g < cur_g: score += 10 # Bonus por acercarse al destino if node_d > cur_d: score += 10 # Bonus por alejarse del inicio return score
Algoritmo A*
astar.py
def astar(start, end, h): """ Búsqueda A* sobre grid de terreno. Retorna set de celdas que forman la ruta. """ closed_set = set() # Nodos ya evaluados open_set = set() # Nodos por evaluar path = set() # Ruta de salida open_set.add(start) while open_set: cur = open_set.pop() if cur == end: return path closed_set.add(cur) path.add(cur) # Obtener vecinos (8 direcciones) y, x = cur neighbors = [] for dy in [-1, 0, 1]: for dx in [-1, 0, 1]: if dy == 0 and dx == 0: continue ny, nx = y + dy, x + dx if 0 <= ny < h.shape[0] and 0 <= nx < h.shape[1]: neighbors.append((ny, nx)) # Encontrar mejor vecino best = None best_score = -1 for node in neighbors: if node in closed_set: continue s = score(cur, node, h, start, end) if s > best_score: best_score = s best = node if best: open_set.add(best) return path
Convertir Ruta a Shapefile
route_to_shapefile.py
import pickle import shapefile def pix2coord(gt, x, y): """Convierte píxel a coordenadas geoespaciales""" c = gt[0] + (x * gt[1]) r = gt[3] + (y * gt[5]) return (c, r) # Cargar ruta guardada con pickle with open("path.p", "rb") as f: path = pickle.load(f) # Convertir a coordenadas geoespaciales coords = [] for y, x in path: coords.append(pix2coord(gt, x, y)) # Crear shapefile de línea with shapefile.Writer("path", shapeType=shapefile.POLYLINE) as w: w.field("NAME") w.record("LeastCostPath") w.line([coords])
MODULE_05
Routing en Redes de Calles
El routing en calles utiliza un grafo conectado de líneas donde cada línea puede tener valores de impedancia (tráfico, velocidad, distancia). Utilizamos NetworkX para calcular la ruta más corta entre dos puntos.
Flujo de Routing
Shapefile Red
→
NetworkX Graph
→
Shortest Path
→
Shapefile Ruta
Implementación con NetworkX
street_routing.py
import networkx as nx import shapefile import math def haversine(p1, p2): """Distancia en metros entre dos puntos lat/lon""" lat1, lon1 = map(math.radians, p1) lat2, lon2 = map(math.radians, p2) dlat = lat2 - lat1 dlon = lon2 - lon1 a = math.sin(dlat/2)**2 + \ math.cos(lat1) * math.cos(lat2) * math.sin(dlon/2)**2 c = 2 * math.asin(math.sqrt(a)) return 6371000 * c # Radio de la Tierra en metros # Cargar red de calles r = shapefile.Reader("roads.shp") shapes = r.shapes() # Crear grafo sg = nx.Graph() # Agregar edges del shapefile for shape in shapes: points = shape.points for i in range(len(points) - 1): p1 = tuple(points[i]) p2 = tuple(points[i + 1]) sg.add_edge(p1, p2) # Calcular distancias para cada edge for n0, n1 in sg.edges(): dist = haversine(n0, n1) sg[n0][n1]["dist"] = dist # Puntos de inicio y fin start = (-89.32, 30.31) end = (-89.28, 30.35) # Encontrar nodos más cercanos nn_start = min(sg.nodes(), key=lambda n: haversine(start, n)) nn_end = min(sg.nodes(), key=lambda n: haversine(end, n)) # Calcular ruta más corta path = nx.shortest_path(sg, source=nn_start, target=nn_end, weight="dist") # Guardar como shapefile with shapefile.Writer("route", shapefile.POLYLINE) as w: w.field("NAME", "C", 40) w.line(parts=[[list(p) for p in path]]) w.record("route")
MODULE_06
Geolocalización de Fotos
Las fotos tomadas con cámaras GPS (incluyendo smartphones) almacenan coordenadas en el header del archivo usando EXIF tags. Podemos extraer esta información para crear shapefiles con las ubicaciones de las fotos.
geolocate_photos.py
import glob import os from PIL import Image from PIL.ExifTags import TAGS import shapefile def exif(img): """Extrae datos EXIF de una imagen""" exif_data = {} try: i = Image.open(img) tags = i._getexif() for tag, value in tags.items(): decoded = TAGS.get(tag, tag) exif_data[decoded] = value except: pass return exif_data def dms2dd(d, m, s, direction): """Convierte Grados/Minutos/Segundos a Decimal""" sec = float((m * 60) + s) dec = float(sec / 3600) deg = float(d + dec) if direction.upper() in ['W', 'S']: deg = deg * -1 return deg def gps(exif_data): """Extrae coordenadas GPS de datos EXIF""" lat, lon = None, None if 'GPSInfo' in exif_data: coords = exif_data['GPSInfo'] # Latitud i = coords[1] d, m, s = coords[2][0][0], coords[2][1][0], coords[2][2][0] lat = dms2dd(d, m, s, i) # Longitud i = coords[3] d, m, s = coords[4][0][0], coords[4][1][0], coords[4][2][0] lon = dms2dd(d, m, s, i) return lat, lon # Procesar directorio de fotos photos = {} photo_dir = "./photos" files = glob.glob(os.path.join(photo_dir, "*.jpg")) for f in files: e = exif(f) lat, lon = gps(e) if lat and lon: photos[f] = [lon, lat] # Crear shapefile de puntos with shapefile.Writer("photos", shapefile.POINT) as w: w.field("NAME", "C", 80) for f, coords in photos.items(): w.point(*coords) w.record(f)
En QGIS puedes configurar acciones para que al hacer clic en un punto se abra la foto correspondiente. Los atributos del shapefile contienen la ruta al archivo de imagen.
MODULE_07
Máscara de Nubes Satelital
Las imágenes satelitales frecuentemente contienen nubes que interfieren con el análisis. Usando datos de quality-assurance de Landsat 8, podemos generar máscaras para filtrar automáticamente las áreas nubladas.
cloud_mask.py
import glob import os import rasterio from l8qa.qa import write_cloud_mask # Directorio con datos Landsat landsat_dir = "l8" # Localizar archivo de quality-assurance src_qa = glob.glob(os.path.join(landsat_dir, '*QA*'))[0] # Generar máscara de nubes with rasterio.open(src_qa) as qa_raster: profile = qa_raster.profile profile.update(nodata=0) write_cloud_mask( qa_raster.read(1), profile, 'cloudmask.tif' )
pip install
pip install rasterio
pip install rio-l8qa
pip install rio-l8qa
Resumen de Algoritmos
01 // NDVI
Índice de vegetación usando bandas roja e infrarroja para evaluar salud de cultivos.
02 // Flood Fill
Algoritmo basado en cola para modelar inundaciones desde un punto de origen.
03 // Color Hillshade
Combinación de shaded relief con clasificación por elevación usando PIL blend.
04 // A* Search
Búsqueda heurística para encontrar rutas de mínimo costo sobre terreno.
05 // Graph Routing
NetworkX para calcular caminos más cortos en redes de calles.
06 // EXIF Extraction
Geolocalización de fotos extrayendo coordenadas GPS de metadatos.