An Introduction to LiDAR Simulation with HELIOS++¶

I: Airborne laser scanning (ALS)

LiDAR simulation can be useful for multiple applications:

• Acquisition planning
• Method development and evaluation
• Training data generation (for machine learning)

In this notebook, we will simulate an airborne LiDAR acquisition, to demonstrate how HELIOS++ can be used for acquisition planning.

You can find a comprehensive documentation in the HELIOS Wiki.

In [1]:

# Imports
import os
import sys
from pathlib import Path
import time
import matplotlib.pyplot as plt
import matplotlib
from ipywidgets import interact
import pyhelios
from pyhelios import SimulationBuilder
from pyhelios.util import flight_planner
import helper_funcs
import numpy as np
from osgeo import gdal
import rasterio as rio
import laspy

# Imports import os import sys from pathlib import Path import time import matplotlib.pyplot as plt import matplotlib from ipywidgets import interact import pyhelios from pyhelios import SimulationBuilder from pyhelios.util import flight_planner import helper_funcs import numpy as np from osgeo import gdal import rasterio as rio import laspy

In [2]:

# making sure we have PyQt5/PyQt6 for interactive plotting - this might take a few minutes, be patient
!{sys.executable} -m pip install PyQt5

# making sure we have PyQt5/PyQt6 for interactive plotting - this might take a few minutes, be patient !{sys.executable} -m pip install PyQt5

Requirement already satisfied: PyQt5 in c:\users\3dgeo\.local\share\mamba\envs\groundfiltering\lib\site-packages (5.15.11)
Requirement already satisfied: PyQt5-sip<13,>=12.15 in c:\users\3dgeo\.local\share\mamba\envs\groundfiltering\lib\site-packages (from PyQt5) (12.17.1)
Requirement already satisfied: PyQt5-Qt5<5.16.0,>=5.15.2 in c:\users\3dgeo\.local\share\mamba\envs\groundfiltering\lib\site-packages (from PyQt5) (5.15.2)

To execute LAStools commands, we need to provide the LAStools directory. This may be a different location on your computer, so you will have to provide the correct path:

In [3]:

# configure LAStools root
lastools_root = "C:/LAStools"

# configure LAStools root lastools_root = "C:/LAStools"

1. The virtual scene¶

The first component of a simulation is the 3D input scene. HELIOS++ supports loading different file format, including Wavefront OBJ triangle meshes, GeoTIFF rasters or XYZ point clouds, which will be converted to voxel models. In this tutorial, we will use GeoTIFF for loading terrain and XYZ files for loading voxelized vegetation.

Extracting vegetation points

Using LAStools "las2txt", we convert the LAZ file of the study size to an ASCII file and extract only vegetation points, which have classification values of 0.

Note: With the exclamation mark (!), we can issue a shell command and hence run the command line tools of LAStools.

In [4]:

Path("data/output").mkdir(parents=True, exist_ok=True)
!$lastools_root/bin/las2txt64.exe -i data/StA_last.laz -o data/output/StA_last_vegetation.xyz -keep_classification 0 -parse xyz

Path("data/output").mkdir(parents=True, exist_ok=True) !$lastools_root/bin/las2txt64.exe -i data/StA_last.laz -o data/output/StA_last_vegetation.xyz -keep_classification 0 -parse xyz

Defining scene part files

We now define the terrain file (GeoTIFF) and the vegetation point cloud (xyz ASCII), as well as a voxel size for voxelizing the vegetation points.

In [5]:

terrain = "data/StA_last_dtm.tiff"
vegetation = "data/output/StA_last_vegetation.xyz"
voxel_size = 0.5

# Optional: Deriving a spatial subset so simulations run faster. You may set to "False" if you have a device with plenty of RAM (>= 64 GB) that can simulate the entire scene.
spatial_subset = True

terrain = "data/StA_last_dtm.tiff" vegetation = "data/output/StA_last_vegetation.xyz" voxel_size = 0.5 # Optional: Deriving a spatial subset so simulations run faster. You may set to "False" if you have a device with plenty of RAM (>= 64 GB) that can simulate the entire scene. spatial_subset = True

In [6]:

