Hallar la altura de la vegetación con archivos LAZ del IGN

 

Cómo obtener la altura de la vegetación (CHM) a partir de LiDAR en QGIS con PDAL

Para calcular la altura real del arbolado o la masa forestal (Canopy Height Model o CHM) no necesitamos complementos externos ni plugins adicionales a partir de QGIS 3.18, ya que el motor de PDAL viene integrado directamente en la Caja de herramientas de procesamiento.

El procedimiento se realiza en tres fases: obtener el Modelo Digital del Terreno (MDT), obtener el Modelo Digital de Vegetación (MDVeg) y restar ambos modelos con la Calculadora ráster.

Todo a partir de archivos LAZ desde el Centro de Descargas del CNIG https://centrodedescargas.cnig.es/CentroDescargas/index.jsp 

Paso 1: Generar el Modelo Digital del Terreno (MDT)

El terreno es una superficie continua que requiere interpolación para evitar vacíos o huecos de información en las zonas sin retornos directo de suelo.

1. Abre la Caja de herramientas de procesamiento (Ctrl + Alt + T).
2. Selecciona la herramienta Filtro (en el bloque de Extracción de nube de puntos).
3. Introduce la expresión para aislar únicamente el suelo: Classification = 2
4. A continuación, ejecuta la herramienta Exportar a ráster (usando triangulación) -en conversión de nube de puntos- utilizando la nube filtrada como capa de entrada:
    • Atributo: Z
    • Resolución: Definir en metros (ej. 1.0).
    • Tile size: 1000 (valor adecuado para teselas del IGN PNOA).
    • Salida: Guardar como GeoTIFF (.tif).

Paso 2: Generar el Modelo Digital de Vegetación (MDVeg)

La vegetación está formada por elementos discontinuos y aislados (copas de árboles). Para evitar los artefactos y rampas falsas que crea la triangulación entre árboles lejanos, utilizaremos una exportación ráster basada en el valor máximo.

1. En la Caja de herramientas, abre la herramienta Filtro
2. Introduce la expresión para filtrar los retornos de vegetación baja, media y alta: Classification = 3 OR Classification = 4 OR Classification = 5
3. Ejecuta la herramienta estándar Exportar a ráster -sin triangulación- sobre la capa filtrada:
    • Atributo: Z
    • Resolución: La misma resolución utilizada para el MDT (ej. 1.0).
    • Método de agregación: Máximo (max).
    • Tile size: 1000.
    • Salida: Guardar como GeoTIFF (.tif).

Paso 3: Calcular la altura de la vegetación con la Calculadora Ráster

Una vez exportados ambos archivos GeoTIFF con alturas absolutas respecto al nivel del mar, abrimos la Calculadora ráster (Ráster > Calculadora ráster...) para calcular la altura relativa del arbolado desde el suelo:

CHM (Altura de Vegetación) = MDVeg - MDT

El ráster resultante mostrará la altura real en metros de cada árbol o masa vegetal de tu zona de estudio.

 

Aquí puedes ver el cálculo de altura de los árboles de El Retiro, Madrid. https://edu.forestry.es/2026/04/arboles-singulares-de-el-retiro.html

 

Automatizar la generación del Canopy Height Model (CHM) en QGIS con un script de Python

Si trabajas habitualmente con nubes de puntos LiDAR (como los datos del PNOA en España) y necesitas extraer la altura de la vegetación, calcular un CHM (Canopy Height Model) mediante la resta entre el modelo de vegetación y el terreno es una tarea recurrente.

Para agilizar este proceso de forma limpia y automatizada, hemos desarrollado un script modular en Python para QGIS que realiza toda la cadena de procesamiento de principio a fin:

• Filtra los puntos de suelo (Clase 2).
• Filtra los puntos de vegetación (Clases 3, 4 y 5).
• Genera el Modelo Digital del Terreno (MDT) mediante triangulación (TIN).
• Genera el Modelo de Vegetación (MDVeg) calculando el valor máximo.
• Calcula el CHM final y borra automáticamente todos los archivos intermedios y temporales del disco, dejando únicamente el ráster limpio.

---

