4D Objects By Change - Creating an analysis

This notebook explains the data preparation step for the extraction of 4D Objects-By-Change (4D-OBCs; Anders et al., 2020). For details about the method, we refer to the articles by Anders et al. (2020; 2021).

import py4dgeo
import numpy as np

The main data structure for the 4D-OBC algorithm is a SpatiotemporalAnalysis object. It is always backed by an archive file, which we specify when instantiating the analysis object. If the given filename already exists, we open it, otherwise an empty archive is created:

analysis = py4dgeo.SpatiotemporalAnalysis("test.zip", force=True)

[2024-05-14 13:00:26][INFO] Creating analysis file test.zip

The analysis object stores all information necessary to (re-)run and extend the analysis. This includes, e.g., the reference epoch, the core points, and the space-time array of change values. However, it does not store all epochs that were used in building up the analysis.

In the following we show how the required components are added to the analysis using the example of the two epochs used for demonstration of the M3C2. To be usable during spatiotemporal segmentation, a timestamp is added to these epochs, respectively:

reference_epoch, epoch = py4dgeo.read_from_xyz(
    "plane_horizontal_t1.xyz", "plane_horizontal_t2.xyz"
)
reference_epoch.timestamp = "March 9th 2022, 16:00"
epoch.timestamp = "March 10th 2022, 16:00"

[2024-05-14 13:00:26][INFO] Reading point cloud from file '/home/docs/.cache/py4dgeo/./plane_horizontal_t1.xyz'
[2024-05-14 13:00:26][INFO] Reading point cloud from file '/home/docs/.cache/py4dgeo/./plane_horizontal_t2.xyz'

The reference epoch needs to be added as such to the analysis:

analysis.reference_epoch = reference_epoch

[2024-05-14 13:00:26][INFO] Building KDTree structure with leaf parameter 10
[2024-05-14 13:00:26][INFO] Saving epoch to file '/tmp/tmpr45mmoci/reference_epoch.zip'
[2024-05-14 13:00:26][INFO] Saving a file without normals.

We also need to define the set of core points in use. In this example, we use all points in the reference epoch as core points:

analysis.corepoints = reference_epoch

[2024-05-14 13:00:26][INFO] Saving epoch to file '/tmp/tmpopr_3hw6/corepoints.zip'
[2024-05-14 13:00:26][INFO] Saving a file without normals.

Next, we want to add epochs to the spatiotemporal analysis. For this, we need to define a method that calculates the point cloud distances, for example using the M3C2 algorithm (Lague et al., 2013). For details of the algorithm setup, you can have a look at the M3C2 notebook.

analysis.m3c2 = py4dgeo.M3C2(cyl_radii=[2.0], normal_radii=[2.0])

Having done all of this, we can add an epoch to the analysis:

analysis.add_epochs(epoch)

[2024-05-14 13:00:26][INFO] Removing intermediate results from the analysis file test.zip
[2024-05-14 13:00:26][INFO] Starting: Adding epoch 1/1 to analysis object
[2024-05-14 13:00:26][INFO] Building KDTree structure with leaf parameter 10
[2024-05-14 13:00:26][INFO] Finished in 0.0054s: Adding epoch 1/1 to analysis object
[2024-05-14 13:00:26][INFO] Starting: Rearranging space-time array in memory
[2024-05-14 13:00:26][INFO] Finished in 0.0034s: Rearranging space-time array in memory
[2024-05-14 13:00:26][INFO] Starting: Updating disk-based analysis archive with new epochs
[2024-05-14 13:00:26][INFO] Finished in 0.0009s: Updating disk-based analysis archive with new epochs

Note that the add_epochs can be passed several epochs at the same time to save some costly memory operations: add_epochs(*list_of_epochs). Also, add_epochs can be called again on existing analysis objects allowing you to update your analysis after additional data acquisition at the same site.

Having done this, we can inspect some of the data contained in the analysis:

np.set_printoptions(threshold=5)

print(f"Space-time distance array: {analysis.distances}")

Space-time distance array: [[-0.10269298]
 [-0.10275566]
 [-0.10072999]
 ...
 [-0.10194934]
 [-0.1026113 ]
 [-0.1022107 ]]

print(f"Uncertainty array of M3C2 calculation: {analysis.uncertainties}")

Uncertainty array of M3C2 calculation: [[(0.00565384, 0.00417673, 5, 0.00491526, 5)]
 [(0.0019427 , 0.00149693, 7, 0.00215317, 7)]
 [(0.00363399, 0.00330837, 8, 0.00406884, 8)]
 ...
 [(0.00299548, 0.00228501, 8, 0.00343241, 7)]
 [(0.00444868, 0.00394557, 8, 0.00473712, 7)]
 [(0.0067821 , 0.00662842, 6, 0.00482219, 5)]]

print(f"Timestamp deltas of analysis: {analysis.timedeltas}")

Timestamp deltas of analysis: [datetime.timedelta(days=1)]

Note that the add_epochs can be passed several epochs at the same time to save some costly memory operations. Also, add_epochs can be called again on existing analysis objects allowing you to update your analysis after additional data acquisition at the same site.

Sometimes, you will want to apply preprocessing in the form of smoothing to the calculated distance array. You can do so by setting the analysis object’s smoothed_distances property. When this was set, the 4D-OBC algorithm will operate on the smoothed data instead of the raw data stored in the analysis object. py4dgeo provides a smoothing implementation using a sliding window approach:

analysis.smoothed_distances = py4dgeo.temporal_averaging(
    analysis.distances, smoothing_window=24
)

[2024-05-14 13:00:26][INFO] Starting: Smoothing temporal data
[2024-05-14 13:00:26][INFO] Finished in 0.0002s: Smoothing temporal data

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.

• Lague, D., Brodu, N., & Leroux, J. (2013). Accurate 3D comparison of complex topography with terrestrial laser scanner: Application to the Rangitikei canyon (N-Z). ISPRS Journal of Photogrammetry and Remote Sensing, 82, pp. 10-26. doi: 10.1016/j.isprsjprs.2013.04.009.