Phase unwrapping

In this tutorial, we demostrate how to unwrap the interferograms on refined measurement points (PS+DS).

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

logger = mc.get_logger()

import holoviews as hv
hv.extension('bokeh')
from bokeh.models import WheelZoomTool
from holoviews import opts
hv.output(widget_location='bottom')

ⓘ

!mkdir phase_unwrapping -p

pc_ph = './pixel_refinement/hix/pc/pc_ph.zarr'
pc_e = './pixel_refinement/hix/pc/pc_e.zarr/'
pc_n = './pixel_refinement/hix/pc/pc_n.zarr/'

pc_unw = './phase_unwrapping/pc_unw.zarr'
pc_ph_zarr = zarr.open(pc_ph,mode='r')
tnet = mr.TempNet.from_bandwidth(pc_ph_zarr.shape[1],bandwidth=1)

mc.gamma_mcf_pt(pc_e, pc_n, pc_ph, pc_unw, tnet.image_pairs)

2025-11-02 17:19:05 - log_args - INFO - running function: gamma_mcf_pt
2025-11-02 17:19:05 - log_args - INFO - fetching args:
2025-11-02 17:19:05 - log_args - INFO - pc_x = './pixel_refinement/hix/pc/pc_e.zarr/'
2025-11-02 17:19:05 - log_args - INFO - pc_y = './pixel_refinement/hix/pc/pc_n.zarr/'
2025-11-02 17:19:05 - log_args - INFO - ph = './pixel_refinement/hix/pc/pc_ph.zarr'
2025-11-02 17:19:05 - log_args - INFO - unw_ph = './phase_unwrapping/pc_unw.zarr'
2025-11-02 17:19:05 - log_args - INFO - image_pairs = array([[ 0,  1],
       [ 1,  2],
       [ 2,  3],
       [ 3,  4],
       [ 4,  5],
       [ 5,  6],
       [ 6,  7],
       [ 7,  8],
       [ 8,  9],
       [ 9, 10],
       [10, 11],
       [11, 12],
       [12, 13],
       [13, 14],
       [14, 15],
       [15, 16]], dtype=int32)
2025-11-02 17:19:05 - log_args - INFO - ref_point = 1
2025-11-02 17:19:05 - log_args - INFO - out_chunks = None
2025-11-02 17:19:05 - log_args - INFO - n_workers = 1
2025-11-02 17:19:05 - log_args - INFO - threads_per_worker = 2
2025-11-02 17:19:05 - log_args - INFO - dask_cluster_arg = {}
2025-11-02 17:19:05 - log_args - INFO - fetching args done.
2025-11-02 17:19:05 - gamma_mcf_pt - INFO - load coordinates
2025-11-02 17:19:05 - gamma_mcf_pt - INFO - Done
2025-11-02 17:19:05 - zarr_info - INFO - ./pixel_refinement/hix/pc/pc_ph.zarr zarray shape, chunks, dtype: (157365, 17), (200000, 1), complex64
2025-11-02 17:19:05 - gamma_mcf_pt - INFO - starting dask local cluster.
2025-11-02 17:19:06 - gamma_mcf_pt - INFO - dask local cluster started.
2025-11-02 17:19:06 - dask_cluster_info - INFO - dask cluster: LocalCluster(dashboard_link='http://127.0.0.1:8787/status', workers=1, threads=2, memory=1.95 TiB)
2025-11-02 17:19:06 - darr_info - INFO - ph dask array shape, chunksize, dtype: (157365, 17), (157365, 1), complex64
2025-11-02 17:19:06 - gamma_mcf_pt - INFO - phase wrapping with mcf.
2025-11-02 17:19:06 - gamma_mcf_pt - INFO - got unwrapped phase.
2025-11-02 17:19:06 - darr_info - INFO - unw_ph dask array shape, chunksize, dtype: (157365, 16), (157365, 1), float32
2025-11-02 17:19:06 - gamma_mcf_pt - INFO - save unw_ph
2025-11-02 17:19:06 - zarr_info - INFO - ./phase_unwrapping/pc_unw.zarr zarray shape, chunks, dtype: (157365, 16), (200000, 1), float32
2025-11-02 17:19:06 - gamma_mcf_pt - INFO - computing graph setted. doing all the computing.
2025-11-02 17:19:14 - gamma_mcf_pt - INFO - computing finished.7.3s

