network.py

"""Utilities for interferogram network selection."""
############################################################
# Program is part of MintPy                                #
# Copyright (c) 2013, Zhang Yunjun, Heresh Fattahi         #
# Author: Zhang Yunjun, 2016                               #
############################################################
# Recommend import:
#   from mintpy.utils import network as pnet


import itertools
import os
import sys

import numpy as np
from matplotlib import pyplot as plt
from matplotlib.tri import Triangulation
from scipy import sparse

from mintpy.objects import ifgramStack, sensor
from mintpy.utils import ptime, readfile

SPEED_OF_LIGHT = 299792458   # m/s, speed of light


##################################################################
BASELINE_LIST_FILE = """
# Date  Bperp    dop0/PRF  dop1/PRF   dop2/PRF   PRF    slcDir
070106     0.0   0.03      0.0000000  0.000000   2155.2 /KyushuT422F650AlosA/SLC/070106/
070709  2631.9   0.07      0.0000000  0.000000   2155.2 /KyushuT422F650AlosA/SLC/070709/
070824  2787.3   0.07      0.0000000  0.000000   2155.2 /KyushuT422F650AlosA/SLC/070824/
...
"""

IFGRAM_LIST_FILE = """
060713-070113
060828-070113
060828-070831
...
"""


##################################################################
def read_pairs_list(date12ListFile, dateList=[]):
    """Read Pairs List file like below:
    070311-070426
    070311-070611
    ...
    """
    # Read date12 list file
    date12List = sorted(list(np.loadtxt(date12ListFile, dtype=bytes).astype(str)))

    # Get dateList from date12List
    if not dateList:
        dateList = []
        for date12 in date12List:
            dates = date12.split('-')
            if not dates[0] in dateList:
                dateList.append(dates[0])
            if not dates[1] in dateList:
                dateList.append(dates[1])
        dateList.sort()
    date6List = ptime.yymmdd(dateList)

    # Get pair index
    pairs_idx = []
    for date12 in date12List:
        dates = date12.split('-')
        pair_idx = [date6List.index(dates[0]), date6List.index(dates[1])]
        pairs_idx.append(pair_idx)

    return pairs_idx


def write_pairs_list(pairs, dateList, outName):
    """Write pairs list file."""
    dateList6 = ptime.yymmdd(dateList)
    fl = open(outName, 'w')
    for idx in pairs:
        date12 = dateList6[idx[0]]+'-'+dateList6[idx[1]]+'\n'
        fl.write(date12)
    fl.close()
    return 1


def read_baseline_file(baselineFile, exDateList=[]):
    """Read bl_list.txt without dates listed in exDateList
    # Date  Bperp    dop0/PRF  dop1/PRF   dop2/PRF      PRF    slcDir
    070106     0.0   0.03      0.0000000  0.00000000000 2155.2 /scratch/KyushuT422F650AlosA/SLC/070106/
    070709  2631.9   0.07      0.0000000  0.00000000000 2155.2 /scratch/KyushuT422F650AlosA/SLC/070709/
    070824  2787.3   0.07      0.0000000  0.00000000000 2155.2 /scratch/KyushuT422F650AlosA/SLC/070824/
    ...

    Examples:
        date8List, perpBaseList, dopList, prfList, slcDirList = read_baseline_file(baselineFile)
        date8List, perpBaseList, dopList, prfList, slcDirList = read_baseline_file(baselineFile,['080520','100726'])
        date8List, perpBaseList = read_baseline_file(baselineFile)[0:2]
    """
    exDateList = ptime.yymmdd(exDateList)
    if not exDateList:
        exDateList = []

    # Read baseline file into lines
    fb = open(baselineFile)
    lines = []
    for line in fb:
        line = str.replace(line, '\n', '').strip()
        lines.append(line)
    fb.close()

    # Read each line and put the values into arrays
    date6List = []
    perpBaseList = []
    dopplerList = []
    slcDirList = []
    for line in lines:
        c = line.split()    # splits on white space
        date = c[0]
        if date not in exDateList:
            date6List.append(date)
            perpBaseList.append(float(c[1]))
            try:
                dop = np.array([float(c[2]), float(c[3]), float(c[4])])
                prf = float(c[5])
                dop *= prf
                dopplerList.append(dop)
            except:
                pass
            try:
                slcDirList.append(c[6])
            except:
                pass

    date8List = ptime.yyyymmdd(date6List)
    return date8List, perpBaseList, dopplerList, slcDirList


def date12_list2index(date12_list, date_list=[]):
    """Convert list of date12 string into list of index"""
    # Get dateList from date12List
    if not date_list:
        m_dates = [date12.split('-')[0] for date12 in date12_list]
        s_dates = [date12.split('-')[1] for date12 in date12_list]
        date_list = list(set(m_dates + s_dates))
    date6_list = ptime.yymmdd(sorted(ptime.yyyymmdd(date_list)))

    # Get pair index
    pairs_idx = []
    for date12 in date12_list:
        dates = date12.split('-')
        pair_idx = [date6_list.index(dates[0]), date6_list.index(dates[1])]
        pairs_idx.append(pair_idx)

    return pairs_idx


