OccPy tutorial notebook for Terrestrial Laser Scanning (TLS) data - Individual OccPy Run per LAZ file¶

This is basically the same notebook as TLS_notebook.ipynb, however, each input laz file is treated individually, instead of passing the input folder to OccPy and let OccPy treat the input. Treating laz files individually can be beneficial from a performance perspective and it can further introduce an increased amount of flexibility for linking scan position to the specific laz point cloud. The process though, is roughly the same.

First, we download the test data and configure paths

In [1]:

import json
import os
import shutil
from pathlib import Path
import numpy as np
import pooch
import glob

import json import os import shutil from pathlib import Path import numpy as np import pooch import glob

In [2]:

# Download test data once to cache and mirror it into the repository for stable relative paths
repo_root = Path('..', '..').resolve()
data_notebooks = repo_root / 'data_notebooks'
tls_data_dir = data_notebooks / 'TLS_demo'
data_notebooks.mkdir(parents=True, exist_ok=True)

p = pooch.create(
    path=pooch.os_cache('occpy_test_data'),
    base_url='https://zenodo.org/records/17750604/files/',
    registry={'TLS_demo.zip': 'md5:ff15ef6b1b6e655e33c50020089dad56'},
)
p.fetch('TLS_demo.zip', processor=pooch.Unzip(members=['TLS_demo']), progressbar=True)

cache_data_dir = Path(p.path) / 'TLS_demo.zip.unzip' / 'TLS_demo'
if not tls_data_dir.exists():
    shutil.copytree(cache_data_dir, tls_data_dir)

print(f'Repository notebook data folder: {tls_data_dir}')
config_file = str(repo_root / 'config' / 'settings_TLS_tutorial_out_indivLAZ.JSON')

# Download test data once to cache and mirror it into the repository for stable relative paths repo_root = Path('..', '..').resolve() data_notebooks = repo_root / 'data_notebooks' tls_data_dir = data_notebooks / 'TLS_demo' data_notebooks.mkdir(parents=True, exist_ok=True) p = pooch.create( path=pooch.os_cache('occpy_test_data'), base_url='https://zenodo.org/records/17750604/files/', registry={'TLS_demo.zip': 'md5:ff15ef6b1b6e655e33c50020089dad56'}, ) p.fetch('TLS_demo.zip', processor=pooch.Unzip(members=['TLS_demo']), progressbar=True) cache_data_dir = Path(p.path) / 'TLS_demo.zip.unzip' / 'TLS_demo' if not tls_data_dir.exists(): shutil.copytree(cache_data_dir, tls_data_dir) print(f'Repository notebook data folder: {tls_data_dir}') config_file = str(repo_root / 'config' / 'settings_TLS_tutorial_out_indivLAZ.JSON')

Repository notebook data folder: C:\Users\Kueken\dev\occpy_dev\OccPy\data_notebooks\TLS_demo

OccPy configurations are saved in JSON files.

There are a few minimal requirements for a run:

• laz_in: input laz file or directory containing laz files. If given multiple las files, will assume multi-position TLS. If given single las file, will assume single-position TLS unless is_mobile is set.
• plot_dim: grid for occlusion mapping: [minX, minY, minZ, maxX, maxY, maxZ]
• vox_dim: voxel size - currently this applies to x, y, and z dimension. Non-cubic voxels are on the list of features to be added

Optional arguments:

• is_mobile: whether the acquisition is mobile (MLS/ULS) or static (TLS) (default: False)
• single_return: whether the data is single return or multi return data (default: False)
• out_dir: output directory (default: ./output)
• verbose: set logging level (default: False)
• debug: set logging level (default: False)
• lower_threshold: lower threshold above ground to exclude from occlusion mapping in voxels (default: 0)
• points_per_iter: number of points read in from laz file in each iteration (default: 10000000)
• delimiter: csv delimiter for scan position file (default: ",")
• root_folder: if given, will assume other paths are relative to this root folder and will prepend it to the paths (default: None)
• str_idxs_ScanPosID: string indices of where the scan position identifier is written in the laz file name. If not given, will use file name as ID (without extension) (default: None)
• output_voxels: whether to export .ply voxel grids (large files, slow) (default: False)

Additionally, we use the config file to save paths to other files, such as the positions file, DSM and DTM.

In [3]:

