4

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?

5

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]]  
  • Yes, I had read about this as well but need to figure out a LAS file's coordinate system information with python code (as in a script format), rather than use the commands or GUI functionalities that tools such as liblas and LAStools provide. So is there any way to read a file's spatial reference within its header or perhaps via using the gdal libraries? That way I could then use that information to export it elsewhere, for example. – spartmar Jul 26 '16 at 18:20
  • Yes, I found out one way. See my Editing Note. – xunilk Jul 28 '16 at 10:19

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.