Extrayendo la altura de los edificios a partir de datos LiDAR

La obtención de modelos virtuales urbanos 3D se ha convertido hoy en día en la base para un cada vez mayor número de aplicaciones. El potencial de estos modelos va más allá la representación tridimensional de escenarios virtuales. Este tipo de información tiene cada vez más significación en campos como la planificación urbana, el control de variables medioambientales, la gestión de desastres, las telecomunicaciones o la administración de instalaciones entre otros.

Prueba de su importancia cada vez mayor son los esfuerzos en poder integrar los modelos 3D urbanos del mundo real que se utilizan en aplicaciones GIS para análisis geoespaciales con los, cada vez más implantados, modelos de diseños BIM de propósitos constructivos del ámbito de la ingeniería y la arquitectura, destinados a la administración de todo el ciclo de vida de los edificios. En ambos casos tanto las representaciones BIM como de objetos geoespaciales poseen sus propios modelos semánticos, siendo los más importantes los modelos de datos CityGML y IFC respectivamente. Del primero ya se habló en este blog de manera suscinta en un artículo anterior.

En este artículo nos centraremos en la metodología GIS para la extracción de datos de altura de edificaciones a partir de nubes de puntos LiDAR destinados a poder genera estos modelos 3D urbanos. Para este ejemplo utilizaremos software libre y datos provenientes de fuentes abiertas:

• Una base de datos PostGIS.
• La suite de herramientas de LAStools, o alternativamente LibLAS.
• Datos vectoriales de OpenStreetMap con la geometría 2D de edificios.
• Datos LiDAR del Centro de Descargas del CNIG.

Obteniendo e importando los datos básicos necesarios

El primero paso fundamental es descarga la cartografía con la que vamos a trabajar. En este ejemplo obtendremos de OpenStreetMap la geometría 2D de la base de los edificios del casco antiguo de la ciudad de Burgos, en España.

Existen diversas formas de descargar cartografía de OSM. Yo aquí voy a acceder directamente a los datos utilizando la API extendida de OSM y descargarlos en formato .osm con Wget desde el terminal de Linux.

wget -O burgos.osm http://api.openstreetmap.fr/xapi?way[building=*][bbox=-3.710,42.339,-3.696,42.344]

A continuación procedemos mediante GDAL/OGR a importar el archivo a nuestra base de datos PostGIS que hemos creado previamente en PostgreSQL. Aprovechamos para reproyectar los datos a EPGS:25830:

ogr2ogr -f PostgreSQL "PG:dbname=burgos user=postgres password=admin host=localhost port=5432" -t_srs EPSG:25830 -lco GEOMETRY_NAME=geom -lco FID=gid burgos.osm

OGR crea varias tablas para los diferentes tipos de geometrías. A nosotros la única que nos interesa es la tabla denominada «multipolygons», la cual renombraremos a «edificaciones» para darle un nombre más descriptivo y eliminaremos las restantes.

ALTER TABLE multipolygons RENAME TO edificaciones; 
DROP TABLE IF EXISTS lines, multilinestrings, other_relations, points;

A continuación procederemos a obtener los ficheros LiDAR para esta zona de Burgos del Centro de Descargas del CNIG. Aceptamos la licencia y descargamos todos los archivos .laz uno por uno, o bien en el caso de lotes más grandes podemos utilizar alguna alternativa que nos agilice el proceso como iMacros.

Procesando los datos LiDAR

Los archivos LiDAR descargados se encuentran en formato LAZ, un tipo de compresión estándar de ficheros LAS. LAS a su vez es un formato de archivo binario de intercambio no muy complejo, que además de la nube de puntos mantiene información específica de la naturaleza de los datos LiDAR que almacena.

Trabajar con datos LiDAR puede llegar a ser duro en cuanto al manejo y tiempo de procesado. Los datos LiDAR son un tipo de Big Data que dado su gran volumen presentan cierta dificultad en su gestión y análisis respecto a las prácticas tradiciones a las que estamos acostumbrados a aplicar en el ámbito GIS. Para un proyecto típico el volumen de este tipo de datos puede alcanzar fácilmente el orden de terabytes. Aunque para este ejemplo manejaremos unos pocos megabytes de información limitada al centro histórico de la ciudad de Burgos, como podremos imaginar el procesamiento de estos datos para el fin que perseguimos puede llegar a requerir un cálculo intensivo en función del dominio al que lo queremos aplicar.

Aquí recurriremos a la suite LASTools, que facilita un conjunto de herramientas para gestionar este tipos de archivos LAZ/LAS, obviando la interfaz gráfica para agilizar el proceso. En mi caso utilizaré en Linux los ejecutables .exe de Windows mediante Wine.

El primer paso será unir todos los ficheros en uno solo, para a continuación exportarlo a un archivo CSV con el que nos será más sencillo trabajar (este proceso nos puede llevar unos minutos):

emilio@glados:~/lidar$ wine lasmerge -i *.laz -o lidar.laz
emilio@glados:~/lidar$ wine las2txt -i lidar.laz -o lidar.csv -parse xyz

Si listamos las primeras líneas del archivo CSV creado observamos que por defecto se guardan los valores XYZ de la nube de puntos en columnas delimitadas por espacios y sin cabecera.

emilio@glados:~/lidar$ head -n 7 lidar.csv 
440000.30 4688395.84 845.52
440000.34 4688395.12 845.45
440000.38 4688394.50 845.53
440000.42 4688393.85 845.46
440000.45 4688393.23 845.50
440000.49 4688392.57 845.42
440000.52 4688391.99 845.43