# Load config and resolve all relative paths once for downstream plotting/normalization calls
with open(config_file, 'r') as file:
    config = json.load(file)

config

# Load config and resolve all relative paths once for downstream plotting/normalization calls with open(config_file, 'r') as file: config = json.load(file) config

Out[3]:

{'root_folder': '../..',
 'laz_in': 'data_notebooks/TLS_demo/LAZ',
 'tif_in': {'DTM': 'data_notebooks/TLS_demo/Grids/Ramerenwald_DTM_20250305.tif',
  'DSM': 'data_notebooks/TLS_demo/Grids/Ramerenwald_DSM_20250305.tif'},
 'out_dir': 'output/TLS_individual',
 'vox_dim': 0.1,
 'lower_threshold': 1,
 'points_per_iter': 1000000,
 'plot_dim': [2676515, 1246063, 545, 2676525, 1246113, 590],
 'ScanPos': 'data_notebooks/TLS_demo/ScanPos/ScanPositions.txt',
 'output_voxels': False,
 'single_return': False,
 'is_mobile': False,
 'str_idxs_ScanPosID': [7, 10]}

Running OccPy¶

In [4]:

from occpy.OccPy import OccPy
from occpy.util import read_sensorpos_file

from occpy.OccPy import OccPy from occpy.util import read_sensorpos_file

RIEGL RDBlib is not available
RIEGL RiVlib is not available
Jupyter environment detected. Enabling Open3D WebVisualizer.
[Open3D INFO] WebRTC GUI backend enabled.
[Open3D INFO] WebRTCWindowSystem: HTTP handshake server disabled.

Let's load the sensor position information. The function is the same as the occpy class function, but returns the positions, so we can link each position to its scan manually below.

In [5]:

sens_pos = read_sensorpos_file(os.path.join(config['root_folder'], config['ScanPos']),  # Path to the sensor position file
                               delimiter=",",                                           # Delimiter used in the file (e.g., comma, tab)
                               hdr_scanpos_id="ID",                                     # Header name for the scan position ID column
                               hdr_x='X',                                               # Header name for the X coordinate column
                               hdr_y='Y',                                               # Header name for the Y coordinate column
                               hdr_z='Z',                                               # Header name for the Z coordinate column
                               sens_pos_id_offset=0)                                    # Offset to be added to the scan position IDs (if needed, i.e. if the IDs in the ScanPos file is not matching the LAZ file names)

