4D Objects-By-Change - Analysis

This notebook explains how the extraction of 4D objects-by-change (4D-OBCs; Anders et al., 2020) is run in py4dgeo. For details about the algorithm, we refer to the articles by Anders et al. (2020; 2021).

Implemented by

Dominic Kempf (@dokempf, Heidelberg University)

Author(s) of the method Katharina Anders, Lukas Winiwarter, Bernhard Höfle (Heidelberg University)

Original publication of the method Anders, K., Winiwarter, L., Mara, H., Lindenbergh, R., Vos, S. E., & Höfle, B. (2021). Fully automatic spatiotemporal segmentation of 3D LiDAR time series for the extraction of natural surface changes. ISPRS Journal of Photogrammetry and Remote Sensing, 173, pp. 297-308. doi: 10.1016/j.isprsjprs.2021.01.015.

The concept and method are explained in this scientific talk:

import py4dgeo

import tempfile
import shutil

from pathlib import Path
from py4dgeo.util import find_file

Numba is not installed.
Use `pip install numba` to install it to enable faster computations in Vapc.
You can run py4dgeo.Vapc without Numba, but it may be slower.

The necessary data for the analysis is stored in an analysis file. There is a dedicated notebook on creating these analysis files. In this notebook, we assume that the file already exists and contains a space-time array with corresponding metadata. You can pass an absolute path to the analysis class or have py4dgeo locate files for relative paths:

# Fetch original test data
test_filename = "synthetic.zip"
original_file = find_file(test_filename)
print(f"File downloaded to: {original_file}")

# Create temporary directory
tmpdir = tempfile.mkdtemp()
zip_file = Path(tmpdir) / test_filename

# Copy original file into temporary directory
shutil.copy(find_file(test_filename), zip_file)
print(f"Working copy: {zip_file}")

Downloading file 'usage_data.zip' from 'doi:10.5281/zenodo.18378272/usage_data.zip' to '/home/docs/.cache/py4dgeo'.
Unzipping contents of '/home/docs/.cache/py4dgeo/usage_data.zip' to '/home/docs/.cache/py4dgeo/extracted'
Downloading file 'synthetic.zip' from 'doi:10.5281/zenodo.18378272/synthetic.zip' to '/home/docs/.cache/py4dgeo'.
Unzipping contents of '/home/docs/.cache/py4dgeo/synthetic.zip' to '/home/docs/.cache/py4dgeo/extracted'
Downloading file 'pbm3c2.zip' from 'doi:10.5281/zenodo.18378272/pbm3c2.zip' to '/home/docs/.cache/py4dgeo'.
Unzipping contents of '/home/docs/.cache/py4dgeo/pbm3c2.zip' to '/home/docs/.cache/py4dgeo/extracted'

File downloaded to: /home/docs/.cache/py4dgeo/synthetic.zip
Working copy: /tmp/tmpfzsph8b_/synthetic.zip

analysis = py4dgeo.SpatiotemporalAnalysis(zip_file)

If needed, the data can be retrieved from the analysis object. Note that this is a lazy operation, meaning that only explicitly requesting such data will trigger the retrieval from disk into memory:

analysis.distances

