May 27

Using Arcpy to convert las files to mulitpoint shapefile

Recently, I have been doing some heavy processing of las files to eventually generate a digital elevation model (DEM). I found myself in a dilemma of having thousands of las files and needing a way to convert them to mulitpoint  shapefiles (mps). I have used lastools in the past, but typically steer away from most of the utilities due to the licensing and the introduction of artificial noise if you have not paid for a full license – I don’t want to take any chances with the end product. ArcGIS also has a tool to convert las to multipoint, but I did not want to go this route for two reasons: 1) The number of points are much too large for a single mps and 2) I don’t want a mps for each las file. This leads to the following scenario; I need to combine, say 30 las files to a single mps, which is a tedious task to by hand in ArcGIS. However, it is easily done with a little python script.

The code below combines 40 las files to a single mps.

 

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
# Import system modules
import arcpy
from arcpy import env
import os
import glob
 
workspace = 'Input Directory'
env.workspace = workspace
odir = 'Output Directory'
arcpy.CheckOutExtension("3D")
ptspacing = 1.8         # Average point spacing
classcode = "2"         # Ground
sr = arcpy.SpatialReference("NAD 1983 UTM Zone 17N") # Coodinate System
 
numFiles = 0
numshpFiles = 0
for file in glob.glob('*.las'):
    numFiles = numFiles + 1
 
totalnumFiles = numFiles
print 'Starting to convert ' + str(totalnumFiles) + ' files...'
 
numFiles = 0
flist = ''
for file in glob.glob('*.las'):
    numFiles = numFiles + 1
    if numFiles != 1:
        flist = flist + ';' + file
    else:
        flist = file
    if numFiles % 40 == 0  or numFiles == totalnumFiles::
        numshpFiles = numshpFiles + 1
        ofile = odir + "/shpGroup_" + str(numshpFiles)
        print ofile
        arcpy.LASToMultipoint_3d(flist, ofile, ptspacing, classcode)
        arcpy.DefineProjection_management(ofile+'.shp', sr)
        flist = ''  # Clear file list

 

Permanent link to this article: http://www.mattbilskie.com/using-arcpy-to-convert-las-files-to-mulitpoint-shapefile/