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?

• Aaron: One approach would be to extract the WKT spatial reference using liblas (see link) and then convert the WKT to an arcpy spatial reference object. Check out this related post: gis.stackexchange.com/q/203392
• Scott4993: Thanks, and Liblas will work with the ArcGIS distribution of Python? I am having trouble importing the Liblas library.
• Scott4993: I used "pip install" to successfully install Liblas, but I'm having trouble importing it in a python script. I am using Python 2.7.8 with ArcGIS 10.3, is there a compatibility issue?
• Aaron: I like to use anaconda to manage the installs. Check out the following: developers.arcgis.com/python/guide/install-and-set-up/…
• mikewatt: Use laspy, it works fine with arcpy and it's pure python so it's easy to install. liblas is no longer maintained.

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 
   

• Howard Butler: The PDAL website has examples of how to create rasters in its workshop at pdal.io/workshop/exercises/analysis/dtm/dtm.html and pdal.io/workshop/exercises/analysis/rasterize/rasterize.html if you're up for skipping Arc altogether too...
• Scott4993: Thanks I will try this out. Does the script save the WKT to a file or does it output it into the terminal? Also, is LibLAS off the table as far as an option goes?
• Howard Butler: libLAS does not support LAS 1.4 files. It would need a significant update to do so.
• Scott4993: Thanks. And would it possible to call this PDAL script you gave me embedded within my ArcPy script and use the output in other processes?
• Howard Butler: Yes. You should be able to adapt the code as necessary.