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?