Python y Datos de Elevación
Procesamiento de terreno, LIDAR y visualización 3D con NumPy y GDAL
Los datos de elevación representan uno de los tipos más fascinantes de datos geoespaciales. Pueden mostrar propiedades tanto de datos vectoriales como raster, resultando en productos de datos únicos. Se utilizan para visualización de terreno, clasificación de cobertura del suelo, modelado hidrológico, enrutamiento de transporte, extracción de características y muchos otros propósitos.
En este capítulo aprenderás a procesar archivos ASCII Grid, crear mapas de relieve sombreado, generar contornos de elevación, trabajar con datos LIDAR y crear mallas 3D (TIN).
Archivos ASCII Grid
Los archivos ASCII Grid (ASCIIGRID) son un tipo de datos raster usualmente asociado con datos de elevación. Este formato de cuadrícula almacena datos como texto en filas y columnas de igual tamaño con un encabezado simple. Cada celda almacena un único valor numérico que puede representar alguna característica del terreno como elevación, pendiente o dirección de flujo.
Estructura del Encabezado
| Línea | Parámetro | Descripción |
|---|---|---|
| 1 | ncols | Número de columnas (eje X) |
| 2 | nrows | Número de filas (eje Y) |
| 3 | xllcorner | Coordenada X de la esquina inferior izquierda |
| 4 | yllcorner | Coordenada Y de la esquina inferior izquierda |
| 5 | cellsize | Tamaño de celda (resolución) |
| 6 | NODATA_value | Valor para celdas sin datos |
Lectura de Grids con NumPy
NumPy puede leer el formato ASCII Grid directamente usando el método loadtxt(), que está diseñado para leer arrays desde archivos de texto. Las primeras seis líneas consisten en el encabezado.
import numpy as np from linecache import getline # Archivo ASCII Grid de elevación source = "dem.asc" # Parsear el encabezado usando linecache hdr = [getline(source, i) for i in range(1, 7)] values = [float(h.split(" ")[-1].strip()) for h in hdr] # Desempaquetar valores del encabezado cols, rows, lx, ly, cell, nd = values # Cargar el DEM en un array de NumPy # skiprows=6 omite las líneas del encabezado arr = np.loadtxt(source, skiprows=6) print(f"Dimensiones: {arr.shape}") print(f"Elevación mín: {arr.min()} m") print(f"Elevación máx: {arr.max()} m")
Elevación mín: 356.0 m
Elevación máx: 2700.0 m
Escritura de Grids
Escribir grids es igual de fácil. Usamos numpy.savetxt() y construimos el encabezado como una cadena:
import numpy as np # Crear el encabezado usando format strings header = "ncols {}\n".format(myArray.shape[1]) header += "nrows {}\n".format(myArray.shape[0]) header += "xllcorner 277750.0\n" header += "yllcorner 6122250.0\n" header += "cellsize 1.0\n" header += "NODATA_value -9999" # Guardar el array con el encabezado with open("output.asc", "wb") as f: f.write(bytes(header, 'UTF-8')) np.savetxt(f, myArray, fmt="%1.2f")
linecache permite leer líneas específicas de un archivo de manera eficiente sin cargar todo el archivo en memoria.
Creación de Relieve Sombreado
Los mapas de relieve sombreado colorean la elevación de tal manera que parece que el terreno está iluminado por una luz en ángulo bajo, creando puntos brillantes y sombras. Este estilo estético crea una ilusión casi fotográfica que facilita comprender la variación del terreno.
Conceptos Clave
Implementación del Algoritmo
El script usa un método de ventana de 3x3 para escanear la imagen y procesar eficientemente cada celda considerando sus vecinos:
from linecache import getline import numpy as np # Archivos de entrada y salida source = "dem.asc" slopegrid = "slope.asc" aspectgrid = "aspect.asc" shadegrid = "relief.asc" # Parámetros de iluminación azimuth = 315.0 # Dirección del sol altitude = 45.0 # Ángulo del sol z = 1.0 # Exageración de elevación scale = 1.0 # Resolución NODATA = -9999 # Conversiones de grados a radianes deg2rad = 3.141592653589793 / 180.0 rad2deg = 180.0 / 3.141592653589793 # Parsear el encabezado 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 xres = cell yres = cell * -1 # Cargar el DEM arr = np.loadtxt(source, skiprows=6) # Crear ventanas 3x3 para procesamiento window = [] for row in range(3): for col in range(3): window.append(arr[row:(row + arr.shape[0] - 2), col:(col + arr.shape[1] - 2)])
Cálculo de Slope y Aspect
# Procesar ventanas 3x3 en direcciones X e Y x = ((z * window[0] + z * window[3] + z * window[3] + z * window[6]) - (z * window[2] + z * window[5] + z * window[5] + z * window[8])) / \ (8.0 * xres * scale) y = ((z * window[6] + z * window[7] + z * window[7] + z * window[8]) - (z * window[0] + z * window[1] + z * window[1] + z * window[2])) / \ (8.0 * yres * scale) # Calcular pendiente (slope) slope = 90.0 - np.arctan(np.sqrt(x * x + y * y)) * rad2deg # Calcular orientación (aspect) aspect = np.arctan2(x, y) # Calcular el relieve sombreado shaded = np.sin(altitude * deg2rad) * np.sin(slope * deg2rad) + \ np.cos(altitude * deg2rad) * np.cos(slope * deg2rad) * \ np.cos((azimuth - 90.0) * deg2rad - aspect) # Escalar valores de 0-1 a 0-255 shaded = shaded * 255
C = Celda central siendo procesada
Creación de Contornos de Elevación
Un contorno es una isolínea a lo largo de la misma elevación en un dataset. Los contornos generalmente se escalonan a intervalos para crear una forma intuitiva de representar datos de elevación, tanto visual como numéricamente, usando un dataset vectorial eficiente en recursos.
Parámetros de gdal.ContourGenerate()
Implementación con GDAL
from osgeo import gdal, ogr, osr # Archivo DEM de entrada source = "dem.asc" # Crear el shapefile de salida ogr_ds = ogr.GetDriverByName("ESRI Shapefile").CreateDataSource("contour.shp") # Definir la referencia espacial (EPSG:3157) srs = osr.SpatialReference() srs.ImportFromEPSG(3157) # Crear la capa de líneas ogr_lyr = ogr_ds.CreateLayer("contour", srs, ogr.wkbLineString) # Añadir campos para ID y elevación field_defn = ogr.FieldDefn("ID", ogr.OFTInteger) ogr_lyr.CreateField(field_defn) field_defn = ogr.FieldDefn("ELEV", ogr.OFTReal) ogr_lyr.CreateField(field_defn) # Abrir el DEM ds = gdal.Open(source) # Generar contornos cada 400 metros gdal.ContourGenerate( ds.GetRasterBand(1), # Banda fuente 400, # Intervalo de contorno 10, # Base del contorno [], # Niveles fijos (vacío) 0, # Usar NODATA: no 0, # Valor NODATA ogr_lyr, # Capa de destino 0, # Índice campo ID 1 # Índice campo ELEV ) ogr_ds = None
Visualización con PNGCanvas
import shapefile import pngcanvas # Abrir los contornos r = shapefile.Reader("contour.shp") # Configurar conversión mundo a píxeles xdist = r.bbox[2] - r.bbox[0] ydist = r.bbox[3] - r.bbox[1] iwidth = 800 iheight = 600 xratio = iwidth / xdist yratio = iheight / ydist contours = [] for shape in r.shapes(): for i in range(len(shape.parts)): pixels = [] if i < len(shape.parts) - 1: pt = shape.points[shape.parts[i]:shape.parts[i+1]] else: pt = shape.points[shape.parts[i]:] for x, y in pt: px = int(iwidth - ((r.bbox[2] - x) * xratio)) py = int((r.bbox[3] - y) * yratio) pixels.append([px, py]) contours.append(pixels) # Crear el canvas y dibujar canvas = pngcanvas.PNGCanvas(iwidth, iheight) canvas.color = [0xff, 0, 0, 0xff] # Rojo for c in contours: canvas.polyline(c) with open("contours.png", "wb") as f: f.write(canvas.dump())
Marching Squares, que es complejo de implementar con álgebra lineal de NumPy. GDAL proporciona esta funcionalidad a través de su API Python.
Trabajando con Datos LIDAR
LIDAR (Light Detection and Ranging) utiliza haces láser finitos que golpean el suelo cientos de miles de veces por segundo para recolectar una enorme cantidad de ubicaciones (x, y, z) muy precisas, así como tiempo e intensidad. Los valores de intensidad permiten clasificar y colorizar los datos LIDAR.
Formato LAS
El formato más común para datos LIDAR es LAS (LIDAR Exchange Format). Para leerlo usamos la biblioteca laspy:
pip install laspy
Creación de Grid desde LIDAR
Este script proyecta los puntos LIDAR a una cuadrícula regular con tamaño de celda de un metro, promediando puntos múltiples en la misma celda e interpolando huecos:
from laspy.file import File import numpy as np # Archivos de entrada y salida source = "lidar.las" target = "lidar.asc" cell = 1.0 # Tamaño de celda en metros NODATA = 0 # Abrir archivo LAS las = File(source, mode="r") # Obtener límites del dataset min_coords = las.header.min max_coords = las.header.max # Calcular dimensiones del grid xdist = max_coords[0] - min_coords[0] ydist = max_coords[1] - min_coords[1] cols = int(xdist / cell) + 1 rows = int(ydist / cell) + 1 # Arrays para conteo y suma de elevaciones count = np.zeros((rows, cols)).astype(np.float32) zsum = np.zeros((rows, cols)).astype(np.float32) # Proyectar puntos al grid ycell = -1 * cell projx = (las.x - min_coords[0]) / cell projy = (las.y - min_coords[1]) / ycell # Convertir a índices enteros ix = projx.astype(np.int32) iy = projy.astype(np.int32) # Agregar valores de elevación for x, y, z in np.nditer([ix, iy, las.z]): count[y, x] += 1 zsum[y, x] += z
Interpolación de Huecos
# Evitar división por cero nonzero = np.where(count > 0, count, 1) # Promediar valores Z zavg = zsum / nonzero # Interpolación simple: promedio de vecinos mean = np.ones((rows, cols)) * np.mean(zavg) left = np.roll(zavg, -1, 1) lavg = np.where(left > 0, left, mean) right = np.roll(zavg, 1, 1) ravg = np.where(right > 0, right, mean) interpolate = (lavg + ravg) / 2 fill = np.where(zavg > 0, zavg, interpolate) # Crear encabezado ASCII Grid header = "ncols {}\n".format(fill.shape[1]) header += "nrows {}\n".format(fill.shape[0]) header += "xllcorner {}\n".format(min_coords[0]) header += "yllcorner {}\n".format(min_coords[1]) header += "cellsize {}\n".format(cell) header += "NODATA_value {}\n".format(NODATA) # Guardar el grid with open(target, "wb") as f: f.write(bytes(header, 'UTF-8')) np.savetxt(f, fill, fmt="%1.2f")
Visualización de LIDAR con PIL
Podemos crear imágenes de salida en Python usando la Python Imaging Library (PIL). El módulo PIL.ImageOps tiene funciones para ecualización de histograma y mejora automática de contraste.
Imagen Básica Mejorada
import numpy as np from PIL import Image, ImageOps # Archivo DEM LIDAR en ASCII Grid source = "lidar.asc" target = "lidar.bmp" # Cargar el ASCII DEM arr = np.loadtxt(source, skiprows=6) # Convertir array a imagen PIL im = Image.fromarray(arr).convert("RGB") # Mejorar la imagen im = ImageOps.equalize(im) # Ecualizar histograma im = ImageOps.autocontrast(im) # Auto contraste # Guardar la imagen im.save(target)
Colorización con Rampa de Calor (HSV)
Usamos el módulo colorsys para conversión de espacio de color. HSV permite crear rampas de color usando un solo valor de matiz:
import numpy as np from PIL import Image, ImageOps import colorsys source = "lidar.asc" target = "lidar_color.bmp" # Cargar y preparar imagen arr = np.loadtxt(source, skiprows=6) im = Image.fromarray(arr).convert('L') im = ImageOps.equalize(im) im = ImageOps.autocontrast(im) # Construir rampa de color HSV # De azul (h=0.67) a rojo (h=0) palette = [] h = 0.67 # Matiz inicial (azul) s = 1 # Saturación v = 1 # Valor step = h / 256.0 for i in range(256): rp, gp, bp = colorsys.hsv_to_rgb(h, s, v) r = int(rp * 255) g = int(gp * 255) b = int(bp * 255) palette.extend([r, g, b]) h -= step # Aplicar paleta y guardar im.putpalette(palette) im.save(target)
Red Irregular Triangulada (TIN)
Una Red Irregular Triangulada (TIN) es una representación vectorial de un dataset de puntos como una superficie vectorial de puntos conectados como triángulos. Un algoritmo determina qué puntos son absolutamente necesarios para representar el terreno con precisión.
Beneficios del TIN
Implementación del TIN
import pickle import math import time import numpy as np import shapefile from laspy.file import File import voronoi # pip install http://git.io/vOEuJ # Archivos source = "clippedLAS.las" target = "mesh" # Clase Point requerida por voronoi class Point: def __init__(self, x, y): self.px = x self.py = y def x(self): return self.px def y(self): return self.py # Abrir archivo LAS las = File(source, mode="r") # Extraer puntos points = [] print("Ensamblando puntos...") for x, y in np.nditer((las.x, las.y)): points.append(Point(x, y)) # Triangulación de Delaunay print("Componiendo triángulos...") triangles = voronoi.computeDelaunayTriangulation(points)
Creación del Shapefile PolygonZ
print("Creando shapefile...") # PolygonZ permite vértices con valores Z w = shapefile.Writer(target, shapefile.POLYGONZ) w.field("X1", "C", "40") w.field("X2", "C", "40") w.field("X3", "C", "40") w.field("Y1", "C", "40") w.field("Y2", "C", "40") w.field("Y3", "C", "40") w.field("Z1", "C", "40") w.field("Z2", "C", "40") w.field("Z3", "C", "40") tris = len(triangles) max_segment = 3 # Filtrar triángulos muy grandes for i in range(tris): t = triangles[i] x1, y1, z1 = las.x[t[0]], las.y[t[0]], las.z[t[0]] x2, y2, z2 = las.x[t[1]], las.y[t[1]], las.z[t[1]] x3, y3, z3 = las.x[t[2]], las.y[t[2]], las.z[t[2]] # Filtrar triángulos grandes del convex hull if math.sqrt((x2-x1)**2+(y2-y1)**2) > max_segment: continue if math.sqrt((x3-x2)**2+(y3-y2)**2) > max_segment: continue if math.sqrt((x3-x1)**2+(y3-y1)**2) > max_segment: continue part = [[x1,y1,z1,0], [x2,y2,z2,0], [x3,y3,z3,0]] w.poly(parts=[part]) w.record(x1, x2, x3, y1, y2, y3, z1, z2, z3) w.close() print("Completado.")
voronoi.py de Bill Simons basada en el algoritmo de Steve Fortune. Filtramos triángulos muy grandes que son artefactos comunes del convex hull.
Resumen del Capítulo
01 // ASCII Grid
Lectura y escritura de rasters de elevación usando NumPy y linecache para procesamiento eficiente.
02 // Relieve Sombreado
Cálculo de slope, aspect y shaded relief usando ventanas 3x3 y operaciones matriciales.
03 // Contornos
Generación de isolíneas vectoriales con gdal.ContourGenerate() usando el algoritmo Marching Squares.
04 // LIDAR
Transformación de nubes de puntos LAS a ASCII Grid con interpolación de huecos.
05 // Visualización
Mejora de imágenes con PIL y colorización usando rampas HSV.
06 // TIN
Creación de superficies 3D eficientes con triangulación de Delaunay en shapefiles PolygonZ.
Casos de Uso
Las técnicas de este capítulo son herramientas fundamentales para análisis de terreno que se utilizan en:
- Planificación de transporte — Diseño de carreteras y análisis de rutas óptimas
- Planificación de construcción — Evaluación de sitios y movimiento de tierras
- Modelado hidrológico — Análisis de drenaje y cuencas
- Exploración geológica — Identificación de formaciones y recursos
- Gestión de emergencias — Modelado de inundaciones y deslizamientos
- Agricultura de precisión — Análisis de pendientes para irrigación