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