if spatial_subset:
    # bounding box of subset
    bbox = [20103.0, 312815.0, 20596.0, 312410.0]

    # spatial subset of terrain using gdal_translate
    terrain_sub = "data/output/StA_last_filtered_sub.tiff"
    !gdal_translate -projwin 20103.0 312815.0 20596.0 312410.0 $terrain $terrain_sub

    # spatial subset of vegetation
    vegetation_sub = "data/output/StA_last_vegetation_sub.xyz"
    vegetation_arr = np.loadtxt(vegetation)
    idx_sub = (vegetation_arr[:, 0] > bbox[0]) & (vegetation_arr[:, 0] < bbox[2]) & (vegetation_arr[:, 1] > bbox[3]) & (vegetation_arr[:, 1] < bbox[1])
    vegetation_arr_sub = vegetation_arr[idx_sub, :]
    np.savetxt(vegetation_sub, vegetation_arr_sub, fmt="%.6f")

    terrain = terrain_sub
    vegetation = vegetation_sub

if spatial_subset: # bounding box of subset bbox = [20103.0, 312815.0, 20596.0, 312410.0] # spatial subset of terrain using gdal_translate terrain_sub = "data/output/StA_last_filtered_sub.tiff" !gdal_translate -projwin 20103.0 312815.0 20596.0 312410.0 $terrain $terrain_sub # spatial subset of vegetation vegetation_sub = "data/output/StA_last_vegetation_sub.xyz" vegetation_arr = np.loadtxt(vegetation) idx_sub = (vegetation_arr[:, 0] > bbox[0]) & (vegetation_arr[:, 0] < bbox[2]) & (vegetation_arr[:, 1] > bbox[3]) & (vegetation_arr[:, 1] < bbox[1]) vegetation_arr_sub = vegetation_arr[idx_sub, :] np.savetxt(vegetation_sub, vegetation_arr_sub, fmt="%.6f") terrain = terrain_sub vegetation = vegetation_sub

Input file size is 801, 751
0...10...20...30...40...50...60...70...80...90...100 - done.

Writing the scene

HELIOS++ uses XML files for configuring scanners, platforms, surveys and scenes. Below, we write such an XML file and fill in the values we defined before (terrain, vegetation, voxel_size). Note how each scene part is embedded in a <part> tag and uses a different "loader" filter depending on the file format.

In [7]:

# Writing the scene XML file
scene_content = f"""<?xml version="1.0" encoding="UTF-8"?>
<document>
    <scene id="helios_scene" name="HELIOS scene">
        
        <part>
            <filter type="geotiffloader">
                <param type="string" key="filepath" value="{terrain}" />
            </filter>
        </part>
        <part>
            <filter type="xyzloader">
                <param type="string" key="filepath" value="{vegetation}" />
                <param type="string" key="separator" value=" " />
                <param type="double" key="voxelSize" value="{voxel_size}" />
            </filter>
        </part>
    </scene>
</document>
"""

# Writing the scene XML file scene_content = f""" """

In [8]:

# save scene file to the current directory
scene_file = "StA_helios_scene.xml"
with open(scene_file, "w") as f:
    f.write(scene_content)

# save scene file to the current directory scene_file = "StA_helios_scene.xml" with open(scene_file, "w") as f: f.write(scene_content)

2. Selection of platform and scanner¶

In this scenario, we will select an airborne platform ("Cirrus SR-22"), similar to the original acquisition. As scanner, we choose the "RIEGL LMS-Q560".

Feel free to try out other scanners in python/pyhelios/data/scanners_als.xml but please notice that different scanner types require different survey specifications.

In [9]:

platform = "sr22"
scanner = "riegl_lms-q560"

platform = "sr22" scanner = "riegl_lms-q560"

3. Specification of the survey¶

First, we want to mimic the original acquisition.

Waypoints (trajectory)

For simplification, we will include only two flight passes, which you can define interactively below. A plot will open in a new window, where you can select the waypoints by clicking on the DEM image.

In [10]:

n_pos = 4
dtm, tf, bounds, origin_left, origin_bottom, origin_right, origin_top = helper_funcs.read_raster(terrain)

n_pos = 4 dtm, tf, bounds, origin_left, origin_bottom, origin_right, origin_top = helper_funcs.read_raster(terrain)

In [11]:

%matplotlib qt
waypoints = helper_funcs.interactive_flight_trajectory(dtm, tf, n_pos=n_pos)
waypoints = np.array(waypoints)

%matplotlib qt waypoints = helper_funcs.interactive_flight_trajectory(dtm, tf, n_pos=n_pos) waypoints = np.array(waypoints)

Ignoring fixed x limits to fulfill fixed data aspect with adjustable data limits.