def get_date12_list(fname, dropIfgram=False):
    """Read date12 info from input file: Pairs.list or multi-group hdf5 file
    Parameters: fname       - string, path/name of input multi-group hdf5 file or text file
                dropIfgram  - bool, check the "dropIfgram" dataset in ifgramStack hdf5 file
    Returns:    date12_list - list of string in YYYYMMDD_YYYYMMDD format
    Example:
        date12List = get_date12_list('ifgramStack.h5')
        date12List = get_date12_list('ifgramStack.h5', dropIfgram=True)
        date12List = get_date12_list('coherenceSpatialAvg.txt')
    """
    date12_list = []
    ext = os.path.splitext(fname)[1].lower()
    if ext == '.h5':
        ftype = readfile.read_attribute(fname)['FILE_TYPE']
        if ftype == 'ifgramStack':
            date12_list = ifgramStack(fname).get_date12_list(dropIfgram=dropIfgram)
        else:
            return None
    else:
        date12_list = np.loadtxt(fname, dtype=bytes, usecols=0).astype(str).tolist()
        # for txt file with only one interferogram
        if isinstance(date12_list, str):
            date12_list = [date12_list]

    date12_list = sorted(date12_list)
    date12_list = ptime.yyyymmdd_date12(date12_list)
    return date12_list


def critical_perp_baseline(sensor_name, inc_angle, print_msg=False):
    """Critical Perpendicular Baseline for each satellite"""
    # Jers: 5.712e3 m (near_range=688849.0551m)
    # Alos: 6.331e3 m (near_range=842663.2917m)
    # Tsx : 8.053e3 m (near_range=634509.1271m)
    sensor_dict = sensor.SENSOR_DICT[sensor_name.lower()]
    wvl = SPEED_OF_LIGHT / sensor_dict['carrier_frequency']
    near_range = 688849  # Yunjun 5/2016, case for Jers, need a automatic way to get this number
    rg_bandwidth = sensor_dict['chirp_bandwidth']
    bperp_c = wvl * (rg_bandwidth / SPEED_OF_LIGHT) * near_range * np.tan(inc_angle * np.pi / 180.0)
    if print_msg:
        print(f'Critical Perpendicular Baseline: {bperp_c} m')
    return bperp_c


def calculate_doppler_overlap(dop_a, dop_b, bandwidth_az):
    """Calculate Overlap Percentage of Doppler frequency in azimuth direction
    Inputs:
        dop_a/b      : np.array of 3 floats, doppler frequency
        bandwidth_az : float, azimuth bandwidth
    Output:
        dop_overlap  : float, doppler frequency overlap between a & b.
    """
    # Calculate mean Doppler difference between a and b
    no_of_rangepix = 5000
    ddiff = np.zeros(10)
    for i in range(10):
        rangepix = (i-1)*no_of_rangepix/10 + 1
        da = dop_a[0]+(rangepix-1)*dop_a[1]+(rangepix-1)**2*dop_a[2]
        db = dop_b[0]+(rangepix-1)*dop_b[1]+(rangepix-1)**2*dop_b[2]
        ddiff[i] = abs(da - db)
    ddiff_mean = np.mean(ddiff)

    #dopOverlap_prf = bandwidth_az - ddiff
    #dopOverlap_percent = np.mean(dopOverlap_prf / bandwidth_az * 100)
    ##overlapLessthanZero = np.nonzero(dopOverlap_percent < 0)
    #dopOverlap_percent = np.mean(dopOverlap_percent)
    # return dopOverlap_percent

    # Doppler overlap
    dop_overlap = (bandwidth_az - ddiff_mean) / bandwidth_az
    return dop_overlap


