An Easy Way to Work and Visualize Lidar Data in Python

Ingest, process, and Visualize 3D Point Cloud Data in Python

Working with 3D point cloud data is challenging. The collection of the points in one Lidar image is usually extensive and might reach millions of points depending on the image size. Handling the height dimension also poses other complexities.

A point cloud is a set of data points in space. The points may represent a 3D shape or object. Each point position has its set of Cartesian coordinates (X, Y, Z).

However, we have different tools that make working with Lidar data more straightforward than ever before. This blog post will teach you an easy way to ingest, manipulate, and visualize Lidar data in Python.

Reading & Accessing Lidar Data

In this tutorial, we use Laspy, a Python library for lidar LAS/LAZ IO, to ingest the point cloud data. Later, we will use open3D, a modern library for 3D data processing, to visualize the 3D lidar data. So let us import these libraries first.

import laspy
import open3d as o3d
import numpy as np

Let us read the lidar data with Laspy. Reading is done using laspy.read() function. You can also use laspy.open() if you only want the metadata but not the points.

las = laspy.read(“data/lidar.las”)
las

And once we read with laspy.read() , we get the metadata as well as the point cloud data. The las header for our lidar data looks like this.

<LasData(1.1, point fmt: <PointFormat(1, 0 bytes of extra dims)>, 277573 points, 1 vlrs)>

It contains the las header, the point format, point count and vlrs. We can also call them individually to get a separate output for each of them like this.

las.header
las.header.point_format
las.header.point_count
las.vlrs

The ingested file also contains the point data. First, let us examine the available features for the lidar file we have read.

list(las.point_format.dimension_names)

In this case, we have the following features available. Note this depends on the LAS version.

['X',
 'Y',
 'Z',
 'intensity',
 'return_number',
 'number_of_returns',
 'scan_direction_flag',
 'edge_of_flight_line',
 'classification',
 'synthetic',
 'key_point',
 'withheld',
 'scan_angle_rank',
 'user_data',
 'point_source_id',
 'gps_time']

We have X, Y, and Z for the point data, intensity, classification, GPS time, and some other essential dimensions. Let us, for example, see some of these dimensions.

las.X
las.intensity
las.gps_time

The output is an array of numbers representing the X, Intensity and GPS time.

# X
array([27799997, 27799997, 27799952, ..., 27775004, 27775002, 27775001])

#Intensity
array([10, 15, 12, ..., 30, 41, 35], dtype=uint16)

# gps_time
array([5880.963028, 5880.963032, 5880.978038, ..., 5886.739728,
       5886.739733, 5886.739738])

Classification is another essential dimension that you can call with las.classification that will provide an array of numbers. However, you can call the following code if you want to get the unique classification code numbers from your lidar data.

set(list(las.classification))

In our case, we have only four classes.

1 = Unassigned

2 = Ground

5 = High Vegetation

6 = Building

The classification scheme for Lidar data is more than the classes we have. Consult with the ArcGIS documentation page for full details of classes and corresponding LAS format.

Creating, Filtering, and Writing Point Cloud Data

To create 3D point cloud data, we can stack together with the X, Y, and Z dimensions, using Numpy like this.

point_data = np.stack([las.X, las.Y, las.Z], axis=0).transpose((1, 0))

point_data

The output is a 3-dimensional array and looks like this.

array([[ 68506260, 390272760,     31571],
       [ 68506057, 390272764,     31572],
       [ 68505659, 390272774,     31599],
       ...,
       [ 68747636, 389918893,     31996],
       [ 68747762, 389918889,     32296],
       [ 68747965, 389918898,     32223]])

We know the classification codes for our data, and therefore we can filter out any classes we want. For example, we can run the following code to filter out only the building class.

buildings = laspy.create(point_format=las.header.point_format, file_version=las.header.version)
buildings.points = las.points[las.classification == 6]

You can also write the building class points we have filtered out to a .las file quickly with the.write()function like this.

buildings.write(‘buildings.las’)

And there you have a new lidar image with only buildings saved in your local machine. Let us move on to point cloud data visualization.

3D Point Cloud Visualization

Laspy has no visualization methods so that we will use the open3d library. We first create the open3D geometries and pass the point data we have created earlier. Finally, we use the open3d visualization to draw geometries.

geom = o3d.geometry.PointCloud()
geom.points = o3d.utility.Vector3dVector(point_data)
o3d.visualization.draw_geometries([geom])

The nice thing about this is that it also works with Jupyter Notebook. It can render the 3D visualization in a separate output window. Here is the 3D visualization for the data with Open3D.

Conclusion

Great job. You know how to read lidar point cloud data and manipulate the data. You have also seen how to visualize the Lidar data in Python. For a basic understanding of Lidar data, where to access it freely, and its different formats, you can read my previous article on Lidar Data.