You have chosen your flight lines.
[[20175.89, 312521.69], [20471.51, 312780.97], [20568.74, 312686.69], [20256.42, 312450.97]]

Let's plot the chosen configured flight paths on top of the terrain.

In [12]:

%matplotlib inline
plt.imshow(dtm, cmap="terrain", extent=[origin_left, origin_right, origin_bottom, origin_top])
plt.plot(waypoints[:, 0], waypoints[:, 1])
plt.show()

%matplotlib inline plt.imshow(dtm, cmap="terrain", extent=[origin_left, origin_right, origin_bottom, origin_top]) plt.plot(waypoints[:, 0], waypoints[:, 1]) plt.show()

Acquisition settings

Now we configure the acquisition settings.

In [13]:

# compute mean DTM height
dtm_height = np.mean(dtm)

pulse_freq = 100_000
scan_freq = 80
flight_v = 60
alt = dtm_height + 1000  # summing up mean DTM height and desired height above ground level (AGL)
scan_angle = 45

# compute mean DTM height dtm_height = np.mean(dtm) pulse_freq = 100_000 scan_freq = 80 flight_v = 60 alt = dtm_height + 1000 # summing up mean DTM height and desired height above ground level (AGL) scan_angle = 45

Writing the survey file

Finally, we can write the survey file, which is also an XML file. You can find further examples in the HELIOS++ demo data folder.

In [14]:

legs = ""
active = True
for i in range(n_pos):
    if i % 2 == 0:
        active = "true"
    else:
        active = "false"
    
    x, y = waypoints[i]
    
    legs += f'''
    <leg>
            <platformSettings x="{x}" y="{y}" template="platform_als"/>
            <scannerSettings template="scanner_als" active="{active}"/>
    </leg>'''

legs = "" active = True for i in range(n_pos): if i % 2 == 0: active = "true" else: active = "false" x, y = waypoints[i] legs += f''' '''

Note how the survey below links to the scanner and scene that we defined earlier.

In [15]:

survey_content = f'''<?xml version="1.0" encoding="UTF-8"?>
<document>
    <platformSettings id="platform_als" movePerSec_m="{flight_v}" z="{alt}" />
    <scannerSettings active="true" id="scanner_als" pulseFreq_hz="{pulse_freq}" scanAngle_deg="{scan_angle}" scanFreq_hz="{scan_freq}" trajectoryTimeInterval_s="0.01"/>
    <survey name="als_survey" platform="data/platforms.xml#{platform}" scanner="data/scanners_als.xml#{scanner}" scene="{scene_file}#helios_scene">
    <FWFSettings beamSampleQuality="3" winSize_ns="1.5"/>
    {legs}
    </survey>
</document>
''' 

survey_content = f''' {legs} '''

In [16]:

survey_file = "als_StA_survey.xml"
with open(survey_file, 'w') as f:
    f.write(survey_content)

survey_file = "als_StA_survey.xml" with open(survey_file, 'w') as f: f.write(survey_content)

4. Running the simulation¶

In [17]:

pyhelios.loggingDefault()
simB = SimulationBuilder(str(survey_file), ['assets/'], 'helios_output/')
simB.setLasOutput(True)
simB.setZipOutput(True)

# build the simulation
sim = simB.build()

# Start the simulation.
start_time = time.time()
sim.start()

if sim.isStarted():
    print(f'Simulation has started!\nSurvey Name: {sim.sim.getSurvey().name}\n{sim.sim.getScanner().toString()}')

while sim.isRunning():
    duration = time.time()-start_time
    mins = duration // 60
    secs = duration % 60
    print('\r'+f'Simulation has been running for {int(mins)} min and {int(secs)} sec. Please wait.', end='')
    time.sleep(1)

if sim.isFinished():
    print('\n'+'Simulation has finished!')

output = sim.join()
meas, traj = pyhelios.outputToNumpy(output)

pyhelios.loggingDefault() simB = SimulationBuilder(str(survey_file), ['assets/'], 'helios_output/') simB.setLasOutput(True) simB.setZipOutput(True) # build the simulation sim = simB.build() # Start the simulation. start_time = time.time() sim.start() if sim.isStarted(): print(f'Simulation has started!\nSurvey Name: {sim.sim.getSurvey().name}\n{sim.sim.getScanner().toString()}') while sim.isRunning(): duration = time.time()-start_time mins = duration // 60 secs = duration % 60 print('\r'+f'Simulation has been running for {int(mins)} min and {int(secs)} sec. Please wait.', end='') time.sleep(1) if sim.isFinished(): print('\n'+'Simulation has finished!') output = sim.join() meas, traj = pyhelios.outputToNumpy(output)

