Librería liblas para manejar datos LIDAR con Python

En un post pasado se hizo referencia, principalmente, al uso de la librería liblas para hacer un análisis exploratorio, a través de líneas de comando en consola de bash (Debian), de archivos *.las que contenían datos LIDAR. Se comenzó por la línea de comandos porque un desarrollador de liblas señalaba la dificultad de acceder a parámetros importantes del archivo usando Python y recomendaba más bien el uso de la librería laspy.

La documentación de la API de Python para liblas se encuentra aquí y obviamente lo fundamental es abrir el archivo y acceder a la información del header y a los datos. Para ello vamos a usar el mismo archivo anterior (Autzen_Stadium.zip); pero en su versión reproyectada a EPSG:4326 (lidar_rp.las).

El módulo liblas tiene varios sub módulos (‘color’, ‘core’, ‘file’, ‘get_version’, ‘guid’, ‘header’, ‘las’, ‘point’, ‘schema’, ‘srs’, ‘sys’, ‘version’, ‘vlr’) y, al igual que en PyQGIS, existen métodos dentro de ellos que permiten inter convertirlos entre si. Por ejemplo, en el código siguiente, se usa el método ‘File’ de ‘file’ para obtener un objeto de liblas.file.File. Este objeto se convierte en un objeto de liblas.header.Header mediante el método ‘header’ de liblas.file.File y en uno de liblas.srs.SRS mediante el método ‘srs’ del primero.

Por tanto, no es tan difícil extraer la proyección del archivo; a menos que no esté allí registrada o esté en un formato no compatible. La ejecución del script a continuación en la Python Console de QGIS:

from liblas import file

input_file = '/home/zeito/pyqgis_data/lidar_rp.las'

f = file.File(input_file, mode='r')

header_file = f.header
srs_file = header_file.srs

wkt = srs_file.wkt
proj4 = srs_file.proj4

print wkt
print proj4

crs1 = QgsCoordinateReferenceSystem()
crs2 = QgsCoordinateReferenceSystem()

crs1.createFromWkt(wkt)
crs2.createFromProj4(proj4)

print crs1.postgisSrid()
print crs2.postgisSrid()

f.close()

produce el resultado siguiente; donde se observa que el ‘srs’ extraido en formato WKT o proj4 permitió crear dos objetos de la clase QgsCoordinateReferenceSystem de PyQGIS que exhibieron correctamente el valor del código EPSG:4326 correspondiente al archivo *.las.

execfile(u'/home/zeito/pyqgis_scripts/lidar2.py'.encode('UTF-8'))
GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],
PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433],AUTHORITY["EPSG","4326"]]
+proj=longlat +datum=WGS84 +no_defs 
4326
4326

La inclusión del código a continuación:

.
.
.
points = [ QgsPoint(point.x, point.y) for point in f]

epsg = crs1

uri = "Point?crs=epsg:" + str(epsg) + "&field=id:integer""&index=yes"

mem_layer = QgsVectorLayer(uri,
                           'point',
                           'memory')

prov = mem_layer.dataProvider()

feats = [ QgsFeature() for i in range(len(points)) ]

for i, feat in enumerate(feats):
    feat.setAttributes([i])
    feat.setGeometry(QgsGeometry.fromPoint(points[i]))

prov.addFeatures(feats)

QgsMapLayerRegistry.instance().addMapLayer(mem_layer)

producirá una memory layer con la ubicación de los puntos tomados con el vuelo LIDAR. Modificando el archivo lidar_rp.las por el original lidar.las permitirá también obtener su memory_layer y observar que coinciden totalmente. Se puede comprobar que la herramienta las2las funciona adecuadamente.

lidar1

Ahora sólo resta probar que esta nube de puntos está ubicada en la zona correspondiente al Autzen Stadium de Oregon. Seleccionamos un punto arbitrario mediante el plugin de captura de coordenadas y lo empleamos en Google Maps; con las coordenadas invertidas (44.0565982143,-123.068589286) porque esa es la convención allí empleada. La imagen siguiente permite comprobar que efectivamente los puntos corresponden a la zona previamente señalada.

lidar2

Esta entrada fue publicada en Código Python, LIDAR, PyQGIS, SIG. Guarda el enlace permanente.

Una respuesta a Librería liblas para manejar datos LIDAR con Python

  1. Pingback: Caracterización de archivos *.las (datos LIDAR) con módulo liblas de Python | El Blog de José Guerrero

Responder

Por favor, inicia sesión con uno de estos métodos para publicar tu comentario:

Logo de WordPress.com

Estás comentando usando tu cuenta de WordPress.com. Cerrar sesión / Cambiar )

Imagen de Twitter

Estás comentando usando tu cuenta de Twitter. Cerrar sesión / Cambiar )

Foto de Facebook

Estás comentando usando tu cuenta de Facebook. Cerrar sesión / Cambiar )

Google+ photo

Estás comentando usando tu cuenta de Google+. Cerrar sesión / Cambiar )

Conectando a %s