def simulate_coherence_v2(date12_list, decor_time=200.0, coh_resid=0.2, inc_angle=40, sensor_name='Sen',
                          display=False):
    """Simulate coherence version 2 (without using bl_list.txt file).
    Parameters: date12_list - list of string in YYYYMMDD_YYYYMMDD format, indicating pairs configuration
                decor_time  - float, decorrelation rate in days, time for coherence to drop to 1/e of its initial value
                coh_resid   - float, long-term coherence, minimum attainable coherence value
                inc_angle   - float, incidence angle in degrees
                sensor_name - string, SAR sensor name
                display     - bool, display result as matrix or not
    Returns:    coh         - 2D np.array in size of (ifgram_num)
    """
    num_pair = len(date12_list)
    date1s = [x.split('_')[0] for x in date12_list]
    date2s = [x.split('_')[1] for x in date12_list]
    date_list = sorted(list(set(date1s + date2s)))
    tbase_list = ptime.date_list2tbase(date_list)[0]

    SNR = 22  # NESZ = -22 dB from Table 1 in https://sentinels.copernicus.eu/web/sentinel/
    coh_thermal = 1. / (1. + 1./SNR)

    # bperp
    rng = np.random.default_rng(2)
    pbase_list = rng.normal(0, 50, num_pair).tolist()
    pbase_c = critical_perp_baseline(sensor_name, inc_angle)

    coh = np.zeros(num_pair, dtype=np.float32)
    for i in range(num_pair):
        date1, date2 = date12_list[i].split('_')
        ind1, ind2 = date_list.index(date1), date_list.index(date2)
        tbase = tbase_list[ind2] - tbase_list[ind1]
        pbase = pbase_list[ind2] - pbase_list[ind1]

        coh_geom = (pbase_c - abs(pbase)) / pbase_c
        coh_temp = np.multiply((coh_thermal - coh_resid), np.exp(-1*abs(tbase)/decor_time)) + coh_resid
        coh[i] = coh_geom * coh_temp

    if display:
        print('critical perp baseline: %.f m' % pbase_c)
        coh_mat = coherence_matrix(date12_list, coh)
        plt.figure()
        plt.imshow(coh_mat, vmin=0.0, vmax=1.0, cmap='jet')
        plt.xlabel('Image number')
        plt.ylabel('Image number')
        cbar = plt.colorbar()
        cbar.set_label('Coherence')
        plt.title('Coherence matrix')
        plt.show()

    return coh


def simulate_coherence(date12_list, baseline_file='bl_list.txt', sensor_name='Env', inc_angle=22.8,
                       decor_time=200.0, coh_resid=0.2, display=False):
    """Simulate coherence for a given set of interferograms
    Inputs:
        date12_list  - list of string in YYMMDD-YYMMDD format, indicating pairs configuration
        baseline_file - string, path of baseline list text file
        sensor_name     - string, SAR sensor name
        inc_angle  - float, incidence angle
        decor_time - float / 2D np.array in size of (1, pixel_num)
                     decorrelation rate in days, time for coherence to drop to 1/e of its initial value
        coh_resid  - float / 2D np.array in size of (1, pixel_num)
                     long-term coherence, minimum attainable coherence value
        display    - bool, display result as matrix or not
    Output:
        cohs       - 2D np.array in size of (ifgram_num, pixel_num)
    Example:
        date12_list = pnet.get_date12_list('ifgram_list.txt')
        cohs = simulate_coherences(date12_list, 'bl_list.txt', sensor_name='Tsx')

    References:
        Guarnieri, A. M. (2013), Introduction to RADAR, Politecnico di Milano Dipartimento di Elettronica
            e Informazione, Milano.
        Zebker, H. A., & Villasenor, J. (1992). Decorrelation in interferometric radar echoes.
            IEEE-TGRS, 30(5), 950-959.
        Hanssen, R. F. (2001). Radar interferometry: data interpretation and error analysis
            (Vol. 2). Dordrecht, Netherlands: Kluwer Academic Pub.
        Morishita, Y., & Hanssen, R. F. (2015). Temporal decorrelation in L-, C-, and X-band satellite
            radar interferometry for pasture on drained peat soils. IEEE-TGRS, 53(2), 1096-1104.
        Parizzi, A., Cong, X., & Eineder, M. (2009). First Results from Multifrequency Interferometry.
            A comparison of different decorrelation time constants at L, C, and X Band. ESA Scientific
            Publications(SP-677), 1-5.
    """
    date_list, pbase_list, dop_list = read_baseline_file(baseline_file)[0:3]
    tbase_list = ptime.date_list2tbase(date_list)[0]

    # Thermal decorrelation (Zebker and Villasenor, 1992, Eq.4)
    SNR = 19.5  # hardwired for Envisat (Guarnieri, 2013)
    coh_thermal = 1. / (1. + 1./SNR)

    pbase_c = critical_perp_baseline(sensor_name, inc_angle)
    bandwidth_az = sensor.SENSOR_DICT[sensor_name.lower()]['doppler_bandwidth']

    date12_list = ptime.yyyymmdd_date12(date12_list)
    ifgram_num = len(date12_list)

    if isinstance(decor_time, (int, np.int16, np.int32, float, np.float32, np.float64)):
        pixel_num = 1
        decor_time = float(decor_time)
    else:
        pixel_num = decor_time.shape[1]
    if decor_time == 0.:
        decor_time = 0.01
    cohs = np.zeros((ifgram_num, pixel_num), np.float32)
    for i in range(ifgram_num):
        if display:
            sys.stdout.write('\rinterferogram = %4d/%4d' % (i, ifgram_num))
            sys.stdout.flush()
        m_date, s_date = date12_list[i].split('_')
        m_idx = date_list.index(m_date)
        s_idx = date_list.index(s_date)

        pbase = pbase_list[s_idx] - pbase_list[m_idx]
        tbase = tbase_list[s_idx] - tbase_list[m_idx]

        # Geometric decorrelation (Hanssen, 2001, Eq. 4.4.12)
        coh_geom = (pbase_c - abs(pbase)) / pbase_c
        if coh_geom < 0.:
            coh_geom = 0.

        # Doppler centroid decorrelation (Hanssen, 2001, Eq. 4.4.13)
        if not dop_list:
            coh_dc = 1.
        else:
            coh_dc = calculate_doppler_overlap(dop_list[m_idx],
                                               dop_list[s_idx],
                                               bandwidth_az)
            if coh_dc < 0.:
                coh_dc = 0.

        # Option 1: Temporal decorrelation - exponential delay model (Parizzi et al., 2009; Morishita and Hanssen, 2015)
        coh_temp = np.multiply((coh_thermal - coh_resid), np.exp(-1*abs(tbase)/decor_time)) + coh_resid

        coh = coh_geom * coh_dc * coh_temp
        cohs[i, :] = coh
    #epsilon = 1e-3
    #cohs[cohs < epsilon] = epsilon
    if display:
        print('')

    if display:
        print('critical perp baseline: %.f m' % pbase_c)
        cohs_mat = coherence_matrix(date12_list, cohs)
        plt.figure()
        plt.imshow(cohs_mat, vmin=0.0, vmax=1.0, cmap='jet')
        plt.xlabel('Image number')
        plt.ylabel('Image number')
        cbar = plt.colorbar()
        cbar.set_label('Coherence')
        plt.title('Coherence matrix')
        plt.show()
    return cohs


