Demo of a full DePSI processing workflow¶

This notebook uses a coregistered and georeferenced Sentinel-1 interferogram stack over Amsterdam as an example dataset. Please download the data and unzip it locally.

In [ ]:

import numpy as np
from matplotlib import pyplot as plt
import matplotlib.colors as plc
from datetime import datetime

from pathlib import Path

import depsi
from depsi.transformations import latlonh_to_xyz

import sarxarray

import numpy as np from matplotlib import pyplot as plt import matplotlib.colors as plc from datetime import datetime from pathlib import Path import depsi from depsi.transformations import latlonh_to_xyz import sarxarray

Create an interferometric stack in Zarr format¶

In [ ]:

WAVELENGTH = 0.055465763  # Sentinel-1 wavelength in meters
MOTHER_IDX = 0  # Reference acquisition index

WAVELENGTH = 0.055465763 # Sentinel-1 wavelength in meters MOTHER_IDX = 0 # Reference acquisition index

In [ ]:

# Path to the interferogram dataset
path = Path("./nl_amsterdam_s1_asc_t088")

# Make a list of SLCs to read
f_slc = "cint_srd.raw"
list_slcs = [p for p in path.rglob("*_cint_srd.raw")]
list_slcs.sort()

# Geo-referenced coordinates
f_lat = [path / "lat.raw"]
f_lon = [path / "lon.raw"]

# Get time coordinates from filenames
# Note: this assumes filenames start with date in YYYYMMDD format
# It is not a robust way to get time coordinates
# We recommend using metadata files for real applications
time_coords = []
for f in list_slcs:
    date_str = f.name.split("_")[0]
    time_coords.append(np.datetime64(datetime.strptime(date_str, "%Y%m%d")))
time_coords = np.array(time_coords)

# Metadata of the stack, assume known.
shape = (2000, 4000)

list_slcs

# Path to the interferogram dataset path = Path("./nl_amsterdam_s1_asc_t088") # Make a list of SLCs to read f_slc = "cint_srd.raw" list_slcs = [p for p in path.rglob("*_cint_srd.raw")] list_slcs.sort() # Geo-referenced coordinates f_lat = [path / "lat.raw"] f_lon = [path / "lon.raw"] # Get time coordinates from filenames # Note: this assumes filenames start with date in YYYYMMDD format # It is not a robust way to get time coordinates # We recommend using metadata files for real applications time_coords = [] for f in list_slcs: date_str = f.name.split("_")[0] time_coords.append(np.datetime64(datetime.strptime(date_str, "%Y%m%d"))) time_coords = np.array(time_coords) # Metadata of the stack, assume known. shape = (2000, 4000) list_slcs

In [ ]:

# Load complex data
reading_chunks = (500, 500)
stack = sarxarray.from_binary(
    list_slcs, shape, dtype=np.complex64, chunks=reading_chunks
)

# Assign time coordinates
stack["time"] = (("time"), time_coords)

# Load space coordinates
lat = sarxarray.from_binary(
    f_lat, shape, vlabel="lat", dtype=np.float32, chunks=reading_chunks
)
lon = sarxarray.from_binary(
    f_lon, shape, vlabel="lon", dtype=np.float32, chunks=reading_chunks
)
stack = stack.assign_coords(
    lat=(("azimuth", "range"), lat.squeeze().lat.data),
    lon=(("azimuth", "range"), lon.squeeze().lon.data),
)

# DEBUG ONLY: simulate h2ph data
# In real applications, h2ph should be loaded from actual data, e.g. a binary file or a Zarr store
rng = np.random.default_rng(42)
h2ph = rng.uniform(1e-4, 1e-3, (2000, 4000, 10))
stack = stack.assign(h2ph=(("azimuth", "range", "time"), h2ph))

stack