El siguiente paso será importar este archivo CSV a PostgreSQL. Para ello creamos una nueva tabla en nuestra base de datos PostGIS en la que copiamos la nube de puntos.

CREATE TABLE lidar ( x DOUBLE PRECISION, y DOUBLE PRECISION, z DOUBLE PRECISION );
COPY lidar FROM '/tmp/lidar.csv' DELIMITER ' ';

A continuación añadimos un nuevo campo geométrico en la tabla «lidar» y lo actualizamos con las coordenadas de las columnas XY (este proceso nos volverá a llevar algo de tiempo). Finalmente creamos un índice espacial que nos permitirá acelerar las consultas que realicemos con posterioridad:

ALTER TABLE lidar ADD geom GEOMETRY(Point, 25830);
UPDATE lidar SET geom = ST_PointFromText('POINT(' || x || ' ' || y ||')', 25830);
CREATE INDEX lidar_idx ON lidar USING GIST (geom);

Calculando la altura de las edificaciones

Una vez tenemos los ingredientes, es hora de preparar la pizza. Vamos a determinar la altura del edificio discriminando lógicamente la elevación del terreno sobre el que se asienta. Antes de empezar debemos tener en cuenta que la altitud registrada por el LiDAR no permite distinguir entre un edificio situado en una vaguada o si este se encuentra a una cota elevada, por lo que debemos calcular la altura de este respecto a la elevación del terreno en su rasante.

Lo primero que haremos es crear nuevas columnas en la tabla «edificaciones» para almacenar los datos de elevaciones:

ALTER TABLE edificaciones ADD COLUMN elevacion_cubierta NUMERIC(6,2) NULL;
ALTER TABLE edificaciones ADD COLUMN elevacion_rasante NUMERIC(6,2) NULL;
ALTER TABLE edificaciones ADD COLUMN altura NUMERIC(6,2) NULL;

A continuación comenzamos calculando la elevación en lo alto del edificio. Para ello obtendremos la altura media de todos lo puntos que caen dentro de lo que es la planta del edificio. Aunque en este caso únicamente utilizaremos el promedio para no complicar el ejemplo, particularmente suelo también calcular la moda, con el fin de detectar cubiertas singulares o muy variables con patios interiores que pueden distorsionar la medida.

UPDATE edificaciones SET elevacion_cubierta = elevacion_media 
FROM (SELECT e.gid, AVG(l.z) AS elevacion_media FROM lidar l, edificaciones e WHERE ST_Within(l.geom, e.geom) GROUP BY e.gid) t 
WHERE edificaciones.gid = t.gid;

A continuación pasamos a calcular la elevación pero esta vez en la base de la edificación. La idea es tomar los valores que caen en el interior de un buffer de 1 metro de ancho situado a nivel de calle, pero no realizar este área de influencia próxima a la base, sino desplazarla 1 metro más hacia el exterior de la rasante del edificio. Con ello nos aseguramos, en la medida de los posible, que la estimación sea más precisa y evitar posibles errores derivados de obstáculos (terrazas, vehículos aparcados, mobiliario urbano, etc.).

Para ello creamos una copia de la tabla «edificaciones» sobre la que realizaremos este buffer.

CREATE TABLE edificaciones_buffer AS 
SELECT gid, elevacion_rasante, geom FROM edificaciones;

Y pasamos a calcular la diferencia de superficie entre un buffer de 2 metros en torno a los edificios y de otro de 1 metro también alrededor de la base de estos.

UPDATE edificaciones_buffer SET geom = buffer 
FROM ( SELECT e.gid, ST_Multi(ST_Difference(ST_Multi(ST_Buffer(e.geom, 2)), ST_Multi(ST_Buffer(e.geom, 1)))) AS buffer FROM edificaciones_buffer e WHERE e.gid > 0 ) t 
WHERE edificaciones_buffer.gid = t.gid;

Ya tenemos nuestra área de influencia de un metro separada otra más de la base del edificio. Procedemos a continuación a actualiza la tabla «edificaciones_buffer» con el valor mínimo de la elevación que hay dentro de este área de influencia.

UPDATE edificaciones_buffer SET elevacion_rasante = elevacion_min 
FROM ( SELECT e.gid, MIN(l.z) AS elevacion_min FROM lidar l, edificaciones_buffer e WHERE ST_within(l.geom, e.geom) GROUP BY e.gid ) t 
WHERE edificaciones_buffer.gid = t.gid;

Actualizamos ahora la columna en la tabla ‘edificaciones’ con ese valor:

UPDATE edificaciones e SET elevacion_rasante = e2.elevacion_rasante 
FROM edificaciones_buffer e2 
WHERE e.gid = e2.gid;

Y por último ya solo nos resta calcular la altura definitiva con la diferencia entre la medida obtenida en la cubierta y la capturada en la base del edificio:

UPDATE edificaciones SET altura = (elevacion_cubierta - elevacion_rasante);

Como conclusión podemos ver el resultado de nuestro trabajo realizando una extrusión prismática de la planta de los edificios en función del valor de altura que hemos calculado. La escena será un pequeño modelo urbano 3D con un nivel de abstracción LoD1 de bloques simple con cubiertas planas como el que se muestra a continuación: Haz clic aquí para ver la escena 3D en el navegador