SimulationBuilder is building simulation ...
SimulationBuilder built simulation in 31.579064099998504 seconds
Simulation has started!
Survey Name: als_survey
Scanner: riegl_lms-q560
Device[0]: riegl_lms-q560
	Average Power: 4 W
	Beam Divergence: 0.5 mrad
	Wavelength: 1064 nm
	Visibility: 23 km

Simulation has been running for 0 min and 1 sec. Please wait.
Simulation has finished!

5. Comparison to the real data¶

Using side views and top views, we visualise the simulated and real-world data to compare them.

In [18]:

# Load the real-world data
las = laspy.read("data/StA_last.laz")
ref = np.array([las.x, las.y, las.z]).T

# Load the real-world data las = laspy.read("data/StA_last.laz") ref = np.array([las.x, las.y, las.z]).T

In [19]:

%matplotlib inline
# plotting a cross section
fig, (ax1, ax2) = plt.subplots(nrows=2, ncols=1, figsize=(20, 8))

for pc, ax, title in zip([meas, ref], [ax1, ax2], ["Simulated", "Real-world"]):
    if spatial_subset:
        idx = (pc[:, 0] > bbox[0]) & (pc[:, 0] < bbox[2])
        pc = pc[idx, :]
    idx = (pc[:, 1] > 312580) & (pc[:, 1] < 312590)
    ax.scatter(pc[idx, 0], pc[idx, 2], c=pc[idx, 2], s=0.02, cmap="terrain")
    ax.set_xlabel("X [m]")
    ax.set_ylabel("Z [m]")
    ax.set_title(title)
    if spatial_subset:
        ax.set_xlim(xmin=20100, xmax=20600)

plt.tight_layout()
plt.show()

%matplotlib inline # plotting a cross section fig, (ax1, ax2) = plt.subplots(nrows=2, ncols=1, figsize=(20, 8)) for pc, ax, title in zip([meas, ref], [ax1, ax2], ["Simulated", "Real-world"]): if spatial_subset: idx = (pc[:, 0] > bbox[0]) & (pc[:, 0] < bbox[2]) pc = pc[idx, :] idx = (pc[:, 1] > 312580) & (pc[:, 1] < 312590) ax.scatter(pc[idx, 0], pc[idx, 2], c=pc[idx, 2], s=0.02, cmap="terrain") ax.set_xlabel("X [m]") ax.set_ylabel("Z [m]") ax.set_title(title) if spatial_subset: ax.set_xlim(xmin=20100, xmax=20600) plt.tight_layout() plt.show()

In [20]:

fig, (ax1, ax2) = plt.subplots(nrows=1, ncols=2, figsize=(20, 10))
meas_ground = meas[meas[:, 14] == 0, :]
ref_ground = ref[las["classification"] == 2, :]

for pc, ax in zip([meas_ground, ref_ground], [ax1, ax2]):
    ax.scatter(pc[::2, 0], pc[::2, 1], c=pc[::2, 2], s=0.02, cmap="terrain")  # only plot every 2nd point so it is faster
    ax.set_xlabel("X [m]")
    ax.set_ylabel("Y [m]")
    if spatial_subset:
        ax.set_xlim(bbox[0], bbox[2])
        ax.set_ylim(bbox[1], bbox[3])
    ax.set_aspect("equal")
    
ax1.set_title("Simulated")
ax2.set_title("Real-World")
plt.tight_layout()
plt.show()

fig, (ax1, ax2) = plt.subplots(nrows=1, ncols=2, figsize=(20, 10)) meas_ground = meas[meas[:, 14] == 0, :] ref_ground = ref[las["classification"] == 2, :] for pc, ax in zip([meas_ground, ref_ground], [ax1, ax2]): ax.scatter(pc[::2, 0], pc[::2, 1], c=pc[::2, 2], s=0.02, cmap="terrain") # only plot every 2nd point so it is faster ax.set_xlabel("X [m]") ax.set_ylabel("Y [m]") if spatial_subset: ax.set_xlim(bbox[0], bbox[2]) ax.set_ylim(bbox[1], bbox[3]) ax.set_aspect("equal") ax1.set_title("Simulated") ax2.set_title("Real-World") plt.tight_layout() plt.show()

You should see clear differences in point density and ground coverage. You may also see where flight strips overlap and where they don't.

In [ ]: