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?
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.