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')

Coherence matrix estimation

The coherence matrix estimation is kind of complicated as it involves both raster and point cloud data. We also want to save the output point cloud data in the hilbert order.

First, we generate DS candidate point cloud is_shp by ras chunks (which is a directory):

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 we estimate the coherence matrix. n_looks_dir also writes the effective number of independent looks of the SHPs of each point: neighbouring pixels of an SLC are correlated (oversampling, interpolation in the coregistration), so n SHPs are fewer independent looks (decision 0025). The weighted temporal coherence and the connectivity of the coherent image pairs below need it:

ds_can_coh_dir = './ds_processing/gix/ds_can_coh'
ds_can_n_looks_dir = './ds_processing/gix/ds_can_n_looks'

mc.emperical_co_pc(rslc,ds_can_is_shp_dir,ds_can_gix,ds_can_coh_dir,n_looks_dir=ds_can_n_looks_dir,cuda=True)

The output of emperical_co_pc is point cloud data in the raster chunk order. So we prepare the key to sort the point cloud data to hilbert order (the hillbert index of DS candidate can be generated by the way):

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 and sort ds_can_coh and ds_can_n_looks to hilbert order:

ds_can_coh = './ds_processing/hix/ds_can/ds_can_coh.zarr'
ds_can_coh_ave = './ds_processing/hix/ds_can/ds_can_coh_ave.zarr'
ds_can_n_looks = './ds_processing/hix/ds_can/ds_can_n_looks.zarr'

chunks = zarr.open(ds_can_gix,mode='r').chunks[0]
mc.pc_concat([ds_can_coh_dir,ds_can_n_looks_dir],[ds_can_coh,ds_can_n_looks],key=[ds_can_key,gix2hix_key],chunks=chunks)

n_points = zarr.open(ds_can_coh,mode='r').shape[0]
mc.data_reduce(ds_can_coh,ds_can_coh_ave,map_func=np.abs,reduce_func=np.sum,post_map_func=lambda x: x/n_points)

ds_can_coh_ave_zarr = zarr.open(ds_can_coh_ave,mode='r')
# average coherence matrix of the DS candidates: x is the secondary image, y the reference image
mc.view(np.abs(mr.uncompress_coh(ds_can_coh_ave_zarr[:])), label='average coherence', clim=(0, 1))

Visualize adpative multi-looked interferograms

Visualize interferograms after adaptive multi-looking and compare with the raw 1 by 1 look interferograms:

mc.pc_pyramid(
    './ds_processing/hix/ds_can/ds_can_coh.zarr',
    './ds_processing/hix/ds_can/ds_can_coh_pyramid',
    yx = './ds_processing/hix/ds_can/ds_can_gix.zarr',
    ras_resolution=2,
)

rslc_zarr = zarr.open(rslc,mode='r')
tnet = mr.TempNet.from_bandwidth(rslc_zarr.shape[2])

dates = './load_data/meta.toml'
raw_intf = mc.view('./load_data/rslc_pyramid', show='intf_all', dates=dates, label='Raw interferograms')
# the phase of the compressed coherence matrix, with sliders for the reference and secondary image
ds_can_intf = mc.view('./ds_processing/hix/ds_can/ds_can_coh_pyramid', show='coh',
                      label='Adaptively multilooked interferograms')
raw_intf + ds_can_intf

Phase linking

To get rid of the closure phase, phase linking is applied here. 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:

ds_can_coh = './ds_processing/hix/ds_can/ds_can_coh.zarr'
ds_can_ph = './ds_processing/hix/ds_can/ds_can_ph.zarr'

mc.emi(ds_can_coh,ds_can_ph)

We can also estimate a temporal coherence for DS. With n_looks it also 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 it also gives n_components, 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, so that such a pure noise point has all images connected with probability alpha (default 1e-3):

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.ds_temp_coh(ds_can_coh,ds_can_ph,ds_can_t_coh,n_looks=ds_can_n_looks,
               t_coh_w=ds_can_t_coh_w,n_components=ds_can_n_components)

Plot it:

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'))

The full coherence matrix may be too big to be saved. It is also possible to do the coherence matrix estimation, phase linking and temporal coherence together; t_coh_w_dir and n_components_dir give the same weighted temporal coherence and connected components as above:

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'

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)
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)

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_coh = './ds_processing/hix/ds/ds_coh.zarr'
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_coh, ds_can_ph, './ds_processing/hix/ds_can/ds_can_gix.zarr/' ],
    pc=[ds_coh,ds_ph, ds_gix],threads_per_worker=32)

Then plot them:

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

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,
)

raw_intf = mc.view('./load_data/rslc_pyramid', show='intf_all', dates=dates, label='Raw interferograms')
ds_intf = mc.view('./ds_processing/hix/ds/ds_coh_pyramid', show='coh', label='Adaptively multilooked interferograms')
ds_pl_intf = mc.view('./ds_processing/hix/ds/ds_ph_pyramid', show='intf_all', label='Phase linked interferograms')
raw_intf + ds_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')