sens_pos = read_sensorpos_file(os.path.join(config['root_folder'], config['ScanPos']), # Path to the sensor position file delimiter=",", # Delimiter used in the file (e.g., comma, tab) hdr_scanpos_id="ID", # Header name for the scan position ID column hdr_x='X', # Header name for the X coordinate column hdr_y='Y', # Header name for the Y coordinate column hdr_z='Z', # Header name for the Z coordinate column sens_pos_id_offset=0) # Offset to be added to the scan position IDs (if needed, i.e. if the IDs in the ScanPos file is not matching the LAZ file names)

Now we iterate over all laz files, matching each one to its sensor position and running raytracing.

In [6]:

# get all laz files in the input folder
laz_files = sorted(glob.glob(os.path.join(config["root_folder"], config['laz_in'], '*.laz')))

for laz_in in laz_files:
    scan_name = os.path.basename(laz_in)
    print(f'Processing {scan_name}...')

    # Build a per-scan config dict (keeps canonical config untouched).
    scan_cfg = dict(config)
    scan_cfg['laz_in'] = os.path.abspath(laz_in)
    scan_cfg['out_dir'] = os.path.join(config["root_folder"], config['out_dir'], scan_name[:-4])
    scan_cfg["root_folder"] = None  # Avoid confusion with per-scan absolute paths

    occpy = OccPy(config=scan_cfg)

    # Match current LAZ to sensor position entry.
    scan_id = int(scan_name[7:10])  # Adapt if LAZ naming convention differs.
    scanpos_X = sens_pos.loc[sens_pos['ScanPos'] == scan_id, 'sensor_x'].values[0]
    scanpos_Y = sens_pos.loc[sens_pos['ScanPos'] == scan_id, 'sensor_y'].values[0]
    scanpos_Z = sens_pos.loc[sens_pos['ScanPos'] == scan_id, 'sensor_z'].values[0]

    occpy.define_sensor_pos_singlePos(scan_pos_id=scan_id,
                                      x=scanpos_X,
                                      y=scanpos_Y,
                                      z=scanpos_Z)

    occpy.do_raytracing()

# get all laz files in the input folder laz_files = sorted(glob.glob(os.path.join(config["root_folder"], config['laz_in'], '*.laz'))) for laz_in in laz_files: scan_name = os.path.basename(laz_in) print(f'Processing {scan_name}...') # Build a per-scan config dict (keeps canonical config untouched). scan_cfg = dict(config) scan_cfg['laz_in'] = os.path.abspath(laz_in) scan_cfg['out_dir'] = os.path.join(config["root_folder"], config['out_dir'], scan_name[:-4]) scan_cfg["root_folder"] = None # Avoid confusion with per-scan absolute paths occpy = OccPy(config=scan_cfg) # Match current LAZ to sensor position entry. scan_id = int(scan_name[7:10]) # Adapt if LAZ naming convention differs. scanpos_X = sens_pos.loc[sens_pos['ScanPos'] == scan_id, 'sensor_x'].values[0] scanpos_Y = sens_pos.loc[sens_pos['ScanPos'] == scan_id, 'sensor_y'].values[0] scanpos_Z = sens_pos.loc[sens_pos['ScanPos'] == scan_id, 'sensor_z'].values[0] occpy.define_sensor_pos_singlePos(scan_pos_id=scan_id, x=scanpos_X, y=scanpos_Y, z=scanpos_Z) occpy.do_raytracing()

Processing ScanPos033_sort_filtered.laz...
INFO: optional arguments: ['out_dir', 'output_voxels', 'verbose', 'debug', 'lower_threshold', 'points_per_iter', 'delimiter', 'root_folder', 'single_return', 'str_idxs_ScanPosID']

Tracing Pulses...: 100%|██████████| 4623932/4623932 [00:42<00:00, 109437.65pulses/s]

b'#### 3062939 Pulses were traversed of possible 3062939 Pulses\n#### 1527078 Returns have been registered by the algorithm \n#### 2468919 Returns were found outside the voxel grid \n#### 603 Returns were missed during the traversal!\n#### 0 Pulses did not intersect voxel grid!\n'
Processing ScanPos035_sort_filtered.laz...
INFO: optional arguments: ['out_dir', 'output_voxels', 'verbose', 'debug', 'lower_threshold', 'points_per_iter', 'delimiter', 'root_folder', 'single_return', 'str_idxs_ScanPosID']

Tracing Pulses...: 100%|██████████| 4713530/4713530 [00:27<00:00, 172596.78pulses/s]

b'#### 2225942 Pulses were traversed of possible 3143921 Pulses\n#### 527592 Returns have been registered by the algorithm \n#### 2365359 Returns were found outside the voxel grid \n#### 178 Returns were missed during the traversal!\n#### 917979 Pulses did not intersect voxel grid!\n'
Processing ScanPos045_sort_filtered.laz...
INFO: optional arguments: ['out_dir', 'output_voxels', 'verbose', 'debug', 'lower_threshold', 'points_per_iter', 'delimiter', 'root_folder', 'single_return', 'str_idxs_ScanPosID']

Tracing Pulses...: 100%|██████████| 5394449/5394449 [00:28<00:00, 188823.13pulses/s]

b'#### 3030796 Pulses were traversed of possible 3065216 Pulses\n#### 1512084 Returns have been registered by the algorithm \n#### 2952204 Returns were found outside the voxel grid \n#### 1119 Returns were missed during the traversal!\n#### 34420 Pulses did not intersect voxel grid!\n'
Processing ScanPos047_sort_filtered.laz...
INFO: optional arguments: ['out_dir', 'output_voxels', 'verbose', 'debug', 'lower_threshold', 'points_per_iter', 'delimiter', 'root_folder', 'single_return', 'str_idxs_ScanPosID']

Tracing Pulses...: 100%|██████████| 4829061/4829061 [00:30<00:00, 160188.51pulses/s]

b'#### 2987270 Pulses were traversed of possible 2987270 Pulses\n#### 1461202 Returns have been registered by the algorithm \n#### 2600354 Returns were found outside the voxel grid \n#### 939 Returns were missed during the traversal!\n#### 0 Pulses did not intersect voxel grid!\n'
Processing ScanPos061_sort_filtered.laz...
INFO: optional arguments: ['out_dir', 'output_voxels', 'verbose', 'debug', 'lower_threshold', 'points_per_iter', 'delimiter', 'root_folder', 'single_return', 'str_idxs_ScanPosID']

Tracing Pulses...: 100%|██████████| 4981696/4981696 [00:27<00:00, 179090.79pulses/s]

b'#### 2913516 Pulses were traversed of possible 2913516 Pulses\n#### 2440049 Returns have been registered by the algorithm \n#### 1654218 Returns were found outside the voxel grid \n#### 2452 Returns were missed during the traversal!\n#### 0 Pulses did not intersect voxel grid!\n'

Now we have an output for each individual scan position. This would allow us to analyze how each scan position contributes to the overall oclusion.
However, here we just combine the individual outputs, which should give us the same output as in the joint TLS notebook.

In [7]:

# Combine outputs
comb_dir = os.path.join(config["root_folder"], config['out_dir'], "merged_voxgrids")


# create a new output directory for the combined outputs
os.makedirs(comb_dir, exist_ok=True)

fCont = glob.glob(os.path.join(config["root_folder"], config['out_dir'], "ScanPos*"))  

idx = 0
for scan in fCont:
    print(f"Processing {os.path.basename(scan)}")
    if idx==0:
        Nhit = np.load(os.path.join(scan, "Nhit.npy"))
        Nmiss = np.load(os.path.join(scan, "Nmiss.npy"))
        Nocc = np.load(os.path.join(scan, "Nocc.npy"))
    else:
        Nhit = Nhit + np.load(os.path.join(scan, "Nhit.npy"))
        Nmiss = Nmiss + np.load(os.path.join(scan, "Nmiss.npy"))
        Nocc = Nocc + np.load(os.path.join(scan,"Nocc.npy"))
    idx = idx+1

# Calculate Classification according to Bienert et al. and save outputs
Classification = np.zeros(Nhit.shape, dtype=int)

Classification[np.logical_and.reduce((Nhit > 0, Nmiss >= 0, Nocc >= 0))] = 1 # voxels that were observed
Classification[np.logical_and.reduce((Nhit == 0, Nmiss > 0, Nocc >= 0))] = 2  # voxels that are empty
Classification[np.logical_and.reduce((Nhit == 0, Nmiss == 0, Nocc > 0))] = 3  # voxels that are hidden (occluded)
Classification[np.logical_and.reduce((Nhit == 0, Nmiss == 0, Nocc == 0))] = 4  # voxels that were not observed #

print("Saving combined grids")
np.save(os.path.join(comb_dir, "Nhit.npy"), Nhit)
np.save(os.path.join(comb_dir, "Nmiss.npy"), Nmiss)
np.save(os.path.join(comb_dir, "Nocc.npy"), Nocc)
np.save(os.path.join(comb_dir,"Classification.npy"), Classification)

# Combine outputs comb_dir = os.path.join(config["root_folder"], config['out_dir'], "merged_voxgrids") # create a new output directory for the combined outputs os.makedirs(comb_dir, exist_ok=True) fCont = glob.glob(os.path.join(config["root_folder"], config['out_dir'], "ScanPos*")) idx = 0 for scan in fCont: print(f"Processing {os.path.basename(scan)}") if idx==0: Nhit = np.load(os.path.join(scan, "Nhit.npy")) Nmiss = np.load(os.path.join(scan, "Nmiss.npy")) Nocc = np.load(os.path.join(scan, "Nocc.npy")) else: Nhit = Nhit + np.load(os.path.join(scan, "Nhit.npy")) Nmiss = Nmiss + np.load(os.path.join(scan, "Nmiss.npy")) Nocc = Nocc + np.load(os.path.join(scan,"Nocc.npy")) idx = idx+1 # Calculate Classification according to Bienert et al. and save outputs Classification = np.zeros(Nhit.shape, dtype=int) Classification[np.logical_and.reduce((Nhit > 0, Nmiss >= 0, Nocc >= 0))] = 1 # voxels that were observed Classification[np.logical_and.reduce((Nhit == 0, Nmiss > 0, Nocc >= 0))] = 2 # voxels that are empty Classification[np.logical_and.reduce((Nhit == 0, Nmiss == 0, Nocc > 0))] = 3 # voxels that are hidden (occluded) Classification[np.logical_and.reduce((Nhit == 0, Nmiss == 0, Nocc == 0))] = 4 # voxels that were not observed # print("Saving combined grids") np.save(os.path.join(comb_dir, "Nhit.npy"), Nhit) np.save(os.path.join(comb_dir, "Nmiss.npy"), Nmiss) np.save(os.path.join(comb_dir, "Nocc.npy"), Nocc) np.save(os.path.join(comb_dir,"Classification.npy"), Classification)

Processing ScanPos033_sort_filtered
Processing ScanPos035_sort_filtered
Processing ScanPos045_sort_filtered
Processing ScanPos047_sort_filtered
Processing ScanPos061_sort_filtered
Saving combined grids

Now we can height normalize the merged output grid

In [8]:

from occpy.util import normalize_occlusion_output

from occpy.util import normalize_occlusion_output

In [9]:

# normalize outputs
Nhit_norm, Nmiss_norm, Nocc_norm, Classification_norm, chm = normalize_occlusion_output(input_folder=comb_dir,
                           PlotDim=config['plot_dim'],
                           vox_dim=config['vox_dim'],
                           dtm_file=os.path.join(config["root_folder"], config['tif_in']['DTM']),
                           dsm_file=os.path.join(config["root_folder"], config['tif_in']['DSM']),
                           lower_threshold=1,
                           output_voxels=False)

# normalize outputs Nhit_norm, Nmiss_norm, Nocc_norm, Classification_norm, chm = normalize_occlusion_output(input_folder=comb_dir, PlotDim=config['plot_dim'], vox_dim=config['vox_dim'], dtm_file=os.path.join(config["root_folder"], config['tif_in']['DTM']), dsm_file=os.path.join(config["root_folder"], config['tif_in']['DSM']), lower_threshold=1, output_voxels=False)

Saving normalized output files into directory as .npy...

And finally we can visualize the occlusion grid.

In [10]:

from occpy.visualization import get_Occl_TransectFigure_BinaryOcclusion, get_Occlusion_ProfileFigure

from occpy.visualization import get_Occl_TransectFigure_BinaryOcclusion, get_Occlusion_ProfileFigure

In [11]:

# define figure properties
fig_prop = dict(fig_size=(3.5, 3.2),
                label_size=8,
                label_size_ticks=6,
                label_size_tiny=5,
                out_format='png',)

# get transect figure
%matplotlib inline
get_Occl_TransectFigure_BinaryOcclusion(Nhit_norm, Classification_norm, plot_dim=config['plot_dim'], vox_dim=config['vox_dim'],
                                        out_dir=os.path.join(config["root_folder"], config['out_dir']), axis=0, start_ind=0, end_ind=100, chm=chm, vertBuffer=10, fig_prop=fig_prop, show_plots=True)
# define figure properties for occlusion profile
fig_prop = dict(fig_size=(1.75, 3.2),
                label_size=8,
                label_size_ticks=6,
                label_size_tiny=5,
                out_format='png', )
get_Occlusion_ProfileFigure(Classification_norm, plot_dim=config['plot_dim'], vox_dim=config['vox_dim'], out_dir=os.path.join(config["root_folder"], config['out_dir']), low_thresh=0, vertBuffer=10, max_percentage=100, fig_prop=fig_prop, show_plots=True)

# define figure properties fig_prop = dict(fig_size=(3.5, 3.2), label_size=8, label_size_ticks=6, label_size_tiny=5, out_format='png',) # get transect figure %matplotlib inline get_Occl_TransectFigure_BinaryOcclusion(Nhit_norm, Classification_norm, plot_dim=config['plot_dim'], vox_dim=config['vox_dim'], out_dir=os.path.join(config["root_folder"], config['out_dir']), axis=0, start_ind=0, end_ind=100, chm=chm, vertBuffer=10, fig_prop=fig_prop, show_plots=True) # define figure properties for occlusion profile fig_prop = dict(fig_size=(1.75, 3.2), label_size=8, label_size_ticks=6, label_size_tiny=5, out_format='png', ) get_Occlusion_ProfileFigure(Classification_norm, plot_dim=config['plot_dim'], vox_dim=config['vox_dim'], out_dir=os.path.join(config["root_folder"], config['out_dir']), low_thresh=0, vertBuffer=10, max_percentage=100, fig_prop=fig_prop, show_plots=True)