Ground Classification, without any classified data? (Laspy, Python)

Question

I am trying to classify ground point in .las file. Points are never classified which is 0. I read something about, if the number of return is 1, we can declare it as ground. Is it true? las.number_of_returns == las.return_number. After that with that points I will create a KDTree so I can determine the all ground points.

Is this a correct approach?

I found the following code and I'm trying to change it.

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

in_file = laspy.read(las_filepath)
output_file = 'out.las'

point_records = in_file.points.copy()

single_veg = np.where(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['X'], 
                                  single_veg_points['Y'],
                                  single_veg_points['Z'])).transpose()

ground_only = np.where(in_file.number_of_returns == in_file.return_number)

ground_points = in_file.points[ground_only]

ground_points_xyz = np.array((ground_points['X'],
                              ground_points['Y'],
                              ground_points['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]['Z']
        record_z = single_veg_points[idx]['Z']
        if neighbor_z >= record_z:
            single_veg_points[idx]['raw_classification'] = 2

#update points just once
point_records[single_veg] = single_veg_points

out_file = laspy.open(output_file, mode='w', header=in_file.header)
out_file.points = point_records
out_file.close()

Answer 1

If you apply the filter "Number of Returns" == "Return Number" it means that the laser beam has encountered a flat surface which includes ground and also other flat surfaces like building rooftops, side-walls of the buildings, terraces, bridges, elevated-roads etc.

So to answer your question, if you filter all points with "Number of returns == Return Number" you will still have to filter the other flat surfaces from your data, maybe you can filter them by height so that all other points depicting the other flat surfaces like building rooftops, side-walls of the buildings, terraces, bridges, elevated-roads etc. gets eliminated.

Let's understand this further. Look at the diagram below. The Aerial vehicle with LiDAR sensor is shooting the lasers, but as the distance increases the laser ray spreads out, this phenomenon is called beam divergence. Due to the beam divergence, the spread out laser beam can encounter multiple surfaces and multiple returns are registered due to this. In this example the beam is going through a deciduous tree and it is registering returns from different branches that it encounters and also a last return (Return#5) is registered from the ground.

So for e.g. if Number of Returns==5 and Return Number ==5 as in the example above it should return a ground point.

However, if the Number of Returns ==1 and Return Number ==1 the point may have been registered from a flat surface like terrace, roof-top, elevated road or ground. So you will have to apply that height filter.

Alternate approach.

If you can assume that the majority of the points in your region of interest represents ground, you can do a Planar RANSAC and identify inliers and outliers of the RANSAC Plane. The inliers would be your ground points. Refer to this implementation https://github.com/abhilshit/aerial-lidar-classification/blob/main/app/PlanarRansac.py

Or you can also create an ensemble of the both methods to get better results.

Source of the image: My blog on understanding aerial lidar, check out for further detailed explanation https://medium.com/ai-ml-cv-in-enriching-digital-maps-navigation/analysing-an-aerial-lidar-point-cloud-dataset-6a31df9e1573