Notebook

Interferograms¶

On July 5th, 2019, an earthquake with a magnitude of 7.1 mainshock struck eastern California, near the city of Ridgecrest. The seismic event produced a surface rupture spanning more than 50 kilometers with a complex vertical and horizontal offset pattern along the main fault line. SAR imagery can be employed for accurately measuring and describing ground motion through a well-established technique called SAR Interferometry. In this framework, the phase information contained in Synthetic Aperture Radar (SAR) data is employed. In this notebook, we will dive into the main interferometric SAR processing operations which involves retrieving the difference between the phase signals of repeated SAR acquisitions to analyze the shape and deformation of the Earth's surface. In our case, we will use a pair of Single Look Complex (SLC) Sentinel-1 images to obtain an interferogram of the Ridgecrest earthquake.

This notebook will outline the process of working with interferograms and the steps needed to extract valuable information. Here, we will focus on displaying products generated by the Sentinel Application Platform (SNAP) software from the European Space Agency (ESA).

Photo by Brian Olson / California Geological Survey

In [ ]:

import base64
from io import BytesIO

import folium
import holoviews as hv
import hvplot.xarray  # noqa: F401
import intake
import matplotlib.pyplot as plt
import seaborn as sns

hv.extension("bokeh")

Single Look Complex (SLC) Data¶

We introduce now another level-1 radar product type, which is called Single Look Complex (SLC). Interferometry can only be performed with SLC data. What are the main differences between SLC and GRD (the other level-1 radar product)?

