0

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()
5
  • 1
    You can also use lidR: r-lidar.github.io/lidRbook/gnd.html
    – BERA
    May 20, 2022 at 11:10
  • @BERA Im working on python and notebook :/ May 20, 2022 at 11:11
  • 1
    if number of return is 1, we can declare it as ground. Is it true? No
    – JRR
    May 20, 2022 at 14:11
  • @JRR what is the mathematical expression for to define ground then by usin this parametres, I couldnt find anything can help to solve this May 20, 2022 at 14:14
  • 1
    classifying ground point is complex and requires dedicated algorithm. I can't help in python but you can read the book chapter linked by @BERA. It is R but the overall idea is the same
    – JRR
    May 20, 2022 at 15:40

1 Answer 1

3

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.

Laser beam divergence resulting into multiple returns.  Source: https://medium.com/ai-ml-cv-in-enriching-digital-maps-navigation/analysing-an-aerial-lidar-point-cloud-dataset-6a31df9e1573

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

1
  • 1
    Vov, thats a lot. Thank you very much this was so important for me. I learn a lot from you. May 23, 2022 at 22:56

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.