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:
Find all single-return vegetation points within my data set.
Iterate through each of the above points and get a list of all its neighbors within an 80 meter radius.
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?