# Load complex data reading_chunks = (500, 500) stack = sarxarray.from_binary( list_slcs, shape, dtype=np.complex64, chunks=reading_chunks ) # Assign time coordinates stack["time"] = (("time"), time_coords) # Load space coordinates lat = sarxarray.from_binary( f_lat, shape, vlabel="lat", dtype=np.float32, chunks=reading_chunks ) lon = sarxarray.from_binary( f_lon, shape, vlabel="lon", dtype=np.float32, chunks=reading_chunks ) stack = stack.assign_coords( lat=(("azimuth", "range"), lat.squeeze().lat.data), lon=(("azimuth", "range"), lon.squeeze().lon.data), ) # DEBUG ONLY: simulate h2ph data # In real applications, h2ph should be loaded from actual data, e.g. a binary file or a Zarr store rng = np.random.default_rng(42) h2ph = rng.uniform(1e-4, 1e-3, (2000, 4000, 10)) stack = stack.assign(h2ph=(("azimuth", "range", "time"), h2ph)) stack

In [ ]:

mrm = stack.slcstack.mrm()
mrm = mrm.compute()

fig, ax = plt.subplots(figsize=(5, 5))
ax.imshow(mrm, cmap="gray")
ax.set_aspect(2)
im = mrm.plot(ax=ax, cmap="gray")
im.set_clim([0, 40000])
im.colorbar.remove()
ax.axis("off")

mrm = stack.slcstack.mrm() mrm = mrm.compute() fig, ax = plt.subplots(figsize=(5, 5)) ax.imshow(mrm, cmap="gray") ax.set_aspect(2) im = mrm.plot(ax=ax, cmap="gray") im.set_clim([0, 40000]) im.colorbar.remove() ax.axis("off")

PS selection¶

In [ ]:

stm_preselention_pnts = depsi.ps_selection(
    stack,
    method="nmad",
    threshold=0.05,  # Points with NMAD below this threshold are selected
    output_chunks=5000,
    mem_persist=False,
)

# Assign metadata attributes
stm_preselention_pnts = stm_preselention_pnts.assign_attrs({"wavelength": WAVELENGTH})
stm_preselention_pnts = stm_preselention_pnts.assign_attrs({"mother_idx": MOTHER_IDX})

# Compute single differential phase w.r.t. the mother acquisition
decimal_years = (
    stm_preselention_pnts["time"]
    - stm_preselention_pnts["time"].isel(time=stm_preselention_pnts.attrs["mother_idx"])
).dt.days / 365.25
stm_preselention_pnts = stm_preselention_pnts.assign_coords(
    {"decimal_years": (("time"), decimal_years.data)}
)
sd_phase = stm_preselention_pnts["phase"] - stm_preselention_pnts["phase"].isel(
    time=stm_preselention_pnts.attrs["mother_idx"]
)
stm_preselention_pnts = stm_preselention_pnts.assign(sd_phase=sd_phase)


# Create local coordinates with meter units
latlonh = np.array(
    [
        stm_preselention_pnts["lat"].values,
        stm_preselention_pnts["lon"].values,
        np.zeros_like(stm_preselention_pnts["lat"].values),
    ]
).T
xyz = latlonh_to_xyz(latlonh).T
stm_preselention_pnts = stm_preselention_pnts.assign_coords(
    {
        "local_x": (("space"), xyz[:, 0]),
        "local_y": (("space"), xyz[:, 1]),
    }
)

stm_preselention_pnts

