[GIS] Classifying LiDAR ground points using laspy

classificationlaspylidarnumpyscipy.spatial

I have a classified LiDAR point cloud (.las) with misclassifications. The walls of steep surfaces that should be ground (class 2) have been misclassified as vegetation (classes 3, 4, 5). I am attempting to reclassify these erroneous vegetation points as ground.

My approach is to:

  1. Find all single-return vegetation points within my data set.

  2. Iterate through each of the above points and get a list of all its neighbors within an 80 meter radius.

  3. If one of those neighbors is a ground point with a higher elevation, then I reclassify the original point from vegetation (3, 4, or 5) to ground (2).

I have tried to implement this workflow using a combination of laspy, numpy, and scipy.spatial:

import laspy
import numpy as np
from scipy.spatial import cKDTree
import matplotlib as plt

#read in my original, misclassified LAS file. Create a point_records array
in_file = laspy.file.File(r"misclassified.las", mode = "r")

point_records = in_file.points

#find all points that and single-return vegetation points (class 3, 4, or 5)

#note: num_returns is not included in the point_records above. So I am 
#creating new arrays here.
class_arr = in_file.raw_classification
return_arr = in_file.num_returns

single_veg_pts = np.where(np.logical_and(np.logical_or(class_arr == 3, 
                                                       class_arr == 4, 
                                                       class_arr == 5), return_arr == 1))


#get the filtered points' records
single_veg_pt_records = in_file.points[single_veg_pts]


#create an x, y, z array of all points to build a 3D KDTree
point_records_xyz = np.array((point_records['point']['X'],
point_records['point']['Y'], point_records['point']['Z'])).transpose()

ctree = cKDTree(point_records_xyz)

#loop through my single vegetation point records to compare their X, Y, Z 
#against my 3-D KDTree matrix. When a nearest neighbor is found...
for idx in single_veg_pt_records:
    neighbors = ctree.query_ball_point(np.array((idx['point']['X'],
    idx['point']['Y'],idx['point']['Z'])).transpose(), 8200)

for neighbor in neighbors:
   #if the neighbor point is class 2 and above idx, change the class of idx 
   #in point_records.

The final for loop is where I am hitting a wall.

How do I compare idx to all of its neighbors. How can I modify the original set of points (point_records) using the index value of idx?

Best Answer

I was able to come up with a solution. The following code works, but it feels like there is a much simpler way to do this operation. If anyone has a cleaner or more efficiect way to filter and reclassify lidar with laspy I'd love to accept that answer.

import laspy
import numpy as np
from scipy.spatial import cKDTree

in_file = laspy.file.File(input_file, mode = "r")
point_records = in_file.points.copy()

#filter all single return points of classes 3, 4, or 5 (vegetation)
single_veg = np.where(np.logical_and(np.logical_or(np.logical_or(
                                                   in_file.raw_classification == 3, 
                                                   in_file.raw_classification == 4), 
                                                   in_file.raw_classification == 5), in_file.num_returns == 1))

# pull out full point records of filtered points, and create an XYZ array for KDTree
single_veg_points = in_file.points[single_veg]     
single_veg_points_xyz = np.array((single_veg_points['point']['X'], 
                                  single_veg_points['point']['Y'],
                                  single_veg_points['point']['Z'])).transpose()

#filter class 2 points (ground) create XYZ array for KDTree           
ground_only = np.where(in_file.raw_classification == 2)
ground_points = in_file.points[ground_only]

ground_points_xyz = np.array((ground_points['point']['X'],
                              ground_points['point']['Y'],
                              ground_points['point']['Z'])).transpose()

#create a KDTree to query against
ctree = cKDTree(ground_points_xyz)

#For every single return veg point query against all points within 20 meters.
#If a higher elevation ground point is nearby, then change the classification 
#from vegetation to ground in the original point array.
for idx, record in enumerate(single_veg_points_xyz):
    neighbors = ctree.query_ball_point(record, 2000)
    for neighbor in neighbors:
        neighbor_z = ground_points[neighbor]['point']['Z']
        record_z = single_veg_points[idx]['point']['Z']
        if neighbor_z >= record_z:
            single_veg_points[idx]['point']['raw_classification'] = 2

#update points just once
point_records[single_veg] = single_veg_points

out_file = laspy.file.File(output_file, mode = "w", header=in_file.header)
out_file.points = point_records
out_file.close()
Related Question