2

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?

  • 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 – Aaron Jan 23 at 14:35
  • Thanks, and Liblas will work with the ArcGIS distribution of Python? I am having trouble importing the Liblas library. – Scott4993 Jan 23 at 15:29
  • 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? – Scott4993 Jan 23 at 17:09
  • I like to use anaconda to manage the installs. Check out the following: developers.arcgis.com/python/guide/install-and-set-up/… – Aaron Jan 23 at 17:23
  • Use laspy, it works fine with arcpy and it's pure python so it's easy to install. liblas is no longer maintained. – mikewatt Jan 23 at 17:47
3

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 
    
  • 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... – Howard Butler Jan 23 at 20:48
  • 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? – Scott4993 Jan 24 at 12:09
  • libLAS does not support LAS 1.4 files. It would need a significant update to do so. – Howard Butler Jan 24 at 13:20
  • 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? – Scott4993 Jan 25 at 12:44
  • Yes. You should be able to adapt the code as necessary. – Howard Butler Jan 25 at 16:32

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.