stm_preselention_pnts = depsi.ps_selection( stack, method="nmad", threshold=0.05, # Points with NMAD below this threshold are selected output_chunks=5000, mem_persist=False, ) # Assign metadata attributes stm_preselention_pnts = stm_preselention_pnts.assign_attrs({"wavelength": WAVELENGTH}) stm_preselention_pnts = stm_preselention_pnts.assign_attrs({"mother_idx": MOTHER_IDX}) # Compute single differential phase w.r.t. the mother acquisition decimal_years = ( stm_preselention_pnts["time"] - stm_preselention_pnts["time"].isel(time=stm_preselention_pnts.attrs["mother_idx"]) ).dt.days / 365.25 stm_preselention_pnts = stm_preselention_pnts.assign_coords( {"decimal_years": (("time"), decimal_years.data)} ) sd_phase = stm_preselention_pnts["phase"] - stm_preselention_pnts["phase"].isel( time=stm_preselention_pnts.attrs["mother_idx"] ) stm_preselention_pnts = stm_preselention_pnts.assign(sd_phase=sd_phase) # Create local coordinates with meter units latlonh = np.array( [ stm_preselention_pnts["lat"].values, stm_preselention_pnts["lon"].values, np.zeros_like(stm_preselention_pnts["lat"].values), ] ).T xyz = latlonh_to_xyz(latlonh).T stm_preselention_pnts = stm_preselention_pnts.assign_coords( { "local_x": (("space"), xyz[:, 0]), "local_y": (("space"), xyz[:, 1]), } ) stm_preselention_pnts

In [ ]:

# Plot preselected points over mrm
fig, ax = plt.subplots(figsize=(5, 5))
ax.imshow(mrm, cmap="gray")
ax.set_aspect(2)
im = mrm.plot(ax=ax, cmap="gray")
im.set_clim([0, 40000])
im.colorbar.remove()
ax.axis("off")
ax.scatter(
    stm_preselention_pnts["range"],
    stm_preselention_pnts["azimuth"],
    s=0.3,
    c=stm_preselention_pnts["full_ts_nmad"],
    marker="s",
    cmap="jet_r",
    norm=plc.LogNorm(vmin=0.02, vmax=0.04),
)

# Plot preselected points over mrm fig, ax = plt.subplots(figsize=(5, 5)) ax.imshow(mrm, cmap="gray") ax.set_aspect(2) im = mrm.plot(ax=ax, cmap="gray") im.set_clim([0, 40000]) im.colorbar.remove() ax.axis("off") ax.scatter( stm_preselention_pnts["range"], stm_preselention_pnts["azimuth"], s=0.3, c=stm_preselention_pnts["full_ts_nmad"], marker="s", cmap="jet_r", norm=plc.LogNorm(vmin=0.02, vmax=0.04), )

Network estimation¶

In [ ]:

# Select network points
# First selcect stm_preselention_pnts with full_ts_nmad <0.03
mask = stm_preselention_pnts["full_ts_nmad"] < 0.03
stm_network_pnts = stm_preselention_pnts.where(mask.compute(), drop=True)
# Then Make points sparser
# Assuming setinel-1 pixels size
stm_network_pnts = depsi.network_stm_selection(
    stm_network_pnts,
    min_dist=300,
    azimuth_spacing=10,
    range_spacing=5,
    sortby_var="full_ts_nmad",  # The sort is in ascending order so that points with smaller NMAD are prioritized
)

stm_network_pnts

# Select network points # First selcect stm_preselention_pnts with full_ts_nmad <0.03 mask = stm_preselention_pnts["full_ts_nmad"] < 0.03 stm_network_pnts = stm_preselention_pnts.where(mask.compute(), drop=True) # Then Make points sparser # Assuming setinel-1 pixels size stm_network_pnts = depsi.network_stm_selection( stm_network_pnts, min_dist=300, azimuth_spacing=10, range_spacing=5, sortby_var="full_ts_nmad", # The sort is in ascending order so that points with smaller NMAD are prioritized ) stm_network_pnts

In [ ]:

# Plot preselected points over mrm
fig, ax = plt.subplots(figsize=(5, 5))
ax.imshow(mrm, cmap="gray")
ax.set_aspect(2)
im = mrm.plot(ax=ax, cmap="gray")
im.set_clim([0, 40000])
im.colorbar.remove()
ax.axis("off")
ax.scatter(
    stm_network_pnts["range"],
    stm_network_pnts["azimuth"],
    s=2,
    c=stm_network_pnts["full_ts_nmad"],
    marker="s",
    cmap="jet_r",
    norm=plc.LogNorm(vmin=0.01, vmax=0.05),
)

# Plot preselected points over mrm fig, ax = plt.subplots(figsize=(5, 5)) ax.imshow(mrm, cmap="gray") ax.set_aspect(2) im = mrm.plot(ax=ax, cmap="gray") im.set_clim([0, 40000]) im.colorbar.remove() ax.axis("off") ax.scatter( stm_network_pnts["range"], stm_network_pnts["azimuth"], s=2, c=stm_network_pnts["full_ts_nmad"], marker="s", cmap="jet_r", norm=plc.LogNorm(vmin=0.01, vmax=0.05), )

In [ ]:

# Form network arcs
stm_network_arcs = depsi.form_network(
    stm_network_pnts,
    key_phase="sd_phase",
    key_h2ph="h2ph",
    key_Btemporal="decimal_years",
    key_xcrds="local_x",
    key_ycrds="local_y",
    network_method="delaunay",
    max_length=800,  # in meters
)
stm_network_arcs

# Form network arcs stm_network_arcs = depsi.form_network( stm_network_pnts, key_phase="sd_phase", key_h2ph="h2ph", key_Btemporal="decimal_years", key_xcrds="local_x", key_ycrds="local_y", network_method="delaunay", max_length=800, # in meters ) stm_network_arcs

In [ ]:

# arc estimation using periodogram
_, ambiguities, _, _, ens_coh = depsi.periodogram(
    stm_network_arcs,
    key_dphase="d_phase",
    key_h2ph="h2ph",
    key_Btemporal="Btemp",
)
stm_network_arcs["ambiguities"] = ambiguities
stm_network_arcs["temp_coh"] = ens_coh

stm_network_arcs = stm_network_arcs.compute()
stm_network_arcs

# arc estimation using periodogram _, ambiguities, _, _, ens_coh = depsi.periodogram( stm_network_arcs, key_dphase="d_phase", key_h2ph="h2ph", key_Btemporal="Btemp", ) stm_network_arcs["ambiguities"] = ambiguities stm_network_arcs["temp_coh"] = ens_coh stm_network_arcs = stm_network_arcs.compute() stm_network_arcs

In [ ]:

# Plot preselected points over mrm
fig, ax = plt.subplots(figsize=(5, 5))
ax.imshow(mrm, cmap="gray")
ax.set_aspect(2)
im = mrm.plot(ax=ax, cmap="gray")
im.set_clim([0, 40000])
im.colorbar.remove()
ax.axis("off")

# Plot network arcs
# x of sources and targets
xx = np.stack(
    [
        stm_network_pnts.isel(space=stm_network_arcs["source"])["range"].values,
        stm_network_pnts.isel(space=stm_network_arcs["target"])["range"].values,
    ]
).T
# y of sources and targets
yy = np.stack(
    [
        stm_network_pnts.isel(space=stm_network_arcs["source"])["azimuth"].values,
        stm_network_pnts.isel(space=stm_network_arcs["target"])["azimuth"].values,
    ]
).T
# Visualize created arcs
cmap = plt.cm.jet_r
norm = plc.Normalize(vmin=0, vmax=1.0)
temp_coh = np.abs(stm_network_arcs["temp_coh"].data)
for i in range(stm_network_arcs.sizes["space"]):
    ax.plot(xx[i], yy[i], color=cmap(norm(temp_coh[i])), linewidth=0.5)

# Plot preselected points over mrm fig, ax = plt.subplots(figsize=(5, 5)) ax.imshow(mrm, cmap="gray") ax.set_aspect(2) im = mrm.plot(ax=ax, cmap="gray") im.set_clim([0, 40000]) im.colorbar.remove() ax.axis("off") # Plot network arcs # x of sources and targets xx = np.stack( [ stm_network_pnts.isel(space=stm_network_arcs["source"])["range"].values, stm_network_pnts.isel(space=stm_network_arcs["target"])["range"].values, ] ).T # y of sources and targets yy = np.stack( [ stm_network_pnts.isel(space=stm_network_arcs["source"])["azimuth"].values, stm_network_pnts.isel(space=stm_network_arcs["target"])["azimuth"].values, ] ).T # Visualize created arcs cmap = plt.cm.jet_r norm = plc.Normalize(vmin=0, vmax=1.0) temp_coh = np.abs(stm_network_arcs["temp_coh"].data) for i in range(stm_network_arcs.sizes["space"]): ax.plot(xx[i], yy[i], color=cmap(norm(temp_coh[i])), linewidth=0.5)

In [ ]:

# Integrate ambiguities from arcs to points
stm_network_arcs_solved, stm_network_pnts_solved = depsi.spatial_integration(
    stm_network_pnts,
    stm_network_arcs,
    key_arc_quality="temp_coh",
    threshold_arc_quality=0.45,
)
stm_network_pnts_solved

# Integrate ambiguities from arcs to points stm_network_arcs_solved, stm_network_pnts_solved = depsi.spatial_integration( stm_network_pnts, stm_network_arcs, key_arc_quality="temp_coh", threshold_arc_quality=0.45, ) stm_network_pnts_solved

In [ ]:

stm_network_pnts_solved, param_list = depsi.estimate_model_params(
    stm_network_pnts_solved,
    key_observations="unwrapped_phase",
    key_time="decimal_years",
)
stm_network_pnts_solved = stm_network_pnts_solved.compute()
stm_network_pnts_solved

stm_network_pnts_solved, param_list = depsi.estimate_model_params( stm_network_pnts_solved, key_observations="unwrapped_phase", key_time="decimal_years", ) stm_network_pnts_solved = stm_network_pnts_solved.compute() stm_network_pnts_solved

In [ ]:

# Plot network points over mrm
fig, ax = plt.subplots(figsize=(5, 5))
ax.imshow(mrm, cmap="gray")
ax.set_aspect(2)
im = mrm.plot(ax=ax, cmap="gray")
im.set_clim([0, 40000])
im.colorbar.remove()
ax.axis("off")
ax.scatter(
    stm_network_pnts_solved["range"].values,
    stm_network_pnts_solved["azimuth"].values,
    s=2,
    c=stm_network_pnts_solved["pnt_velocity"].values,
    marker="s",
    cmap="jet_r",
)

# Plot network points over mrm fig, ax = plt.subplots(figsize=(5, 5)) ax.imshow(mrm, cmap="gray") ax.set_aspect(2) im = mrm.plot(ax=ax, cmap="gray") im.set_clim([0, 40000]) im.colorbar.remove() ax.axis("off") ax.scatter( stm_network_pnts_solved["range"].values, stm_network_pnts_solved["azimuth"].values, s=2, c=stm_network_pnts_solved["pnt_velocity"].values, marker="s", cmap="jet_r", )

Densification¶

In [ ]:

stm_densified = depsi.densification(
    stm_preselention_pnts,
    stm_network_pnts_solved,
    key_h2ph="h2ph",
    key_Btemporal="decimal_years",
)
stm_densified = stm_densified.compute()
stm_densified

stm_densified = depsi.densification( stm_preselention_pnts, stm_network_pnts_solved, key_h2ph="h2ph", key_Btemporal="decimal_years", ) stm_densified = stm_densified.compute() stm_densified

In [ ]:

stm_densified, param_list = depsi.estimate_model_params(
    stm_densified, key_observations="unwrapped_phase", key_time="decimal_years"
)
stm_densified = stm_densified.compute()
stm_densified

stm_densified, param_list = depsi.estimate_model_params( stm_densified, key_observations="unwrapped_phase", key_time="decimal_years" ) stm_densified = stm_densified.compute() stm_densified

In [ ]:

# Plot network points over mrm
fig, ax = plt.subplots(figsize=(5, 5))
ax.imshow(mrm, cmap="gray")
ax.set_aspect(2)
im = mrm.plot(ax=ax, cmap="gray")
im.set_clim([0, 40000])
im.colorbar.remove()
ax.axis("off")
ax.scatter(
    stm_densified["range"].values,
    stm_densified["azimuth"].values,
    s=0.3,
    c=stm_densified["pnt_velocity"].values,
    marker="s",
    cmap="jet_r",
)

# Plot network points over mrm fig, ax = plt.subplots(figsize=(5, 5)) ax.imshow(mrm, cmap="gray") ax.set_aspect(2) im = mrm.plot(ax=ax, cmap="gray") im.set_clim([0, 40000]) im.colorbar.remove() ax.axis("off") ax.scatter( stm_densified["range"].values, stm_densified["azimuth"].values, s=0.3, c=stm_densified["pnt_velocity"].values, marker="s", cmap="jet_r", )