SLC vs GRD:
- SLC contains complex-value data (amplitude and phase) vs GRD contains intensity only (amplitude)
- SLC geometry is Slant Range (radar's line of sight) vs GRD data are projected onto ground range
- SLC resolution is full vs GRD has lower resolution (it is multi-looked)
- SLC supports phase-based applications (Interferometry) vs GRD supports only amplitude-based ones
- SLC has larger file sizes compared to GRD

In [ ]:

url = "https://huggingface.co/datasets/martinschobben/microwave-remote-sensing/resolve/main/microwave-remote-sensing.yml"
cat = intake.open_catalog(url)
iw1_ds = cat.iw1.read()
iw2_ds = cat.iw2.read()
iw3_ds = cat.iw3.read()

Let's plot all three sub-swaths to view the full scene acquired by the satellite. The acquisition times for each swath on July 10th, 2019 are the following:

IW1 at 01:50:01 - 01:50:26
IW2 at 01:49:59 - 01:50:24
IW3 at 01:50:00 - 01:50:25

In [ ]:

fig, ax = plt.subplots(1, 3, figsize=(15, 7), sharey=True)

datasets = [iw1_ds, iw2_ds, iw3_ds]
val_range = dict(vmin=0, vmax=255, cmap="gray")

for i, ds in enumerate(datasets):
    im = ds.intensity.plot(ax=ax[i], add_colorbar=False, **val_range)
    ax[i].tick_params(axis="both", which="major")

cbar = fig.colorbar(im, ax=ax, orientation="horizontal", shrink=0.9, pad=0.2)

plt.show()

We don’t need all three of the subswaths for our notebook, so we will focus on IW2 and display its intensity and phase measurements.

In [ ]:

# Compute the intensity and phase from complex data
cmap_hls = sns.hls_palette(as_cmap=True)

fig, axes = plt.subplots(1, 2, figsize=(15, 6))

ds.intensity.plot(ax=axes[0], cmap="gray", robust=True)
axes[0].set_title("Intensity Measurement of IW2")

ds.phase.plot(ax=axes[1], cmap=cmap_hls)
axes[1].set_title("Phase Measurement of IW2")

plt.tight_layout()

Intensity is represented in an 8-bit format (ranging from 0 to 255), while phase measurements range from $- \pi$ to $\pi$ . At first glance, phase does not correspond to easily observable physical properties of the ground. However, the phase becomes incredibly valuable when, for example, it is used comparatively between two successive phase measurements (two Sentinel-1 images acquired at different times over the same area). Here are the processing steps needed to retrieve a difference between the phases of two radar acquisitions:

Coregistering¶

Before creating an interferogram, measurements from two different dates need to be coregistered. This means that each pixel from the two acquisitions must be precisely aligned so that they are representing the same ground object. Accurate and successful co-registration of the two (or more) images is vital for interferometry processing. We call the "master" image the reference image (typically the earliest acquisition in time) to which we coregister the "slave" image (typically acquired later in time).

In [ ]:

coregistered_ds = cat.coreg.read()

fig, axes = plt.subplots(1, 2, figsize=(18, 8))
coregistered_ds.band_data.sel(band=1).plot(ax=axes[0], cmap="gray", robust=True)
axes[0].set_title("Master Phase Measurement - 28 Jun 2019")

coregistered_ds.band_data.sel(band=2).plot(ax=axes[1], cmap="gray", robust=True)
axes[1].set_title("Slave Phase Measurement - 10 Jul 2019")

plt.tight_layout()

Interferogram Formation and Coherence Estimation¶

The interferogram formation process combines the amplitudes of both images and calculates the difference between their respective phases at each SAR image pixel (cross-multiplication of the master image with the complex conjugate of the slave image).

After building up the interferogram, we have to take into account the presence of other contributing terms that could hinder our goal of measuring the surface deformation due to the earthquake. For example, we need to subtract from the interferogram the flat-earth phase contribution, which is a signal contribution due to the curvature of the Earth's surface. This is here done automatically through the SNAP software operators.

In general, the accuracy of interferometric measurements are influenced by many contributors that could result in a loss of coherence. But what is coherence? It is a measure of phase correlation between the master and slave image. Interferometric coherence (γ) can be expressed as:

$γ = γ_{proc}*γ_{geom}*γ_{vol}*γ_{SNR}*γ_{temp}$

where $γ_{proc}$ refers to inaccuracies in the processing (e.g., coregistration errors), $γ_{geom}$ refers to the baseline decorrelation (different position of satellites during the two acquisitions), $γ_{vol}$ refers to volume decorrelation (vegetation related), $γ_{SNR}$ refers to the radar instrument thermal noise and $γ_{temp}$ refers to the decorrelation caused by change of position of the objects in the scene during the time interval of the images acquisitions (e.g., plant growth, wind-induced movements or ground deformation due to earthquakes, landslides).

Therefore, we can conclude that interferometric accuracy is sensitive to many processes, hence isolating the ground deformation signal involves several operations. On the other hand, interferometric coherence sensitivity could be exploited to track and map phenomena that cause its degradation (e.g., vegetation features, and water content).

In [ ]:

interferogram_ds = cat.inter.read()

cmap_hls_hex = sns.color_palette("hls", n_colors=256).as_hex()

interferogram_ds = interferogram_ds.where(interferogram_ds != 0)
igf_da = interferogram_ds.sel(band=1).band_data
coh_da = interferogram_ds.sel(band=2).band_data

# Invert y axis
igf_da["y"] = igf_da.y[::-1]
coh_da["y"] = coh_da.y[::-1]

igf_plot = igf_da.hvplot.image(
    x="x",
    y="y",
    cmap=cmap_hls_hex,
    width=600,
    height=600,
    dynamic=False,
)

coh_plot = coh_da.hvplot.image(
    x="x",
    y="y",
    cmap="viridis",
    width=600,
    height=600,
    dynamic=False,
).opts(clim=(0, 1))

(igf_plot + coh_plot).opts(shared_axes=True)

Now we can observe patterns that emerged between the two acquisitions. If you look at the data range in the interferogram (left plot), you’ll notice it spans approximately one wavelength, from $-\pi$ to $\pi$ . On the right, you find a plot of the interferometric coherence (values ranging between 0 and 1), where low coherence is found along the ground surface ruptures caused by the earthquake. Please note, that the interferogram has already undergone a deburst operation (all bursts merged into a single image).

Topographic Phase Removal¶

Since the local topography is an additional phase term constituting the interferogram that we built up so far, we need to make an estimate of its impact in order to further remove it to keep only the ground deformation-related phase. For this purpose, we use a reference known DEM to simulate an interferogram and to subtract it from the original interferogram.

In [ ]:

topo_ds = cat.topo.read()

fig, axes = plt.subplots(1, 3, figsize=(18, 5))

igf_da.plot(ax=axes[0], cmap=cmap_hls)
axes[0].set_title("Interferogram With Topographic Phase")

topo_ds.topo.plot(ax=axes[1], cmap="gist_earth")
axes[1].set_title("Topography")

topo_ds.Phase.plot(ax=axes[2], cmap=cmap_hls)
axes[2].set_title("Interferogram Without Topographic Phase")

plt.tight_layout()

Multi-looking, Goldstein Phase Filtering and Geocoding¶

In order to improve the phase signatures contained within our interferogram and get a generally higher signal-to-noise (SNR) ratio, we will perform two additional operations called multi-looking and Goldstein phase filtering. Multi-looking is the process of averaging adjacent pixels using a moving window of the interferogram to reduce noise (at the cost of reducing the spatial resolution). Coherence is involved in this operation to flag and set areas to no data that are considered unreliable (low coherence) and to keep the reliable ones (high coherence).

Finally, to make data interpretable, we geocode the wrapped interferogram. So far we performed the interferometric processing in the radar geometry. The transformation into geographic coordinates will help us to perform further comparisons in a real-world coordinate system.

In [ ]:

geocoded_ds = cat.geocode.read()

step = 4  # if you want to zoom in, suggestion is to make this step smaller

geocoded_ds = geocoded_ds.where(geocoded_ds != 0)
igf_data = geocoded_ds.sel(band=1).band_data
coh_da = geocoded_ds.sel(band=2).band_data

igf_plot = igf_data.isel(x=slice(0, -1, step), y=slice(0, -1, step)).hvplot.image(
    x="x", y="y", cmap=cmap_hls_hex, width=600, height=600, dynamic=False
)

coh_plot = (
    coh_da.isel(x=slice(0, -1, step), y=slice(0, -1, step))
    .hvplot.image(x="x", y="y", cmap="viridis", width=600, height=600, dynamic=False)
    .opts(clim=(0, 1))
)

(igf_plot + coh_plot).opts(shared_axes=True)

In the above plot, we can compare georeferenced data in the form of the interferogram (left) and the coherence (right). Along the earthquake fault line, low coherence between the two phase acquisitions is visible. This occurs due to extreme changes in terrain heights or displacements, which are beyond the sensitivity of the SAR sensor. This area of low coherence indicates higher uncertainty in the interferogram. However, this isn't necessarily a drawback, as it helps to clearly identify the earthquake epicenter.

You can also explore and zoom into regions with "fringe patterns" to observe ground movement. Each fringe cycle (e.g., from red to red or blue to blue) corresponds to ground motion in this case. The fringe patterns indicate motion in the line-of-sight (LOS) of the satellite (Sentinel-1 has a mean incidence angle of 38°) in terms of either uplift (relative motion of the ground towards the satellite) or sinking (relative motion of the ground away from the satellite). If the interferogram phase changes from 0 to -3.14 (cycles in the negative direction), the surface is moving away from the satellite (i.e., sinking movement). Vice versa, if cycles go in the positive direction (from 0 to +3.14), it would mean a relative uplifting movement of the ground. In areas with no ground motion, fringe patterns disappear. The radar's sensitivity to motion depends on its wavelength. For Sentinel-1 (~5.6cm), a full fringe cycle ( $2\pi$ ) represents a displacement of ~2.8 cm in the LOS direction.

Visualisation of the Earthquake Event on July 5th, 2019¶

In [ ]:

step = 4  # Downsample data for visualization
igf_data_subset = igf_data.isel(x=slice(0, -1, step), y=slice(0, -1, step))


def array_to_img(data_array, cmap):
    fig, ax = plt.subplots(figsize=(6, 6), dpi=600)
    data_array.plot(ax=ax, cmap=cmap, add_colorbar=False, add_labels=False)
    ax.set_axis_off()
    buf = BytesIO()
    plt.savefig(buf, format="png", bbox_inches="tight", pad_inches=0, transparent=True)
    plt.close(fig)
    return base64.b64encode(buf.getvalue()).decode("utf-8")


igf_image = array_to_img(igf_data_subset, cmap=cmap_hls)
bounds = [
    [float(igf_data["y"].min()), float(igf_data["x"].min())],
    [float(igf_data["y"].max()), float(igf_data["x"].max())],
]

m = folium.Map(
    location=[float(igf_data["y"].mean()), float(igf_data["x"].mean())],
    zoom_start=10,
)
folium.TileLayer(
    tiles=(
        "https://server.arcgisonline.com/ArcGIS/rest/"
        + "services/World_Imagery/MapServer/tile/{z}/{y}/{x}"
    ),
    attr="Tiles © Esri",
    name="ESRI World Imagery",
).add_to(m)
folium.TileLayer(
    tiles=(
        "https://server.arcgisonline.com/ArcGIS/rest/"
        + "services/Reference/World_Boundaries_and_Places/"
        + "MapServer/tile/{z}/{y}/{x}"
    ),
    attr="Tiles © Esri",
    name="ESRI Labels",
    overlay=True,
).add_to(m)

folium.raster_layers.ImageOverlay(
    image=f"data:image/png;base64,{igf_image}",
    bounds=bounds,
    opacity=0.65,
    name="Interferogram",
).add_to(m)
folium.LayerControl().add_to(m)

m