DS Processing¶

Notebook

Converted from nbs/Tutorials/CLI/CampiFlegrei/03_ds.ipynb (text and code cells; the outputs are not stored in the repository). Run it in Jupyter or VS Code for the interactive maps.

Data: Campi Flegrei. Sentinel-1A IW, descending track 22, VV, 92 dates from 2019-01-02 to 2021-12-29 (every 12 days, reference 2020-07-07), cropped to the caldera (Pozzuoli and the Gulf of Pozzuoli) at full resolution. A bell-shaped uplift of tens of centimetres in the centre; the town gives many PS. The data and how they were made are described in data/CampiFlegrei/README.md.

The thresholds below are starting values; look at the figures and adjust them to the data.

In this tutorial, we demostrate how to do DS processing with the decorrelation CLI interface. The input data for DS processing is prepared in Load Data.

!ls ./load_data

SHP selection¶

Here we use two sample Kolmogorov-Smirnov (KS) test for Spatially Homogeneous Pixel (SHP) identification.

The outputs are is_shp, True for the pixels of a window that are SHPs of its centre pixel, and shp_num, the number of SHPs of every pixel. A pixel is an SHP when the test does not reject that both samples come from the same distribution at the significance level alpha, i.e. its p-value is at least alpha. The p-values themselves are an optional output (pvalue), needed only to select the SHPs again at another level with mc.select_shp.

import numpy as np
import moraine.cli as mc
import moraine as mr

rslc = './load_data/rslc.zarr'
is_shp = './ds_processing/gix/ras_is_shp.zarr'
shp_num = './ds_processing/gix/ras_shp_num.zarr'
az_half_win = 3    # 7 x 27 pixels, about 100 m x 60 m on the ground: 92 images need many SHPs
r_half_win = 13
alpha = 0.05   # significance level: SHPs have a p-value >= alpha
method = 'ks'

logger = mc.get_logger()

mc.shp_test(rslc,is_shp,shp_num,
            az_half_win=az_half_win,r_half_win=r_half_win,alpha=alpha)

import zarr
from matplotlib import pyplot as plt
import colorcet

is_shp indicates whether a pixel in the window is an SHP of its centre pixel, shp_num is the number of SHPs of every pixel.

shp_num_zarr = zarr.open(shp_num,mode='r')
mc.view(shp_num, label='number of SHP')

DS candidates selection¶

ds_can_gix = './ds_processing/gix/ds_can_gix.zarr'
is_ds_can = './ds_processing/gix/ras_is_ds_can.zarr'

nshp_min = 80   # DS candidates: at least this number of SHP (of 189 in the window)
mc.pc_logic_ras(shp_num,ds_can_gix,f'ras>={nshp_min}',chunks=200000)
mc.gix2bool(ds_can_gix, is_ds_can, shape= shp_num_zarr.shape)

ds_can_idx is a int array with shape of (num_of_DS_can, 2). It indicate the position of DS candidate in radar coordinate.

mc.view(is_ds_can, label='DS candidates')

Phase linking¶

For every DS candidate, the coherence matrix is estimated from its SHPs, the phase history is estimated from it by phase linking (EMI), and the temporal coherence is computed, all in one step: the coherence matrices are not saved. With 92 images a coherence matrix has 4186 image pairs (33 KB per point as complex64), tens of GB for all the candidates of this scene.

First, we generate the DS candidate point cloud is_shp by raster chunks (a directory with one zarr per chunk); ds_can_key sorts the points of the chunks back to the order of ds_can_gix:

ds_can_is_shp_dir = './ds_processing/gix/ds_can_is_shp'
ds_can_key = './ds_processing/gix/ds_can_key.zarr'

mc.ras2pc_ras_chunk(ds_can_gix,is_shp,ds_can_is_shp_dir,ds_can_key,chunks=(1000,1000))

