In this notebook we will perform time series-based surface change analysis on a time series of permanent TLS point clouds of the sandy beach at Kijkduin for a timespan of around 6 months (Vos et al., 2022). An introduction to the case study and dataset can be found here.
The objective is to assess surface dynamics with 4D objects-by-change (following Anders et al., 2021). Look into the related article for comparison of possible surface dynamics of the use case and help for deciding on suitable parameters, etc.
The workflow is introduced throughout this notebook. You can also make use of the software documentation!
This task is solved using Python with the py4dgeo library.
You can use CloudCompare or GIS Software (e.g. QGIS) to check the data and visualize your results.
The dataset is a subsampled version of the original time series, using 12-hourly epochs of point clouds and spatial subsampling to 50 cm. The dataset can be downloaded (module3.zip) from the E-learning course E-TRAINEE. In the data directory kijkduin, you find the prepared input point clouds and a core points point cloud, which is manually cleaned from noise.
E-TRAINEE is an e-learning course on Time Series Analysis in Remote Sensing for Understanding Human-Environment Interactions. This course has been developed by research groups from four partner universities – Charles University, Heidelberg University, University of Innsbruck, and University of Warsaw.
Prepare the analysis by compiling the list of files (epochs) and read the timestamps from the file names (format YYMMDD_hhmmss) into datetime objects. Use the point cloud files and timestamps to create a py4dgeo SpatiotemporalAnalysis object. For this you need to instantiate the M3C2 algorithm. You can use the point cloud file 170115_150816_aoi_50cm.laz as core points. Explore the point cloud properties in CloudCompare:
Hint: In this flat topography and predominant provess of sand deposition and erosion, it can be suitable to orient the normals purely vertically. In this case, they do not need to be computed, and you can customize the py4dgeo algorithm accordingly.
Use the first point cloud in the time series (list of files) as reference epoch. You can assume a registration error of 1.9 cm for the M3C2 distance calculation (cf. Vos et al., 2022).
Explore the spatiotemporal change information by visualizing the changes at a selected epoch and visualizing the time series at a selected location.
First, we start by setting up the Python environment and data:
# Import required modules
import py4dgeo
import os
import numpy as np
from datetime import datetime
import pooch
# Download data from zenodo and set path to point cloud folder
p = pooch.Pooch(base_url="doi:10.5281/zenodo.10003574/", path=pooch.os_cache("py4dgeo"))
p.load_registry_from_doi()
p.fetch("module3.zip", processor=pooch.Unzip(members=["module3"]))
pc_dir = p.path / "module3.zip.unzip" / "module3" / "kijkduin" / "pointclouds"
print(pc_dir)
# List of point clouds (time series)
pc_list = os.listdir(pc_dir)
pc_list[:5] # print the first elements
In the list of point cloud files you can see that we have one laz file per epoch available. The file name contains the timestamp of the epoch, respectively, in format YYMMDD_hhmmss. To use this information for our analysis, we read the timestamp information from the file names into datetime objects.
# Read the timestamps from file names
timestamps = []
for f in pc_list:
if not f.endswith(".laz"):
continue
# Get the timestamp from the file name
timestamp_str = "_".join(f.split(".")[0].split("_")[1:]) # yields YYMMDD_hhmmss
# Convert string to datetime object
timestamp = datetime.strptime(timestamp_str, "%y%m%d_%H%M%S")
timestamps.append(timestamp)
timestamps[:5]
Now we use the point cloud files and timestamp information to create a SpatiotemporalAnalysis object.
data_path = p.path / "module3.zip.unzip" / "module3" / "kijkduin" / "kijkduin.zip"
analysis = py4dgeo.SpatiotemporalAnalysis(data_path, force=True)
As reference epoch, we use the first epoch in our time series:
# Specify the reference epoch
reference_epoch_file = os.path.join(pc_dir, pc_list[0])
# Read the reference epoch and set the timestamp
reference_epoch = py4dgeo.read_from_las(reference_epoch_file)
reference_epoch.timestamp = timestamps[0]
# Set the reference epoch in the spatiotemporal analysis object
analysis.reference_epoch = reference_epoch
For epochs to be added, we now configure the M3C2 algorithm to derive the change values. We would like to set the normals purely vertically, so we define a customized computation of cylinder directions:
# Inherit from the M3C2 algorithm class to define a custom direction algorithm
class M3C2_Vertical(py4dgeo.M3C2):
def directions(self):
return np.array([0, 0, 1]) # vertical vector orientation
# Specify corepoints, here all points of the reference epoch
analysis.corepoints = reference_epoch.cloud[::]
# Specify M3C2 parameters for our custom algorithm class
analysis.m3c2 = M3C2_Vertical(
cyl_radii=(1.0,), max_distance=10.0, registration_error=0.019
)
Now we add all the other epochs with their timestamps:
# Create a list to collect epoch objects
epochs = []
for e, pc_file in enumerate(pc_list[1:]):
epoch_file = pc_dir / pc_file
epoch = py4dgeo.read_from_las(epoch_file)
epoch.timestamp = timestamps[e]
epochs.append(epoch)
# Add epoch objects to the spatiotemporal analysis object
analysis.add_epochs(*epochs)
# Print the spatiotemporal analysis data for 3 corepoints and 5 epochs, respectively
print(f"Space-time distance array:\n{analysis.distances[:3,:5]}")
print(
f"Uncertainties of M3C2 distance calculation:\n{analysis.uncertainties['lodetection'][:3, :5]}"
)
print(f"Timestamp deltas of analysis:\n{analysis.timedeltas[:5]}")
We visualize the changes in the scene for a selected epoch, together with the time series of surface changes at a selected location. The location here was selected separately in CloudCompare (as the corepoint id).
cp_idx_sel = 15162 # selected core point index
epoch_idx_sel = 28 # selected epoch index
# Import plotting module
import matplotlib.pyplot as plt
# Allow interactive rotation in notebook
%matplotlib inline
# Create the figure
fig = plt.figure(figsize=(15, 5))
ax1 = fig.add_subplot(1, 2, 1)
ax2 = fig.add_subplot(1, 2, 2)
# Get the corepoints
corepoints = analysis.corepoints.cloud
# Get change values of last epoch for all corepoints
distances = analysis.distances
distances_epoch = [d[epoch_idx_sel] for d in distances]
# Get the time series of changes at a specific core point locations
coord_sel = analysis.corepoints.cloud[cp_idx_sel]
timeseries_sel = distances[cp_idx_sel]
# Get the list of timestamps from the reference epoch timestamp and timedeltas
timestamps = [t + analysis.reference_epoch.timestamp for t in analysis.timedeltas]
# Plot the scene
d = ax1.scatter(
corepoints[:, 0],
corepoints[:, 1],
c=distances_epoch[:],
cmap="seismic_r",
vmin=-1.5,
vmax=1.5,
s=1,
zorder=1,
)
plt.colorbar(d, format=("%.2f"), label="Distance [m]", ax=ax1, pad=0.15)
# Add the location of the selected coordinate
ax1.scatter(
coord_sel[0],
coord_sel[1],
facecolor="yellow",
edgecolor="black",
s=100,
zorder=2,
label="Selected location",
marker="*",
)
ax1.legend()
# Configure the plot layout
ax1.set_xlabel("X [m]")
ax1.set_ylabel("Y [m]")
ax1.set_aspect("equal")
ax1.set_title(
"Changes at %s"
% (analysis.reference_epoch.timestamp + analysis.timedeltas[epoch_idx_sel])
)
# Plot the time series
ax2.scatter(timestamps, timeseries_sel, s=5, color="black", label="Raw")
ax2.plot(timestamps, timeseries_sel, color="blue", label="Smoothed")
ax2.set_xlabel("Date")
# Add the epoch of the plotted scene
ax2.scatter(
timestamps[epoch_idx_sel],
timeseries_sel[epoch_idx_sel],
facecolor="yellow",
marker="*",
edgecolor="black",
s=100,
color="red",
label="Selected epoch",
)
ax2.legend()
# Format the date labels
import matplotlib.dates as mdates
dtFmt = mdates.DateFormatter("%b-%d")
plt.gca().xaxis.set_major_formatter(dtFmt)
# Configure the plot layout
ax2.set_ylabel("Distance [m]")
ax2.grid()
ax2.set_ylim(-0.3, 1.6)
ax2.set_title("Time series at selected location")
plt.tight_layout()
plt.show()
The map of changes in the scene shows us linear structures of sand deposition near the coast which can be interpreted as sand bars (with knowledge about coastal processes). This is confirmed by the surface behavior over time, expressed in the time series plot. However, the time series is quite noisy especially in this part of the beach, which is regularly covered by water during high tides (leading to missing data) and also varies strongly in surface moisture (influencing the LiDAR range measurement and causing noise). We therefore continue with temporal smoothing of the time series.
We are dealing with a temporal subset of the original hourly time series. The effect of temporal measurement variability may therefore be less pronounced (compared to the assessment in, e.g., Anders et al., 2019). Nonetheless, we may apply temporal smoothing to reduce the influence of noise on your change analysis using a rolling median averaging of one week. This will also fill possible gaps in your data, e.g., lower ranges during poor atmospheric conditions or no data due to water coverage during tides on the lower beach part.
We can visualize the raw and smoothed change values in the time series of your selected location.
We apply a rolling median with a defined temporal window of 14 epochs (corresponding to one week of 12-hourly point clouds) using the temporal_averaging() function in py4dgeo.
analysis.smoothed_distances = py4dgeo.temporal_averaging(
analysis.distances, smoothing_window=14
)
Now we can compare the raw and smoothed time series at our selected location:
# Create the figure
fig, ax = plt.subplots(1, 1, figsize=(7, 5))
# Plot the raw time series
ax.scatter(timestamps, timeseries_sel, color="blue", label="raw", s=5)
# Plot the smoothed time series
timeseries_sel_smooth = analysis.smoothed_distances[cp_idx_sel]
ax.plot(timestamps, timeseries_sel_smooth, color="red", label="smooth")
# Format the date labels
import matplotlib.dates as mdates
dtFmt = mdates.DateFormatter("%b-%d")
plt.gca().xaxis.set_major_formatter(dtFmt)
# Add plot elements
ax.set_xlabel("Date")
ax.set_ylabel("Distance [m]")
ax.grid()
ax.set_ylim(-0.25, 1.25)
plt.tight_layout()
plt.show()
From the smoothed time series at the selected location, we can now much better interpret the surface behavior. In fact, we can distinctly observe that there were two consecutive occurrences of temporary deposition of several weeks. These represent two phases where sand bars are present. They can be extracted as individual objects by the 4D objects-by-change method. Before, we continue with time series clustering and the assessment of overall change patterns in the following.
Now we use the 4D objects-by-change (4D-OBC) method to identify individual surface activities in the beach scene. The objective is to extract temporary occurrences of accumulation or erosion, as occurs when a sand bar is formed during some timespan, or when sand is deposited by anthropogenic works. This type of surface activity is implemented with the original seed detection in py4dgeo, so you do not need to customize the algorithm. Decide for suitable parameters following Anders et al. (2021) - but bear in mind that we are using a different temporal resolution, so you may need to adjust the temporal change detection.
We perform the extraction for selected seed locations, e.g. considering interesting clusters of change patterns identified in the previous step. In principle, the spatiotemporal segmentation can also be applied to the full dataset (all time series at all core point locations are used as potential seed candidates), but long processing time needs to be expected.
In this solution, we will use the selected location from above to extract the sand bars as 4D object-by-change. The strength of the method is that the occurrences are identified individually, even though they have spatial overlap, as they can be separated in the time series information.
We instantiate the RegionGrowingAlgorithm class using the parameters of Anders et al, 2021, and run the object extraction:
# Parametrize the 4D-OBC extraction
algo = py4dgeo.RegionGrowingAlgorithm(
window_width=14,
minperiod=2,
height_threshold=0.05,
neighborhood_radius=1.0,
min_segments=10,
thresholds=[0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9],
seed_candidates=[cp_idx_sel],
)
# Run the algorithm
analysis.invalidate_results(
seeds=True, objects=True, smoothed_distances=False
) # Only required if you want to re-run the algorithm
objects = algo.run(analysis)
seed_timeseries = analysis.smoothed_distances[cp_idx_sel]
plt.plot(
timestamps,
seed_timeseries,
c="black",
linestyle="--",
linewidth=0.5,
label="Seed timeseries",
)
for sid, example_seed in enumerate(analysis.seeds):
seed_end = example_seed.end_epoch
seed_start = example_seed.start_epoch
seed_cp_idx = example_seed.index
plt.plot(
timestamps[seed_start : seed_end + 1],
seed_timeseries[seed_start : seed_end + 1],
label=f"Seed {sid}",
)
# Format the date labels
dtFmt = mdates.DateFormatter("%b-%d")
plt.gca().xaxis.set_major_formatter(dtFmt)
# Add plot elements
plt.xlabel("Date")
plt.ylabel("Distance [m]")
plt.legend()
plt.show()
At the selected location, three separate surface activities are detected as seed for 4D-OBC extraction. We may not be satisfied with the determination of the timing. In a real analysis, we would now aim to improve the temporal detection - either by using a different approach of seed detection (cf. algorithm customization), or by postprocessing the temporal segments from the current seed detection.
Here, we use the result and look into the spatial object properties:
fig, axs = plt.subplots(1, 2, figsize=(15, 7))
ax1, ax2 = axs
# Get indices of 4D-OBC
sel_seed_idx = 1
idxs = objects[sel_seed_idx].indices
# Get change values at end of object for each location
epoch_of_interest = int(
(objects[sel_seed_idx].end_epoch - objects[sel_seed_idx].start_epoch) / 2
+ objects[sel_seed_idx].start_epoch
)
distances_of_interest = analysis.smoothed_distances[:, epoch_of_interest]
# Get the change magnitude between end and start of object for each location
magnitudes_of_interest = (
analysis.smoothed_distances[:, epoch_of_interest]
- analysis.smoothed_distances[:, int(objects[sel_seed_idx].start_epoch)]
)
# Set the colormap according to magnitude at each location in the object
crange = 1.0
import matplotlib.colors as mcolors
cmap = plt.get_cmap("seismic_r").copy()
norm = mcolors.CenteredNorm(halfrange=crange)
cmapvals = norm(magnitudes_of_interest)
# Plot the timeseries of the segmented locations (colored by time series similarity)
for idx in idxs[::10]:
ax1.plot(
timestamps,
analysis.smoothed_distances[idx],
c=cmap(cmapvals[idx]),
linewidth=0.5,
)
# Plot the seed time series
ax1.plot(
timestamps,
analysis.smoothed_distances[cp_idx_sel],
c="black",
linewidth=1.0,
label="Seed timeseries",
)
# Fill the area of the object
ax1.axvspan(
timestamps[objects[sel_seed_idx].start_epoch],
timestamps[objects[sel_seed_idx].end_epoch],
alpha=0.3,
color="grey",
label="4D-OBC timespan",
)
# Add legend and format the date labels
dtFmt = mdates.DateFormatter("%b-%d")
plt.gca().xaxis.set_major_formatter(dtFmt)
ax1.legend()
# Get subset of core points incorporated in 4D-OBC
cloud = analysis.corepoints.cloud
subset_cloud = cloud[idxs, :2]
# Plot coordinates colored by change values at end magnitude of object
d = ax2.scatter(
cloud[:, 0],
cloud[:, 1],
c=magnitudes_of_interest,
cmap="seismic_r",
vmin=-crange,
vmax=crange,
s=1,
)
plt.colorbar(d, format=("%.2f"), label="Change magnitude [m]", ax=ax2)
ax2.set_aspect("equal")
# Plot convex hull of 4D-OBC
from scipy.spatial import ConvexHull
from matplotlib.patches import Polygon
hull = ConvexHull(subset_cloud)
ax2.add_patch(
Polygon(subset_cloud[hull.vertices, 0:2], label="4D-OBC hull", fill=False)
)
# Plot seed location of 4D-OBC
ax2.scatter(
cloud[cp_idx_sel, 0], cloud[cp_idx_sel, 1], marker="*", c="black", label="Seed"
)
# Add plot elements
ax1.set_title("Time series of segmented 4D-OBC locations")
ax1.set_xlabel("Date")
ax1.set_ylabel("Distance [m]")
ax2.set_title(
f"Magnitudes of change in the 4D-OBC timespan\n({timestamps[epoch_of_interest]-timestamps[analysis.objects[sel_seed_idx].start_epoch]} hours)"
)
ax2.set_xlabel("X [m]")
ax2.set_ylabel("Y [m]")
ax2.legend(loc="upper right")
plt.tight_layout()
plt.show()
You may now check the different spatial extents for the three objects extracted from this seed location. In subsequent analysis, the spatial-temporal extent of each object can be used to describe individual activities (e.g., their change rates, magnitudes, ...) and relating them to one another in their timing and spatial distribution.