From 47c9bccd7d805653bf30f0e697f340e1e86c0ac0 Mon Sep 17 00:00:00 2001 From: Yujie Zheng Date: Tue, 26 Jul 2022 17:04:02 -0700 Subject: [PATCH] non-closure phase bias estimation (#765) This PR adds the bias time-series estimation introduced by the non-closure phase, as described in Zheng et al. (2022, TGRS), via `--action quick_estimate / estimate` options. Detailed changes are as below: + `utils.isce_utils.py` - add `guassian_kernel()` and `convolve()` for the complex interferogram filtering. - add `estimate_coherence()` and `unwrap_snaphu()` based on `isce2/topsStack` to re-unwrap the closure phase interferogram, to mitigate the PU error during closure phase bias estimation, including notes on SNAPHU from Piyush. + `objects.stack.ifgramStack`: - move `split2boxes()` from `ifgram_inversion.py` to here for more convenient re-use, add `dim0_size` option for more flexibility - add `get_closure_phase_index()` to return the indices of interferograms that forms the given connection level of closure loop - add `get_sequential_closure_phase()` for the 3D CP calculation and its temporal average. + `ifgram_inversion.estimate_timeseries()`: support `inv_quality_name = no` to turn OFF the calculation. + closure_phase_bias.py: - add `--action mask` for average CP and mask calculation, by refactoring the existing script. - add `compute_unwrap_closure_phase()` to compute the wrapped sequential closure phase, unwrap them using isce2 modules, and then calculate the cumulative unwrapped seq closure phase time-series via `cum_seq_unw_closure_phase_timeseries()`, with parallel support and auto-skip. - add `estimate_wratio()` and `estimate_wratio_all()` for the `W_r` estimation, a.k.a., the weighted closure phase history. - add `bandwidth2num_ifgram()` - add `estimate_bias_timeseries_approx()` for `--action quick_estimate` - add `estimate_bias_timeseries()` for `--action estimate` with parallel support using dask. - 2X speedup of `estimate_bias_timeseries()` via `ifgram_inversion.estimate_timeseries()` with `scipy.linalg`, instead of `np.linalg.pinv`, and auto-skip the zero/nan values in part of the interferograms among the stack. - 2X speed up via common design matrix operation, same as in `ifgram_inversion.py` Co-authored-by: Zhang Yunjun --- mintpy/closure_phase_bias.py | 1243 ++++++++++++++++++++++++++++++---- mintpy/ifgram_inversion.py | 127 ++-- mintpy/objects/cluster.py | 2 +- mintpy/objects/stack.py | 154 ++++- mintpy/utils/isce_utils.py | 220 +++++- 5 files changed, 1543 insertions(+), 203 deletions(-) diff --git a/mintpy/closure_phase_bias.py b/mintpy/closure_phase_bias.py index 83b420d5f..d92217eae 100755 --- a/mintpy/closure_phase_bias.py +++ b/mintpy/closure_phase_bias.py @@ -2,164 +2,1173 @@ ############################################################ # Program is part of MintPy # # Copyright (c) 2013, Zhang Yunjun, Heresh Fattahi # -# Author: Yujie Zheng, Feb 2022 # +# Author: Yujie Zheng, Zhang Yunjun, Feb 2022 # ############################################################ -# Compute average con-nl closure phase and output mask identifying areas suseptible to closure phase errors. +# Recommend import: +# from mintpy import closure_phase_bias as cpbias + import os import sys +import time import argparse import numpy as np +import glob +from datetime import datetime as dt -from mintpy.objects import ifgramStack -from mintpy.utils import readfile, writefile -from mintpy import ifgram_inversion as ifginv +from mintpy.objects import ifgramStack, cluster +from mintpy.utils import ( + arg_group, + ptime, + readfile, + writefile, + isce_utils, + utils as ut, +) +from mintpy.ifgram_inversion import estimate_timeseries ################################################################################ REFERENCE = """reference: - Zheng, Y., et al., (2022) On Closure Phase and Systematic Bias in Multilooked SAR Interferometry, IEEE TGRS, under review (minor revision) + Y. Zheng, H. Fattahi, P. Agram, M. Simons and P. Rosen, (2022). On Closure Phase + and Systematic Bias in Multi-looked SAR Interferometry, in IEEE Trans. Geosci. + Remote Sens., doi:10.1109/TGRS.2022.3167648. """ + EXAMPLE = """example: - closure_phase_bias.py -i inputs/ifgramStack.h5 --nl 20 --numsigma 2.5 + # Note: ONLY sequential network is supported in this implementation. + + # create mask for areas suseptible to biases + closure_phase_bias.py -i inputs/ifgramStack.h5 --nl 5 -a mask + closure_phase_bias.py -i inputs/ifgramStack.h5 --nl 20 -a mask --num-sigma 2.5 + + # estimate non-closure phase bias time-series [quick and approximate solution] + closure_phase_bias.py -i inputs/ifgramStack.h5 --nl 5 --bw 3 -a quick_estimate --num-worker 6 + + # estimate non-closure phase bias time-series + closure_phase_bias.py -i inputs/ifgramStack.h5 --nl 5 --bw 3 -a estimate --num-worker 6 -c local + closure_phase_bias.py -i inputs/ifgramStack.h5 --nl 20 --bw 10 -a estimate --num-worker 6 -c local """ def create_parser(): - parser = argparse.ArgumentParser(description = 'Create an indication map for closure phase bias.') - parser.add_argument('-i','--ifgramstack',type = str, dest = 'ifgram_stack',help = 'interferogram stack file that contains the unwrapped phases') - parser.add_argument('--nl', dest = 'nl', type = int, default = 20, help = 'connection level that we are correcting to (or consider as no bias)') - parser.add_argument('--numsigma',dest = 'numsigma', type = float, default = 3, help = 'Threshold for phase (number of sigmas,0-infty), default to be 3 sigma of a Gaussian distribution (assumed distribution for the cumulative closure phase) with sigma = pi/sqrt(3*num_cp)') - parser.add_argument('--epi',dest = 'episilon', type = float, default = 0.3, help = 'Threshold for amplitude (0-1), default 0.3') - parser.add_argument('--maxMemory', dest = 'max_memory', type = float, default = 8, help = 'max memory to use in GB') - parser.add_argument('-o', dest = 'outdir', type = str, default = '.', help = 'output file directory') + parser = argparse.ArgumentParser(description = 'Phase non-closure related biases correction', + formatter_class = argparse.RawTextHelpFormatter, + epilog=REFERENCE+'\n'+EXAMPLE) + + parser.add_argument('-i','--ifgramstack', type=str, dest='stack_file', + help='interferogram stack file that contains the unwrapped phases') + parser.add_argument('--wm','--water-mask', dest='water_mask_file', + help='Water mask to skip pixels on water body.\n' + 'Default: waterMask.h5 if exists, otherwise None.') + + # bandwidth and bias free connection level + parser.add_argument('--nl','--conn-level', dest='nl', type=int, default=20, + help='connection level that we are correcting to (or consider as no bias)\n' + '(default: %(default)s)') + parser.add_argument('--bw', dest='bw', type=int, default=10, + help='bandwidth of time-series analysis that you want to correct') + + parser.add_argument('-a','--action', dest='action', type=str, default='mask', + choices={'mask', 'quick_estimate', 'estimate'}, + help='action to take (default: %(default)s):\n' + 'mask - create a mask of areas susceptible to closure phase errors\n' + 'quick_estimate - quick and approximate estimation on how bias decays with time\n' + ' output sequential closure phase files\n' + 'estimate - estimate how bias decays with time\n' + ' processed for each pixel on a pixel by pixel basis [slow]') + + # mask configuration + mask = parser.add_argument_group('Mask', 'Configuration for closure phase bias mask') + mask.add_argument('--num-sigma', dest='num_sigma', type=float, default=3, + help='Threashold for phase, in number of sigmas (default: %(default)s).\n' + 'Assuming a Gaussian distribution for the cumulative closure phase' + ' with sigma = pi / sqrt(3*num_cp)') + mask.add_argument('--eps','--epsilon', dest='epsilon', type=float, default=0.3, + help='Threashold for the normalized amplitude in [0-1] (default: %(default)s).') + + # compute + parser = arg_group.add_parallel_argument(parser) + parser = arg_group.add_memory_argument(parser) + + # output + parser.add_argument('-o', dest='outdir', type=str, default='./', help='output file directory') + return parser + def cmd_line_parse(iargs=None): parser = create_parser() - inps = parser.parse_args(args = iargs) + inps = parser.parse_args(args=iargs) + + # --water-mask option + if not inps.water_mask_file and os.path.isfile('./waterMask.h5'): + inps.water_mask_file = os.path.abspath('./waterMask.h5') + return inps -# Obtain sum of consecutive complex sequential closure phase of connection n -def cum_seq_closurePhase(SLC_list, date12_list_all, ifgram_stack, ref_phase, n, box): + +################################# Mask ####################################### +def calc_closure_phase_mask(stack_file, bias_free_conn, num_sigma=3, threshold_amp=0.3, + outdir='./', max_memory=4.0): + """Calculate a mask for areas suseptible to biases, based on the average closure phase tau. + + Equation: tau = 1 / K * Sigma_{k=1}^K (np.exp(j * Phi_k^{nl})) + where K is the number of closure phase for connection nl, Phi_k^{nl} is the k-th sequential + closure phase for connection nl, as defined in equation (21). + Reference: Section VI in Zheng et al. (2022, TGRS). + + Parameters: stack_file - str, path for ifgramStack.h5 file + bias_free_conn - int, connection level at which we assume is bias-free + num_sigma - float, number of sigmas for computing phase threshold + threshold_amp - float, threshold of ampliutde of the cumulative sequential closure phase + outdir - str, directory of output files + max_mermory - float, maximum memory in GB for each patch processed + Returns: mask - 2D np.ndarray of size (length, width) in boolean, 0 for areas suseptible to biases. + Saved to file: maskClosurePhase.h5 + avg_cp - 2D np.ndarray of size (length, width) in complex64, average cum. seq. closure phase + Saved to file: avgCpxClosurePhase.h5 """ - Input parameters: - SLC_list : list of SLC dates - date12_list_all: date12 of all the interferograms stored in the ifgramstack file - ifgram_stack: stack file - refphase : reference phase - n : connection level of the closure phase - box : bounding box for the patch + + # basic info + stack_obj = ifgramStack(stack_file) + stack_obj.open(print_msg=False) + meta = dict(stack_obj.metadata) + length, width = stack_obj.length, stack_obj.width + date_list = stack_obj.get_date_list(dropIfgram=True) + num_cp = stack_obj.get_closure_phase_index(bias_free_conn).shape[0] + + ## What is a good thredshold? + # Assume that it's pure noise so that the phase is uniform distributed from -pi to pi. + # The standard deviation of phase in each loop is: + # pi/sqrt(3) + # (technically should be smaller because when forming loops there should be a reduction in phase variance) + # The standard deviation of phase in cumulative wrapped closure phase is: + # pi/sqrt(3)/sqrt(num_cp) -- again another simplification assuming no correlation. + # We use 3\delta as threshold -- 99.7% confidence + threshold_pha = np.pi / np.sqrt(3 * num_cp) * num_sigma + + # key info + print('\n'+'-'*80) + print('calculating the mask to flag areas suseptible to non-closure-phase related biases (as zero) ...') + print(f'number of valid acquisitions: {len(date_list)} ({date_list[0]} - {date_list[-1]})') + print(f'average complex closure phase threshold in amplitude/correlation: {threshold_amp}') + print(f'average complex closure phase threshold in phase: {num_sigma} sigma ({threshold_pha:.1f} rad)') + + # calculate the average complex closure phase + # process block-by-block to save memory + print('\ncalculating the average complex closure phase') + print(f'length / width: {length} / {width}') + box_list, num_box = stack_obj.split2boxes(max_memory=max_memory, dim0_size=stack_obj.numIfgram+num_cp*2) + + avg_cp = np.zeros([length,width], dtype=np.complex64) + for i, box in enumerate(box_list): + if num_box > 1: + print('\n------- processing patch {} out of {} --------------'.format(i+1, num_box)) + print(f'box: {box}') + print('box width: {}'.format(box[2] - box[0])) + print('box length: {}'.format(box[3] - box[1])) + + avg_cp[box[1]:box[3], box[0]:box[2]], num_cp = stack_obj.get_sequential_closure_phase( + box=box, + conn=bias_free_conn, + post_proc='mean', + )[1:] + + # mask out no-data pixels + geom_file = 'geometryGeo.h5' if 'Y_FIRST' in meta.keys() else 'geometryRadar.h5' + geom_file = os.path.join(os.path.dirname(stack_file), geom_file) + if os.path.isfile(geom_file): + geom_ds_names = readfile.get_dataset_list(geom_file) + ds_names = [x for x in ['incidenceAngle', 'waterMask'] if x in geom_ds_names] + if len(ds_names) > 0: + print(f'mask out pixels with no-data-value (zero {ds_names[0]} from file: {os.path.basename(geom_file)})') + no_data_mask = readfile.read(geom_file, datasetName=ds_names[0])[0] == 0 + avg_cp[no_data_mask] = np.nan + + # create mask + print('\ncreate mask for areas suseptible to non-closure phase biases') + mask = np.ones([length,width], dtype=bool) + + # mask areas with potential bias + print(f'set pixels with average complex closure phase angle > {num_sigma} sigma ({threshold_pha:.1f} rad) to 0.') + mask[np.abs(np.angle(avg_cp)) > threshold_pha] = 0 + + # unmask areas with low correlation + # where it's hard to know wheter there is bias or not + print(f'set pixels with average complex closure phase amplitude (correlation) < {threshold_amp} to 1.') + mask[np.abs(np.abs(avg_cp) < threshold_amp)] = 1 + + # write file 1 - mask + mask_file = os.path.join(outdir, 'maskClosurePhase.h5') + meta['FILE_TYPE'] = 'mask' + meta['DATA_TYPE'] = 'bool' + writefile.write(mask, out_file=mask_file, metadata=meta) + + # write file 2 - average closure phase + avg_cp_file = os.path.join(outdir, 'avgCpxClosurePhase.h5') + meta['FILE_TYPE'] = 'mask' + meta['DATA_TYPE'] = 'float32' + ds_dict = { + 'amplitude' : np.abs(avg_cp), + 'phase' : np.angle(avg_cp), + } + writefile.write(ds_dict, out_file=avg_cp_file, metadata=meta) + + return mask, avg_cp + + +################################################################################ +def unwrap_closure_phase(int_file, cor_file, unw_file): + """Unwrap the input wrapped sequential closure phase interferogram. """ - cp_idx = [] - NSLC = len(SLC_list) - for i in range(NSLC-n): - ifgram = [] - flag = True - for j in range(n): - ifgram.append('{}_{}'.format(SLC_list[i+j],SLC_list[i+j+1])) - ifgram.append('{}_{}'.format(SLC_list[i],SLC_list[i+n])) - for ifgram_name in ifgram: - if ifgram_name not in date12_list_all: - flag = False # if missing an interferogram, we won't make the corresponding closure phase - if flag: - cp_idx.append([date12_list_all.index(ifgram[j]) for j in range(n+1)]) - - cp_idx = np.array(cp_idx, np.int16) - cp_idx = np.unique(cp_idx, axis = 0) - - num_cp = len(cp_idx) - print('Number of closure measurements expected, ', len(SLC_list)-n) - print('Number of closure measurements found, ', num_cp) - - if num_cp <1: - print('No closure phase measurements found, abort') - raise Exception("No triplets found!") - - box_width = box[2] - box[0] - box_length = box[3] - box[1] - phase = readfile.read(ifgram_stack, box=box,print_msg=False)[0] - cum_cp = np.zeros((box_length, box_width), np.complex64) - for i in range(num_cp): - cp0_w = np.zeros ((box_length, box_width), np.float32) - for j in range(n): - idx = cp_idx[i,j] - cp0_w = cp0_w + phase[idx,:,:] - ref_phase[idx] - idx = cp_idx[i,n] - cp0_w = cp0_w - (phase[idx,:,:]-ref_phase[idx]) - cum_cp = cum_cp + (np.exp(1j*cp0_w)) - - # cum_cp = np.angle(cum_cp) - return cum_cp, num_cp - -def main(iargs = None): - inps = cmd_line_parse(iargs) - stack_obj = ifgramStack(inps.ifgram_stack) + + if os.path.isfile(cor_file) and os.path.isfile(unw_file): + print(f'unwrapped interferogram file exists: {unw_file}, skip re-unwrapping.') + else: + isce_utils.estimate_coherence(int_file, cor_file) + isce_utils.unwrap_snaphu(int_file, cor_file, unw_file) + + return unw_file + + +def cum_seq_unw_closure_phase_timeseries(conn, conn_dir, date_list, meta): + '''Output cumulative conn-n sequential closure phase in time-series format, + which is the weighted phase history of the temporally inconsistent process (f^n / n_l). + + Reference: Equation (25) and (28) in Zheng et al. (2022, TGRS). + + Parameters: conn - int, connection level of closure phases + conn_dir - str, path of sequential closure phases file for connection - n + date_list - list of str, SLC dates + meta - dict, metadata of ifgramStack.h5 + Returns: bias_ts - 3D np.ndarray in size of (num_ifgram, length, width) in float32, + cumulative sequential closure phase time series, + saved to file: cumSeqClosurePhase.h5 + common_mask - 2D np.ndarray in size of (length, width) in bool, + mask based on common connected components, + saved to file: maskConnComp.h5 + ''' + + # output file + cum_cp_file = os.path.join(conn_dir, 'cumSeqClosurePhase.h5') + mask_file = os.path.join(conn_dir, 'maskConnComp.h5') + + # update mode checking + if os.path.isfile(cum_cp_file) and os.path.isfile(mask_file): + msg = 'cumulative seq closure phase time series and mask files exist, skip re-generating.' + msg += f'\n{cum_cp_file}\n{mask_file}' + print(msg) + + # read data + bias_ts = readfile.read(cum_cp_file)[0] + common_mask = readfile.read(mask_file)[0] + + return bias_ts, common_mask + + # basic info + length, width = int(meta['LENGTH']), int(meta['WIDTH']) + ref_y, ref_x = int(meta['REF_Y']), int(meta['REF_X']) + + unw_files = sorted(glob.glob(os.path.join(conn_dir, '*.unw'))) + num_file = len(unw_files) + + print('calculate the cumulative seq closure phase time series ...') + cp_phase = np.zeros((num_file, length, width), dtype=np.float32) + mask = np.zeros((num_file, length, width), dtype=np.float32) + + prog_bar = ptime.progressBar(maxValue=num_file) + for i, unw_file in enumerate(unw_files): + prog_bar.update(i+1, suffix=f'{i+1}/{num_file} {os.path.basename(unw_file)}') + + unw = readfile.read(unw_file, datasetName='phase')[0] + unw -= unw[ref_y, ref_x] + cp_phase[i] = unw + + conn_comp_file = unw_file + '.conncomp' + conn_comp = readfile.read(conn_comp_file)[0] + mask[i] = np.where(conn_comp >= 1, 1, np.nan) + + prog_bar.close() + + # compute cumulative sequential closure phase - f^n + # equation (25) in Zheng et al. (2022, TGRS) + num_date = len(date_list) + bias_ts = np.zeros((num_date, length, width), dtype=np.float32) + bias_ts[1:num_date-conn+1, :, :] = np.cumsum(cp_phase, 0) + for i in range(num_date-conn+1, num_date): + bias_ts[i] = (i - num_date + conn) * cp_phase[-1] + bias_ts[num_date - conn] + + # equation (28) + bias_ts /= conn + + # write bias time series to HDF5 file + ds_dict = { + 'timeseries' : [np.float32, (num_date, length, width), bias_ts], + 'date' : [np.dtype('S8'), (num_date,), np.array(date_list, np.string_)], + } + meta['FILE_TYPE'] = 'timeseries' + writefile.layout_hdf5(cum_cp_file, ds_dict, metadata=meta) + + # write mask to HDF5 file + common_mask = np.where(np.isnan(np.sum(mask,0)), False, True) + meta['FILE_TYPE'] = 'mask' + writefile.write(common_mask, out_file=mask_file, metadata=meta) + + return bias_ts, common_mask + + +def compute_unwrap_closure_phase(stack_file, conn, num_worker=1, outdir='./', max_memory=4.0): + '''Compute the following phase stack & time-series of connection-conn: + + + wrapped seq closure phase stack + + unwrapped seq closure phase stack + + cumulative unwrapped seq closure phase time-series + at directory: outdir/closurePhase/conn{conn} + + Parameters: stack_file - str, path for ifgramStack.h5 + conn - int, connection level + max_mermory - float, maximum memory in GB for each patch processed + outdir - str, path for output files + ''' + # output directory + conn_dir = os.path.join(outdir, f'closurePhase/conn{conn}') + os.makedirs(conn_dir, exist_ok=True) + + # update mode checking + cum_cp_file = os.path.join(conn_dir, 'cumSeqClosurePhase.h5') + if os.path.isfile(cum_cp_file): + print(f'cumulative unwrapped seq closure phase time-series exists at: {cum_cp_file}, skip re-generating.') + return + + print('-'*60) + print('step 1/3: calculate and filter the wrapped sequential closure phase stack ...') + + # basic info + stack_obj = ifgramStack(stack_file) stack_obj.open() length, width = stack_obj.length, stack_obj.width - date12_list = stack_obj.get_date12_list(dropIfgram=True) - date12_list_all = stack_obj.get_date12_list(dropIfgram=False) - print('scene length, width', length, width) - ref_phase = stack_obj.get_reference_phase(unwDatasetName = 'unwrapPhase') - inps.length = length - inps.width = width - # retrieve the list of SLC dates from ifgramStack.h5 - ifgram0 = date12_list[0] - date1, date2 = ifgram0.split('_') - SLC_list = [date1, date2] - for ifgram in date12_list: - date1, date2 = ifgram.split('_') - if date1 not in SLC_list: - SLC_list.append(date1) - if date2 not in SLC_list: - SLC_list.append(date2) - SLC_list.sort() - print('number of SLC found : ', len(SLC_list)) - print('first SLC: ', SLC_list[0]) - print('last SLC: ', SLC_list[-1]) + meta = dict(stack_obj.metadata) + print(f'scene size: {length} x {width}') + date_list = stack_obj.get_date_list(dropIfgram=True) + num_date = len(date_list) + print(f'number of acquisitions found: {num_date}') + print(f'start / end date: {date_list[0]} / {date_list[-1]}') + # number of expected closure phase + num_cp = num_date - conn + num_digit = len(str(num_cp)) - # split igram_file into blocks to save memory - box_list, num_box = ifginv.split2boxes(inps.ifgram_stack,inps.max_memory) - closurephase = np.zeros([length,width],np.complex64) - #process block-by-block - for i, box in enumerate(box_list): - box_width = box[2] - box[0] - box_length = box[3] - box[1] + ## default output binary filenames + fbases = [os.path.join(conn_dir, f'filt_{x+1:0{num_digit}}') for x in range(num_cp)] + int_files = [f'{x}.int' for x in fbases] + cor_files = [f'{x}.cor' for x in fbases] + unw_files = [f'{x}.unw' for x in fbases] + + if all(os.path.isfile(x) for x in int_files): + print('ALL the filtered closure phase file exist, skip re-generation.') + + else: + # process block-by-block + # split igram_file into blocks to save memory + box_list, num_box = stack_obj.split2boxes(max_memory=max_memory, + dim0_size=stack_obj.numIfgram+num_cp*2) + closure_phase = np.zeros([num_cp, length, width],np.float32) + for i, box in enumerate(box_list): print(box) if num_box > 1: print('\n------- processing patch {} out of {} --------------'.format(i+1, num_box)) - print('box width: {}'.format(box_width)) - print('box length: {}'.format(box_length)) + print('box length: {}'.format(box[3] - box[1])) + print('box width : {}'.format(box[2] - box[0])) - closurephase[box[1]:box[3],box[0]:box[2]], numcp = cum_seq_closurePhase(SLC_list, date12_list_all, inps.ifgram_stack, ref_phase,inps.nl,box) + closure_phase[:, box[1]:box[3], box[0]:box[2]] = stack_obj.get_sequential_closure_phase( + box=box, + conn=conn, + )[0] + ## filter the closure phase and re-unwrap + print('-'*80) + print('filter the wrapped closure phase stack with a Gaussian kernel of 5 x 5 ...') + print(f'number of wrapped closure phase: {num_cp}') + kernel = isce_utils.gaussian_kernel(5, 5, 1, 1) + for i, int_file in enumerate(int_files): + if not os.path.isfile(int_file): + # filter the closure phase interferogram + closure_phase_filt = isce_utils.convolve( + data=np.exp(1j*closure_phase[i]), + kernel=kernel, + ).astype(np.complex64) - # What is a good thredshold? - # Assume that it's pure noise so that the phase is uniform distributed from -pi to pi. - # The standard deviation of phase in each loop is pi/sqrt(3) (technically should be smaller because when forming loops there should be a reduction in phase variance) - # The standard deviation of phase in cumulative wrapped closure phase is pi/sqrt(3)/sqrt(num_cp) -- again another simplification assuming no correlation. - # We use 3\delta as default threshold -- 99.7% confidence + # write to binary file in isce2 format + print(f'write file: {int_file}') + with open(int_file, mode='wb') as fid: + closure_phase_filt.tofile(fid) + + # write metadata in isce2/roipac format + meta['FILE_TYPE'] = '.int' + meta['INTERLEAVE'] = 'BIP' + meta['DATA_TYPE'] = 'complex64' + meta['BANDS'] = 1 + writefile.write_isce_xml(meta, int_file) + writefile.write_roipac_rsc(meta, int_file+'.rsc') + del closure_phase + + print('-'*60) + print('step 2/3: unwrap the filtered wrapped closure phase stack ...') + print(f'number of closure phase: {num_cp}') + if all(os.path.isfile(x) for x in unw_files): + print('ALL the unwrapped closure phase file exist, skip re-generation.') + + else: + num_core, run_parallel, Parallel, delayed = ut.check_parallel( + num_cp, + print_msg=False, + maxParallelNum=num_worker) + + if run_parallel and num_core > 1: + print(f'parallel processing using {num_core} cores') + Parallel(n_jobs=num_core)(delayed(unwrap_closure_phase)(x, y, z) + for x, y, z in zip(int_files, cor_files, unw_files)) + + else: + for x, y, z in zip(int_files, cor_files, unw_files): + unwrap_closure_phase(x, y, z) + + ## calc the cumulativev unwrapped closure phase time-series + print('-'*60) + print('step 3/3: calculate the unwrapped cumulative sequential closure phase time-series ...') + cum_seq_unw_closure_phase_timeseries(conn, conn_dir, date_list, meta) + + return + + +################################################################################ +def read_cum_seq_closure_phase4conn(conn, outdir='./', box=None, print_msg=False): + '''Read cumulative sequential closure phase from individual closure phase directory. + + Reference: Eq. (25) in Zheng et al. (2022). + + Parameters: conn - integer, connection level of sequential closure phases + outdir - string, directory of conn{n}_cumSeqClosurePhase.h5 + box - list in size of (4,) in integer, coordinates of bounding box + Returns: bias_ts - 3D array in size of (num_date, box_lengh, box_wid) in float, + cumulative sequential closure phases + ''' + cum_cp_file = os.path.join(outdir, f'closurePhase/conn{conn}/cumSeqClosurePhase.h5') + if print_msg: + print(f'read timeseries from file: {cum_cp_file}') + + bias_ts = readfile.read(cum_cp_file, box=box, print_msg=False)[0] + return bias_ts + + +def estimate_wratio(tbase, conn, bias_free_conn, wvl, box, outdir='./', mask=False): + '''Estimate W_r & velocity bias for the given connection level. + + W_r is the M x M diagonal matrix with W_r(ii) = w(i*delta_t) / w(delta_t), i = 1,2,...,M. + This is defined in Equation (20), to be used for bias estimation; and can be calculated from + the weighted phase history (cumulative seq closure phase) based on Equation (29). + + Parameters: tbase - list(float) or array, in size of (num_date), time in accumulated years + conn - integer, connection-level + bias_free_conn - integer, minimum connection-level that we think is bias-free + wvl - float, wavelength + box - list in size of (4,) in integer, coordinates of bounding box + outdir - str, the working directory + mask - bool, whether to mask out areas with average bias velocity less than 1 mm/year + Returns: wratio_connN - 2D np.ndarray of size (box_len, box_wid) in float32, W_r of the input connection level + vel_bias_connN - 2D np.ndarray of size (box_len, box_wid) in float32, velocity bias of the input connection level + ''' + delta_t = tbase[-1] - tbase[0] + phase2range = -1 * wvl / (4 * np.pi) + + # cum/vel_bias at bias free connection level + cum_bias_connF = read_cum_seq_closure_phase4conn(bias_free_conn, outdir, box)[-1,:,:] + vel_bias_connF = cum_bias_connF / delta_t * phase2range + + # calc wratio at input connection level + box_wid = box[2] - box[0] + box_len = box[3] - box[1] + wratio_connN = np.ones([box_len, box_wid], dtype=np.float32) + + if conn > 1: + cum_bias_connN = read_cum_seq_closure_phase4conn(conn, outdir, box)[-1,:,:] + # skip invalid pixels + flag = np.multiply(~np.isnan(cum_bias_connF), cum_bias_connF != 0) + # Equation (29) + wratio_connN[flag] = 1 - cum_bias_connN[flag] / cum_bias_connF[flag] + + # bound within [0, 1] + wratio_connN[wratio_connN > 1] = 1 + wratio_connN[wratio_connN < 0] = 0 - if inps.numsigma: - threshold_cp = np.pi/np.sqrt(3)/np.sqrt(numcp)*inps.numsigma else: - threshold_cp = np.pi/np.sqrt(3)/np.sqrt(numcp)*3 # 3/sigma, 99.7% confidence + cum_bias_connN = None + + # vel_bias at input connection level + vel_bias_connN = np.multiply(wratio_connN, vel_bias_connF) + if mask: + # if average velocity smaller than 1 mm/year (hardcoded here), mask out for better visual + # this option is only turned on while outputing wratio.h5 file. + wratio_connN[abs(vel_bias_connF) < 0.001] = np.nan + + # debug mode + debug_mode = False + if debug_mode: + from matplotlib import pyplot as plt + + data_list = [wratio_connN, vel_bias_connN, cum_bias_connN, + None, vel_bias_connF, cum_bias_connF] + titles = ['w_ratio', f'bias_vel_{conn}', f'bias_{conn}', + None, 'bias_vel_F', 'bias_F'] + + fig, axs = plt.subplots(nrows=2, ncols=3, figsize=[12, 6]) + for ax, data, title in zip(axs.flatten(), data_list, titles): + if data is not None: + im = ax.imshow(data, cmap='jet', interpolation='nearest') + ax.set_title(title) + fig.colorbar(im, ax=ax) + else: + ax.axis('off') + fig.tight_layout() + plt.show() + + return wratio_connN, vel_bias_connN + + +def estimate_wratio_all(bw, bias_free_conn, outdir, box): + '''Estimate diaginal matrix W_r for all connections levels within the given bandwidth. + + Parameters: bias_free_conn - integer, minimum connection-level that we think is bias-free + bw - integer, bandwidth of given time-sereis analysis + box - list in size of (4,) in integer, coordinates of bounding box + outdir - string, the working directory + Returns: wratio - 3D array in size of (bw+1, length, width) in float32, + the first slice (w[0,:,:]) is a padding + to ensure that wratio[n,:,:] = w(n * delta_t) / w(delta_t). + ''' + box_wid = box[2] - box[0] + box_len = box[3] - box[1] + cum_bias_connF = read_cum_seq_closure_phase4conn(bias_free_conn, outdir, box)[-1,:,:] + # skip invalid pixels + flag = np.multiply(~np.isnan(cum_bias_connF), cum_bias_connF != 0) + + wratio = np.ones([bw+1, box_len, box_wid], dtype=np.float32) + for n in np.arange(2, bw+1): + cum_bias_connN = read_cum_seq_closure_phase4conn(n, outdir, box)[-1,:,:] + wratio[n,flag] = 1 - cum_bias_connN[flag] / cum_bias_connF[flag] + + # bound into [0, 1] + wratio[wratio > 1] = 1 + wratio[wratio < 0] = 0 + + return wratio + + +def get_avg_time_span_within_bandwidth(date_ordinal, bw): + '''Compute average temporal span (days) for all interferogram within the given bandwidth + + Parameters: date_ordinal - list of size (num_date,) in integer, time in days + bw - integer, bandwidth of time-series analysis + Return: avg_time - float, average time-span in days + ''' + avg_time = 0 + num_ifgram = 0 + for level in range(1, bw+1): + slc_date_firstN = date_ordinal[0:level] + slc_date_lastN = date_ordinal[-level:] + for i in range(level): + avg_time += slc_date_lastN[i] - slc_date_firstN[i] + num_ifgram += len(date_ordinal) - level + + avg_time /= num_ifgram + + return avg_time + + +def get_avg_time_span4conn(date_ordinal, conn): + '''Compute the average temporal span (days) for connection-n interferograms - mask = np.ones([length,width],np.float32) - mask[np.abs(np.angle(closurephase))>threshold_cp] = 0 # this masks areas with potential bias - mask[np.abs(np.abs(closurephase)/numcp < inps.episilon)] = 1 # this unmasks areas with low correlation (where it's hard to know wheter there is bias either) + Parameters: date_ordinal - list of size (num_date,) in integer, time in days + conn - int, connection level of interferograms + Return: avg_time - float, average time-span in days + ''' + slc_date_firstN = date_ordinal[0:conn] + slc_date_lastN = date_ordinal[-conn:] + avg_time = 0 + for i in range(conn): + avg_time += slc_date_lastN[i] - slc_date_firstN[i] - # save mask + num_ifgram = len(date_ordinal) - conn + avg_time /= num_ifgram + + return avg_time + + +def estimate_bias_timeseries_approx_patch(bias_free_conn, bw, tbase, date_ordinal, wvl, box, outdir): + '''Quick and approximate estimate of the bias time-series of a certain bandwidth (bw) for a bounding box + + Note: This estimate is not exact, but often close enough. + It is good for a quick estimate to see how big the biases are. + + Parameters: bias_free_conn - integer, connection level that we assume bias-free + bw - integer, bandwidth of the given time-series analysis + tbase - 1D np.ndarray in size of (num_date,) in float, time in accumulated years + date_ordinal - list of size (num_date,) in integer, time in days + wvl - float, wavelength of the SAR system + box - list in size of (4,) in integer, coordinates of bounding box + outdir - string, directory for outputing files + Returns: bias_ts - 3D array in size of (num_date, box_len, box_wid) in float, bias timeseries + ''' + print('\n'+'-'*60) + print(f'quick and approximate estimation of bias time series for bandwidth = {bw}') + # basic info + phase2range = -1 * wvl / (4 * np.pi) + num_date = tbase.size + + # average temporal span for ifgrams of connection-1 to connection-bw + deltat_n = np.asarray([get_avg_time_span4conn(date_ordinal, n) for n in range(1, bw+1)]) + avg_time_span = get_avg_time_span_within_bandwidth(date_ordinal, bw) + + # the bias in a bandwidth-bw analysis is similar to bias in connectoin-p interferograms + p = (np.abs(deltat_n - avg_time_span)).argmin() + 1 + print(f'average connection level within the bandwidth = {p}') + + # get wratio + kwargs1 = dict( + bias_free_conn=bias_free_conn, + wvl=wvl, + box=box, + outdir=outdir, + ) + wratio_p = estimate_wratio(tbase, conn=p, **kwargs1)[0] + wratio_2 = estimate_wratio(tbase, conn=2, **kwargs1)[0] + wratio_p[np.isnan(wratio_p)] = 0 + wratio_2[np.isnan(wratio_2)] = 0 + + wratio_2[abs(wratio_2 - 1) < 0.1] = np.nan + ratio_p2 = wratio_p / (1 - wratio_2) + + # get bias_ts + kwargs2 = dict(outdir=outdir, box=box, print_msg=True) + bias_ts = read_cum_seq_closure_phase4conn(2, **kwargs2) * phase2range + bias_ts_bf = read_cum_seq_closure_phase4conn(bias_free_conn, **kwargs2) * phase2range + for i in range(num_date): + bias_ts[i,:,:] *= ratio_p2 + bias_ts_bf[i,:,:] *= wratio_p + + flag = np.isnan(bias_ts) + bias_ts[flag] = bias_ts_bf[flag] + + return bias_ts + + +def estimate_bias_timeseries_approx(stack_file, bias_free_conn, bw, water_mask_file=None, outdir='./', max_memory=4.0): + '''Quick & approximate estimation of the bias time series and Wr. + + Reference: Eq. (20) in Zheng et al. (2022, TGRS). + + Parameters: stack_file - string, path for ifgramStack.h5 + bias_free_conn - integer, connection level at which we assume is bias-free + bw - integer, bandwidth of the given time-series. + wvl - float, wavelength of the SAR System + outdir - str, directory for output files + max_mermory - float, maximum memory in GB for each patch processed + Returns: bias_ts_file - str, path to the HDF5 file for the approximate bias time series in (num_date, length, width) + wratio_file - str, path to the HDF5 file for the wratio and bias velocity in (bw, length, width) + Shows how fast the bias-inducing signal decays with temporal baseline. + ''' + print('\n'+'-'*80) + print(f'quick estimation of the non-closure phase bias time-series for bandwidth={bw} (Zheng et al., 2022) ...') + + stack_obj = ifgramStack(stack_file) + stack_obj.open() + length, width = stack_obj.length, stack_obj.width meta = dict(stack_obj.metadata) - meta['FILE_TYPE'] = 'mask' - ds_name_dict = {'cpmask': [np.float32, (length, width), mask],} - writefile.layout_hdf5(os.path.join(inps.outdir,'cpmask.h5'), ds_name_dict, meta) + wvl = float(meta['WAVELENGTH']) + + # time info + date_list = stack_obj.get_date_list(dropIfgram=True) + num_date = len(date_list) + + tbase = np.array(ptime.date_list2tbase(date_list)[0], dtype=np.float32) / 365.25 + date_str_fmt = ptime.get_date_str_format(date_list[0]) + date_ordinal = [dt.strptime(x, date_str_fmt).toordinal() for x in date_list] + + # split igram_file into blocks to save memory + box_list, num_box = stack_obj.split2boxes(max_memory=max_memory, dim0_size=num_date) + + # initiate output files + # 1 - wratio file + wratio_file = os.path.join(outdir, 'wratio.h5') + meta['FILE_TYPE'] = 'wratio' + meta['DATA_TYPE'] = 'float32' + meta['UNIT'] = '1' + ds_name_dict = { + 'wratio' : [np.float32, (bw, length, width), None], + 'velocityBias' : [np.float32, (bw, length, width), None], + } + ds_unit_dict = { + 'wratio' : '1', + 'velocityBias' : 'm/year', + } + writefile.layout_hdf5(wratio_file, ds_name_dict, metadata=meta, ds_unit_dict=ds_unit_dict) + + # 2 - time series file + bias_ts_file = os.path.join(outdir, 'timeseriesBiasApprox.h5') + meta['FILE_TYPE'] = 'timeseries' + meta['UNIT'] = 'm' + ds_name_dict = { + 'timeseries' : [np.float32, (num_date, length, width), None], + 'date' : [np.dtype('S8'), (num_date,), np.array(date_list, np.string_)], + } + writefile.layout_hdf5(bias_ts_file, ds_name_dict, meta) + + # process block-by-block + for i, box in enumerate(box_list): + box_wid = box[2] - box[0] + box_len = box[3] - box[1] + print(box) + if num_box > 1: + print('\n------- processing patch {} out of {} --------------'.format(i+1, num_box)) + print('box width: {}'.format(box_wid)) + print('box length: {}'.format(box_len)) + + # read water mask + if water_mask_file: + print(f'skip pixels on water (zero value from file: {os.path.basename(water_mask_file)})') + water_mask = readfile.read(water_mask_file, box=box)[0] + else: + water_mask_file = None + + # 1 - estimate the wratio(_velocity) + wratio = np.zeros([bw, box_len, box_wid], dtype=np.float32) + bias_vel = np.zeros([bw, box_len, box_wid], dtype=np.float32) + for j in range(bw): + print(f'estimation W_ratio for connection level: {j+1}') + wratio[j, :, :], bias_vel[j, :, :] = estimate_wratio( + tbase, + conn=j+1, + bias_free_conn=bias_free_conn, + wvl=wvl, + box=box, + outdir=outdir, + mask=True) - # also save the average closure phase - ds_name_dict2 = {'phase': [np.float32, (length, width), np.angle(closurephase)], - 'amplitude':[np.float32,(length,width),np.abs(closurephase)/numcp],} - writefile.layout_hdf5(os.path.join(inps.outdir,'avgwcp.h5'), ds_name_dict2, meta) + if water_mask_file: + wratio[:, water_mask==0] = np.nan + bias_vel[:, water_mask==0] = np.nan + # write the block to disk + block = [0, bw, box[1], box[3], box[0], box[2]] + + writefile.write_hdf5_block( + wratio_file, + data=wratio, + datasetName='wratio', + block=block) + + writefile.write_hdf5_block( + wratio_file, + data=bias_vel, + datasetName='velocityBias', + block=block) + + # 2 - estimate the bias time series + bias_ts = estimate_bias_timeseries_approx_patch( + bias_free_conn=bias_free_conn, + bw=bw, + tbase=tbase, + date_ordinal=date_ordinal, + wvl=wvl, + box=box, + outdir=outdir) + + if water_mask_file: + bias_ts[:, water_mask==0] = np.nan + + # write the block to disk + block = [0, len(date_list), box[1], box[3], box[0], box[2]] + writefile.write_hdf5_block( + bias_ts_file, + data=bias_ts, + datasetName='timeseries', + block=block) + + return bias_ts_file, wratio_file + + + +################################################################################ +def bandwidth2num_ifgram(bw, num_date): + '''Get the number of interferograms for the network with the given bandwidth. + + Reference: Equation (15) in Zheng et al. (2022) + + Parameters: bw - int, bandwith + num_date - int, number of acquisitions + Returns: num_ifgram - int, number of interferograms + ''' + return int(bw * (num_date * 2 - bw - 1) / 2) + + +def get_design_matrix_Wr(date12_list, bw, box, bias_free_conn, outdir='./'): + '''Computes the matrix W_r for a given bounding box, following Equation (20). + + Parameters: date12_list - list(str), interfeorgram pairs list in YYYYMMDD_YYYYMMDD + bw - integer, bandwidth of time-series analysis + bias_free_conn - integer, minimum connection-level that we think is bias-free + box - list in size of (4,) in integer, coordinates of bounding box + outdir - string, the working directory + Returns: Wr - 2D np.ndarray in size of (num_ifgram, num_pix) in float32, + each row stores the diagnal component of W (Eq. 16 in Zheng et al., 2022) for one pixel. + A - 2D np.ndarray in size of (num_ifgram, num_date) in float32 + ''' + # get design matrix A + A = ifgramStack.get_design_matrix4timeseries(date12_list=date12_list, refDate='no')[0] + num_ifgram = A.shape[0] + + # get w(delta_t) * phi^x - section VI-A + wratio_all = estimate_wratio_all(bw, bias_free_conn, outdir, box) + + # intial output value + num_pix = (box[2] - box[0]) * (box[3] - box[1]) + Wr = np.zeros((num_ifgram, num_pix), dtype=np.float32) + for i in range(num_ifgram): + # get the connection level + Aline = list(A[i,:]) + idx1 = Aline.index(-1) + idx2 = Aline.index(1) + conn = idx2 - idx1 + if conn > bw: + print('Existing max-conn-level > the input bandwidth, ' + 'use modify_network.py inputs/ifgramStack.h5 to adjust the max-conn-level.') + + # assign to Wr matrix + Wr[i, :] = wratio_all[conn, :, :].reshape(-1) + + return Wr, A + + +def estimate_bias_timeseries_patch(stack_file, bias_free_conn, bw, wvl, box, water_mask_file=None, outdir='./'): + '''Estimate the bias time-series of a certain bandwidth (bw) for a bounding box. + + Reference: Zheng et al. (2022, TGRS). + + Parameters: stack_file - string, path for ifgramStack.h5 + bias_free_conn - integer, connection level at which we assume is bias-free + bw - integer, bandwidth of the given time-series. + wvl - float, wavelength of the SAR System + box - list in size of (4,) in integer, coordinates of bounding box + outdir - string, directory for output files + Returns: bias_ts_bwn - 3D array of size (bw, box_len, box_wid) of float, estimated bias time-series + box - list in size of (4,) in integer, coordinates of bounding box, output for parallel computing + ''' + phase2range = -1 * wvl / (4 * np.pi) + box_wid, box_len = box[2] - box[0], box[3] - box[1] + num_pix = box_wid * box_len + + stack_obj = ifgramStack(stack_file) + stack_obj.open(print_msg=False) + date12_list = stack_obj.get_date12_list(dropIfgram=True) + num_ifgram = len(date12_list) + + # time info + date_list = stack_obj.get_date_list(dropIfgram=True) + num_date = len(date_list) + # tbase in the unit of years + tbase = np.array(ptime.date_list2tbase(date_list)[0], dtype=np.float32) / 365.25 + tbase_diff = np.diff(tbase).reshape(-1,1) + + + ## mask of pixels to invert + mask = np.ones(num_pix, dtype=np.bool_) + + # water mask + if water_mask_file and os.path.isfile(water_mask_file): + print(f'skip pixels (on the water) with zero value in file: {os.path.basename(water_mask_file)}') + water_mask = readfile.read(water_mask_file, box=box)[0].flatten() + mask *= np.array(water_mask, dtype=np.bool_) + del water_mask + else: + water_mask_file = None + + # invert pixels on mask(s) + num_pix2inv = int(np.sum(mask)) + idx_pix2inv = np.where(mask)[0] + print('number of pixels to invert: {} out of {} ({:.1f}%)'.format( + num_pix2inv, num_pix, num_pix2inv/num_pix*100)) + + + ## 1. get bias time-series for bw-1 analysis + kwargs = dict(outdir=outdir, box=box, print_msg=True) + bias_ts_bw1_rough = read_cum_seq_closure_phase4conn(bias_free_conn, **kwargs).reshape(num_date, -1) + bias_ts_bw1_fine = read_cum_seq_closure_phase4conn(2, **kwargs).reshape(num_date, -1) + + bias_vel_bw1 = bias_ts_bw1_fine[-1,:] * phase2range / (tbase[-1] - tbase[0]) + flag = np.where(np.abs(bias_vel_bw1) < 0.001, 0, 1).astype(np.bool_) + num_pix_less, num_pix_more = np.sum(flag[mask]), np.sum(~flag[mask]) + digit = len(str(num_pix)) + msg = 'number of pixels with bandwidth=1 velocity bias ' + msg += f'< | >= 0.1 mm/yr: {num_pix_less:{digit}d} | {num_pix_more:{digit}d} out of {num_pix} ' + msg += f'({num_pix_less/num_pix*100:.0f}% | {num_pix_less/num_pix*100:.0f}%)' + print(msg) + + # scale bias_ts_bw1_fine based on bias_ts_bw1_rough + r2f_flag = np.multiply(~np.isnan(bias_ts_bw1_fine[-1,:]), bias_ts_bw1_fine[-1,:] != 0) + r2f_scale = np.ones((num_pix), dtype=np.float32) + r2f_scale[r2f_flag] = bias_ts_bw1_rough[-1,r2f_flag] / bias_ts_bw1_fine[-1,r2f_flag] + for i in range(num_date): + bias_ts_bw1_fine[i,:] *= r2f_scale + del r2f_flag, r2f_scale + + + # 2. construct bias_stack = W * A * Phi^X = Wr * A * w(delta_t) * Phi^X + # Equation (20) in Zheng et al. (2022, TGRS) + # this matrix is a num_pix by num_ifgram matrix, each row stores the diagnal component of the Wr matrix for that pixel + print('estimating bias_stack = Wr * A * w(delta_t) * Phi^X (Zheng et al., 2022, TGRS) ...') + Wr, A = get_design_matrix_Wr(date12_list, bw, box, bias_free_conn, outdir) + wPhi_x = np.array(bias_ts_bw1_rough, dtype=np.float32) + wPhi_x[:, flag] = bias_ts_bw1_fine[:, flag] + + bias_stack = np.zeros((num_ifgram, num_pix), dtype=np.float32) + prog_bar = ptime.progressBar(maxValue=num_pix2inv) + for i in range(num_pix2inv): + idx = idx_pix2inv[i] + + # calculate the bias_stack = W * A * phi^x = W^r * A * w(delta_t) * phi^x + bias_stack[:, idx] = np.linalg.multi_dot([np.diag(Wr[:, idx]), A, wPhi_x[:, idx]]).flatten() + + prog_bar.update(i+1, every=3000, suffix='{}/{} pixels'.format(i+1, num_pix2inv)) + prog_bar.close() + del bias_ts_bw1_rough, bias_ts_bw1_fine, wPhi_x, Wr + + + # 3. estimate bias time-series from bias stack: bias_ts = A+ * bias_stack + # Equation (20) in Zheng et al. (2022, TGRS) + # perform phase velocity inversion as per the original SBAS paper rather doing direct phase inversion. + print('estimating bias time-series from bias stack via the SBAS approach ...') + A1, B1 = stack_obj.get_design_matrix4timeseries(date12_list=date12_list) + kwargs = { + 'A' : A1, + 'B' : B1, + 'tbase_diff' : tbase_diff, + 'weight_sqrt' : None, + 'min_norm_velocity' : True, + 'inv_quality_name' : 'no', + } + + # a. split mask into mask_all/par_net + # mask for valid (~NaN) observations in ALL ifgrams (share one B in sbas inversion) + mask_all_net = np.all(~np.isnan(bias_stack), axis=0) * np.all(bias_stack != 0, axis=0) + mask_all_net *= mask + mask_par_net = mask ^ mask_all_net + msg = 'estimating time-series for pixels with valid stack values' + + # b. invert once for all pixels with obs in ALL ifgrams + bias_ts = np.zeros((num_date, num_pix), dtype=np.float32) + if np.sum(mask_all_net) > 0: + num_pix_all = int(np.sum(mask_all_net)) + print(f'{msg} in all ifgrams ({num_pix_all} pixels; {num_pix_all/num_pix2inv*100:.0f}%) ...') + + # invert + bias_ts[:, mask_all_net] = estimate_timeseries(y=bias_stack[:, mask_all_net], **kwargs)[0] + + # c. invert pixel-by-pixel for pixels with obs NOT in all ifgrams + if np.sum(mask_par_net) > 0: + num_pix_par = int(np.sum(mask_par_net)) + idx_pix_par = np.where(mask_par_net)[0] + print(f'{msg} in some ifgrams ({num_pix_par} pixels; {num_pix_par/num_pix2inv*100:.0f}%) ...') + + prog_bar = ptime.progressBar(maxValue=num_pix_par) + for i in range(num_pix_par): + idx = idx_pix_par[i] + # invert + bias_ts[:, idx] = estimate_timeseries(y=bias_stack[:, idx], **kwargs)[0].flatten() + + prog_bar.update(i+1, every=200, suffix='{}/{} pixels'.format(i+1, num_pix_par)) + prog_bar.close() + del bias_stack + + bias_ts = bias_ts.reshape(num_date, box_len, box_wid) * phase2range + + return bias_ts, box + + +def estimate_bias_timeseries(stack_file, bias_free_conn, bw, cluster_kwargs, water_mask_file=None, outdir='./', max_memory=4.0): + '''Run the bias time-series estimation. + + Parameters: stack_file - string, path for ifgramStack.h5 + bias_free_conn - integer, connection level at which we assume is bias-free + bw - integer, bandwidth of the given time-series. + cluster_kwargs - dictonary containing settings of parallel computing. To turn off, set parallel['clustertype']='' + outdir - string, directory for output files + max_memory - float, maximum memory in GB for each patch processed + Returns: bias_ts_file - str, path to the bias time series file: timeseriesBias.h5 + ''' + print('\n'+'-'*80) + print(f'estimating the non-closure phase bias time-series for bandwidth={bw} (Zheng et al., 2022) ...') + + stack_obj = ifgramStack(stack_file) + stack_obj.open(print_msg=False) + length, width = stack_obj.length, stack_obj.width + meta = dict(stack_obj.metadata) + wvl = float(meta['WAVELENGTH']) + + date12_list = stack_obj.get_date12_list(dropIfgram=True) + date_list = stack_obj.get_date_list(dropIfgram=True) + num_ifgram = len(date12_list) + num_date = len(date_list) + + # check the bandwidth of ifgramStack + num_ifgram_exp = bandwidth2num_ifgram(bw, num_date) + print(f'number of interferograms expected for bandwidth={bw}: {num_ifgram_exp}') + print(f'number of interferograms kept in ifgramStack.h5 : {num_ifgram}') + if num_ifgram != num_ifgram_exp: + msg = f'number of the kept interferograms ({num_ifgram}) is NOT the same as expected ({num_ifgram_exp})!' + msg += f'\n This indicates the bandwidth between ifgramStack.h5 and the user input ({bw}) are NOT consistent!' + msg += '\n Modify the network of interferograms to be consistent via modify_network.py by:' + msg += f'\n 1) set mintpy.network.connNumMax={bw} and 2) re-run modify_network.py -t smallbaselineApp.cfg' + raise Exception(msg) + + # estimate for bias time-series + bias_ts_file = os.path.join(outdir, 'timeseriesBias.h5') + ds_name_dict = { + 'timeseries' : [np.float32, (num_date, length, width), None], + 'date' : [np.dtype('S8'), (num_date,), np.array(date_list, np.string_)], + } + meta['FILE_TYPE'] = 'timeseries' + meta['DATA_TYPE'] = 'float32' + meta['REF_DATE'] = date_list[0] + meta['UNIT'] = 'm' + meta['BANDS'] = '1' + meta['DATE12'] = f'{date_list[0][2:]}-{date_list[-1][2:]}' + writefile.layout_hdf5(bias_ts_file, ds_name_dict, meta) + + data_kwargs = { + "stack_file" : stack_file, + "bias_free_conn" : bias_free_conn, + "bw" : bw, + "wvl" : wvl, + "water_mask_file" : water_mask_file, + "outdir" : outdir, + } + + # split igram_file into blocks to save memory + box_list, num_box = stack_obj.split2boxes(max_memory=max_memory, dim0_size=num_ifgram*2+num_date*3) + num_threads_dict = cluster.set_num_threads("1") + + for i, box in enumerate(box_list): + box_wid, box_len = box[2] - box[0], box[3] - box[1] + print(box) + if num_box > 1: + print('\n------- processing patch {} out of {} --------------'.format(i+1, num_box)) + print('box width : {}'.format(box_wid)) + print('box length: {}'.format(box_len)) + + #update box argument in the input data + data_kwargs['box'] = box + + if not cluster_kwargs['cluster_type']: + # non-parallel + bias_ts = estimate_bias_timeseries_patch(**data_kwargs)[:-1] + + else: + # parallel + print('\n\n------- start parallel processing using Dask -------') + + # initiate the output data + bias_ts = np.zeros((num_date, box_len, box_wid), dtype=np.float32) + + # initiate dask cluster and client + cluster_obj = cluster.DaskCluster(**cluster_kwargs) + cluster_obj.open() + + # run dask + bias_ts = cluster_obj.run( + func=estimate_bias_timeseries_patch, + func_data=data_kwargs, + results=[bias_ts], + ) + + # close dask cluster and client + cluster_obj.close() + print('------- finished parallel processing -------\n\n') + + writefile.write_hdf5_block( + bias_ts_file, + data=bias_ts, + datasetName='timeseries', + block=[0, num_date, box[1], box[3], box[0], box[2]], + ) + + # roll back to the original number of threads + cluster.roll_back_num_threads(num_threads_dict) + + return bias_ts_file + + +################################################################################ +def main(iargs=None): + inps = cmd_line_parse(iargs) + start_time = time.time() + + # common inputs + kwargs = dict(outdir=inps.outdir, max_memory=inps.maxMemory) + + if inps.action == 'mask': + calc_closure_phase_mask( + stack_file=inps.stack_file, + bias_free_conn=inps.nl, + num_sigma=inps.num_sigma, + threshold_amp=inps.epsilon, + **kwargs) + + elif inps.action.endswith ('estimate'): + # compute the unwrapped closure phase bias time-series + # and re-unwrap to mitigate the impact of phase unwrapping errors + # which can dominate the true non-closure phase. + # max(2, inps.bw) is used to ensure we have conn-2 closure phase processed + conn_list = np.arange(2, max(2, inps.bw) + 1).tolist() + [inps.nl] + conn_list = sorted(list(set(conn_list))) + for conn in conn_list: + print('\n'+'-'*80) + print('calculating the unwrapped closure phase for ' + f'connection level = {conn} out of {conn_list} ...') + compute_unwrap_closure_phase( + stack_file=inps.stack_file, + conn=conn, + num_worker=int(inps.numWorker), + **kwargs) + + if inps.action == 'quick_estimate': + estimate_bias_timeseries_approx( + stack_file=inps.stack_file, + bias_free_conn=inps.nl, + bw=inps.bw, + water_mask_file=inps.water_mask_file, + **kwargs) + + elif inps.action == 'estimate': + cluster_kwargs = { + "cluster_type" : inps.cluster, + "num_worker" : inps.numWorker, + "config_name" : inps.config} + estimate_bias_timeseries( + stack_file=inps.stack_file, + bias_free_conn=inps.nl, + bw=inps.bw, + cluster_kwargs=cluster_kwargs, + water_mask_file=inps.water_mask_file, + **kwargs) + + # used time + m, s = divmod(time.time() - start_time, 60) + print('time used: {:02.0f} mins {:02.1f} secs.\n'.format(m, s)) + + return + + +################################################################################ if __name__ == '__main__': main(sys.argv[1:]) diff --git a/mintpy/ifgram_inversion.py b/mintpy/ifgram_inversion.py index 29ac0695c..63498166d 100755 --- a/mintpy/ifgram_inversion.py +++ b/mintpy/ifgram_inversion.py @@ -335,7 +335,7 @@ def estimate_timeseries(A, B, y, tbase_diff, weight_sqrt=None, min_norm_velocity y - 2D np.ndarray in size of (num_pair, num_pixel), phase/offset of all interferograms with no-data value: NaN. tbase_diff - 2D np.ndarray in size of (num_date-1, 1), - differential temporal baseline history + differential temporal baseline history, in the unit of years weight_sqrt - 2D np.ndarray in size of (num_pair, num_pixel), square root of weight of all interferograms min_norm_velocity - bool, assume minimum-norm deformation velocity, or not @@ -345,6 +345,7 @@ def estimate_timeseries(A, B, y, tbase_diff, weight_sqrt=None, min_norm_velocity inv_quality_name - str, inversion quality type/name temporalCoherence for phase residual for offset + no to turn OFF the calcualtion Returns: ts - 2D np.ndarray in size of (num_date, num_pixel), phase time-series inv_quality - 1D np.ndarray in size of (num_pixel), temporal coherence (for phase) or residual (for offset) num_inv_obs - 1D np.ndarray in size of (num_pixel), number of observations (ifgrams / offsets) @@ -395,10 +396,11 @@ def estimate_timeseries(A, B, y, tbase_diff, weight_sqrt=None, min_norm_velocity X, e2 = linalg.lstsq(B, y, cond=rcond)[:2] # calc inversion quality - inv_quality = calc_inv_quality(B, X, y, e2, - inv_quality_name=inv_quality_name, - weight_sqrt=weight_sqrt, - print_msg=print_msg) + if inv_quality_name != 'no': + inv_quality = calc_inv_quality(B, X, y, e2, + inv_quality_name=inv_quality_name, + weight_sqrt=weight_sqrt, + print_msg=print_msg) # assemble time-series ts_diff = X * np.tile(tbase_diff, (1, num_pixel)) @@ -414,10 +416,11 @@ def estimate_timeseries(A, B, y, tbase_diff, weight_sqrt=None, min_norm_velocity X, e2 = linalg.lstsq(A, y, cond=rcond)[:2] # calc inversion quality - inv_quality = calc_inv_quality(A, X, y, e2, - inv_quality_name=inv_quality_name, - weight_sqrt=weight_sqrt, - print_msg=print_msg) + if inv_quality_name != 'no': + inv_quality = calc_inv_quality(A, X, y, e2, + inv_quality_name=inv_quality_name, + weight_sqrt=weight_sqrt, + print_msg=print_msg) # assemble time-series ts[1: ,:] = X @@ -569,41 +572,6 @@ def calc_inv_quality(G, X, y, e2, inv_quality_name='temporalCoherence', weight_s ###################################### File IO ############################################ -def split2boxes(ifgram_file, max_memory=4, print_msg=True): - """Split into chunks in rows to reduce memory usage - Parameters: dataset_shape - tuple of 3 int - max_memory - float, max memory to use in GB - print_msg - bool - Returns: box_list - list of tuple of 4 int - num_box - int, number of boxes - """ - ifg_obj = ifgramStack(ifgram_file) - ifg_obj.open(print_msg=False) - - # dataset size: defo obs (phase / offset) + weight + time-series - length = ifg_obj.length - width = ifg_obj.width - ds_size = (ifg_obj.numIfgram * 2 + ifg_obj.numDate + 5) * length * width * 4 - - num_box = int(np.ceil(ds_size * 1.5 / (max_memory * 1024**3))) - y_step = int(np.ceil((length / num_box) / 10) * 10) - num_box = int(np.ceil(length / y_step)) - if print_msg and num_box > 1: - print('maximum memory size: %.1E GB' % max_memory) - print('split %d lines into %d patches for processing' % (length, num_box)) - print(' with each patch up to %d lines' % y_step) - - # y_step / num_box --> box_list - box_list = [] - for i in range(num_box): - y0 = i * y_step - y1 = min([length, y0 + y_step]) - box = (0, y0, width, y1) - box_list.append(box) - - return box_list, num_box - - def check_design_matrix(ifgram_file, weight_func='var'): """ Check Rank of Design matrix for weighted inversion @@ -1029,8 +997,19 @@ def ifgram_inversion_patch(ifgram_file, box=None, ref_phase=None, obs_ds_name='u num_inv_obs = num_inv_obs.reshape(num_row, num_col) return ts, ts_cov, inv_quality, num_inv_obs, box + # common inversion options + kwargs = { + 'A' : A, + 'B' : B, + 'tbase_diff' : tbase_diff, + 'min_norm_velocity' : min_norm_velocity, + 'min_redundancy' : min_redundancy, + 'inv_quality_name' : inv_quality_name, + } + # 2.2 un-weighted inversion (classic SBAS) if weight_sqrt is None: + msg = f'estimating time-series for pixels with valid {obs_ds_name} in' # a. split mask into mask_all/part_net # mask for valid (~NaN) observations in ALL ifgrams (share one B in sbas inversion) @@ -1041,19 +1020,14 @@ def ifgram_inversion_patch(ifgram_file, box=None, ref_phase=None, obs_ds_name='u # b. invert once for all pixels with obs in all ifgrams if np.sum(mask_all_net) > 0: - print(('estimating time-series for pixels with valid {} in all ifgrams' - ' ({:.0f} pixels; {:.1f}%) ...').format(obs_ds_name, - np.sum(mask_all_net), - np.sum(mask_all_net)/num_pixel2inv*100)) - (tsi, - inv_quali, - num_obsi) = estimate_timeseries(A, B, - y=stack_obs[:, mask_all_net], - tbase_diff=tbase_diff, - weight_sqrt=None, - min_norm_velocity=min_norm_velocity, - min_redundancy=min_redundancy, - inv_quality_name=inv_quality_name) + num_pixel2inv_all = int(np.sum(mask_all_net)) + print(f'{msg} all ifgrams ({num_pixel2inv_all} pixels; {num_pixel2inv_all/num_pixel2inv*100:.1f}%) ...') + + # run + tsi, inv_quali, num_obsi = estimate_timeseries( + y=stack_obs[:, mask_all_net], + weight_sqrt=None, + **kwargs) # save result to output matrices ts[:, mask_all_net] = tsi @@ -1062,29 +1036,25 @@ def ifgram_inversion_patch(ifgram_file, box=None, ref_phase=None, obs_ds_name='u # c. pixel-by-pixel for pixels with obs not in all ifgrams if np.sum(mask_part_net) > 0: - print(('estimating time-series for pixels with valid {} in some ifgrams' - ' ({:.0f} pixels; {:.1f}%) ...').format(obs_ds_name, - np.sum(mask_part_net), - np.sum(mask_all_net)/num_pixel2inv*100)) num_pixel2inv_part = int(np.sum(mask_part_net)) idx_pixel2inv_part = np.where(mask_part_net)[0] + print(f'{msg} some ifgrams ({num_pixel2inv_part} pixels; {num_pixel2inv_part/num_pixel2inv*100:.1f}%) ...') + prog_bar = ptime.progressBar(maxValue=num_pixel2inv_part) for i in range(num_pixel2inv_part): idx = idx_pixel2inv_part[i] - (tsi, - inv_quali, - num_obsi) = estimate_timeseries(A, B, - y=stack_obs[:, idx], - tbase_diff=tbase_diff, - weight_sqrt=None, - min_norm_velocity=min_norm_velocity, - min_redundancy=min_redundancy, - inv_quality_name=inv_quality_name) + + # run + tsi, inv_quali, num_obsi = estimate_timeseries( + y=stack_obs[:, idx], + weight_sqrt=None, + **kwargs) # save result to output matrices ts[:, idx] = tsi.flatten() inv_quality[idx] = inv_quali num_inv_obs[idx] = num_obsi + prog_bar.update(i+1, every=200, suffix='{}/{} pixels'.format(i+1, num_pixel2inv_part)) prog_bar.close() @@ -1094,15 +1064,12 @@ def ifgram_inversion_patch(ifgram_file, box=None, ref_phase=None, obs_ds_name='u prog_bar = ptime.progressBar(maxValue=num_pixel2inv) for i in range(num_pixel2inv): idx = idx_pixel2inv[i] - (tsi, - inv_quali, - num_obsi) = estimate_timeseries(A, B, - y=stack_obs[:, idx], - tbase_diff=tbase_diff, - weight_sqrt=weight_sqrt[:, idx], - min_norm_velocity=min_norm_velocity, - min_redundancy=min_redundancy, - inv_quality_name=inv_quality_name) + + # run + tsi, inv_quali, num_obsi = estimate_timeseries( + y=stack_obs[:, idx], + weight_sqrt=weight_sqrt[:, idx], + **kwargs) # save result to output matrices ts[:, idx] = tsi.flatten() @@ -1312,7 +1279,7 @@ def ifgram_inversion(inps=None): ## 3. run the inversion / estimation and write to disk # 3.1 split ifgram_file into blocks to save memory - box_list, num_box = split2boxes(inps.ifgramStackFile, max_memory=inps.maxMemory) + box_list, num_box = stack_obj.split2boxes(max_memory=inps.maxMemory) # 3.2 prepare the input arguments for *_patch() data_kwargs = { diff --git a/mintpy/objects/cluster.py b/mintpy/objects/cluster.py index f56c0b65f..dfdc57071 100644 --- a/mintpy/objects/cluster.py +++ b/mintpy/objects/cluster.py @@ -302,7 +302,7 @@ def collect_result(self, futures, results, box, submission_time): # message num_future += 1 sub_t = time.time() - submission_time - print("FUTURE #{} complete. Time used: {:.0f} seconds".format(num_future, sub_t)) + print(f"\nFUTURE #{num_future} complete. Time used: {sub_t:.0f} seconds") # catch result - sub_box # and convert the abosulte sub_box into local col/row start/end relative to the primary box diff --git a/mintpy/objects/stack.py b/mintpy/objects/stack.py index c55d519c6..b3fb449c8 100644 --- a/mintpy/objects/stack.py +++ b/mintpy/objects/stack.py @@ -966,10 +966,158 @@ def get_max_connection_number(self): num_conn[i] = np.where(Ai == 1)[0] - np.where(Ai == -1)[0] return np.max(num_conn) - # Functions for Unwrap error correction + + def split2boxes(self, max_memory=4, dim0_size=None, print_msg=True): + """Split into chunks in rows to reduce memory usage. + + Parameters: max_memory - float, max memory to use in GB + dim0_size - the 1st dimension size of all used datasets + e.g., dim0_size = num_pair * 2 + num_date + print_msg - bool + Returns: box_list - list of tuple of 4 int + num_box - int, number of boxes + """ + self.open(print_msg=False) + length = self.length + width = self.width + + # dimension in time: phase/offset, weight, timeseries, etc. + if not dim0_size: + # for time series estimation + dim0_size = self.numIfgram * 2 + self.numDate + ds_size = dim0_size * length * width * 4 + + num_box = int(np.ceil(ds_size * 1.5 / (max_memory * 1024**3))) + y_step = int(np.ceil((length / num_box) / 10) * 10) + num_box = int(np.ceil(length / y_step)) + if print_msg and num_box > 1: + print('maximum memory size: %.1E GB' % max_memory) + print('split %d lines into %d patches for processing' % (length, num_box)) + print(' with each patch up to %d lines' % y_step) + + # y_step / num_box --> box_list + box_list = [] + for i in range(num_box): + y0 = i * y_step + y1 = min([length, y0 + y_step]) + box = (0, y0, width, y1) + box_list.append(box) + + return box_list, num_box + + + # Functions for closure phase bias + def get_closure_phase_index(self, conn, dropIfgram=True): + """Get the indices of interferograms that forms the given connection level closure loop. + + Parameters: conn - int, connection level + dropIfgram - bool, exclude the dropped interferograms. + Returns: cp_idx - 2D np.ndarray in int16 in size of (num_cp, conn + 1) + Each row for the indices of interferograms for one closure loop. + num_cp <= num_date - conn + """ + date12_list = self.get_date12_list(dropIfgram=False) + date_list = self.get_date_list(dropIfgram=dropIfgram) + num_date = len(date_list) + + # get the closure index + cp_idx = [] + for i in range(num_date - conn): + # compose the connection-n pairs + cp_date12_list = [] + for j in range(conn): + cp_date12_list.append('{}_{}'.format(date_list[i+j], date_list[i+j+1])) + cp_date12_list.append('{}_{}'.format(date_list[i], date_list[i+conn])) + + # add to cp_idx, ONLY IF all pairs exist for this closure loop + if all(x in date12_list for x in cp_date12_list): + cp_idx.append([date12_list.index(x) for x in cp_date12_list]) + + # list(list) to 2D array + cp_idx = np.array(cp_idx, dtype=np.int16) + cp_idx = np.unique(cp_idx, axis=0) + + return cp_idx + + + def get_sequential_closure_phase(self, box, conn, post_proc=None): + """Computes wrapped sequential closure phases for a given conneciton level. + + Reference: Equation (21) in Zheng et al. (2022, TGRS) + For conn = 5, seq_closure_phase = p12 + p23 + p34 + p45 + p56 - p16. + + Parameters: box - tuple of 4 int, bounding box in (x0, y0, x1, y1) + conn - int, connection level of the closure phase + post_proc - str, post processing of the closure phase: + None - 3D array in float32, seq closure phase + sum - 2D array in complex64, sum in time of the complex seq closure phase + mean - 2D array in complex64, mean in time of the complex seq closure phase + Returns: cp_w - 3D np.ndarray in float32 in size of (num_cp, box_len, box_wid) + wrapped sequential closure phase for the given connection level. + sum_cp - None or 2D np.ndarray in complex64 in size of (box_len, box_width) + wrapped average seq closure phase for the given connection level, + controlled by post_proc. + num_cp - int, number of seq closure phase for the given connection level. + """ + # basic info + num_date = len(self.get_date_list(dropIfgram=True)) + box_wid = box[2] - box[0] + box_len = box[3] - box[1] + + ## get the closure index + cp_idx = self.get_closure_phase_index(conn=conn, dropIfgram=True) + num_cp = cp_idx.shape[0] + print(f'number of closure measurements expected: {num_date - conn}') + print(f'number of closure measurements found : {num_cp}') + + if not post_proc: + if num_cp < num_date - conn: + msg = f'num_cp ({num_cp}) < num_date - conn ({num_date - conn})' + msg += ' --> some interferograms are missing!' + raise Exception(msg) + else: + if num_cp < 1: + raise Exception(f"No triplets found at connection level: {conn}!") + + ## read data + phase = self.read(box=box, print_msg=False) + ref_phase = self.get_reference_phase(dropIfgram=False) + for i in range(phase.shape[0]): + mask = phase[i] != 0. + phase[i][mask] -= ref_phase[i] + + ## calculate the 3D complex seq closure phase + cp_w = np.zeros((num_cp, box_len, box_wid), dtype=np.complex64) + for i in range(num_cp): + + # calculate closure phase - cp0_w + idx_plus, idx_minor = cp_idx[i, :-1], cp_idx[i, -1] + cp0_w = np.sum(phase[idx_plus], axis=0) - phase[idx_minor] + + # get the wrapped closure phase + cp_w[i] = np.exp(1j * cp0_w) + + ## post-processing + if not post_proc: + sum_cp = None + + elif post_proc == 'sum': + sum_cp = np.sum(cp_w, axis=0) + + elif post_proc == 'mean': + sum_cp = np.mean(cp_w, axis=0) + + else: + raise ValueError(f'un-recognized post_proc={post_proc}! Available choices: sum, mean.') + + return np.angle(cp_w), sum_cp, num_cp + + + # Functions for unwrapping error correction @staticmethod def get_design_matrix4triplet(date12_list): """Generate the design matrix of ifgram triangle for unwrap error correction using phase closure + Parameters: date12_list : list of string in YYYYMMDD_YYYYMMDD format Returns: C : 2D np.array in size of (num_tri, num_ifgram) consisting 0, 1, -1 for 3 SAR acquisition in t1, t2 and t3 in time order, @@ -1021,10 +1169,12 @@ def get_design_matrix4triplet(date12_list): return np.stack(C_list).astype(np.float32) - # Functions for Network Inversion + + # Functions for network inversion / time series estimation @staticmethod def get_design_matrix4timeseries(date12_list, refDate=None): """Return design matrix of the input ifgramStack for timeseries estimation + Parameters: date12_list - list of string in YYYYMMDD_YYYYMMDD format refDate - str, date in YYYYMMDD format set to None for the 1st date diff --git a/mintpy/utils/isce_utils.py b/mintpy/utils/isce_utils.py index 534fc418f..09f80a8b6 100644 --- a/mintpy/utils/isce_utils.py +++ b/mintpy/utils/isce_utils.py @@ -5,6 +5,7 @@ ############################################################ # 2020-07: Talib Oliver-Cabrera, add UAVSAR support w/in stripmapStack # 2020-10: Cunren Liang, add alosStack support +# 2022-06: Yujie Zheng, add standard processing from isce2 # Group contents: # metadata # geometry @@ -16,11 +17,15 @@ import os -import re +import datetime import glob import shelve -import datetime +import re +import time + import numpy as np +from scipy import ndimage + from mintpy.objects.constants import SPEED_OF_LIGHT, EARTH_RADIUS from mintpy.objects import sensor from mintpy.utils import ( @@ -585,7 +590,7 @@ def read_tops_baseline(baseline_file): def read_stripmap_baseline(baseline_file): """Read baseline file generated by ISCE/stripmapStack processor. - Example: + Example: baselines/20200111_20200125.txt PERP_BASELINE_BOTTOM 173.97914535263297 PERP_BASELINE_TOP 174.05612879066618 @@ -882,3 +887,212 @@ def get_sensing_datetime_list(proj_dir, date_list=None): return sensingMid, sensingStart, sensingStop + + +############################## Standard Processing ########################################### + +def gaussian_kernel(sx, sy, sig_x, sig_y): + '''Generate a guassian kernal (with all elements sum to 1). + + Parameters: sx/y - int, dimensions of kernal + sig_x/y - float, standard deviation of the guassian distribution + ''' + # ensure sx/y are odd number + sx += 1 if np.mod(sx, 2) == 0 else 0 + sy += 1 if np.mod(sy, 2) == 0 else 0 + + x, y = np.meshgrid(np.arange(sx), np.arange(sy)) + x += 1 + y += 1 + + xc = (sx + 1) / 2 + yc = (sy + 1) / 2 + fx = ((x-xc)**2.) / (2.*sig_x**2.) + fy = ((y-yc)**2.) / (2.*sig_y**2.) + + k = np.exp(-1.0 * (fx+fy)) + a = 1./np.sum(k) + k = a * k + + return k + + +def convolve(data, kernel): + '''Convolve / filter the complex data based on the given kernel. + + Parameters: data - 2D np.ndarray in complex + kernel - 2D np.ndarray in float, convolution kernel + ''' + R = ndimage.convolve(data.real, kernel, mode='constant', cval=0.0) + I = ndimage.convolve(data.imag, kernel, mode='constant', cval=0.0) + return R + 1J*I + + +def estimate_coherence(intfile, corfile): + '''Estimate the spatial coherence (phase sigma) of the wrapped interferogram. + + Parameters: intfile - str, path to the *.int file + corfile - str, path to the output correlation file + ''' + import isce + import isceobj + from mroipac.icu.Icu import Icu + + # create filt interferogram file object + filtImage = isceobj.createIntImage() + filtImage.load(intfile + '.xml') + filtImage.setAccessMode('read') + filtImage.createImage() + + # create phase sigma correlation file object + phsigImage = isceobj.createImage() + phsigImage.dataType='FLOAT' + phsigImage.bands = 1 + phsigImage.setWidth(filtImage.getWidth()) + phsigImage.setFilename(corfile) + phsigImage.setAccessMode('write') + phsigImage.createImage() + + # setup Icu() object + icuObj = Icu(name='sentinel_filter_icu') + icuObj.configure() + icuObj.unwrappingFlag = False + icuObj.useAmplitudeFlag = False + #icuObj.correlationType = 'NOSLOPE' + + # run + icuObj.icu(intImage=filtImage, phsigImage=phsigImage) + phsigImage.renderHdr() + + # close + filtImage.finalizeImage() + phsigImage.finalizeImage() + + return + + +def unwrap_snaphu(int_file, cor_file, unw_file, defo_max=2.0, max_comp=32, + init_only=True, init_method='MCF', cost_mode='SMOOTH'): + '''Unwrap interferograms using SNAPHU via isce2. + + Modified from ISCE-2/topsStack/unwrap.py + Notes from Piyush: + SNAPHU is an iterative solver, starting from the initial solution. It can get + stuck in an infinite loop. + The initial solution is created using MCF or MST method. The MST initial solution + typically require lots of iterations and may not be a good starting point. + DEFO cost mode requires geometry info for the program to interpret the coherence + correctly and setup costs based on that. DEFO always sounds more theoretical + to me, but I haven not fully explored it. TOPO cost mode requires spatial baseline. + SMOOTH cost mode is purely data driven. + Amplitude of zero is a mask in all cost modes. For TOPO mode, amplitude is used to find + layover; for SMOOTH mode, only non-zero amplitude matters. + + Default configurations in ISCE-2/topsStack: + init_only = True + init_method = 'MCF' + cost_mode = 'SMOOTH' + Default configurations in FRInGE: + init_only = False + init_method = 'MST' + cost_mode = 'DEFO' + + Parameters: int_file - str, path to the wrapped interferogram file + cor_file - str, path to the correlation file + unw_file - str, path to the output unwrapped interferogram file + defo_max - float, maximum number of cycles for the deformation phase + max_comp - int, maximum number of connected components + init_only - bool, initlize-only mode + init_method - str, algo used for initialization: MCF, MST + cost_mode - str, statistical-cost mode: TOPO, DEFO, SMOOTH, NOSTATCOSTS + Returns: unw_file - str, path to the output unwrapped interferogram file + ''' + import isce + from contrib.Snaphu.Snaphu import Snaphu + + start_time = time.time() + + # configurations - atr + atr = readfile.read_attribute(int_file) + width = int(atr['WIDTH']) + length = int(atr['LENGTH']) + altitude = float(atr['HEIGHT']) + earth_radius = float(atr['EARTH_RADIUS']) + wavelength = float(atr['WAVELENGTH']) + rg_looks = int(atr['RLOOKS']) + az_looks = int(atr['ALOOKS']) + corr_looks = float(atr.get('NCORRLOOKS', rg_looks * az_looks / 1.94)) + + ## setup SNAPHU + # https://web.stanford.edu/group/radar/softwareandlinks/sw/snaphu/snaphu.conf.full + # https://github.com/isce-framework/isce2/blob/main/contrib/Snaphu/Snaphu.py + print('phase unwrapping with SNAPHU ...') + print('SNAPHU cost mode: {}'.format(cost_mode)) + print('SNAPHU init only: {}'.format(init_only)) + print('SNAPHU init method: {}'.format(init_method)) + print('SNAPHU max number of connected components: {}'.format(max_comp)) + + snp = Snaphu() + + # file IO + snp.setInput(int_file) + snp.setOutput(unw_file) + snp.setCorrfile(cor_file) + snp.setCorFileFormat('FLOAT_DATA') + snp.setWidth(width) + + # runtime options + snp.setCostMode(cost_mode) + snp.setInitOnly(init_only) + snp.setInitMethod(init_method) + + # geometry parameters + # baseline info is not used in deformation mode, but is very important in topography mode + snp.setAltitude(altitude) + snp.setEarthRadius(earth_radius) + snp.setWavelength(wavelength) + snp.setRangeLooks(rg_looks) + snp.setAzimuthLooks(az_looks) + snp.setCorrLooks(corr_looks) + + # deformation mode parameters + snp.setDefoMaxCycles(defo_max) + + # connected component control + # grow connectedc components if init_only is True + # https://github.com/isce-framework/isce2/blob/main/contrib/Snaphu/Snaphu.py#L413 + snp.setMaxComponents(max_comp) + + ## run SNAPHU + snp.prepare() + snp.unwrap() + print('finished SNAPHU running') + + # mask out wired values from SNAPHU + # based on https://github.com/isce-framework/isce2/pull/326 + flag = np.fromfile(int_file, dtype=np.complex64).reshape(length, width) + data = np.memmap(unw_file, dtype='float32', mode='r+', shape=(length*2, width)) + data[0:length*2:2, :][np.nonzero(flag == 0)] = 0 + data[1:length*2:2, :][np.nonzero(flag == 0)] = 0 + + ## render metadata + print('write metadata file: {}.xml'.format(unw_file)) + atr['FILE_TYPE'] = '.unw' + atr['DATA_TYPE'] = 'float32' + atr['INTERLEAVE'] = 'BIL' + atr['BANDS'] = '2' + writefile.write_isce_xml(atr, unw_file) + + if snp.dumpConnectedComponents: + print('write metadata file: {}.conncomp.xml'.format(unw_file)) + atr['FILE_TYPE'] = '.conncomp' + atr['DATA_TYPE'] = 'uint8' + atr['INTERLEAVE'] = 'BIP' + atr['BANDS'] = '1' + writefile.write_isce_xml(atr, f'{unw_file}.conncomp') + + # time usage + m, s = divmod(time.time() - start_time, 60) + print('time used: {:02.0f} mins {:02.1f} secs.'.format(m, s)) + + return unw_file