##################################################################
def threshold_doppler_overlap(date12_list, date_list, dop_list, bandwidth_az, dop_overlap_min=0.15):
    """Remove pairs/interoferogram with doppler overlap larger than critical value
    Inputs:
        date12_list : list of string, for date12 in YYMMDD-YYMMDD format
        date_list   : list of string, for date in YYMMDD/YYYYMMDD format, optional
        dop_list    : list of list of 3 float, for centroid Doppler frequency
        bandwidth_az    : float, bandwidth in azimuth direction
        dop_overlap_min : float, minimum overlap of azimuth Doppler frequency
    Outputs:
        date12_list : list of string, for date12 in YYMMDD-YYMMDD format
    """
    if not date12_list:
        return []
    # Get date6_list
    if not date_list:
        m_dates = [date12.split('-')[0] for date12 in date12_list]
        s_dates = [date12.split('-')[1] for date12 in date12_list]
        date_list = sorted(ptime.yyyymmdd(list(set(m_dates + s_dates))))
    date6_list = ptime.yymmdd(date_list)

    # Threshold
    date12_list_out = []
    for date12 in date12_list:
        date1, date2 = date12.split('-')
        idx1 = date6_list.index(date1)
        idx2 = date6_list.index(date2)
        dop_overlap = calculate_doppler_overlap(dop_list[idx1],
                                                dop_list[idx2],
                                                bandwidth_az)
        if dop_overlap >= dop_overlap_min:
            date12_list_out.append(date12)
    return date12_list_out


def threshold_perp_baseline(date12_list, date_list, pbase_list, pbase_max, pbase_min=0.0):
    """Remove pairs/interoferogram out of [pbase_min, pbase_max]
    Inputs:
        date12_list : list of string for date12 in YYMMDD-YYMMDD format
        date_list   : list of string for date in YYMMDD/YYYYMMDD format, optional
        pbase_list  : list of float for perpendicular spatial baseline
        pbase_max   : float, maximum perpendicular baseline
        pbase_min   : float, minimum perpendicular baseline
    Output:
        date12_list_out : list of string for date12 in YYMMDD-YYMMDD format
    Example:
        date12_list = threshold_perp_baseline(date12_list, date_list, pbase_list, 500)
    """
    if not date12_list:
        return []

    # Get date_list
    delimiter = [i for i in ['-', '_'] if i in date12_list[0]][0]
    if not date_list:
        m_dates = [date12.split(delimiter)[0] for date12 in date12_list]
        s_dates = [date12.split(delimiter)[1] for date12 in date12_list]
        date_list = sorted(list(set(m_dates + s_dates)))

    if not len(date_list) == len(pbase_list):
        print('ERROR: number of existing dates is not equal to number of perp baseline!')
        print('date list is needed for threshold filtering!')
        print('skip filtering.')
        return date12_list

    # Threshold
    date12_list_out = []
    for date12 in date12_list:
        date1, date2 = date12.split(delimiter)
        pbase1 = pbase_list[date_list.index(date1)]
        pbase2 = pbase_list[date_list.index(date2)]
        pbase = abs(pbase2 - pbase1)

        if pbase_min <= pbase <= pbase_max:
            date12_list_out.append(date12)

    return date12_list_out


