Cómo escribir en el header de un archivo *.las con el módulo liblas de Python

En el post pasado se consideró la manera de crear un archivo *.las a partir de una nube de puntos con el módulo liblas de Python. El script usado empleaba las opciones por defecto y, por tanto, mucha información útil no era grabada al header.

Una de ellas venía representada por la escala (scale) para las variables principales: X, Y, Z. Por defecto se asumía como la lista [1.0, 1.0, 1.0] y, por consiguiente, las referidas variables eran truncadas y almacenadas como enteros. Para resolver este problema se considera la asignación de la escala [0.01, 0.01, 0.01] y la precisión a sólo dos dígitos decimales (ver código) para los valores mínimos y máximos. Estos se determinan en tiempo de ejecución del script (porque la librería liblas no lo gerencia automáticamente).

Un aspecto verdaderamente problemático es la asignación de la proyección. Lo lógico es que cualquiera de esta dos asignaciones funcionase:

.
.
.
#first method
hout.srs.set_wkt(crs.toWkt())
hout.srs.set_proj4(crs.toProj4())

#second method
hout.srs.wkt = crs.toWkt()
hout.srs.proj4 = crs.toProj4()
.
.
.

pero fue imposible; a pesar de que la impresión de cualquiera de esos comandos se evaluaba a ‘True’ señalando que aparentemente no existía error de sintaxis.

Para resolver el problema de la proyección, que en este caso tiene un EPSG:32612, se usó el comando las2las con la sintaxis siguiente:

las2las -i /home/zeito/pyqgis_data/test.las -o /home/zeito/pyqgis_data/test_new.las --a_srs EPSG:32612 --t_srs EPSG:32612

o con la lectura del header de otro archivo que si la contiene; tal como se encuentra en el script a continuación:

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()
crs = layers[1].crs()

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

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

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

my_srs = f.header.srs

hout = header.Header()

#copying useful information to the header
hout.scale = [0.01, 0.01, 0.01]
hout.min = [round(min(coor_x), 2), round(min(coor_y), 2), round(min(values), 2)]
hout.max = [round(max(coor_x), 2), round(max(coor_y), 2), round(max(values), 2)]

#first method
#print hout.srs.set_wkt(crs.toWkt())
#hout.srs.set_proj4(crs.toProj4())

#second method
#hout.srs.wkt = crs.toWkt()
#hout.srs.proj4 = crs.toProj4()

hout.srs = my_srs

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!"

Después de ejecutado en la Python Console de QGIS el script anterior, la caracterización del archivo test.las correspondiente produjo los resultados esperados:

scale:  [0.01, 0.01, 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:  [479263.19, 4473425.22, 3385.0]
479263.19 4473425.22 3385.0
header file min:  [354975.4, 4414910.91, 1357.0]
354975.4 4414910.91 1357.0
>>>header_file.srs.wkt
'PROJCS["WGS 84 / UTM zone 12N",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"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],
PARAMETER["central_meridian",-111],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],
PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],
AUTHORITY["EPSG","32612"]]'
>>>header_file.srs.proj4
'+proj=utm +zone=12 +datum=WGS84 +units=m +no_defs '
Esta entrada fue publicada en Código Python, SIG, Software Libre. Guarda el enlace permanente.

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