Appendix I — Homework Exercise 9: Phase Unwrapping Exercise

Your task in this exercise is to read the provided data, plot it, unwrap and plot the phase, and calculate the displacement. Then interpret your results. The following data were acquired from Sentinel-1 on April 23, 2018, from 16:15:24 to 16:15:51 over the region of Hawaii (same data as 08_homework_exercise).

import cmcrameri as cmc  # noqa: F401
import intake
import matplotlib.pyplot as plt
import numpy as np
import seaborn as sns
import xarray as xr

from mrs.catalog import get_intake_url
from mrs.processing import unwrap_array  # noqa: F401

url = get_intake_url()
cat = intake.open_catalog(url)
ds = cat.ex9.read().compute()


# Set cyclic and linear colormaps
cmap_cyc = sns.color_palette("hls", as_cmap=True)
cmap_lin = "cmc.roma_r"
cmap_disp = "cmc.vik"

# Create a mask for the areas which have no data
mask = (ds.phase != 0).astype(bool)

https://git.geo.tuwien.ac.at/public_projects/microwave-remote-sensing/-/raw/main/microwave-remote-sensing.yml

/home/runner/work/microwave-remote-sensing/microwave-remote-sensing/.conda_envs/microwave-remote-sensing/lib/python3.13/site-packages/xarray/coding/variables.py:432: ComplexWarning: Casting complex values to real discards the imaginary part
  data = data.astype(dtype=dtype, copy=True)

def compute_complex_band(phase_band, intensity_band):
    """
    Computes the complex-valued band for SAR interferograms from phase and intensity bands.
    
    Parameters:
    - phase_band (numpy.ndarray): Phase band in radians.
    - intensity_band (numpy.ndarray): Intensity band (magnitude).
    
    Returns:
    - numpy.ndarray: Complex-valued band with dtype complex128.
    """
    real_part = intensity_band * np.cos(phase_band)
    imag_part = intensity_band * np.sin(phase_band)
    complex_band = real_part + 1j * imag_part
    return complex_band


ds["cmplx"] = compute_complex_band(ds["phase"], ds["intensity"])
ds

<xarray.Dataset> Size: 588MB
Dimensions:      (y: 6110, x: 3435)
Coordinates:
  * y            (y) float64 49kB -0.0001007 -0.000302 ... -1.23 -1.23
  * x            (x) float64 27kB 0.0001525 0.0004576 0.0007626 ... 1.047 1.048
    band         int64 8B 1
    spatial_ref  int64 8B 0
Data variables:
    coh          (y, x) float32 84MB 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0
    complex      (y, x) float64 168MB 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0
    intensity    (y, x) float32 84MB 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0
    phase        (y, x) float32 84MB 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0
    cmplx        (y, x) complex64 168MB 0j 0j 0j 0j 0j 0j ... 0j 0j 0j 0j 0j 0j

xarray.Dataset

• Dimensions:
  • y: 6110
  • x: 3435
• Coordinates: (4)
  • y
    (y)
    float64
    -0.0001007 -0.000302 ... -1.23

    array([-1.006741e-04, -3.020223e-04, -5.033705e-04, ..., -1.229734e+00,
           -1.229935e+00, -1.230137e+00], shape=(6110,))

  • x
    (x)
    float64
    0.0001525 0.0004576 ... 1.047 1.048

    array([1.525203e-04, 4.575610e-04, 7.626016e-04, ..., 1.047052e+00,
           1.047357e+00, 1.047662e+00], shape=(3435,))

  • band
    ()
    int64
    1

    array(1)

  • spatial_ref
    ()
    int64
    0

    GeoTransform :

    0.0 0.0003050406542399742 0.0 0.0 0.0 -0.0002013481828520014

    crs_wkt :

    GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]

    geographic_crs_name :

    WGS 84

    grid_mapping_name :

    latitude_longitude

    horizontal_datum_name :

    World Geodetic System 1984

    inverse_flattening :

    298.257223563

    longitude_of_prime_meridian :

    0.0

    prime_meridian_name :

    Greenwich

    reference_ellipsoid_name :

    WGS 84

    semi_major_axis :

    6378137.0

    semi_minor_axis :

    6356752.314245179

    spatial_ref :

    GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]

    array(0)