def threshold_temporal_baseline(date12_list, btemp_max, keep_seasonal=True, btemp_min=0.0):
    """Remove pairs/interferograms out of min/max/seasonal temporal baseline limits
    Inputs:
        date12_list : list of string for date12 in YYMMDD-YYMMDD or YYYYMMDD_YYYYMMDD format
        btemp_max   : float, maximum temporal baseline
        btemp_min   : float, minimum temporal baseline
        keep_seasonal : keep interferograms with seasonal temporal baseline
    Output:
        date12_list_out : list of string for date12
    Example:
        date12_list = threshold_temporal_baseline(date12_list, 200)
        date12_list = threshold_temporal_baseline(date12_list, 200, False)
    """
    if not date12_list:
        return []

    # Get date list and tbase list
    delimiter = [i for i in ['-', '_'] if i in date12_list[0]][0]
    m_dates = [date12.split(delimiter)[0] for date12 in date12_list]
    s_dates = [date12.split(delimiter)[1] for date12 in date12_list]
    date_list = sorted(list(set(m_dates + s_dates)))
    tbase_list = ptime.date_list2tbase(date_list)[0]

    # Threshold
    date12_list_out = []
    for date12 in date12_list:
        date1, date2 = date12.split(delimiter)
        tbase1 = tbase_list[date_list.index(date1)]
        tbase2 = tbase_list[date_list.index(date2)]
        tbase = int(abs(tbase2 - tbase1))

        if btemp_min <= tbase <= btemp_max:
            date12_list_out.append(date12)

        elif keep_seasonal and tbase/30 in [11, 12]:
            date12_list_out.append(date12)

    return date12_list_out


def coherence_matrix(date12_list, coh_list, diag_value=np.nan, fill_triangle='both', date_list=None):
    """Return coherence matrix based on input date12 list and its coherence
    Inputs:
        date12_list - list of string in YYMMDD-YYMMDD format
        coh_list    - list of float, average coherence for each interferograms
        diag_value  - number, value to be filled in the diagonal
        fill_triangle - str, 'both', 'upper', 'lower'
    Output:
        coh_matrix  - 2D np.array with dimension length = date num
                      np.nan value for interferograms non-existed.
                      1.0 for diagonal elements
    """
    # Get date list
    date12_list = ptime.yymmdd_date12(date12_list)
    if not date_list:
        m_dates = [date12.split('-')[0] for date12 in date12_list]
        s_dates = [date12.split('-')[1] for date12 in date12_list]
        date_list = sorted(list(set(m_dates + s_dates)))
    date_list = ptime.yymmdd(date_list)
    date_num = len(date_list)

    coh_mat = np.zeros([date_num, date_num])
    coh_mat[:] = np.nan
    for date12 in date12_list:
        date1, date2 = date12.split('-')
        idx1 = date_list.index(date1)
        idx2 = date_list.index(date2)
        coh = coh_list[date12_list.index(date12)]
        if fill_triangle in ['upper', 'both']:
            coh_mat[idx1, idx2] = coh  # symmetric
        if fill_triangle in ['lower', 'both']:
            coh_mat[idx2, idx1] = coh

    if diag_value is not np.nan:
        for i in range(date_num):    # diagonal value
            coh_mat[i, i] = diag_value
    return coh_mat


def threshold_coherence_based_mst(date12_list, coh_list):
    """Return a minimum spanning tree of network based on the coherence inverse.
    Inputs:
        date12_list - list of string in YYMMDD-YYMMDD format
        coh_list    - list of float, average coherence for each interferogram
    Output:
        mst_date12_list - list of string in YYMMDD-YYMMDD format, for MST network of interferograms
    """
    # coh_list --> coh_mat --> weight_mat
    coh_mat = coherence_matrix(date12_list, coh_list)
    mask = ~np.isnan(coh_mat)
    wei_mat = np.zeros(coh_mat.shape)
    wei_mat[:] = np.inf
    wei_mat[mask] = 1/coh_mat[mask]

    # MST path based on weight matrix
    wei_mat_csr = sparse.csr_matrix(wei_mat)
    mst_mat_csr = sparse.csgraph.minimum_spanning_tree(wei_mat_csr)

    # Get date6_list
    date12_list = ptime.yymmdd_date12(date12_list)
    m_dates = [date12.split('-')[0] for date12 in date12_list]
    s_dates = [date12.split('-')[1] for date12 in date12_list]
    date6_list = ptime.yymmdd(sorted(ptime.yyyymmdd(list(set(m_dates + s_dates)))))

    # Convert MST index matrix into date12 list
    [s_idx_list, m_idx_list] = [date_idx_array.tolist()
                                for date_idx_array in sparse.find(mst_mat_csr)[0:2]]
    mst_date12_list = []
    for m_idx, s_idx in zip(m_idx_list, s_idx_list):
        idx = sorted([m_idx, s_idx])
        date12 = date6_list[idx[0]]+'-'+date6_list[idx[1]]
        mst_date12_list.append(date12)
    return mst_date12_list


