Cómo crear un archivo *.las a partir de una nube de puntos con el módulo liblas de Python

Explorando archivos *.las de muestra, he podido detectar que algunos no presentan en el header parte de la información que deberían. Por esta razón, ya que en el post anterior se produjo un código para la caracterización, voy a crear un script plantilla para corroborar que información se graba por defecto al convertir un shapefile, con 10653 puntos aleatorios y los valores de un ráster dem que subyacen sobre ellos, en un archivo *.las.

Las capas a usar lucen en la Map View de QGIS de la manera siguiente:

copy1

La idea es emplear métodos de PyQGIS para obtener las coordenadas de los puntos y la elevación correspondiente a las capas involucradas y luego, mediante metodos de liblas, crear un header y escribir al archivo *.las de salida toda la información pertinente a éste y la de los puntos 3D. El código es el siguiente:

from liblas import file, header, point

mapcanvas = iface.mapCanvas()

layers = mapcanvas.layers()

my_points = [ feat.geometry().asPoint() for feat in layers[0].getFeatures()]
coor_x = [ feat.geometry().asPoint()[0] for feat in layers[0].getFeatures() ]
coor_y = [ feat.geometry().asPoint()[1] for feat in layers[0].getFeatures() ]

#raster layer
rprovider = layers[1].dataProvider()

values = [ rprovider.identify(my_point, QgsRaster.IdentifyFormatValue).results()[1]
           for my_point in my_points ]

hout = header.Header()

outfil = "/home/zeito/pyqgis_data/test.las"

lout = file.File(outfil, mode = 'w', header = hout)

print "Writing  to output las file..."

n = layers[0].featureCount()

for i in range(n):
    
    pt = point.Point()
    pt.set_header(hout)

    pt.x = float(coor_x[i])
    pt.y = float(coor_y[i])
    pt.z = float(values[i])
    
    lout.write(pt)

lout.close()

print "Done!"

Cuando se ejecuta el código anterior en la Python Console de QGIS se produce un archivo de salida, test.las, que contiene la nube de puntos. Ahora se va a caracterizar con el código que se encuentra aquí. El resultado, en este caso, es el siguiente:

scale:  [1.0, 1.0, 1.0]
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:  [0.0, 0.0, 0.0]
479263.0 4473425.0 3385.0
header file min:  [0.0, 0.0, 0.0]
354975.0 4414911.0 1357.0

Se observa que la escala es [1.0, 1.0, 1.0] por lo que la precisión de los valores almacenadas en las raw_x, raw_y y raw_z será baja. Se puede comprobar imprimiendo un pequeño slice para cada una de las variables almacenadas:

>>>coor_x[0:5]
[445435.0, 429074.0, 363100.0, 405446.0, 406693.0]
>>>coor_y[0:5]
[4456461.0, 4417073.0, 4442217.0, 4433562.0, 4428798.0]
>>>coor_z[0:5]
[1455.0, 1526.0, 1756.0, 1684.0, 1907.0]

y comprobando con los valores originales (antes de la caracterización):

>>>coor_x[0:5]
[445434.5227452647, 429074.22135222325, 363100.4648394152, 405446.2243630681, 406693.3111839685]
>>>coor_y[0:5]
[4456460.640777081, 4417072.924192545, 4442217.153909652, 4433561.757064614, 4428797.662633395]
>>>coor_z[0:5]
[1455.0, 1526.0, 1756.0, 1684.0, 1907.0]

Por otra parte, se observa que el header no contiene valores mínimos ni máximos por lo que esta tarea no es asumida automáticamente por la librería. Además, como es lógico suponer, tampoco existe información de la proyección.

En el post siguiente se tratará de añadir esta información al header, así como, cualquier otra que se crea pertinente.

Esta entrada fue publicada en PyQGIS, SIG, Software Libre. Guarda el enlace permanente.

Una respuesta a Cómo crear un archivo *.las a partir de una nube de puntos con el módulo liblas de Python

  1. 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