Then the phase linking. EMI regularizes the coherence matrices whose magnitude is not positive definite (regularize, default true; decisions 0023 and 0026), which is the case for most candidates when the number of images approaches the number of independent looks of the SHPs. Neighbouring pixels of an SLC are correlated (oversampling, interpolation in the coregistration), so n SHPs are fewer independent looks: the effective number of looks of every point comes from the positions of its SHPs and the speckle correlation of the rslc (decision 0025).

Besides the temporal coherence t_coh, t_coh_w_dir gives the weighted temporal coherence: the image pairs weighted by their squared coherence without the noise bias (decision 0024), so that incoherent pairs do not lower it. It only checks the pairs that are coherent, so n_components_dir gives the number of connected components of the graph whose nodes are the images and whose edges are the coherent image pairs: with 1 every image is linked to every other by coherent pairs, with more the phase between the groups of images is not constrained by any pair, whatever the weighted temporal coherence is (decision 0028). An image pair is coherent when its squared coherence is above the level that a point without any coherent pair reaches with a probability chosen from alpha and the number of images.

ds_can_ph_dir = './ds_processing/gix/ds_can_ph'
ds_can_t_coh_dir = './ds_processing/gix/ds_can_t_coh'
ds_can_t_coh_w_dir = './ds_processing/gix/ds_can_t_coh_w'
ds_can_n_components_dir = './ds_processing/gix/ds_can_n_components'

ds_can_ph = './ds_processing/hix/ds_can/ds_can_ph.zarr'
ds_can_t_coh = './ds_processing/hix/ds_can/ds_can_t_coh.zarr'
ds_can_t_coh_w = './ds_processing/hix/ds_can/ds_can_t_coh_w.zarr'
ds_can_n_components = './ds_processing/hix/ds_can/ds_can_n_components.zarr'

mc.emperical_co_emi_temp_coh_pc(rslc,ds_can_is_shp_dir,ds_can_gix,ds_can_ph_dir,ds_can_t_coh_dir,
                                t_coh_w_dir=ds_can_t_coh_w_dir,n_components_dir=ds_can_n_components_dir,
                                chunks=(1000,1000),cuda=True)

The outputs are point clouds in the raster chunk order. We compute the hilbert index of the DS candidates and the key that sorts them to hilbert order:

ds_can_hix_unsorted = './ds_processing/gix/ds_can_hix.zarr'
ds_can_hix = './ds_processing/hix/ds_can/ds_can_hix.zarr'
gix2hix_key = './ds_processing/gix/gix2hix.zarr'

mc.pc_hix(ds_can_gix, ds_can_hix_unsorted, shape=zarr.open(rslc,mode='r').shape[:2])
mc.pc_sort(
    ds_can_hix_unsorted, ds_can_hix,
    shape=zarr.open(rslc,mode='r').shape[:2],
    pc_in = './ds_processing/gix/ds_can_gix.zarr',
    pc = './ds_processing/hix/ds_can/ds_can_gix.zarr',
    key=gix2hix_key,
)

Then concat the outputs and sort them to hilbert order:

mc.pc_concat(
    [ds_can_ph_dir,ds_can_t_coh_dir,ds_can_t_coh_w_dir,ds_can_n_components_dir],
    [ds_can_ph,ds_can_t_coh,ds_can_t_coh_w,ds_can_n_components],
    key=[ds_can_key,gix2hix_key],
    chunks=200000)

Plot the temporal coherence, the weighted temporal coherence and the connected components:

t_coh = './ds_processing/hix/ds_can/ras_ds_can_t_coh.zarr'
t_coh_w = './ds_processing/hix/ds_can/ras_ds_can_t_coh_w.zarr'
rslc_zarr = zarr.open(rslc,mode='r')
mc.pc2ras(ds_can_hix, ds_can_t_coh,t_coh,rslc_zarr.shape[:2])
mc.pc2ras(ds_can_hix, ds_can_t_coh_w,t_coh_w,rslc_zarr.shape[:2])
# n_components is an integer array, `pc2ras` fills with nan: a pyramid instead
n_components_pyramid = './ds_processing/hix/ds_can/ds_can_n_components_pyramid'
mc.pc_pyramid(ds_can_n_components, n_components_pyramid, yx='./ds_processing/hix/ds_can/ds_can_gix.zarr', ras_resolution=2)

