PS Processing

Notebook

Converted from nbs/Tutorials/CLI/CampiFlegrei/02_ps.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 PS processing with the decorrelation CLI interface. The input data for PS processing is prepared in Load Data.

!ls ./load_data

PS candidates selection

Amplitude dispersion index calculation

Here were calculate the amplitude dispersion index (ADI) which is used to identify persistent scatterer candidates.

import zarr
import numpy as np
from matplotlib import pyplot as plt
import colorcet
import moraine.cli as mc
import moraine as mr

rslc = './load_data/rslc.zarr'
adi = './ps_processing/gix/ras_adi.zarr'
rslc_zarr = zarr.open(rslc, mode='r')

logger = mc.get_logger()

mc.amp_disp(rslc,adi)

Then visualize the amplitude dispersion index with mc.view. Small arrays are shown directly; for large ones, build a pyramid first (mc.ras_pyramid) and view the pyramid.

mc.view(adi, label='amplitude dispersion index', clim=(0, 1))

Then we can select PS candidates by setting a threshold to the amplitude dispersion index:

adi_max = 0.4   # PS candidates: amplitude dispersion index below this
ps_can_1_gix = './ps_processing/gix/ps_can_1_gix.zarr'
mc.pc_logic_ras(adi, ps_can_1_gix, f'ras<{adi_max}', chunks=200000)

Noise2Fringe

Here we use Noise2Fringe to select another set of PS candidates.

tnet = mr.TempNet.from_bandwidth(rslc_zarr.shape[-1],bandwidth=1)

tnet.image_pairs

n2f_intf = './ps_processing/gix/ras_n2f_intf.zarr'
mc.n2f(rslc, n2f_intf, tnet.image_pairs, cuda=True)

Then we estimate a temporal coherence:

temp_coh = './ps_processing/gix/ras_temp_coh.zarr'
tnet = mr.TempNet.from_bandwidth(rslc_zarr.shape[-1],1)

mc.temp_coh(n2f_intf, rslc, temp_coh, tnet.image_pairs)

And select another set of PS candidates:

t_coh_min = 0.58   # PS candidates: temporal coherence of Noise2Fringe above this
ps_can_2_gix = './ps_processing/gix/ps_can_2_gix.zarr'
mc.pc_logic_ras(temp_coh, ps_can_2_gix, f'ras>{t_coh_min}', chunks=20000)

Merge two candidates set

ps_can_gix = './ps_processing/gix/ps_can_gix.zarr'
mc.pc_union(ps_can_1_gix, ps_can_2_gix, ps_can_gix)

Then we convert the ps candiates to hillbert coordinates:

mc.pc_hix('./ps_processing/gix/ps_can_gix.zarr/', './ps_processing/gix/ps_can_hix.zarr', shape=rslc_zarr.shape[:2])
ps_can_hix = './ps_processing/hix/ps_can_hix.zarr'
mc.pc_sort('./ps_processing/gix/ps_can_hix.zarr', ps_can_hix)

ps_can_rslc = './ps_processing/hix/ps_can_rslc.zarr'
mc.ras2pc(ps_can_hix,
          [rslc,'./load_data/e.zarr','./load_data/n.zarr',
           './load_data/lon.zarr','./load_data/lat.zarr',
           './ps_processing/gix/ras_adi.zarr'],
          [ps_can_rslc, './ps_processing/hix/ps_can_e.zarr','./ps_processing/hix/ps_can_n.zarr',
           './ps_processing/hix/ps_can_lon.zarr', './ps_processing/hix/ps_can_lat.zarr',
           './ps_processing/hix/ps_can_adi.zarr',],
)

The point cloud data is sorted according to the hillbert index so they close points in the list are closed to each other in map.

#| hide
e = zarr.open('./ps_processing/hix/ps_can_e.zarr',mode='r')[:]
n = zarr.open('./ps_processing/hix/ps_can_n.zarr',mode='r')[:]
plt.scatter(n,e,c=np.arange(e.shape[0]))
plt.show()

visualize the PS candidates

mc.pc_pyramid(
    './ps_processing/hix/ps_can_rslc.zarr',
    './ps_processing/hix/ps_can_rslc_pyramid',
    x = './ps_processing/hix/ps_can_e.zarr',
    y = './ps_processing/hix/ps_can_n.zarr/',
    ras_resolution = 20,
)

mc.view('./ps_processing/hix/ps_can_rslc_pyramid', show='intf_all', dates='./load_data/meta.toml',
        label='PS candidates')