Tuesday, 18 June 2019

gdal - Finding the coordinate system info from a LAS/LiDAR file with Python?


I found several examples of how one can query and manually set a coordinate system on a tif/raster file (GDAL/Python: How do I get coordinate system name from SpatialReference? or Is it possible to get the EPSG value from an OSR SpatialReference class using the OGR Python API? for example). But if I have a LAS file with Lidar point data and want to read its coordinate system programmatically with python, how could I go about doing that? I tried the two scripts below with no luck:


# first:

import gdal, ogr, os, osr
from laspy.file import File
driver = ogr.GetDriverByName('LOSLAS') ## maybe this driver would read LAS info, but hasn't worked
dataset = driver.Open(r'file.las') ## this could be the problem
layer = dataset.GetLayer()

spatialRef = layer.GetSpatialRef()
feature = layer.GetNextFeature()
geom = feature.GetGeometryRef()
spatialRef = geom.GetSpatialReference()
print spatialRef, geom

# second:

from osgeo import gdal,osr
dataset = gdal.Open(r'file.las') ## again this command doesn't work

spatialRef = osr.SpatialReference(wkt = dataset.GetProjection()) ## get projection function also throws an error
print 'projection = ', spatialRef
if spatialRef.IsProjected:
print spatialRef.GetAttrValue('projcs')
print spatialRef.GetAttrValue('geogcs')

Any thoughts on what I am doing wrong or what other libraries I could try? Again I want to be able to read the spatial referencing of a LAS file using python, without having to convert it to a raster first or use LAStools functionalities. Maybe there's a way to query this within the header?



Answer



In this question: Getting the projection of a LAS file using liblas?


Howard Butler (liblas developer) said that:




To be honest, I think the easiest way [Getting the projection of a LAS file using liblas python API] would be to shell out to lasinfo with the --xml argument, and then use ElementTree or some such to pluck out what you need from the XML.



You can also use lasinfo with the --no-check argument to filter ESPG code. I tried out it (bash console) with this Autzen_Stadium.zip las file and it works well.


lasinfo --no-check /home/zeito/pyqgis_data/lidar.las|grep -i epsg
AUTHORITY["EPSG","2994"],

Editing Note:


By using liblas in a script (Python Console of QGIS), I found out one way to read a file's spatial reference within its header's file.


from liblas import file


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

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

header_file = f.header
srs_file = header_file.srs

print srs_file.wkt


whose result was:


>>>execfile(u'/home/zeito/pyqgis_scripts/lidar2.py'.encode('UTF-8'))
LOCAL_CS["NAD83(HARN) / Oregon Lambert (ft)",GEOGCS["NAD83(HARN)",DATUM["unknown",SPHEROID["unretrievable - using WGS84",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]],AUTHORITY["EPSG","2994"],UNIT["foot",0.3048]]

No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...