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?

• til_b: What is the error message (if any)?
• Barbarossa: No error. I just didn't know how to filter, and was unsure if there was a better way to get LAS into 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
)

• Barbarossa: The filter seems to work but only writes 5 records. I tried to filter only classes 1 and 2, and then tried to filter all but 1 and 2, both giving me only 5 results. Any ideas?
• Barbarossa: These 5 records are in a 1-d array.
• om_henners: Sorry, updated the code above as it requires specification of the axis to do the any calculation along (without it it performs the or across all dimensions of the array).

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

You can use pdal. Take a look at the example: Reading using Numpy Arrays

import pdal
data = "https://github.com/PDAL/PDAL/blob/master/test/data/las/1.2-with-color.las?raw=true"
pipeline = pdal.Reader.las(filename=data).pipeline()
count = pipeline.execute()
arrays = pipeline.arrays
arrays[0]
""" array([(637012.24, 849028.31, 431.66, 143, 1, 1, 1, 0, 1, 0, 0, 0, 0, -9., 132, 7326, 245380.78254963, 68, 77, 88),
 (636896.33, 849087.7 , 446.39, 18, 1, 2, 1, 0, 1, 0, 0, 0, 0, -11., 128, 7326, 245381.45279924, 54, 66, 68),
 (636784.74, 849106.66, 426.71, 118, 1, 1, 0, 0, 1, 0, 0, 0, 0, -10., 122, 7326, 245382.13595007, 112, 97, 114),
 ...,
 (637501.67, 853375.75, 417.52, 43, 1, 1, 0, 0, 1, 0, 0, 0, 0, 11., 124, 7334, 249772.21013494, 100, 96, 120),
 (637433.27, 853230.84, 424.08, 31, 1, 1, 0, 0, 1, 0, 0, 0, 0, 11., 125, 7334, 249772.70733372, 176, 138, 164),
 (637342.85, 853240.32, 423.92, 116, 1, 1, 1, 0, 1, 0, 0, 0, 0, 9., 124, 7334, 249773.20172407, 138, 107, 136)],
 dtype=[('X', '<f8'), ('Y', '<f8'), ('Z', '<f8'), ('Intensity', '<u2'), ('ReturnNumber', 'u1'), ('NumberOfReturns', 'u1'), 
 ('ScanDirectionFlag', 'u1'), ('EdgeOfFlightLine', 'u1'), ('Classification', 'u1'), ('Synthetic', 'u1'), 
 ('KeyPoint', 'u1'), ('Withheld', 'u1'), ('Overlap', 'u1'), ('ScanAngleRank', '<f4'), ('UserData', 'u1'), 
 ('PointSourceId', '<u2'), ('GpsTime', '<f8'), ('Red', '<u2'), ('Green', '<u2'), ('Blue', '<u2')]) """
#You can create a pandas dataframe like this, to make the data more readable etc.
import pandas as pd
df = pd.DataFrame(arrays[0])
df.head()
""" X Y Z Intensity ReturnNumber NumberOfReturns ScanDirectionFlag ... ScanAngleRank UserData PointSourceId GpsTime Red Green Blue
0 637012.24 849028.31 431.66 143 1 1 1 ... -9.0 132 7326 245380.782550 68 77 88
1 636896.33 849087.70 446.39 18 1 2 1 ... -11.0 128 7326 245381.452799 54 66 68
2 636784.74 849106.66 426.71 118 1 1 0 ... -10.0 122 7326 245382.135950 112 97 114
3 636699.38 848991.01 425.39 100 1 1 0 ... -6.0 124 7326 245382.825041 178 138 162
4 636601.87 849018.60 425.10 124 1 1 1 ... -4.0 126 7326 245383.388080 134 104 134
 """

Answer 4

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.

• Barbarossa: Thanks @Nicholas, I am using the liblas python library, which is closely tied to LASTools