Caracterización de archivos *.las (datos LIDAR) con módulo liblas de Python

En el post anterior se comenzó a usar la librería liblas de Python para manejar datos LIDAR con el objetivo de crear un shapefile de puntos con proyección que permitiese observar espacialmente la distribución de la nube de puntos de un archivo *.las. En este post se va a ahondar en la información que se puede obtener a partir del header y cómo está asociada a la nube de puntos.

Para encontrar las constantes, métodos y procedimientos del header basta con un dir(liblas.header.Header); tal como se presenta a continuación:

['DeleteVLR', 'GetVLR', '__class__', '__del__', '__delattr__', 
'__dict__', '__doc__', '__format__', '__getattribute__', 
'__hash__', '__init__', '__len__', '__module__', '__new__', 
'__reduce__', '__reduce_ex__', '__repr__', '__setattr__', 
'__sizeof__', '__str__', '__subclasshook__', '__weakref__', 
'add_vlr', 'compressed', 'count', 'data_format_id', 'data_offset', 
'data_record_length', 'dataformat_id', 'date', 'delete_vlr', 'doc', 
'encoding', 'file_signature', 'file_source_id', 'filesource_id', 
'get_compressed', 'get_count', 'get_dataformatid', 'get_dataoffset', 
'get_datarecordlength', 'get_date', 'get_filesignature', 
'get_filesourceid', 'get_global_encoding', 'get_guid', 
'get_headersize', 'get_majorversion', 'get_max', 'get_min', 
'get_minorversion', 'get_offset', 'get_padding', 
'get_pointrecordsbyreturncount', 'get_pointrecordscount', 
'get_projectid', 'get_recordscount', 'get_scale', 'get_schema', 
'get_softwareid', 'get_srs', 'get_systemid', 'get_version', 
'get_vlr', 'get_vlrs', 'get_xml', 'global_encoding', 'guid', 
'handle', 'header_length', 'header_size', 'major', 'major_version', 
'max', 'min', 'minor', 'minor_version', 'num_vlrs', 'offset', 'owned', 
'padding', 'point_records_count', 'point_return_count', 'project_id', 
'records_count', 'return_count', 'scale', 'schema', 'set_compressed', 
'set_count', 'set_dataformatid', 'set_dataoffset', 'set_date', 
'set_filesourceid', 'set_global_encoding', 'set_guid', 'set_majorversion', 
'set_max', 'set_min', 'set_minorversion', 'set_offset', 'set_padding', 
'set_pointrecordsbyreturncount', 'set_pointrecordscount', 'set_scale', 
'set_schema', 'set_softwareid', 'set_srs', 'set_systemid', 'set_version', 
'set_vlrs', 'software_id', 'srs', 'system_id', 'version', 'version_major', 
'version_minor', 'vlrs', 'xml']

El script siguiente imprime algunas de las características más relevantes del archivo *.las considerado (lidar_rp.las); así como prepara las propiedades asociadas a la nube de puntos (Dimensions). El significado de cada una de ellas puede averiguarse aquí.

from liblas import file

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

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

header_file = f.header

print "scale: ", header_file.scale
print "compressed: ", header_file.compressed
print "count: ", header_file.count
print "file signature: ", header_file.file_signature
print "filesource_id: ", header_file.filesource_id
print "file.header length: ", header_file.header_length
print "file.header size: ", header_file.header_size
print "major: ", header_file.major
print "file offset:", header_file.offset

points = [ point for point in f ]

#Dimensions
#1
coor_x = [ point.x for point in f]
#2
coor_y = [ point.y for point in f]
#3
coor_z = [ point.z for point in f]
#4
intensity = [ point.intensity for point in f ]
#5
return_number = [ point.return_number for point in f ]
#6
number_of_returns = [ point.number_of_returns for point in f ]
#7
scan_direction = [ point.scan_direction for point in f ]
#8
flightline_edge = [ point.flightline_edge for point in f ]
#9
classification = [ point.classification for point in f ]
#10
scan_angle = [ point.scan_angle for point in f ]
#11
user_data = [ point.user_data for point in f ]
#12
point_source_id = [ point.point_source_id for point in f ]
#13
time = [ point.time for point in f ]

print "header file max: ", header_file.max
print max(coor_x), max(coor_y), max(coor_z)

print "header file min: ", header_file.min
print min(coor_x), min(coor_y), min(coor_z)

f.close()

La ejecución del archivo anterior en la Python Console de QGIS produce el resultado siguiente:

scale:  [1e-06, 1e-06, 0.01]
compressed:  False
count:  10653
file signature:  LASF
filesource_id:  0
file.header length:  227
file.header size:  227
major:  1
file offset: [-0.0, -0.0, -0.0]
header file max:  [-123.0625, 44.0625, 180.97]
-123.0625 44.0625 180.97
header file max:  [-123.075, 44.05, 123.93]
-123.075 44.05 123.93

Puede observarse que el conteo de puntos arroja 10653 valores y que por esta razón no se imprime ninguna de las listas correspondientes a los parámetros “Dimensions” (X, Y, Z, intensity, etc). No obstante, en las últimas cuatro líneas anteriores se imprimen los valores mínimos y máximos que registra el header las cuales se comparan con los valores obtenidos mediante las métodos ‘min’ y ‘max’ de las listas Python correspondientes (coor_x, coor_y, coor_z).

Por otra parte, para ejemplificar la importancia de la caracterización, vamos a tomar la lista ‘classification’. En http://www.liblas.org/python.html encontramos que:

lidar

En la Python Console exploramos lo siguiente:

>>>min(classification)
1
>>>max(classification)
2

por lo que inmediatamente se tiene que los puntos de este archivo caen en la categoría 1 (Unclassified) o 2 (Ground).

Para finalizar este post, a modo ilustrativo (porque se pueden incluir de una sola vez todos los parámetros de Dimension), voy a introducir la coordenada Z (coor_x) en el campo del mismo nombre de una memory layer de puntos LIDAR utilizando el código a continuación. El archivo de partida (lidar.las) tiene una proyección en pies y la coordenada z se convierte en metros.

from liblas import file
#lidar_newnorm.las
input_file = '/home/zeito/pyqgis_data/lidar.las'

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

header_file = f.header
srs_file = header_file.srs

wkt = srs_file.wkt

crs1 = QgsCoordinateReferenceSystem()

crs1.createFromWkt(wkt)

points = [ QgsPoint(point.x, point.y) for point in f]
coor_z = [ (point.z)*0.3048 for point in f] #change feet in meters 

epsg = 2994

uri = "Point?crs=epsg:" + str(epsg) + "&field=id:integer&field=Z:double""&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, coor_z[i]])
    feat.setGeometry(QgsGeometry.fromPoint(points[i]))

prov.addFeatures(feats)

QgsMapLayerRegistry.instance().addMapLayer(mem_layer)

f.close()

Después de grabar la memory layer como autzen_2994.shp se tiene:

autzen_2994

Ahora se podría generar un dem a partir de la nube de puntos.

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

2 respuestas a Caracterización de archivos *.las (datos LIDAR) con módulo liblas de Python

  1. Pingback: Cómo crear un archivo *.las a partir de una nube de puntos con el módulo liblas de Python | El Blog de José Guerrero

  2. Pingback: Cómo escribir en el header de un archivo *.las con el 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