Extracting projection from each LAS file in directory using ArcPy?

Question

I am using ArcPy to automate the conversion of LAS files to Raster in ArcGIS. I'd like to extract the spatial reference from each LAS file along the way.

Is it possible to do this using ArcPy tools?

I was able to install and import Laspy, but is it possible to convert a Laspy object to an ArcGIS spatial reference object?

Answer 1

Laspy isn't going to give you convenient access to the SRS in a form you can easily consume. LAS files can have either WKT or GeoTIFF keys as the coordinate system description. For consumption in Esri tools (and elsewhere), you always want the WKT.

The most convenient way to get the WKT from an LAS file is to use PDAL. The following script will read a filename from standard input and emit the WKT (or ESRI-morphed WKT) for the file to standard output. You should be able to use this to script for a directory of files using the batching tool of your choice.

1. Install PDAL via Conda

   conda install -c conda-forge pdal python-pdal

2. Save the following script to fetch-wkt.py or something:

   import subprocess
   import json
   import sys
   filename = sys.argv[1]
   args = ['pdal', 'info', filename, '--metadata']
   p = subprocess.Popen(args, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
   ret = p.communicate()
   if (ret[1]):
       # we are in error
       error = ret[1].decode('utf-8')
       sys.stderr.write('%s' % error)
       sys.exit(1)
   j = json.loads(ret[0].decode('utf-8'))
   srs = j['metadata']['srs']
   
   '''compound coordinate system will have horizontal and vertical'''
   comp_cs = srs['compoundwkt']
   horiz_cs = srs['horizontal']
   vert_cs = srs['vertical']
   
   '''Do this if you need it morphed to ESRI'''
   try:
       esri = sys.argv[2]
       import osr
       reference = osr.SpatialReference()
       success = reference.ImportFromWkt(comp_cs)
       if (success):
           print ("unable to import SRS for file '%s'" % filename)
       reference.MorphToESRI()
       esri = reference.ExportToWkt()
       sys.stdout.write('%s' % esri)
   except IndexError:
       sys.stdout.write('%s' % comp_cs)
   

3. Invoke it on a file

   ESRI-morphed version:

   python fetch-wkt.py myfile.las esri
   

   OGC WKT1 version:

   python fetch-wkt.py myfile.las