Phase Unwrapping

phase unwrapping
import zarr
import moraine as mr

source

gamma_mcf_pt

 gamma_mcf_pt (pc_x:numpy.ndarray, pc_y:numpy.ndarray, ph:numpy.ndarray,
               ph_weight:numpy.ndarray=None, ref_point:int=0)

A simple wrapper for mcf_pt in GAMMA software, only work if you have access to mcf_pt.

Type Default Details
pc_x ndarray x coordinate, shape of (N,)
pc_y ndarray y coordinate, shape of (N,)
ph ndarray wrapped phase, shape of (N,) or (N,M)
ph_weight ndarray None point weight, shape of (N,) or (N,M), optional
ref_point int 0 reference point, the first point by default
Returns ndarray unwrapped phase, shape of (N,) or (N,M)

source

mcf_pc

 mcf_pc (pc_x:numpy.ndarray, pc_y:numpy.ndarray, ph:numpy.ndarray)

Minimum cost flow phase unwrapping solver. Note that the coordinates (pc_x, pc_y) must be unique.

Type Details
pc_x ndarray x coordinate, shape of (N,)
pc_y ndarray y coordinate, shape of (N,)
ph ndarray wrapped phase, shape of (N,), np.complex64
Returns ndarray unwrapped phase, shape of (N,)

usage:

rslc = zarr.open('../../data/rslc.zarr/',mode='r')[:]

# SHP selection
az_half_win = 5; r_half_win = 5
az_win = 2*az_half_win+1; r_win = 2*r_half_win+1

rmli = np.abs(rslc)**2
p = mr.ks_test(rmli,az_half_win=az_half_win,r_half_win=r_half_win)
is_shp = p < 0.05

# Select DS candidate
shp_num = np.count_nonzero(is_shp,axis=(-2,-1))
is_ds_can = shp_num >= 50

ds_can_is_shp = is_shp[is_ds_can]
ds_can_gix = np.stack(np.where(is_ds_can),axis=-1)
ds_can_ph, ds_can_emi_quality, ds_can_t_coh = mr.emperical_co_emi_temp_coh_pc(rslc,ds_can_gix,ds_can_is_shp,batch_size=1000)

_is_ds_can_refined = (ds_can_emi_quality>=1.0) & (ds_can_emi_quality <1.2) & (ds_can_t_coh > 0.7) & (ds_can_t_coh <= 1.0)

ds_gix = ds_can_gix[_is_ds_can_refined]
ds_ph = ds_can_ph[_is_ds_can_refined]

ds_hix = mr.pc_hix(ds_gix,shape=rslc.shape[:2])
hix_sort_key = mr.pc_sort(ds_hix)

ds_hix = ds_hix[hix_sort_key]
ds_gix = ds_gix[hix_sort_key]
ds_ph = ds_ph[hix_sort_key]

e = zarr.open('../../data/e.zarr/',mode='r')[:]
n = zarr.open('../../data/n.zarr/',mode='r')[:]

ds_e = mr.ras2pc(e,ds_gix)
ds_n = mr.ras2pc(n,ds_gix)

ds_ph = ds_ph*ds_ph[:,0:1].conj() # set the first image as secondary

Unwrap one interferogram:

ds_ph_ = ds_ph[:,14]*ds_ph[:,9].conj()
unw = mcf_pc(ds_gix[:,1],ds_gix[:,0], ds_ph_)
unw_gm = gamma_mcf_pt(ds_gix[:,1],ds_gix[:,0], ds_ph_) # if you have access to gamma

# This may produce error as ds_e and ds_n may not be unique
# unw = mcf_pc(ds_e, ds_n, ds_ph_)
# unw_gm = gamma_mcf_pt(ds_e, ds_n, ds_ph_)
(459472,) (918706, 3)
import holoviews as hv
import numpy as np
from holoviews import opts
from bokeh.models import WheelZoomTool
hv.extension('bokeh')
unw_plot = mr.pc_plot(unw,ds_n, ds_e)
unw_plot = unw_plot[0]*unw_plot[1]
ph_plot = mr.pc_plot(np.angle(ds_ph_), ds_n, ds_e)
ph_plot = ph_plot[0]*ph_plot[1]
unw_plot = unw_plot.redim(
    x=hv.Dimension('lon', label='Longitude'), y=hv.Dimension('lat',label='Latitude'),
    z=hv.Dimension('unw',label='Unwrapped phase')
)
unw_plot.opts(
    opts.Image(cmap='colorwheel',width=600, height=400, colorbar=True,
               default_tools=['pan',WheelZoomTool(zoom_on_axis=False),'save','reset','hover'],
               active_tools=['wheel_zoom']
              ),
    opts.Points(color='unw', cmap='colorwheel',width=600, height=400, colorbar=True,
                default_tools=['pan',WheelZoomTool(zoom_on_axis=False),'save','reset','hover'],
                active_tools=['wheel_zoom']
               )
)
ph_plot = ph_plot.redim(
    x=hv.Dimension('lon', label='Longitude'), y=hv.Dimension('lat',label='Latitude'),
    z=hv.Dimension('ph',label='Wrapped phase')
)
ph_plot.opts(
    opts.Image(cmap='colorwheel',width=600, height=400, colorbar=True,
               default_tools=['pan',WheelZoomTool(zoom_on_axis=False),'save','reset','hover'],
               active_tools=['wheel_zoom']
              ),
    opts.Points(color='ph', cmap='colorwheel',width=600, height=400, colorbar=True,
                default_tools=['pan',WheelZoomTool(zoom_on_axis=False),'save','reset','hover'],
                active_tools=['wheel_zoom']
               )
)