array([[0.00000000e+00, 2.58459742e-07, 5.09761587e-07, ...,
        0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
       [0.00000000e+00, 2.58459742e-07, 5.09761587e-07, ...,
        0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
       [0.00000000e+00, 2.58459742e-07, 5.09761587e-07, ...,
        0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
       ...,
       [0.00000000e+00, 2.58459742e-07, 5.09761587e-07, ...,
        0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
       [0.00000000e+00, 2.58459742e-07, 5.09761587e-07, ...,
        0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
       [0.00000000e+00, 2.58459742e-07, 5.09761587e-07, ...,
        0.00000000e+00, 0.00000000e+00, 0.00000000e+00]],
      shape=(22500, 95))

In case (part of) the analysis was run before, in this case we do not want to reuse it. We remove previous results via invalidate_results(). By default, all results are invalidated but we may adjust this to keeping the objects or seeds by setting seeds=False or objects=False.

analysis.invalidate_results(seeds=True, objects=True)

[2026-02-26 08:09:28][INFO] Removing intermediate results from the analysis file /tmp/tmpfzsph8b_/synthetic.zip

Next, we will construct an algorithm object. You can do so by simply instatiating the RegionGrowingAlgorithm:

algo = py4dgeo.RegionGrowingAlgorithm(
    neighborhood_radius=2.0,
    seed_subsampling=30,
    window_width=6,
    minperiod=3,
    height_threshold=0.05,
)

The neighborhood_radius parameter for this algorithm is of paramount importance: At each core point, a radius search with this radius is performed to determine which other core points are considered to be local neighbors during region growing. This allows us to perform region growing on a fully unstructured set of core points. The seed_subsampling parameter is used to speed up the generation of seed candidates by only investigating every n-th core point (here n=30) for changes in its timeseries. This is only done to improve the usability of this notebook - for real, full analysis you should omit that parameter (using its default of 1).

Further parameters can be set for the change point detection in the time series (cf. Anders et al., 2020). The window_width defines the width of the sliding temporal window that moves along the signal and determines the discrepancy between the first and the second half of the window (i.e. subsequent time series segments). Change point detection is performed using a C++ reimplementation of the Python library ruptures (Truong et al., 2018). The parameter minperiod controls the minimum period of a detected change feature to be considered as seed candidate. The height_threshold (default: 0.0) specifies the required magnitude of a dectected change (i.e. unsigned difference between magnitude at start epoch and peak magnitude).

Next, we execute the algorithm on our analysis object:

objects = algo.run(analysis)

[2026-02-26 08:09:28][INFO] Restoring epoch from file '/tmp/tmpaqriz_o1/corepoints.zip'
[2026-02-26 08:09:28][INFO] Starting: Find seed candidates in time series
[2026-02-26 08:09:28][INFO] Finished in 0.0939s: Find seed candidates in time series
[2026-02-26 08:09:28][INFO] Starting: Sort seed candidates by priority
[2026-02-26 08:09:28][INFO] Finished in 0.1278s: Sort seed candidates by priority
[2026-02-26 08:09:28][INFO] Starting: Performing region growing on seed candidate 1/70
[2026-02-26 08:09:28][INFO] Finished in 0.0405s: Performing region growing on seed candidate 1/70
[2026-02-26 08:09:28][INFO] Starting: Performing region growing on seed candidate 2/70
[2026-02-26 08:09:28][INFO] Finished in 0.0473s: Performing region growing on seed candidate 2/70

The algorithm returns a list of 4D objects-by-change (segments in the space time domain):

print(f"The segmentation extracted {len(objects)} 4D objects-by-change.")

The segmentation extracted 1 4D objects-by-change.

We can plot single objects interactively:

objects[0].plot()

We can check the timespan and location of the seed by accessing the information stored in the object seed properties:

# Seed coordinate
seed_coord = analysis.corepoints.cloud[objects[0].seed.index]
print(f"Coordinate of seed (XY): {seed_coord[0]:.1f}, {seed_coord[1]:.1f}")

# Seed start and end epoch
print(f"Start epoch: {objects[0].seed.start_epoch}")
print(f"End epoch: {objects[0].seed.end_epoch}")

Coordinate of seed (XY): 59.6, 74.5
Start epoch: 31
End epoch: 71

In this synthetic example there is no reference epoch set. If it were available, we could derive the timestamp of the start and end epochs via the timestamp of the reference epoch in analysis.reference_epoch.timespamp and the temporal offset using the start and end epoch id on analysis.timedeltas.

Note: Similarly to how the M3C2 class can be customized by subclassing from it, it is possible to create subclasses of RegionGrowingAlgorithm to customize some aspects of the region growing algorithm (not covered in detail in this notebook).

References

• Anders, K., Winiwarter, L., Lindenbergh, R., Williams, J. G., Vos, S. E., & Höfle, B. (2020). 4D objects-by-change: Spatiotemporal segmentation of geomorphic surface change from LiDAR time series. ISPRS Journal of Photogrammetry and Remote Sensing, 159, pp. 352-363. doi: 10.1016/j.isprsjprs.2019.11.025.

• Anders, K., Winiwarter, L., Mara, H., Lindenbergh, R., Vos, S. E., & Höfle, B. (2021). Fully automatic spatiotemporal segmentation of 3D LiDAR time series for the extraction of natural surface changes. ISPRS Journal of Photogrammetry and Remote Sensing, 173, pp. 297-308. doi: 10.1016/j.isprsjprs.2021.01.015.

• Truong, C., Oudre, L., Vayatis, N. (2018): ruptures: Change point detection in python. arXiv preprint: abs/1801.00826.