🛠️ Cómo ejecutar este script en QGIS (Paso a Paso para principiantes)

Si nunca has utilizado la Consola de Python en QGIS, no te preocupes. Sigue estos sencillos pasos:

1. Carga el archivo LiDAR: Abre tu proyecto en QGIS y arrastra o carga tu fichero .laz o .las base en el panel de capas.
2. Abre la consola de Python: En el menú superior de QGIS, ve a Complementos > Consola de Python (o utiliza el atajo de teclado Ctrl + Alt + P). Se abrirá una pequeña ventana en la parte inferior o lateral de la pantalla.
3. Abre el Editor de Código: Dentro de la ventana de la consola, haz clic en el botón de Mostrar editor (tiene el icono de un trozo de código con un lápiz). Esto abrirá un panel donde puedes escribir y editar scripts largos.
4. Configura tus parámetros: Copia el código que encontrarás más abajo y pégalo en el editor. Modifica únicamente las dos primeras líneas con el nombre exacto de tu capa y la ruta de la carpeta de salida.
5. Ejecuta el script: Haz clic en el botón verde de Ejecutar script (o pulsa F5 dentro del editor). La consola procesará todo automáticamente y te avisará cuando haya terminado.

---

--- COPIAR DESDE LA LÍNEA DE ABAJO ---

import os from qgis.core import QgsProject import processing # ============================================================================== # CONFIGURACIÓN DEL USUARIO (MODIFICAR ESTOS DOS PARÁMETROS) # ============================================================================== # 1. Nombre (o parte del nombre) de la capa LAZ cargada en el panel de capas de QGIS nombre_capa_laz = "NOMBRE_DE_TU_CAPA_LAZ_AQUI" # 2. Carpeta del disco donde se guardará el ráster CHM final resultante directorio_salida = r"RUTAS\DE\TU\CARPETA\DE\DESTINO" # ============================================================================== # Comprobación y creación automática del directorio de salida os.makedirs(directorio_salida, exist_ok=True) # Rutas base para los archivos intermedios en disco ruta_laz_suelo = os.path.join(directorio_salida, "temp_suelo_clas_2.las") ruta_laz_veg = os.path.join(directorio_salida, "temp_veg_clas_3_4_5.las") ruta_mdt = os.path.join(directorio_salida, "temp_mdt_suelo.tif") ruta_mdveg = os.path.join(directorio_salida, "temp_mdveg_copas.tif") # Nombre del archivo ráster final de alturas (CHM) ruta_salida_chm = os.path.join(directorio_salida, "Altura_arboles_CHM_Final.tif") # Parámetros espaciales resolucion = 1.0 tile_size = 1000.0 # Buscar la capa LAZ de forma dinámica en el proyecto actual de QGIS capa_laz = None for layer in QgsProject.instance().mapLayers().values(): if nombre_capa_laz.lower() in layer.name().lower(): capa_laz = layer break if not capa_laz: print("❌ ERROR: No se encontró ninguna capa con ese nombre en el panel de capas.") print("👉 Asegúrate de cargar el archivo LAZ base en QGIS y de escribir bien su nombre.") else: print(f"✅ Capa LAZ detectada correctamente: {capa_laz.name()}") print("🚀 Iniciando procesamiento automático del CHM...") try: # 1. Filtrar Suelo (Clase 2) print(" -> [1/5] Filtrando puntos de suelo...") processing.run("pdal:filter", { 'INPUT': capa_laz, 'FILTER_EXPRESSION': 'Classification = 2', 'OUTPUT': ruta_laz_suelo }) # 2. Filtrar Vegetación (Clases 3, 4 y 5) print(" -> [2/5] Filtrando puntos de vegetación...") processing.run("pdal:filter", { 'INPUT': capa_laz, 'FILTER_EXPRESSION': 'Classification = 3 OR Classification = 4 OR Classification = 5', 'OUTPUT': ruta_laz_veg }) # 3. Generar MDT (TIN) print(" -> [3/5] Generando Modelo Digital del Terreno (MDT)...") processing.run("pdal:exportrastertin", { 'INPUT': ruta_laz_suelo, 'DIMENSION': 'Z', 'RESOLUTION': resolucion, 'TILE_SIZE': tile_size, 'OUTPUT': ruta_mdt }) # 4. Generar MDVeg (Valor Máximo) print(" -> [4/5] Generando Modelo de Vegetación (MDVeg)...") processing.run("pdal:exportraster", { 'INPUT': ruta_laz_veg, 'DIMENSION': 'Z', 'METHOD': 3, 'RESOLUTION': resolucion, 'TILE_SIZE': tile_size, 'OUTPUT': ruta_mdveg }) # 5. Calcular CHM (MDVeg - MDT) print(" -> [5/5] Calculando Canopy Height Model (CHM final)...") processing.run("gdal:rastercalculator", { 'INPUT_A': ruta_mdveg, 'BAND_A': 1, 'INPUT_B': ruta_mdt, 'BAND_B': 1, 'INPUT_C': None, 'BAND_C': None, 'INPUT_D': None, 'BAND_D': None, 'FORMULA': 'A - B', 'OUTPUT_CREATION_OPTIONS': '', 'OUTPUT': ruta_salida_chm }) print("✨ ¡Proceso completado con éxito!") print(f"📂 El ráster final se ha guardado en: {ruta_salida_chm}") except Exception as e: print(f"❌ Ocurrió un error durante la ejecución del proceso: {e}") finally: # Limpieza profunda automática de ficheros intermedios y auxiliares print("🧹 Limpiando archivos temporales y auxiliares de disco...") archivos_a_borrar = [ ruta_laz_suelo, ruta_laz_veg, ruta_mdt, ruta_mdveg, ruta_mdt + ".aux", ruta_mdt + ".aux.xml", ruta_mdveg + ".aux", ruta_mdveg + ".aux.xml" ] for archivo in archivos_a_borrar: if os.path.exists(archivo): try: os.remove(archivo) except: pass print("🔒 Limpieza profunda finalizada.")

