14

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?

  • What is the error message (if any)? – til_b Nov 5 '13 at 7:39
  • No error. I just didn't know how to filter, and was unsure if there was a better way to get LAS into array. – Barbarossa Nov 5 '13 at 15:30
13

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
)
  • 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 Nov 5 '13 at 22:06
  • These 5 records are in a 1-d array. – Barbarossa Nov 5 '13 at 22:48
  • 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). – om_henners Nov 6 '13 at 0:33
4

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.

1

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.

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.