3

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?

  • Any chance of you posting some data which would allow running your code? – Stephen Rauch Sep 19 '17 at 4:17
  • @StephenRauch here is a small subset of my dataset in LAS format. link – rosswin Sep 19 '17 at 18:14
4

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()

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.