Converting LAS file to numpy array?

Question

I have begun learning how to manipulate LAS data in python and wanted to see how others handle LAS files. I would like to read the points (I am using a numpy array), and filter out classes 1 and 2 (unclassified and ground) to a separate array. I have the following code but cannot seem to get the points filtered.

# Import modules
from liblas import file
import numpy as np

if __name__=="__main__":
    '''Read LAS file and create an array to hold X, Y, Z values'''
    # Get file
    las_file = r"E:\Testing\ground_filtered.las"
    # Read file
    f = file.File(las_file, mode='r')
    # Get number of points from header
    num_points = int(f.__len__())
    # Create empty numpy array
    PointsXYZIC = np.empty(shape=(num_points, 5))
    # Load all LAS points into numpy array
    counter = 0
    for p in f:
        newrow = [p.x, p.y, p.z, p.intensity, p.classification]
        PointsXYZIC[counter] = newrow
        counter += 1

I have seen arcpy.da.featureClassToNumpyArray, but I did not want to import arcpy nor have to convert to shapefile.

How else can I filter/read LAS data into a numpy array?

Answer 1

Your PointsXYZIC is now a numpy array. Which means you can use numpy indexing to filter the data you're interested in. For example you can use an index of booleans to determine which points to grab.

#the values we're classifying against
unclassified = 1
ground = 2

#create an array of booleans
filter_array = np.any(
    [
        PointsXYZIC[:, 4] == unclassified, #The final column to index against
        PointsXYZIC[:, 4] == ground,
    ],
    axis=0
)

#use the booleans to index the original array
filtered_rows = PointsXYZIC[filter_array]

You should now have a numpy array with all the values where the data is unclassified or ground. To get the values that have been classified you could use:

filter_array = np.all(
    [
        PointsXYZIC[:, 4] != unclassified, #The final column to index against
        PointsXYZIC[:, 4] != ground,
    ],
    axis=0
)

Answer 2

Use laspy to read LAS files and easily return the data as numpy arrays you can interact with. laspy is pure python, is almost as fast as libLAS, has more features than the libLAS Python bindings, and is much easier to deploy.

Answer 3

I apologise if you already know of this, but LASTools is a fantastic Open Source tool which now integrates with both ArcGIS and QGIS 2.0 - This has options to filter the data in the manner you are looking at.