(mc.view(t_coh, label='t_coh', clim=(0, 1))
 + mc.view(t_coh_w, label='weighted t_coh', clim=(0, 1))
 + mc.view(n_components_pyramid, label='connected components of the coherent image pairs'))

DS selection¶

Then we select the DS: all images are linked to each other by coherent image pairs (n_components is 1), and the weighted temporal coherence is high. The weighted temporal coherence is higher than the plain one, so its threshold is higher (0.85 here; the quality of the DS rises with it all the way to 1, so the threshold trades number for quality):

t_coh_w_min = 0.85  # minimum weighted temporal coherence of the phase linked DS
ds_connected_hix = './ds_processing/hix/ds/ds_connected_hix.zarr'
ds_connected_t_coh_w = './ds_processing/hix/ds/ds_connected_t_coh_w.zarr'
ds_hix = './ds_processing/hix/ds/ds_hix.zarr'
mc.pc_logic_pc(ds_can_hix, ds_can_n_components, ds_connected_hix, 'pc_in==1')
mc.pc_select_data(ds_can_hix, ds_connected_hix, pc_in=[ds_can_t_coh_w], pc=[ds_connected_t_coh_w])
mc.pc_logic_pc(ds_connected_hix, ds_connected_t_coh_w, ds_hix, f'(pc_in>={t_coh_w_min})&(pc_in<=1.0)')

Note

If too few points are selected, lower t_coh_w_min or increase alpha of ds_temp_coh.

Visualize phase linked interferograms at refined ds points¶

ds_ph = './ds_processing/hix/ds/ds_ph.zarr'
ds_gix = './ds_processing/hix/ds/ds_gix.zarr'

mc.pc_select_data(
    ds_can_hix,ds_hix,
    pc_in=[ds_can_ph, './ds_processing/hix/ds_can/ds_can_gix.zarr/' ],
    pc=[ds_ph, ds_gix],threads_per_worker=32)

Then plot them:

mc.pc_pyramid(
    './ds_processing/hix/ds/ds_ph.zarr',
    './ds_processing/hix/ds/ds_ph_pyramid',
    yx = './ds_processing/hix/ds/ds_gix.zarr',
    ras_resolution=2,
)

dates = './load_data/meta.toml'
raw_intf = mc.view('./load_data/rslc_pyramid', show='intf_all', dates=dates, label='Raw interferograms')
ds_pl_intf = mc.view('./ds_processing/hix/ds/ds_ph_pyramid', show='intf_all', label='Phase linked interferograms')
raw_intf + ds_pl_intf

Visualize phase linked interferograms at refined ds points in earth coordinates¶

Get the needed data for ds:

ds_ph = './ds_processing/hix/ds/ds_ph.zarr'
ds_e = './ds_processing/hix/ds/ds_e.zarr'
ds_n = './ds_processing/hix/ds/ds_n.zarr'

mc.ras2pc(ds_hix,
          ras=['./load_data/e.zarr/','./load_data/n.zarr/'],
          pc=[ds_e, ds_n],)

mc.pc_pyramid(
    './ds_processing/hix/ds/ds_ph.zarr',
    './ds_processing/hix/ds/ds_ph_geo_pyramid',
    x = ds_e,
    y = ds_n,
    ras_resolution=20,
)

# web mercator coordinates: shown north up over a satellite base map
mc.view('./ds_processing/hix/ds/ds_ph_geo_pyramid', show='intf_all', dates='./load_data/meta.toml',
        label='Phase linked interferograms')