• Data variables: (5)
  • coh
    (y, x)
    float32
    0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0

    AREA_OR_POINT :

    Area

    TIFFTAG_RESOLUTIONUNIT :

    1 (unitless)

    TIFFTAG_XRESOLUTION :

    1

    TIFFTAG_YRESOLUTION :

    1

    array([[0., 0., 0., ..., 0., 0., 0.],
           [0., 0., 0., ..., 0., 0., 0.],
           [0., 0., 0., ..., 0., 0., 0.],
           ...,
           [0., 0., 0., ..., 0., 0., 0.],
           [0., 0., 0., ..., 0., 0., 0.],
           [0., 0., 0., ..., 0., 0., 0.]], shape=(6110, 3435), dtype=float32)

  • complex
    (y, x)
    float64
    0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0

    AREA_OR_POINT :

    Area

    TIFFTAG_RESOLUTIONUNIT :

    1 (unitless)

    TIFFTAG_XRESOLUTION :

    1

    TIFFTAG_YRESOLUTION :

    1

    array([[0., 0., 0., ..., 0., 0., 0.],
           [0., 0., 0., ..., 0., 0., 0.],
           [0., 0., 0., ..., 0., 0., 0.],
           ...,
           [0., 0., 0., ..., 0., 0., 0.],
           [0., 0., 0., ..., 0., 0., 0.],
           [0., 0., 0., ..., 0., 0., 0.]], shape=(6110, 3435))

  • intensity
    (y, x)
    float32
    0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0

    AREA_OR_POINT :

    Area

    TIFFTAG_RESOLUTIONUNIT :

    1 (unitless)

    TIFFTAG_XRESOLUTION :

    1

    TIFFTAG_YRESOLUTION :

    1

    array([[0., 0., 0., ..., 0., 0., 0.],
           [0., 0., 0., ..., 0., 0., 0.],
           [0., 0., 0., ..., 0., 0., 0.],
           ...,
           [0., 0., 0., ..., 0., 0., 0.],
           [0., 0., 0., ..., 0., 0., 0.],
           [0., 0., 0., ..., 0., 0., 0.]], shape=(6110, 3435), dtype=float32)

  • phase
    (y, x)
    float32
    0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0

    AREA_OR_POINT :

    Area

    TIFFTAG_RESOLUTIONUNIT :

    1 (unitless)

    TIFFTAG_XRESOLUTION :

    1

    TIFFTAG_YRESOLUTION :

    1

    array([[0., 0., 0., ..., 0., 0., 0.],
           [0., 0., 0., ..., 0., 0., 0.],
           [0., 0., 0., ..., 0., 0., 0.],
           ...,
           [0., 0., 0., ..., 0., 0., 0.],
           [0., 0., 0., ..., 0., 0., 0.],
           [0., 0., 0., ..., 0., 0., 0.]], shape=(6110, 3435), dtype=float32)

  • cmplx
    (y, x)
    complex64
    0j 0j 0j 0j 0j ... 0j 0j 0j 0j 0j

    AREA_OR_POINT :

    Area

    TIFFTAG_RESOLUTIONUNIT :

    1 (unitless)

    TIFFTAG_XRESOLUTION :

    1

    TIFFTAG_YRESOLUTION :

    1

    array([[0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],
           [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],
           [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],
           ...,
           [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],
           [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j],
           [0.+0.j, 0.+0.j, 0.+0.j, ..., 0.+0.j, 0.+0.j, 0.+0.j]],
          shape=(6110, 3435), dtype=complex64)

I.1 Question 1

Plot the above provided data. Please note that we are plotting images in Slant Range geometry (East-West flipped due to the descending acquisition pass).

fig, axes = plt.subplots(1, 3, figsize=(15, 5))
ds.phase.where(mask).plot.imshow(cmap=cmap_cyc, ax=axes[0]).axes.set_title("Phase")
# YOUR CODE HERE -> (Plot intensity map)
# YOUR CODE HERE -> (Plot coherence map)
plt.tight_layout()

I.2 Question 2

Use the above function to unwrap the given interferogram and to plot your result. Use a coherence threshold of 0.2 for your unwrapping and use the provided lines of code below to coarsen the input (to fasten the process). Kernels with values higher than 3 might lead the unwrapping function to crash, therefore please do not change it. It might take up to 1 or 2 minutes to execute the function, let the code run, and the output be printed (it will be cleared out at the end).

# coarsen
kernel = 3
ds = ds.coarsen(x=kernel, y=kernel, boundary="trim").mean()

unwrapped = ...  # YOUR CODE HERE -> perform unwrapping function

unwrapped.plot.imshow(cmap=cmap_cyc)

def displacement(unwrapped: xr.DataArray, lambda_val: float = 0.056) -> xr.DataArray:
    """Calculate the displacement from the unwrapped phase.

    Parameters
    ----------
    unwrapped : xr.DataArray
        Array of the unwrapped phase.
    lambda_val : float, optional
        Wavelength of the radar signal, by default 0.056

    Returns
    -------
    xr.DataArray
        Array of the displacement.

    """
    return unwrapped * -lambda_val / (4 * np.pi)

I.3 Question 3

Use the above function to calculate the displacement from the unwrapped interferogram, then plot your result with the correct colormap. Can you guess which kind of phenomenon is responsible for the ground surface displacement that you observe? Elaborate on the direction of the displacement movement (is it towards or away from satellite?)

# YOUR CODE HERE -> perform displacement function
# YOUR CODE HERE -> plot displacement map