2025-11-02 17:19:14,156 - distributed.worker - ERROR - Failed to communicate with scheduler during heartbeat.
Traceback (most recent call last):
  File "/users/kangl/miniforge3/envs/work2/lib/python3.12/site-packages/distributed/comm/tcp.py", line 226, in read
    frames_nosplit_nbytes_bin = await stream.read_bytes(fmt_size)
                                ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
tornado.iostream.StreamClosedError: Stream is closed

The above exception was the direct cause of the following exception:

Traceback (most recent call last):
  File "/users/kangl/miniforge3/envs/work2/lib/python3.12/site-packages/distributed/worker.py", line 1267, in heartbeat
    response = await retry_operation(
               ^^^^^^^^^^^^^^^^^^^^^^
  File "/users/kangl/miniforge3/envs/work2/lib/python3.12/site-packages/distributed/utils_comm.py", line 416, in retry_operation
    return await retry(
           ^^^^^^^^^^^^
  File "/users/kangl/miniforge3/envs/work2/lib/python3.12/site-packages/distributed/utils_comm.py", line 395, in retry
    return await coro()
           ^^^^^^^^^^^^
  File "/users/kangl/miniforge3/envs/work2/lib/python3.12/site-packages/distributed/core.py", line 1259, in send_recv_from_rpc
    return await send_recv(comm=comm, op=key, **kwargs)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/users/kangl/miniforge3/envs/work2/lib/python3.12/site-packages/distributed/core.py", line 1018, in send_recv
    response = await comm.read(deserializers=deserializers)
               ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/users/kangl/miniforge3/envs/work2/lib/python3.12/site-packages/distributed/comm/tcp.py", line 237, in read
    convert_stream_closed_error(self, e)
  File "/users/kangl/miniforge3/envs/work2/lib/python3.12/site-packages/distributed/comm/tcp.py", line 137, in convert_stream_closed_error
    raise CommClosedError(f"in {obj}: {exc}") from exc
distributed.comm.core.CommClosedError: in <TCP (closed) ConnectionPool.heartbeat_worker local=tcp://127.0.0.1:47038 remote=tcp://127.0.0.1:38921>: Stream is closed

2025-11-02 17:19:14 - gamma_mcf_pt - INFO - dask cluster closed.

mc.pc_pyramid(
    './phase_unwrapping/pc_unw.zarr',
    './phase_unwrapping/pc_unw_geo_pyramid',
    x = pc_e,
    y = pc_n,
    ras_resolution=20,
)

2025-11-02 17:19:17 - log_args - INFO - running function: pc_pyramid
2025-11-02 17:19:17 - log_args - INFO - fetching args:
2025-11-02 17:19:17 - log_args - INFO - pc = './phase_unwrapping/pc_unw.zarr'
2025-11-02 17:19:17 - log_args - INFO - out_dir = './phase_unwrapping/pc_unw_geo_pyramid'
2025-11-02 17:19:17 - log_args - INFO - x = './pixel_refinement/hix/pc/pc_e.zarr/'
2025-11-02 17:19:17 - log_args - INFO - y = './pixel_refinement/hix/pc/pc_n.zarr/'
2025-11-02 17:19:17 - log_args - INFO - yx = None
2025-11-02 17:19:17 - log_args - INFO - ras_resolution = 20
2025-11-02 17:19:17 - log_args - INFO - ras_chunks = (256, 256)
2025-11-02 17:19:17 - log_args - INFO - pc_chunks = 65536
2025-11-02 17:19:17 - log_args - INFO - processes = False
2025-11-02 17:19:17 - log_args - INFO - n_workers = 1
2025-11-02 17:19:17 - log_args - INFO - threads_per_worker = 2
2025-11-02 17:19:17 - log_args - INFO - dask_cluster_arg = {}
2025-11-02 17:19:17 - log_args - INFO - fetching args done.
2025-11-02 17:19:17 - pc_pyramid - INFO - clean out dir
2025-11-02 17:19:17 - zarr_info - INFO - ./phase_unwrapping/pc_unw.zarr zarray shape, chunks, dtype: (157365, 16), (200000, 1), float32
2025-11-02 17:19:17 - pc_pyramid - INFO - rendering point cloud data coordinates:
2025-11-02 17:19:17 - pc_pyramid - INFO - rasterizing point cloud data to grid with bounds: [np.float64(-16498472.188209932), np.float64(8649654.067951033), np.float64(-16470212.188209932), np.float64(8674674.067951033)].
2025-11-02 17:19:18 - zarr_info - INFO - phase_unwrapping/pc_unw_geo_pyramid/x.zarr zarray shape, chunks, dtype: (157365,), (65536,), float64
2025-11-02 17:19:18 - zarr_info - INFO - phase_unwrapping/pc_unw_geo_pyramid/y.zarr zarray shape, chunks, dtype: (157365,), (65536,), float64
2025-11-02 17:19:18 - pc_pyramid - INFO - pc data coordinates rendering ends.
2025-11-02 17:19:18 - zarr_info - INFO - phase_unwrapping/pc_unw_geo_pyramid/idx_0.zarr zarray shape, chunks, dtype: (1252, 1414), (256, 256), int64
2025-11-02 17:19:19 - zarr_info - INFO - phase_unwrapping/pc_unw_geo_pyramid/idx_1.zarr zarray shape, chunks, dtype: (626, 707), (256, 256), int64
2025-11-02 17:19:19 - zarr_info - INFO - phase_unwrapping/pc_unw_geo_pyramid/idx_2.zarr zarray shape, chunks, dtype: (313, 354), (256, 256), int64
2025-11-02 17:19:19 - zarr_info - INFO - phase_unwrapping/pc_unw_geo_pyramid/idx_3.zarr zarray shape, chunks, dtype: (157, 177), (256, 256), int64
2025-11-02 17:19:20 - zarr_info - INFO - phase_unwrapping/pc_unw_geo_pyramid/idx_4.zarr zarray shape, chunks, dtype: (79, 89), (256, 256), int64
2025-11-02 17:19:20 - zarr_info - INFO - phase_unwrapping/pc_unw_geo_pyramid/idx_5.zarr zarray shape, chunks, dtype: (40, 45), (256, 256), int64
2025-11-02 17:19:20 - zarr_info - INFO - phase_unwrapping/pc_unw_geo_pyramid/idx_6.zarr zarray shape, chunks, dtype: (20, 23), (256, 256), int64
2025-11-02 17:19:20 - zarr_info - INFO - phase_unwrapping/pc_unw_geo_pyramid/idx_7.zarr zarray shape, chunks, dtype: (10, 12), (256, 256), int64
2025-11-02 17:19:20 - zarr_info - INFO - phase_unwrapping/pc_unw_geo_pyramid/idx_8.zarr zarray shape, chunks, dtype: (5, 6), (256, 256), int64
2025-11-02 17:19:20 - zarr_info - INFO - phase_unwrapping/pc_unw_geo_pyramid/idx_9.zarr zarray shape, chunks, dtype: (3, 3), (256, 256), int64
2025-11-02 17:19:20 - zarr_info - INFO - phase_unwrapping/pc_unw_geo_pyramid/idx_10.zarr zarray shape, chunks, dtype: (2, 2), (256, 256), int64
2025-11-02 17:19:20 - pc_pyramid - INFO - rasterized idx rendering ends
2025-11-02 17:19:20 - pc_pyramid - INFO - dask local cluster started to render pc data.
2025-11-02 17:19:20 - dask_cluster_info - INFO - dask cluster: LocalCluster(dashboard_link='http://10.211.48.5:8787/status', workers=1, threads=2, memory=1.95 TiB)
2025-11-02 17:19:20 - zarr_info - INFO - phase_unwrapping/pc_unw_geo_pyramid/pc.zarr zarray shape, chunks, dtype: (157365, 16), (65536, 1), float32
2025-11-02 17:19:20 - zarr_info - INFO - phase_unwrapping/pc_unw_geo_pyramid/0.zarr zarray shape, chunks, dtype: (1252, 1414, 16), (256, 256, 1), float32
2025-11-02 17:19:20 - zarr_info - INFO - phase_unwrapping/pc_unw_geo_pyramid/1.zarr zarray shape, chunks, dtype: (626, 707, 16), (256, 256, 1), float32
2025-11-02 17:19:20 - zarr_info - INFO - phase_unwrapping/pc_unw_geo_pyramid/2.zarr zarray shape, chunks, dtype: (313, 354, 16), (256, 256, 1), float32
2025-11-02 17:19:20 - zarr_info - INFO - phase_unwrapping/pc_unw_geo_pyramid/3.zarr zarray shape, chunks, dtype: (157, 177, 16), (256, 256, 1), float32
2025-11-02 17:19:20 - zarr_info - INFO - phase_unwrapping/pc_unw_geo_pyramid/4.zarr zarray shape, chunks, dtype: (79, 89, 16), (256, 256, 1), float32
2025-11-02 17:19:20 - zarr_info - INFO - phase_unwrapping/pc_unw_geo_pyramid/5.zarr zarray shape, chunks, dtype: (40, 45, 16), (256, 256, 1), float32
2025-11-02 17:19:20 - zarr_info - INFO - phase_unwrapping/pc_unw_geo_pyramid/6.zarr zarray shape, chunks, dtype: (20, 23, 16), (256, 256, 1), float32
2025-11-02 17:19:20 - zarr_info - INFO - phase_unwrapping/pc_unw_geo_pyramid/7.zarr zarray shape, chunks, dtype: (10, 12, 16), (256, 256, 1), float32
2025-11-02 17:19:20 - zarr_info - INFO - phase_unwrapping/pc_unw_geo_pyramid/8.zarr zarray shape, chunks, dtype: (5, 6, 16), (256, 256, 1), float32
2025-11-02 17:19:20 - zarr_info - INFO - phase_unwrapping/pc_unw_geo_pyramid/9.zarr zarray shape, chunks, dtype: (3, 3, 16), (256, 256, 1), float32
2025-11-02 17:19:20 - zarr_info - INFO - phase_unwrapping/pc_unw_geo_pyramid/10.zarr zarray shape, chunks, dtype: (2, 2, 16), (256, 256, 1), float32
2025-11-02 17:19:20 - pc_pyramid - INFO - computing graph setted. doing all the computing.
2025-11-02 17:19:22 - pc_pyramid - INFO - computing finished.  1.4s
2025-11-02 17:19:22 - pc_pyramid - INFO - dask cluster closed.

with open('./load_data/meta.toml','r') as f:
    dates = toml.load(f)['dates']
pc_geo_intf_plot = mc.pc_plot('./pixel_refinement/hix/pc/pc_ph_pyramid',post_proc_ras='intf_seq', post_proc_pc='intf_seq',level_increase=0)
pc_geo_intf_plot = pc_geo_intf_plot[0]*pc_geo_intf_plot[1]
pc_geo_intf_plot = pc_geo_intf_plot.redim(
    i=hv.Dimension('i', label='Intf index', range=(0,len(dates)-2), value_format=(lambda i: dates[i]+'_'+dates[i+1])),
    x=hv.Dimension('lon', label='Longitude'),
    y=hv.Dimension('lat',label='Latitude'),
    z=hv.Dimension('Wrapped Phase',range=(-np.pi,np.pi))
)
pc_geo_intf_plot = pc_geo_intf_plot.opts(
    hv.opts.Image(
        cmap='colorwheel',frame_width=500, frame_height=400, colorbar=True,
        default_tools=['pan',WheelZoomTool(zoom_on_axis=False),'save','reset','hover'],
        active_tools=['wheel_zoom'],
        title="Wrapped Phase",
    ),
    hv.opts.Points(
        color='Wrapped Phase', cmap='colorwheel',frame_width=500, frame_height=400, colorbar=True,
        default_tools=['pan',WheelZoomTool(zoom_on_axis=False),'save','reset','hover'],
        active_tools=['wheel_zoom'],
        title="Wrapped Phase",
    ),
)

pc_geo_unw_plot = mc.pc_plot('./phase_unwrapping/pc_unw_geo_pyramid',level_increase=0)
pc_geo_unw_plot = pc_geo_unw_plot[0]*pc_geo_unw_plot[1]
pc_geo_unw_plot = pc_geo_unw_plot.redim(
    i=hv.Dimension('i', label='Intf index', range=(0,len(dates)-2), value_format=(lambda i: dates[i]+'_'+dates[i+1])),
    x=hv.Dimension('lon', label='Longitude'),
    y=hv.Dimension('lat',label='Latitude'),
    z=hv.Dimension('Unwrapped Phase',range=(-10,10))
)
pc_geo_unw_plot = pc_geo_unw_plot.opts(
    hv.opts.Image(
        cmap='colorwheel',frame_width=500, frame_height=400, colorbar=True,
        default_tools=['pan',WheelZoomTool(zoom_on_axis=False),'save','reset','hover'],
        active_tools=['wheel_zoom'],
        #title="Unwrapped Phase",
    ),
    hv.opts.Points(
        color='Unwrapped Phase', cmap='colorwheel',frame_width=500, frame_height=400, colorbar=True,
        default_tools=['pan',WheelZoomTool(zoom_on_axis=False),'save','reset','hover'],
        active_tools=['wheel_zoom'],
        #title="Unwrapped Phase",
    ),
)

hv.element.tiles.EsriImagery()*(pc_geo_intf_plot+pc_geo_unw_plot)