Python – Finding Coordinate System Information from LAS/LiDAR File

coordinate systemgdallaslidarpython

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?

Best 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]]
Related Question