Deformation products¶
Notebook
Converted from nbs/Tutorials/CLI/CampiFlegrei/06_products.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.
In this tutorial, we turn the unwrapped phase of every image from Phase unwrapping into line of sight displacements: relative to a reference area, in millimeters, positive towards the satellite (decision 0040). More details and the expected results: docs/workflows/06_products.md.
Inputs: the hilbert index and the map coordinates of the refined points, the unwrapped phase of every image relative to the first image and the regions of the phase closure correction, and the radar wavelength (in meters) from the metadata.
pc_hix = './pixel_refinement/hix/pc/pc_hix.zarr'
pc_e = './pixel_refinement/hix/pc/pc_e.zarr'
pc_n = './pixel_refinement/hix/pc/pc_n.zarr'
pc_ts = './phase_unwrapping/pc_ts.zarr'
pc_region = './phase_unwrapping/pc_region.zarr'
pc_lon = './products/pc_lon.zarr'
pc_lat = './products/pc_lat.zarr'
ref_polygons = './products/reference.geojson'
pc_ref_mask = './products/pc_ref_mask.zarr'
pc_ts_ref = './products/pc_ts_ref.zarr'
ref_value = './products/ref_value.zarr'
pc_disp = './products/pc_disp.zarr'
wavelength = toml.load('./load_data/meta.toml')['radar_wavelength']
Longitude and latitude of the refined points, to select the reference area:
The reference area should be stable, coherent and in one region of the phase closure correction (points of other regions can be off by whole cycles in some images). The scene is small: every part of it is within about 10 km of Pozzuoli, so no area is far from the uplift; we take a box east of the caldera, in Naples. A polygon drawn on a map works the same way: mc.view(..., polygons=ref_polygons) saves the polygons drawn with the polygon button to the file.
mr.write_polygons(ref_polygons, [[[14.225, 40.82], [14.26, 40.82], [14.26, 40.89], [14.225, 40.89]]])
mc.polygon_mask(ref_polygons, pc_ref_mask, x=pc_lon, y=pc_lat)
The regions of the reference points (-1: points of small regions, corrected one by one):
region = zarr.open(pc_region, mode='r')[:]
ref = zarr.open(pc_ref_mask, mode='r')[:]
ids, counts = np.unique(region[ref], return_counts=True)
dict(zip(ids.tolist(), counts.tolist()))
Spatial reference: the median of the reference points is subtracted from every image, so that it is 0 in every image (method='mean' takes their mean, ref_point one point instead):
Line of sight displacement in millimeters, positive towards the satellite (a shorter range); one radian is 1000 * wavelength / (4 pi) millimeters, 4.41 mm for Sentinel-1:
mc.pc_pyramid(
pc_disp,
'./products/pc_disp_pyramid',
x = pc_e,
y = pc_n,
ras_resolution=20,
method='mean',
)
The displacement of the last image (2021-12-29) on the map; click a point to plot its displacement time series, double click one to make it the reference of the time series. Without show, the slider image chooses the date.