def pair_sort(pairs):
    for i, pair in enumerate(pairs):
        if pair[0] > pair[1]:
            pairs[i][0] = pair[1]
            pairs[i][1] = pair[0]
    pairs = sorted(pairs)
    return pairs


def pair_merge(pairs1, pairs2):
    pairs = pairs1
    for pair in pairs2:
        if pair not in pairs:
            pairs.append(pair)

    pairs = sorted(pairs)
    return pairs


def select_pairs_all(date_list, date_format='YYMMDD'):
    """Select All Possible Pairs/Interferograms
    Input : date_list   - list of date in YYMMDD/YYYYMMDD format
    Output: date12_list - list date12 in YYMMDD-YYMMDD format
    Reference:
        Berardino, P., G. Fornaro, R. Lanari, and E. Sansosti (2002), A new algorithm for surface deformation monitoring
        based on small baseline differential SAR interferograms, IEEE TGRS, 40(11), 2375-2383.
    """
    date8_list = sorted(ptime.yyyymmdd(date_list))
    date6_list = ptime.yymmdd(date8_list)
    date12_list = list(itertools.combinations(date6_list, 2))
    date12_list = [date12[0]+'-'+date12[1] for date12 in date12_list]
    if date_format == 'YYYYMMDD':
        date12_list = ptime.yyyymmdd_date12(date12_list)
    return date12_list


def select_pairs_sequential(date_list, num_conn=2, date_format=None):
    """Select Pairs in a Sequential way:
        For each acquisition, find its num_connection nearest acquisitions in the past time.

    Parameters: date_list   - list of str for date
                num_conn    - int, number of sequential connections
                date_format - str / None, output date format
    Returns:    date12_list - list of str for date12
    """

    date_list = sorted(date_list)
    date_inds = list(range(len(date_list)))

    # Get pairs index list
    date12_inds = []
    for date_ind in date_inds:
        for i in range(num_conn):
            if date_ind-i-1 >= 0:
                date12_inds.append([date_ind-i-1, date_ind])
    date12_inds = [sorted(i) for i in sorted(date12_inds)]

    # Convert index into date12
    date12_list = [f'{date_list[ind12[0]]}_{date_list[ind12[1]]}'
                  for ind12 in date12_inds]

    # adjust output date format
    if date_format is not None:
        if date_format == 'YYYYMMDD':
            date12_list = ptime.yyyymmdd_date12(date12_list)
        elif date_format == 'YYMMDD':
            date12_list = ptime.yymmdd_date12(date12_list)
        else:
            raise ValueError(f'un-supported date format: {date_format}!')

    return date12_list


def select_pairs_hierarchical(date_list, pbase_list, temp_perp_list, date_format='YYMMDD'):
    """Select Pairs in a hierarchical way using list of temporal and perpendicular baseline thresholds
        For each temporal/perpendicular combination, select all possible pairs; and then merge all combination results
        together for the final output (Zhao, 2015).
    Inputs:
        date_list  : list of date in YYMMDD/YYYYMMDD format
        pbase_list : list of float, perpendicular spatial baseline
        temp_perp_list : list of list of 2 floats, for list of temporal/perp baseline, e.g.
                         [[32.0, 800.0], [48.0, 600.0], [64.0, 200.0]]
    Examples:
        pairs = select_pairs_hierarchical(date_list, pbase_list, [[32.0, 800.0], [48.0, 600.0], [64.0, 200.0]])
    Reference:
        Zhao, W., (2015), Small deformation detected from InSAR time-series and their applications in geophysics, Doctoral
        dissertation, Univ. of Miami, Section 6.3.
    """
    # Get all date12
    date12_list_all = select_pairs_all(date_list)

    # Loop of Threshold
    print('List of temporal and perpendicular spatial baseline thresholds:')
    print(temp_perp_list)
    date12_list = []
    for temp_perp in temp_perp_list:
        tbase_max = temp_perp[0]
        pbase_max = temp_perp[1]
        date12_list_tmp = threshold_temporal_baseline(date12_list_all,
                                                      tbase_max,
                                                      keep_seasonal=False)
        date12_list_tmp = threshold_perp_baseline(date12_list_tmp,
                                                  date_list,
                                                  pbase_list,
                                                  pbase_max)
        date12_list += date12_list_tmp
    date12_list = sorted(list(set(date12_list)))
    if date_format == 'YYYYMMDD':
        date12_list = ptime.yyyymmdd_date12(date12_list)
    return date12_list