--- COPIAR HASTA LA LÍNEA DE ARRIBA ---

 

 

VERSIÓN SENCILLA

--- COPIAR DESDE LA LÍNEA DE ABAJO ---

import os from qgis.core import QgsProject import processing # CONFIGURACIÓN laz_name = "NOMBRE_DE_TU_CAPA_LAZ_AQUI" out_dir = r"RUTAS\DE\TU\CARPETA\DE\DESTINO" os.makedirs(out_dir, exist_ok=True) layer = next((l for l in QgsProject.instance().mapLayers().values() if laz_name.lower() in l.name().lower()), None) if layer: s, v, mdt, mdv, chm = [os.path.join(out_dir, x) for x in ["t_suelo.las", "t_veg.las", "t_mdt.tif", "t_mdveg.tif", "CHM_Final.tif"]] processing.run("pdal:filter", {'INPUT': layer, 'FILTER_EXPRESSION': 'Classification = 2', 'OUTPUT': s}) processing.run("pdal:filter", {'INPUT': layer, 'FILTER_EXPRESSION': 'Classification = 3 OR Classification = 4 OR Classification = 5', 'OUTPUT': v}) processing.run("pdal:exportrastertin", {'INPUT': s, 'DIMENSION': 'Z', 'RESOLUTION': 1.0, 'TILE_SIZE': 1000.0, 'OUTPUT': mdt}) processing.run("pdal:exportraster", {'INPUT': v, 'DIMENSION': 'Z', 'METHOD': 3, 'RESOLUTION': 1.0, 'TILE_SIZE': 1000.0, 'OUTPUT': mdv}) processing.run("gdal:rastercalculator", {'INPUT_A': mdv, 'BAND_A': 1, 'INPUT_B': mdt, 'BAND_B': 1, 'FORMULA': 'A - B', 'OUTPUT': chm}) for f in [s, v, mdt, mdv, mdt+".aux", mdt+".aux.xml", mdv+".aux", mdv+".aux.xml"]: if os.path.exists(f): os.remove(f) print(f"¡Completado! Guardado en: {chm}") else: print("❌ No se encontró la capa LAZ.")

--- COPIAR HASTA LA LÍNEA DE ARRIBA ---