def select_pairs_delaunay(date_list, pbase_list, norm=True, date_format='YYMMDD'):
    """Select Pairs using Delaunay Triangulation based on temporal/perpendicular baselines
    Inputs:
        date_list  : list of date in YYMMDD/YYYYMMDD format
        pbase_list : list of float, perpendicular spatial baseline
        norm       : normalize temporal baseline to perpendicular baseline
    Key points
        1. Define a ratio between perpendicular and temporal baseline axis units (Pepe and Lanari, 2006, TGRS).
        2. Pairs with too large perpendicular / temporal baseline or Doppler centroid difference should be removed
           after this, using a threshold, to avoid strong decorrelations (Zebker and Villasenor, 1992, TGRS).
    Reference:
        Pepe, A., and R. Lanari (2006), On the extension of the minimum cost flow algorithm for phase unwrapping
        of multitemporal differential SAR interferograms, IEEE TGRS, 44(9), 2374-2383.
        Zebker, H. A., and J. Villasenor (1992), Decorrelation in interferometric radar echoes, IEEE TGRS, 30(5), 950-959.
    """
    # Get temporal baseline in days
    date6_list = ptime.yymmdd(date_list)
    date8_list = ptime.yyyymmdd(date_list)
    tbase_list = ptime.date_list2tbase(date8_list)[0]

    # Normalization (Pepe and Lanari, 2006, TGRS)
    if norm:
        temp2perp_scale = (max(pbase_list)-min(pbase_list)) / (max(tbase_list)-min(tbase_list))
        tbase_list = [tbase*temp2perp_scale for tbase in tbase_list]

    # Generate Delaunay Triangulation
    date12_idx_list = Triangulation(tbase_list, pbase_list).edges.tolist()
    date12_idx_list = [sorted(idx) for idx in sorted(date12_idx_list)]

    # Convert index into date12
    date12_list = [date6_list[idx[0]]+'-'+date6_list[idx[1]]
                   for idx in date12_idx_list]
    if date_format == 'YYYYMMDD_YYYYMMDD':
        date12_list = ptime.yyyymmdd_date12(date12_list)
    return date12_list


def select_pairs_mst(date_list, pbase_list, date_format='YYMMDD'):
    """Select Pairs using Minimum Spanning Tree technique
        Connection Cost is calculated using the baseline distance in perp and scaled temporal baseline (Pepe and Lanari,
        2006, TGRS) plane.
    Inputs:
        date_list  : list of date in YYMMDD/YYYYMMDD format
        pbase_list : list of float, perpendicular spatial baseline
    References:
        Pepe, A., and R. Lanari (2006), On the extension of the minimum cost flow algorithm for phase unwrapping
        of multitemporal differential SAR interferograms, IEEE TGRS, 44(9), 2374-2383.
        Perissin D., Wang T. (2012), Repeat-pass SAR interferometry with partially coherent targets. IEEE TGRS. 271-280
    """
    # Get temporal baseline in days
    date6_list = ptime.yymmdd(date_list)
    date8_list = ptime.yyyymmdd(date_list)
    tbase_list = ptime.date_list2tbase(date8_list)[0]
    # Normalization (Pepe and Lanari, 2006, TGRS)
    temp2perp_scale = (max(pbase_list)-min(pbase_list)) / (max(tbase_list)-min(tbase_list))
    tbase_list = [tbase*temp2perp_scale for tbase in tbase_list]

    # Get weight matrix
    ttMat1, ttMat2 = np.meshgrid(np.array(tbase_list), np.array(tbase_list))
    ppMat1, ppMat2 = np.meshgrid(np.array(pbase_list), np.array(pbase_list))
    ttMat = np.abs(ttMat1 - ttMat2)  # temporal distance matrix
    ppMat = np.abs(ppMat1 - ppMat2)  # spatial distance matrix

    # 2D distance matrix in temp/perp domain
    weightMat = np.sqrt(np.square(ttMat) + np.square(ppMat))
    weightMat = sparse.csr_matrix(weightMat)  # compress sparse row matrix

    # MST path based on weight matrix
    mstMat = sparse.csgraph.minimum_spanning_tree(weightMat)

    # Convert MST index matrix into date12 list
    [s_idx_list, m_idx_list] = [date_idx_array.tolist()
                                for date_idx_array in sparse.find(mstMat)[0:2]]
    date12_list = []
    for i in range(len(m_idx_list)):
        idx = sorted([m_idx_list[i], s_idx_list[i]])
        date12 = date6_list[idx[0]]+'-'+date6_list[idx[1]]
        date12_list.append(date12)
    if date_format == 'YYYYMMDD':
        date12_list = ptime.yyyymmdd_date12(date12_list)
    return date12_list


def select_pairs_star(date_list, m_date=None, pbase_list=[], date_format='YYMMDD'):
    """Select Star-like network/interferograms/pairs, it's a single reference network, similar to PS approach.
    Usage:
        m_date : reference date, choose it based on the following cretiria:
                 1) near the center in temporal and spatial baseline
                 2) prefer winter season than summer season for less temporal decorrelation
    Reference:
        Ferretti, A., C. Prati, and F. Rocca (2001), Permanent scatterers in SAR interferometry, IEEE TGRS, 39(1), 8-20.
    """
    date8_list = sorted(ptime.yyyymmdd(date_list))
    date6_list = ptime.yymmdd(date8_list)

    # Select reference date if not chosen
    if not m_date:
        m_date = select_reference_date(date8_list, pbase_list)
        print('auto select reference date: '+m_date)

    # Check input reference date
    m_date8 = ptime.yyyymmdd(m_date)
    if m_date8 not in date8_list:
        print('Input reference date does not exist in date list!')
        print(f'Input reference date: {m_date8}')
        print(f'Input date list: {date8_list}')
        m_date8 = None

    # Generate star/ps network
    m_idx = date8_list.index(m_date8)
    date12_idx_list = [sorted([m_idx, s_idx]) for s_idx in range(len(date8_list))
                       if s_idx is not m_idx]
    date12_list = [date6_list[idx[0]]+'-'+date6_list[idx[1]]
                   for idx in date12_idx_list]
    if date_format == 'YYYYMMDD':
        date12_list = ptime.yyyymmdd_date12(date12_list)
    return date12_list


def select_reference_date(date_list, pbase_list=[]):
    """Select super reference date based on input temporal and/or perpendicular baseline info.
    Return reference date in YYYYMMDD format.
    """
    date8_list = ptime.yyyymmdd(date_list)
    if not pbase_list:
        # Choose date in the middle
        m_date8 = date8_list[int(len(date8_list)/2)]
    else:
        # Get temporal baseline list
        tbase_list = ptime.date_list2tbase(date8_list)[0]
        # Normalization (Pepe and Lanari, 2006, TGRS)
        temp2perp_scale = (max(pbase_list)-min(pbase_list)) / (max(tbase_list)-min(tbase_list))
        tbase_list = [tbase*temp2perp_scale for tbase in tbase_list]
        # Get distance matrix
        ttMat1, ttMat2 = np.meshgrid(np.array(tbase_list),
                                     np.array(tbase_list))
        ppMat1, ppMat2 = np.meshgrid(np.array(pbase_list),
                                     np.array(pbase_list))
        ttMat = np.abs(ttMat1 - ttMat2)  # temporal distance matrix
        ppMat = np.abs(ppMat1 - ppMat2)  # spatial distance matrix
        # 2D distance matrix in temp/perp domain
        disMat = np.sqrt(np.square(ttMat) + np.square(ppMat))

        # Choose date minimize the total distance of temp/perp baseline
        disMean = np.mean(disMat, 0)
        m_idx = np.argmin(disMean)
        m_date8 = date8_list[m_idx]
    return m_date8


def select_reference_interferogram(date12_list, date_list, pbase_list, m_date=None):
    """Select reference interferogram based on input temp/perp baseline info
    If m_date is specified, select its closest s_date, which is newer than m_date;
        otherwise, choose the closest pair among all pairs as reference interferogram.
    Example:
        m_date12   = pnet.select_reference_ifgram(date12_list, date_list, pbase_list)
        '080211-080326' = pnet.select_reference_ifgram(date12_list, date_list, pbase_list, m_date='080211')
    """
    pbase_array = np.array(pbase_list, dtype='float64')
    # Get temporal baseline
    date8_list = ptime.yyyymmdd(date_list)
    date6_list = ptime.yymmdd(date8_list)
    tbase_array = np.array(ptime.date_list2tbase(date8_list)[0], dtype='float64')
    # Normalization (Pepe and Lanari, 2006, TGRS)
    temp2perp_scale = (max(pbase_array)-min(pbase_array)) / (max(tbase_array)-min(tbase_array))
    tbase_array *= temp2perp_scale

    # Calculate sqrt of temp/perp baseline for input pairs
    idx1 = np.array([date6_list.index(date12.split('-')[0]) for date12 in date12_list])
    idx2 = np.array([date6_list.index(date12.split('-')[1]) for date12 in date12_list])
    base_distance = np.sqrt((tbase_array[idx2] - tbase_array[idx1])**2 +
                            (pbase_array[idx2] - pbase_array[idx1])**2)

    # Get reference interferogram index
    if not m_date:
        # Choose pair with shortest temp/perp baseline
        m_date12_idx = np.argmin(base_distance)
    else:
        m_date = ptime.yymmdd(m_date)
        # Choose pair contains m_date with shortest temp/perp baseline
        m_date12_idx_array = np.array([date12_list.index(date12) for date12 in date12_list
                                       if m_date+'-' in date12])
        min_base_distance = np.min(base_distance[m_date12_idx_array])
        m_date12_idx = np.where(base_distance == min_base_distance)[0][0]

    m_date12 = date12_list[m_date12_idx]
    return m_date12