From b7b5cd7fd5c20599d7c67b563160d71e6085fc23 Mon Sep 17 00:00:00 2001
From: Yohcy <865533709@qq.com>
Date: Tue, 22 Oct 2024 15:40:00 +0800
Subject: [PATCH 01/12] tropo_ver1

---
 src/mintpy/cli/tropo_HTC.py    |  82 +++++++++
 src/mintpy/smallbaselineApp.py |  22 ++-
 src/mintpy/tropo_HTC.py        | 293 +++++++++++++++++++++++++++++++++
 3 files changed, 395 insertions(+), 2 deletions(-)
 create mode 100644 src/mintpy/cli/tropo_HTC.py
 create mode 100644 src/mintpy/tropo_HTC.py

diff --git a/src/mintpy/cli/tropo_HTC.py b/src/mintpy/cli/tropo_HTC.py
new file mode 100644
index 000000000..61ab2ad51
--- /dev/null
+++ b/src/mintpy/cli/tropo_HTC.py
@@ -0,0 +1,82 @@
+#!/usr/bin/env python3
+############################################################
+# Program is part of MintPy                                #
+# Copyright (c) 2013, Zhang Yunjun, Heresh Fattahi         #
+# Author: Antonio Valentino, Zhang Yunjun, Aug 2022        #
+############################################################
+
+import argparse
+import sys
+
+from mintpy.utils.arg_utils import create_argument_parser
+
+############################################################################
+REFERENCE = """reference:
+    Yang, Q., Yunjun, Z., Wang, R. Y. Heterogeneous InSAR Tropospheric Correction 
+    Based on Local Texture Correlation, IEEE Transactions on Geoscience and Remote Sensing, 62,
+      doi:10.1109/TGRS.2024.3356749.
+"""
+#TODO
+EXAMPLE = """example:
+  tropo_HTC.py  timeseries_ramp_demErr.h5  -v velocity.h5  -g inputs/geometryRadar.h5  -m maskTempCoh.h5
+  tropo_HTC .py  geo_timeseries_demErr.h5  -g geo_geometryRadar.h5     -m geo_maskTempCoh.h5
+"""
+
+def create_parser(subparsers=None):
+    #TODO
+    synopsis = 'Correct Topo-correlated Stratified tropospheric delay'
+    epilog = REFERENCE + '\n' + EXAMPLE
+    name = __name__.split('.')[-1]
+    parser = create_argument_parser(
+        name, synopsis=synopsis, description=synopsis, epilog=epilog, subparsers=subparsers)
+    
+    parser.add_argument('timeseries_file', help='time-series file to be corrected')
+    parser.add_argument('-g', '--geometry', dest='geom_file', required=True,
+                        help='DEM file used for correlation calculation.')
+    parser.add_argument('-m', '--mask', dest='mask_file', required=True,
+                        help='mask file for pixels used for correlation calculation')
+    parser.add_argument('-v', '--velocity', dest='velo_file', required=True,
+                        help='velocity file for generation of deformation mask')
+    
+    parser.add_argument('-w', '--windowsize', type=int, default=141,
+                        help='window size (square window, must be odd number).')
+    parser.add_argument('-r', '--overlapratio', type=float, default=0.4,
+                        help='overlap ratio for window filtering')
+    
+    parser.add_argument('-o', '--outfile', help='output corrected timeseries file name')
+    return parser
+
+def cmd_line_parse(iargs=None):
+    # parse
+    parser = create_parser()
+    inps = parser.parse_args(args=iargs)
+
+    # check: -r / --overlapration option (must be within [0,1])
+    if inps.overlapratio and (not 0.0 <= inps.overlapratio <= 1.0):
+        msg = f'overlap ratio {inps.overlapratio} is NOT within [0.1, 1.0]'
+        raise argparse.ArgumentTypeError(msg)
+    
+    #check: -w / --windowsize option (must be odd number)
+    if inps.windowsize and (inps.windowsize % 2 == 0):
+        msg = f'window size {inps.windowsize} is NOT odd number'
+        raise argparse.ArgumentTypeError(msg)
+    #TODO
+
+    return inps
+
+
+############################################################################
+def main(iargs=None):
+    # parse
+    inps = cmd_line_parse(iargs)
+
+    # import
+    from mintpy.tropo_HTC import run_tropo_HTC
+
+    # run
+    run_tropo_HTC(inps)
+
+
+############################################################################
+if __name__ == '__main__':
+    main(sys.argv[1:])
diff --git a/src/mintpy/smallbaselineApp.py b/src/mintpy/smallbaselineApp.py
index 9e607093a..265688003 100644
--- a/src/mintpy/smallbaselineApp.py
+++ b/src/mintpy/smallbaselineApp.py
@@ -505,6 +505,9 @@ def get_timeseries_filename(template, work_dir='./'):
                     elif method == 'pyaps':
                         fname1 = f'{os.path.splitext(fname0)[0]}_{model}.h5'
 
+                    elif method == 'HTC':
+                        fname1 = f'{os.path.splitext(fname0)[0]}_tropHTC.h5'
+                    
                     else:
                         msg = f'un-recognized tropospheric correction method: {method}'
                         raise ValueError(msg)
@@ -615,8 +618,11 @@ def run_ionospheric_delay_correction(self, step_name):
 
     def run_tropospheric_delay_correction(self, step_name):
         """Correct tropospheric delays."""
-        geom_file = ut.check_loaded_dataset(self.workDir, print_msg=False)[1]
-        mask_file = os.path.join(self.workDir, 'maskTempCoh.h5')
+        #TODO
+        #geom_file = ut.check_loaded_dataset(self.workDir, print_msg=False)[1]
+        geom_file = os.path.join(self.workDir, './inputs/geometry.Geo.h5')
+        #mask_file = os.path.join(self.workDir, 'maskTempCoh.h5')
+        mask_file = os.path.join(self.workDir, 'temporalCoherence.h5')
 
         fnames = self.get_timeseries_filename(self.template, self.workDir)[step_name]
         in_file = fnames['input']
@@ -688,6 +694,18 @@ def get_dataset_size(fname):
 
                         else:
                             raise ValueError(f'un-recognized dataset name: {tropo_model}.')
+            # High-frequency Texture Correlation (Yang et al., 2024)         
+            elif method == 'HTC': 
+                #TODO
+                velocity_dir = self.template['mintpy.tropospheric.velocityDir']
+                window_size = self.template['mintpy.tropospheric.windowSize']
+                overlap_ratio = self.template['mintpy.tropospheric.overlapRatio']
+                iargs = [in_file, '-g', geom_file, '-m', mask_file, '-o', out_file, '-v', velocity_dir, '-w' , window_size, '-r', overlap_ratio]
+                print('tropospheric delay correction with high frequency texture correlation approach')
+                print('\ntropo_HTC.py', ' '.join(iargs))
+                if ut.run_or_skip(out_file=out_file, in_file=in_file) == 'run':
+                    import mintpy.cli.tropo_HTC
+                    mintpy.cli.tropo_HTC.main(iargs)
 
         else:
             print('No tropospheric delay correction.')
diff --git a/src/mintpy/tropo_HTC.py b/src/mintpy/tropo_HTC.py
new file mode 100644
index 000000000..b00dcf338
--- /dev/null
+++ b/src/mintpy/tropo_HTC.py
@@ -0,0 +1,293 @@
+############################################################
+# Program is part of MintPy                                #
+# Copyright (c) 2013, Zhang Yunjun, Heresh Fattahi         #
+# Author: Zhang Yunjun, Heresh Fattahi, 2013               #
+############################################################
+
+
+import os
+import cv2
+import numpy as np
+from scipy.interpolate import interp2d
+
+from mintpy.mask import mask_matrix
+from mintpy.objects import timeseries
+from mintpy.utils import readfile, writefile
+
+############################################################################
+
+def read_topographic_data(geom_file, meta):
+    print('read height & incidenceAngle from file: '+geom_file)
+    dem = readfile.read(geom_file, datasetName='height', print_msg=False)[0]
+    inc_angle = readfile.read(geom_file, datasetName='incidenceAngle', print_msg=False)[0]
+    dem *= 1.0/np.cos(inc_angle*np.pi/180.0)
+
+    ref_y = int(meta['REF_Y'])
+    ref_x = int(meta['REF_X'])
+    dem -= dem[ref_y, ref_x]
+
+    # Design matrix for elevation v.s. phase
+    # dem = dem.flatten()
+    return dem
+
+def read_velocity_data(velo_file, meta):
+    print('read velocity from file: '+velo_file)
+    velocity = readfile.read(velo_file, datasetName='velocity', print_msg=False)[0] #TODO
+
+    return velocity
+
+
+def estimate_local_slope(dem, ts_data, inps, n_ref):
+    """Estimate local slope based on texture correlation for each acquisition of timeseries
+    Parameters: dem     : 2D array in size of (          length, width)
+                ts_data : 3D array in size of (num_date, length, width)
+                inps    : Namespace
+                n_ref   : Index of reference SAR image in timeseries
+    Returns:    X       : 2D array in size of (poly_num+1, num_date)  
+    """#TODO
+
+    """Filtering parameters for obtaining high-frequency texture correlation"""
+    w1 = 9
+    w2 = 13
+    res_step = 0.5
+    Iteration = 10
+    rg = 40
+    step = 0.0001
+
+    velocity = readfile.read(inps.velo_file, datasetName='velocity')[0]
+
+    print('reading mask from file: '+inps.mask_file)
+    mask_coh = readfile.read(inps.mask_file, datasetName='temporalCoherence')[0]
+    maskdef = 0 #TODO
+    mask_def = create_deformation_mask(dem, velocity, maskdef)
+    mask = mask_def * mask_coh
+
+
+    Na, Nr, N = ts_data.shape
+    #Na, Nr = dem.shape
+    W = inps.windowsize
+    w = (W-1)/2
+    overlap = round(inps.overlapratio*(2*w+1))
+    Na_C = np.arange(w + 1, Na - w - overlap + 1, 2 * w - overlap)
+    Nr_C = np.arange(w + 1, Nr - w - overlap + 1, 2 * w - overlap)
+    
+    k_LLF = np.zeros(len(Na_C), len(Nr_C), N)
+    d_LLF = np.zeros(len(Na_C), len(Nr_C), N)
+    k_HTC = np.zeros(len(Na_C), len(Nr_C), N)
+
+    for na in range(0, len(Na_C)):
+        for nr in range(0, len(Nr_C)):
+            # patch
+            nac = Na_C[na]
+            nrc = Nr_C[nr]
+            u = nac - w
+            U = np.where(u < 1, 1, u)
+            d = nac + w
+            D = np.where(d > Na, Na, d)
+            l = nrc - w
+            L = np.where(l < 1, 1, l)
+            r = nrc + w
+            R = np.where(r > Nr, Nr, r)
+
+            if na == len(Na_C):
+                D = Na
+            if nr == len(Nr_C):
+                R = Nr
+            tmp = np.full((Na, Nr), np.nan)
+            tmp[U:D, L:R] = 1
+            mask_process = mask
+            mask_process[np.isnan(tmp)] = np.nan
+
+            # solve
+            result_compare = np.zeros(4, N)
+            for n in range(0, N):
+                if n == n_ref:
+                    continue
+
+                # ----------- initial value ----------- #
+                # mask
+                mask_std = np.ones(Na, Nr) * mask_process
+
+                # iterative linear fitting
+                # res+step = 0.5 # step [rad]
+                # Iteration = 10
+                for iter in range(1, Iteration + 1):
+                    # linear fitting
+                    phase_tmp = ts_data[:, :, n] * mask_process
+                    dem_tmp = dem * mask_std
+                    coe = np.polyfit(dem[~np.isnan(mask_std)], phase_tmp[~np.isnan(mask_std)], 1)
+                    cor_tmp = phase_tmp - coe[0]*dem_tmp - coe[1]
+
+                    # result recording
+                    if iter == 1:
+                        result_compare[0, n] = coe[0]
+                        result_compare[3, n] = coe[1]
+                    
+                    # mask uploading
+                    max_tmp = np.max(np.max(np.abs(cor_tmp))) #TODO
+                    max_tmp = np.where(max_tmp > res_step, max_tmp - res_step, max_tmp)
+                    mask_std[np.abs(cor_tmp) > max_tmp] = np.nan
+
+                    if sum(np.sum(~np.isnan(mask_std))) < sum(np.sum(~np.isnan(mask_process))) / 10: #TODO
+                        break
+
+                # ----------- texture correlation ----------- #        
+                mask_tmp = mask[U:D, L:R]
+                
+                A = dem[U:D, L:R]
+                A_LP = cv2.GaussianBlur(A, (w2, w2), sigmaX=w1, sigmaY=w1)
+                A = A - A_LP
+                A_line = A[~np.isnan(mask_tmp)]
+                A_line = A_line / np.linalg.norm(A_line)
+
+                # range = 40
+                # step = 0.0001
+                # left
+                k_left = -rg * step + coe[0]
+                phase_ts_Scor_tmp = ts_data[:, :, n] - k_left * dem
+                C = phase_ts_Scor_tmp[U:D, L:R]
+                C[np.isnan(C)] = 0
+                C_LP = cv2.GaussianBlur(A, (w2, w2), sigmaX=w1, sigmaY=w1)
+                C = C - C_LP
+                C_line = C[~np.isnan(mask_tmp)]
+                C_line = C_line / np.linalg.norm(C_line)
+                conv_AC = np.sum(A_line * C_line)
+                record_left = np.abs(conv_AC)
+                # right
+                k_right = rg * step + coe[0]
+                phase_ts_Scor_tmp = ts_data[:, :, n] - k_right * dem
+                C = phase_ts_Scor_tmp[U:D, L:R]
+                C[np.isnan(C)] = 0
+                C_LP = cv2.GaussianBlur(A, (w2, w2), sigmaX=w1, sigmaY=w1)
+                C = C - C_LP
+                C_line = C[~np.isnan(mask_tmp)]
+                C_line = C_line / np.linalg.norm(C_line)
+                conv_AC = np.sum(A_line * C_line)
+                record_right = np.abs(conv_AC)
+
+                if record_left >= record_right:
+                    k = np.arange(0, rg + 1) * step + coe[0]
+                    record = np.zeros(len(k))
+                else:
+                    k = np.arange(-rg, 1) * step + coe[0]
+                    record = np.zeros(len(k))
+
+                for i in range(0, len(k)):
+                    phase_ts_Scor_tmp = ts_data[:, :, n] - k[i] * dem
+                    C = phase_ts_Scor_tmp[U:D, L:R]
+                    C[np.isnan(C)] = 0
+                    C_LP = cv2.GaussianBlur(A, (w2, w2), sigmaX=w1, sigmaY=w1)
+                    C = C - C_LP
+                    C_line = C[~np.isnan(mask_tmp)]
+                    C_line = C_line / np.linalg.norm(C_line)
+
+                    conv_AC = np.sum(A_line * C_line)
+                    record[i] = np.abs(conv_AC)
+
+                # slope
+                index = np.where(record == np.min(record))[0]
+                result_compare[1, n] = coe[0]
+                result_compare[2, n] = k[index]
+            
+            # result recording
+            k_LLF[na, nr, :] = result_compare[0, :]
+            d_LLF[na, nr, :] = result_compare[3, :]
+            k_HTC[na, nr, :] = result_compare[2, :]
+    return [k_LLF, d_LLF, k_HTC]
+
+def slope_interpolation(ts_data, inps, k_HTC):
+    """Filtering parameters for obtaining high-frequency texture correlation"""
+    sigma_slope = 7
+    w_slope = 7   
+
+    Na, Nr, N = ts_data.shape
+    #Na, Nr = dem.shape
+    W = inps.windowsize
+    w = (W-1)/2
+    overlap = round(inps.overlapratio*(2*w+1))
+    Na_C = np.arange(w + 1, Na - w - overlap + 1, 2 * w - overlap)
+    Nr_C = np.arange(w + 1, Nr - w - overlap + 1, 2 * w - overlap)
+    # grid
+    Y = Na_C.reshape(-1, 1) * np.ones((1, len(Nr_C)))
+    X = np.ones((len(Na_C), 1)) * Nr_C.reshape(-1, 1)
+
+    Xq, Yq = np.meshgrid(np.arange(1, Nr+1), np.arange(1, Na+1))
+
+    # K_HTC
+    k_HTC_interp = np.zeros(Na, Nr, N)
+    for n in range(0, N):
+        k_HTC_filt = cv2.GaussianBlur(k_HTC[:, :, n], (w_slope, w_slope), sigmaX=sigma_slope, sigmaY=sigma_slope)
+        f_interp = interp2d(X, Y, k_HTC_filt, kind='spline')
+        k_HTC_interp[:, :, n] = f_interp(Xq, Yq).reshape(k_HTC_interp[:, :, n].shape) #TODO 
+    
+    return k_HTC_interp
+
+def intercept_filtering(dem, ts_data, inps, k_HTC_interp, meta):
+    """Filtering parameters for obtaining high-frequency texture correlation"""
+    sigma_intercept = 251 
+    w_intercept = 251 
+
+    ref_y = int(meta['REF_Y'])
+    ref_x = int(meta['REF_X'])
+
+    phase_ts_HTC_low = ts_data
+    intercept = np.zeros(phase_ts_HTC_low.shape)
+    for n in range(0, N):
+        tmp = ts_data[:, :, n] - k_HTC_interp[:, :, n] * dem
+        tmp_filt = cv2.GaussianBlur(tmp, (w_intercept, w_intercept), sigmaX=sigma_intercept, sigmaY=sigma_intercept)
+        tmp = tmp - tmp_filt
+        phase_ts_HTC_low[:, :, n] = tmp
+        intercept[:, :, n] = tmp_filt
+    
+    phase_ts_HTC_low = phase_ts_HTC_low - phase_ts_HTC_low[ref_y, ref_x, :] #TODO
+
+    return phase_ts_HTC_low
+    
+    
+
+def create_deformation_mask(dem, rate, maskdef):
+    if maskdef == 1:
+        mask_def = rate
+        mask_def[np.abs(rate)>1.5] = np.nan
+        mask_def[~np.isnan(mask_def)] = 1
+    else:
+        mask_def = np.ones(dem.shape) 
+    return mask_def
+
+############################################################################
+
+def run_tropo_HTC(inps):
+
+    # read time-series data
+    ts_obj = timeseries(inps.timeseries_file)
+    ts_obj.open()
+    ts_data = ts_obj.read()
+    inps.date_list = list(ts_obj.dateList)
+    n_ref = inps.date_list.index(ts_obj.metadata['REF_DATE'])
+
+    # read topographic data (DEM)
+    dem = read_topographic_data(inps.geom_file, ts_obj.metadata)
+    #TODO
+    # slope estimation
+    k_HTC = estimate_local_slope(dem, ts_data, inps, n_ref)[2]
+    # slope interpolation
+    k_HTC_interp = slope_interpolation(ts_data, inps, k_HTC)
+    # intercept filtering
+    ts_HTC_low = intercept_filtering(dem, ts_data, inps, k_HTC_interp, ts_obj.metadata)
+
+    # write corrected time-series file
+    meta = dict(ts_obj.metadata)
+    meta['mintpy.troposphericDelay.polyOrder'] = str(inps.poly_order)
+    if not inps.outfile:
+        fbase = os.path.splitext(inps.timeseries_file)[0]
+        inps.outfile - f'{fbase}_tropHTC.h5'
+    
+    writefile.write(
+        ts_data,
+        out_file=inps.outfile,
+        metadata=meta,
+        ref_file=inps.timeseries_file
+    )
+
+    return
+

From 1d0510f4f1af2e6b244f4167e8e467f803110bd0 Mon Sep 17 00:00:00 2001
From: Yohcy <865533709@qq.com>
Date: Fri, 25 Oct 2024 16:06:55 +0800
Subject: [PATCH 02/12] tropo_ver2

---
 src/mintpy/cli/{tropo_HTC.py => tropo_htc.py} |   8 +-
 src/mintpy/defaults/smallbaselineApp.cfg      |   5 +
 src/mintpy/defaults/smallbaselineApp_auto.cfg |   5 +
 src/mintpy/smallbaselineApp.py                |  27 ++-
 src/mintpy/{tropo_HTC.py => tropo_htc.py}     | 226 +++++++++++-------
 5 files changed, 162 insertions(+), 109 deletions(-)
 rename src/mintpy/cli/{tropo_HTC.py => tropo_htc.py} (93%)
 rename src/mintpy/{tropo_HTC.py => tropo_htc.py} (52%)

diff --git a/src/mintpy/cli/tropo_HTC.py b/src/mintpy/cli/tropo_htc.py
similarity index 93%
rename from src/mintpy/cli/tropo_HTC.py
rename to src/mintpy/cli/tropo_htc.py
index 61ab2ad51..0b97f907e 100644
--- a/src/mintpy/cli/tropo_HTC.py
+++ b/src/mintpy/cli/tropo_htc.py
@@ -18,8 +18,8 @@
 """
 #TODO
 EXAMPLE = """example:
-  tropo_HTC.py  timeseries_ramp_demErr.h5  -v velocity.h5  -g inputs/geometryRadar.h5  -m maskTempCoh.h5
-  tropo_HTC .py  geo_timeseries_demErr.h5  -g geo_geometryRadar.h5     -m geo_maskTempCoh.h5
+  tropo_htc.py  timeseries_ramp_demErr.h5  -v velocity.h5  -g inputs/geometryRadar.h5  -m maskTempCoh.h5
+  tropo_htc .py  geo_timeseries_demErr.h5  -g geo_geometryRadar.h5     -m geo_maskTempCoh.h5
 """
 
 def create_parser(subparsers=None):
@@ -71,10 +71,10 @@ def main(iargs=None):
     inps = cmd_line_parse(iargs)
 
     # import
-    from mintpy.tropo_HTC import run_tropo_HTC
+    from MintPy.src.mintpy.tropo_htc import run_tropo_htc
 
     # run
-    run_tropo_HTC(inps)
+    run_tropo_htc(inps)
 
 
 ############################################################################
diff --git a/src/mintpy/defaults/smallbaselineApp.cfg b/src/mintpy/defaults/smallbaselineApp.cfg
index 82f56464a..1899742a2 100644
--- a/src/mintpy/defaults/smallbaselineApp.cfg
+++ b/src/mintpy/defaults/smallbaselineApp.cfg
@@ -249,6 +249,11 @@ mintpy.troposphericDelay.minCorrelation = auto  #[0.0-1.0], auto for 0
 ## Set the path below to directory that contains the downloaded *.ztd* files
 mintpy.troposphericDelay.gacosDir = auto # [path2directory], auto for "./GACOS"
 
+## Notes for htc:
+## Set the path below to directory that contains the deformation velocity velocity.h5 files
+mintpy.troposphericDelay.velocityDir = auto # [pah2directory], auto for ".velocity.h5"
+mintpy.troposphericDelay.windowSize = auto
+mintpy.troposphericDelay.overlapRatio = auto
 
 ########## 9. deramp (optional)
 ## Estimate and remove a phase ramp for each acquisition based on the reliable pixels.
diff --git a/src/mintpy/defaults/smallbaselineApp_auto.cfg b/src/mintpy/defaults/smallbaselineApp_auto.cfg
index 8c5bf907f..5313e63c2 100644
--- a/src/mintpy/defaults/smallbaselineApp_auto.cfg
+++ b/src/mintpy/defaults/smallbaselineApp_auto.cfg
@@ -107,6 +107,11 @@ mintpy.troposphericDelay.minCorrelation  = 0
 ## gacos
 mintpy.troposphericDelay.gacosDir        = ./GACOS
 
+## htc
+mintpy.troposphericDelay.velocity        = ./velocity.h5
+mintpy.troposphericDelay.windowSize      = 141
+mintpy.troposphericDelay.overlapRatio    = 0.4
+
 
 ########## deramp
 mintpy.deramp            = no
diff --git a/src/mintpy/smallbaselineApp.py b/src/mintpy/smallbaselineApp.py
index 265688003..509071fdd 100644
--- a/src/mintpy/smallbaselineApp.py
+++ b/src/mintpy/smallbaselineApp.py
@@ -456,8 +456,11 @@ def get_timeseries_filename(template, work_dir='./'):
         Returns:    steps    : dict of dicts, input/output filenames for each step
         """
         work_dir = os.path.abspath(work_dir)
-        fname0 = os.path.join(work_dir, 'timeseries.h5')
-        fname1 = os.path.join(work_dir, 'timeseries.h5')
+        #TODO
+        #fname0 = os.path.join(work_dir, 'timeseries.h5')
+        #fname1 = os.path.join(work_dir, 'timeseries.h5')
+        fname0 = os.path.join(work_dir, 'timeseries_ramp_demErr.h5')
+        fname1 = os.path.join(work_dir, 'timeseries_ramp_demErr.h5')
         atr = readfile.read_attribute(fname0)
 
         phase_correction_steps = [
@@ -505,8 +508,8 @@ def get_timeseries_filename(template, work_dir='./'):
                     elif method == 'pyaps':
                         fname1 = f'{os.path.splitext(fname0)[0]}_{model}.h5'
 
-                    elif method == 'HTC':
-                        fname1 = f'{os.path.splitext(fname0)[0]}_tropHTC.h5'
+                    elif method == 'htc':
+                        fname1 = f'{os.path.splitext(fname0)[0]}_trophtc.h5'
                     
                     else:
                         msg = f'un-recognized tropospheric correction method: {method}'
@@ -620,7 +623,7 @@ def run_tropospheric_delay_correction(self, step_name):
         """Correct tropospheric delays."""
         #TODO
         #geom_file = ut.check_loaded_dataset(self.workDir, print_msg=False)[1]
-        geom_file = os.path.join(self.workDir, './inputs/geometry.Geo.h5')
+        geom_file = os.path.join(self.workDir, 'inputs/geometryGeo.h5')
         #mask_file = os.path.join(self.workDir, 'maskTempCoh.h5')
         mask_file = os.path.join(self.workDir, 'temporalCoherence.h5')
 
@@ -695,17 +698,17 @@ def get_dataset_size(fname):
                         else:
                             raise ValueError(f'un-recognized dataset name: {tropo_model}.')
             # High-frequency Texture Correlation (Yang et al., 2024)         
-            elif method == 'HTC': 
+            elif method == 'htc': 
                 #TODO
-                velocity_dir = self.template['mintpy.tropospheric.velocityDir']
-                window_size = self.template['mintpy.tropospheric.windowSize']
-                overlap_ratio = self.template['mintpy.tropospheric.overlapRatio']
+                velocity_dir = self.template['mintpy.troposphericDelay.velocityDir']
+                window_size = self.template['mintpy.troposphericDelay.windowSize']
+                overlap_ratio = self.template['mintpy.troposphericDelay.overlapRatio']
                 iargs = [in_file, '-g', geom_file, '-m', mask_file, '-o', out_file, '-v', velocity_dir, '-w' , window_size, '-r', overlap_ratio]
                 print('tropospheric delay correction with high frequency texture correlation approach')
-                print('\ntropo_HTC.py', ' '.join(iargs))
+                print('\ntropo_htc.py', ' '.join(iargs))
                 if ut.run_or_skip(out_file=out_file, in_file=in_file) == 'run':
-                    import mintpy.cli.tropo_HTC
-                    mintpy.cli.tropo_HTC.main(iargs)
+                    import mintpy.cli.tropo_htc
+                    mintpy.cli.tropo_htc.main(iargs)
 
         else:
             print('No tropospheric delay correction.')
diff --git a/src/mintpy/tropo_HTC.py b/src/mintpy/tropo_htc.py
similarity index 52%
rename from src/mintpy/tropo_HTC.py
rename to src/mintpy/tropo_htc.py
index b00dcf338..fb96e58ce 100644
--- a/src/mintpy/tropo_HTC.py
+++ b/src/mintpy/tropo_htc.py
@@ -8,7 +8,8 @@
 import os
 import cv2
 import numpy as np
-from scipy.interpolate import interp2d
+from scipy.interpolate import RectBivariateSpline
+from scipy.interpolate import griddata
 
 from mintpy.mask import mask_matrix
 from mintpy.objects import timeseries
@@ -19,15 +20,21 @@
 def read_topographic_data(geom_file, meta):
     print('read height & incidenceAngle from file: '+geom_file)
     dem = readfile.read(geom_file, datasetName='height', print_msg=False)[0]
-    inc_angle = readfile.read(geom_file, datasetName='incidenceAngle', print_msg=False)[0]
-    dem *= 1.0/np.cos(inc_angle*np.pi/180.0)
-
-    ref_y = int(meta['REF_Y'])
-    ref_x = int(meta['REF_X'])
-    dem -= dem[ref_y, ref_x]
-
-    # Design matrix for elevation v.s. phase
-    # dem = dem.flatten()
+    Na, Nr = dem.shape
+    dataLine = dem.flatten()
+    
+    # 找到非NaN值的索引
+    index = np.where(~np.isnan(dataLine))[0]
+    x = np.ceil(index / Na).astype(int)
+    y = np.mod(index, Na).astype(int)
+    
+    # 创建插值网格
+    xi = np.arange(0, Nr)
+    yi = np.arange(0, Na)
+    
+    # 使用三次样条插值
+    dataInp = griddata((x, y), dataLine[index], (xi[None, :], yi[:, None]), method='cubic')
+    dem[np.isnan(dem)] = dataInp[np.isnan(dem)]
     return dem
 
 def read_velocity_data(velo_file, meta):
@@ -37,7 +44,7 @@ def read_velocity_data(velo_file, meta):
     return velocity
 
 
-def estimate_local_slope(dem, ts_data, inps, n_ref):
+def estimate_local_slope(dem, ts_data, inps, n_ref, meta):
     """Estimate local slope based on texture correlation for each acquisition of timeseries
     Parameters: dem     : 2D array in size of (          length, width)
                 ts_data : 3D array in size of (num_date, length, width)
@@ -53,33 +60,45 @@ def estimate_local_slope(dem, ts_data, inps, n_ref):
     Iteration = 10
     rg = 40
     step = 0.0001
+    lamda = 0.3284035 # TODO
+    
+    # TODO
+    #ref_y = int(meta['REF_Y'])
+    #ref_x = int(meta['REF_X'])
+    ref_y = 282
+    ref_x = 204
 
     velocity = readfile.read(inps.velo_file, datasetName='velocity')[0]
 
     print('reading mask from file: '+inps.mask_file)
     mask_coh = readfile.read(inps.mask_file, datasetName='temporalCoherence')[0]
+    mask_coh = np.where(mask_coh < 0.87, np.nan, 1) #TODO
+    
     maskdef = 0 #TODO
     mask_def = create_deformation_mask(dem, velocity, maskdef)
     mask = mask_def * mask_coh
 
-
-    Na, Nr, N = ts_data.shape
+    N, Na, Nr = ts_data.shape
     #Na, Nr = dem.shape
     W = inps.windowsize
     w = (W-1)/2
     overlap = round(inps.overlapratio*(2*w+1))
-    Na_C = np.arange(w + 1, Na - w - overlap + 1, 2 * w - overlap)
-    Nr_C = np.arange(w + 1, Nr - w - overlap + 1, 2 * w - overlap)
+    Na_C = np.arange(w + 1, Na + 1, 2 * w - overlap)
+    Nr_C = np.arange(w + 1, Nr + 1, 2 * w - overlap)
+    
+    k_LLF = np.zeros((N, len(Na_C), len(Nr_C)))
+    d_LLF = np.zeros((N, len(Na_C), len(Nr_C)))
+    k_htc = np.zeros((N, len(Na_C), len(Nr_C)))
+
     
-    k_LLF = np.zeros(len(Na_C), len(Nr_C), N)
-    d_LLF = np.zeros(len(Na_C), len(Nr_C), N)
-    k_HTC = np.zeros(len(Na_C), len(Nr_C), N)
+    ts_data = 4 * np.pi / lamda * ts_data[:N, :, :]
+    reference_value = ts_data[:, ref_y-1, ref_x-1]
+    ts_data = ts_data - reference_value[:, np.newaxis, np.newaxis]
+    ts_data[np.isnan(ts_data)] = 0
 
-    for na in range(0, len(Na_C)):
-        for nr in range(0, len(Nr_C)):
+    for na, nac in enumerate(Na_C):
+        for nr, nrc in enumerate(Nr_C):
             # patch
-            nac = Na_C[na]
-            nrc = Nr_C[nr]
             u = nac - w
             U = np.where(u < 1, 1, u)
             d = nac + w
@@ -89,80 +108,87 @@ def estimate_local_slope(dem, ts_data, inps, n_ref):
             r = nrc + w
             R = np.where(r > Nr, Nr, r)
 
-            if na == len(Na_C):
+            if na == len(Na_C) - 1:
                 D = Na
-            if nr == len(Nr_C):
+            if nr == len(Nr_C) - 1:
                 R = Nr
+            U = int(U)
+            D = int(D)
+            L = int(L)
+            R = int(R)    
             tmp = np.full((Na, Nr), np.nan)
-            tmp[U:D, L:R] = 1
-            mask_process = mask
+            tmp[U-1:D, L-1:R] = 1
+            mask_process = mask.copy()
             mask_process[np.isnan(tmp)] = np.nan
 
             # solve
-            result_compare = np.zeros(4, N)
-            for n in range(0, N):
+            result_compare = np.zeros((N, 4))
+            for n in range(N):
                 if n == n_ref:
                     continue
 
                 # ----------- initial value ----------- #
                 # mask
-                mask_std = np.ones(Na, Nr) * mask_process
+                mask_std = np.ones((Na, Nr)) * mask_process
 
                 # iterative linear fitting
                 # res+step = 0.5 # step [rad]
                 # Iteration = 10
-                for iter in range(1, Iteration + 1):
+                for iter in range(Iteration):
                     # linear fitting
-                    phase_tmp = ts_data[:, :, n] * mask_process
+                    phase_tmp = ts_data[n, :, :] * mask_std
                     dem_tmp = dem * mask_std
-                    coe = np.polyfit(dem[~np.isnan(mask_std)], phase_tmp[~np.isnan(mask_std)], 1)
+                    valid_idx = ~np.isnan(phase_tmp) & ~np.isnan(dem_tmp)
+                    if np.nansum(valid_idx) == 0:
+                        continue
+                    coe = np.polyfit(dem[valid_idx], phase_tmp[valid_idx], 1)
                     cor_tmp = phase_tmp - coe[0]*dem_tmp - coe[1]
 
                     # result recording
-                    if iter == 1:
-                        result_compare[0, n] = coe[0]
-                        result_compare[3, n] = coe[1]
+                    if iter == 0:
+                        result_compare[n, 0] = coe[0]
+                        result_compare[n, 3] = coe[1]
                     
                     # mask uploading
-                    max_tmp = np.max(np.max(np.abs(cor_tmp))) #TODO
+                    max_tmp = np.nanmax(np.abs(cor_tmp)) #TODO
                     max_tmp = np.where(max_tmp > res_step, max_tmp - res_step, max_tmp)
                     mask_std[np.abs(cor_tmp) > max_tmp] = np.nan
 
-                    if sum(np.sum(~np.isnan(mask_std))) < sum(np.sum(~np.isnan(mask_process))) / 10: #TODO
+                    if np.nansum(~np.isnan(mask_std)) < np.nansum(~np.isnan(mask_process)) / 10: #TODO
                         break
 
                 # ----------- texture correlation ----------- #        
-                mask_tmp = mask[U:D, L:R]
+                mask_tmp = mask[U-1:D, L-1:R]
                 
-                A = dem[U:D, L:R]
-                A_LP = cv2.GaussianBlur(A, (w2, w2), sigmaX=w1, sigmaY=w1)
+                A = dem[U-1:D, L-1:R]
+                A_LP = cv2.GaussianBlur(A, (w2, w2), sigmaX=w1, borderType=cv2.BORDER_REPLICATE)
                 A = A - A_LP
                 A_line = A[~np.isnan(mask_tmp)]
                 A_line = A_line / np.linalg.norm(A_line)
-
+                
                 # range = 40
                 # step = 0.0001
                 # left
                 k_left = -rg * step + coe[0]
-                phase_ts_Scor_tmp = ts_data[:, :, n] - k_left * dem
-                C = phase_ts_Scor_tmp[U:D, L:R]
+                phase_ts_Scor_tmp = ts_data[n, :, :] - k_left * dem
+                C = phase_ts_Scor_tmp[U-1:D, L-1:R]
                 C[np.isnan(C)] = 0
-                C_LP = cv2.GaussianBlur(A, (w2, w2), sigmaX=w1, sigmaY=w1)
+                C_LP = cv2.GaussianBlur(C, (w2, w2), sigmaX=w1, borderType=cv2.BORDER_REPLICATE)
                 C = C - C_LP
                 C_line = C[~np.isnan(mask_tmp)]
                 C_line = C_line / np.linalg.norm(C_line)
-                conv_AC = np.sum(A_line * C_line)
+                conv_AC = np.nansum(A_line * C_line)
                 record_left = np.abs(conv_AC)
                 # right
                 k_right = rg * step + coe[0]
-                phase_ts_Scor_tmp = ts_data[:, :, n] - k_right * dem
-                C = phase_ts_Scor_tmp[U:D, L:R]
+                phase_ts_Scor_tmp = ts_data[n, :, :] - k_right * dem
+                C = phase_ts_Scor_tmp[U-1:D, L-1:R]
                 C[np.isnan(C)] = 0
-                C_LP = cv2.GaussianBlur(A, (w2, w2), sigmaX=w1, sigmaY=w1)
+                C_LP = cv2.GaussianBlur(C, (w2, w2), sigmaX=w1, borderType=cv2.BORDER_REPLICATE)
                 C = C - C_LP
                 C_line = C[~np.isnan(mask_tmp)]
                 C_line = C_line / np.linalg.norm(C_line)
-                conv_AC = np.sum(A_line * C_line)
+                conv_AC = np.nansum(A_line * C_line)
                 record_right = np.abs(conv_AC)
 
                 if record_left >= record_right:
@@ -173,75 +199,89 @@ def estimate_local_slope(dem, ts_data, inps, n_ref):
                     record = np.zeros(len(k))
 
                 for i in range(0, len(k)):
-                    phase_ts_Scor_tmp = ts_data[:, :, n] - k[i] * dem
-                    C = phase_ts_Scor_tmp[U:D, L:R]
+                    phase_ts_Scor_tmp = ts_data[n, :, :] - k[i] * dem
+                    C = phase_ts_Scor_tmp[U-1:D, L-1:R]
                     C[np.isnan(C)] = 0
-                    C_LP = cv2.GaussianBlur(A, (w2, w2), sigmaX=w1, sigmaY=w1)
+                    C_LP = cv2.GaussianBlur(C, (w2, w2), sigmaX=w1, borderType=cv2.BORDER_REPLICATE)
                     C = C - C_LP
                     C_line = C[~np.isnan(mask_tmp)]
                     C_line = C_line / np.linalg.norm(C_line)
 
-                    conv_AC = np.sum(A_line * C_line)
+                    conv_AC = np.nansum(A_line * C_line)
                     record[i] = np.abs(conv_AC)
 
                 # slope
-                index = np.where(record == np.min(record))[0]
-                result_compare[1, n] = coe[0]
-                result_compare[2, n] = k[index]
+                #index = np.where(record == np.min(record))[0]
+                index = np.argmin(record)
+                result_compare[n, 1] = coe[0]
+                result_compare[n, 2] = k[index]
             
             # result recording
-            k_LLF[na, nr, :] = result_compare[0, :]
-            d_LLF[na, nr, :] = result_compare[3, :]
-            k_HTC[na, nr, :] = result_compare[2, :]
-    return [k_LLF, d_LLF, k_HTC]
+            k_LLF[:, na, nr] = result_compare[:, 0]
+            d_LLF[:, na, nr] = result_compare[:, 3]
+            k_htc[:, na, nr] = result_compare[:, 2]
+    return [k_LLF, d_LLF, k_htc]
 
-def slope_interpolation(ts_data, inps, k_HTC):
+def slope_interpolation(ts_data, inps, k_htc):
     """Filtering parameters for obtaining high-frequency texture correlation"""
     sigma_slope = 7
     w_slope = 7   
 
-    Na, Nr, N = ts_data.shape
+    N, Na, Nr = ts_data.shape
     #Na, Nr = dem.shape
     W = inps.windowsize
     w = (W-1)/2
     overlap = round(inps.overlapratio*(2*w+1))
-    Na_C = np.arange(w + 1, Na - w - overlap + 1, 2 * w - overlap)
-    Nr_C = np.arange(w + 1, Nr - w - overlap + 1, 2 * w - overlap)
+    Na_C = np.arange(w + 1, Na + 1, 2 * w - overlap)
+    Nr_C = np.arange(w + 1, Na + 1, 2 * w - overlap)
     # grid
     Y = Na_C.reshape(-1, 1) * np.ones((1, len(Nr_C)))
-    X = np.ones((len(Na_C), 1)) * Nr_C.reshape(-1, 1)
+    X = np.ones((len(Na_C), 1)) * Nr_C
 
-    Xq, Yq = np.meshgrid(np.arange(1, Nr+1), np.arange(1, Na+1))
-
-    # K_HTC
-    k_HTC_interp = np.zeros(Na, Nr, N)
+    Xq, Yq = np.meshgrid(np.arange(0, Nr), np.arange(0, Na))
+    
+    # K_htc
+    k_htc_interp = np.zeros((N, Na, Nr))
     for n in range(0, N):
-        k_HTC_filt = cv2.GaussianBlur(k_HTC[:, :, n], (w_slope, w_slope), sigmaX=sigma_slope, sigmaY=sigma_slope)
-        f_interp = interp2d(X, Y, k_HTC_filt, kind='spline')
-        k_HTC_interp[:, :, n] = f_interp(Xq, Yq).reshape(k_HTC_interp[:, :, n].shape) #TODO 
+        k_htc_filt = cv2.GaussianBlur(k_htc[n, :, :], (w_slope, w_slope), sigmaX=sigma_slope, borderType=cv2.BORDER_REPLICATE)
+        spline = RectBivariateSpline(Na_C, Nr_C, k_htc_filt) #TODO
+        k_htc_interp[n, :, :] = spline.ev(Xq, Yq)
+        #f_interp = interp2d(X.flatten(), Y.flatten(), k_htc_filt, kind='spline')
+        #k_htc_interp[n, :, :] = f_interp(Xq, Yq).reshape(Xq.shape) #TODO 
     
-    return k_HTC_interp
+    return k_htc_interp
 
-def intercept_filtering(dem, ts_data, inps, k_HTC_interp, meta):
+def intercept_filtering(dem, ts_data, inps, k_htc_interp, meta):
     """Filtering parameters for obtaining high-frequency texture correlation"""
     sigma_intercept = 251 
     w_intercept = 251 
 
-    ref_y = int(meta['REF_Y'])
-    ref_x = int(meta['REF_X'])
+    # TODO
+    #ref_y = int(meta['REF_Y'])
+    #ref_x = int(meta['REF_X'])
+    ref_y = 282
+    ref_x = 204
+    N, Na, Nr = ts_data.shape
+    lamda = 0.3284035
 
-    phase_ts_HTC_low = ts_data
-    intercept = np.zeros(phase_ts_HTC_low.shape)
+    ts_data = 4 * np.pi / lamda * ts_data[:N, :, :]
+    reference_value = ts_data[:, ref_y-1, ref_x-1]
+    ts_data = ts_data - reference_value[:, np.newaxis, np.newaxis]
+    ts_data[np.isnan(ts_data)] = 0
+
+    phase_ts_htc_low = ts_data
+    
+    intercept = np.zeros(phase_ts_htc_low.shape)
     for n in range(0, N):
-        tmp = ts_data[:, :, n] - k_HTC_interp[:, :, n] * dem
-        tmp_filt = cv2.GaussianBlur(tmp, (w_intercept, w_intercept), sigmaX=sigma_intercept, sigmaY=sigma_intercept)
+        tmp = ts_data[n, :, :] - k_htc_interp[n, :, :] * dem
+        tmp_filt = cv2.GaussianBlur(tmp, (w_intercept, w_intercept), sigmaX=sigma_intercept, borderType=cv2.BORDER_REPLICATE)
         tmp = tmp - tmp_filt
-        phase_ts_HTC_low[:, :, n] = tmp
-        intercept[:, :, n] = tmp_filt
-    
-    phase_ts_HTC_low = phase_ts_HTC_low - phase_ts_HTC_low[ref_y, ref_x, :] #TODO
+        phase_ts_htc_low[n, :, :] = tmp
+        intercept[n, :, :] = tmp_filt
+    reference_phase = phase_ts_htc_low[:, ref_y - 1, ref_x - 1]   
+    phase_ts_htc_low = phase_ts_htc_low - reference_phase[:, np.newaxis, np.newaxis] #TODO
 
-    return phase_ts_HTC_low
+    return phase_ts_htc_low
     
     
 
@@ -256,7 +296,7 @@ def create_deformation_mask(dem, rate, maskdef):
 
 ############################################################################
 
-def run_tropo_HTC(inps):
+def run_tropo_htc(inps):
 
     # read time-series data
     ts_obj = timeseries(inps.timeseries_file)
@@ -269,21 +309,21 @@ def run_tropo_HTC(inps):
     dem = read_topographic_data(inps.geom_file, ts_obj.metadata)
     #TODO
     # slope estimation
-    k_HTC = estimate_local_slope(dem, ts_data, inps, n_ref)[2]
+    k_htc = estimate_local_slope(dem, ts_data, inps, n_ref, ts_obj.metadata)[2]
     # slope interpolation
-    k_HTC_interp = slope_interpolation(ts_data, inps, k_HTC)
+    k_htc_interp = slope_interpolation(ts_data, inps, k_htc)
     # intercept filtering
-    ts_HTC_low = intercept_filtering(dem, ts_data, inps, k_HTC_interp, ts_obj.metadata)
-
+    ts_htc_low = intercept_filtering(dem, ts_data, inps, k_htc_interp, ts_obj.metadata)
+    lamda = 0.3284035 #TODO
+    ts_htc_data = lamda / 4 /np.pi * ts_htc_low
     # write corrected time-series file
     meta = dict(ts_obj.metadata)
-    meta['mintpy.troposphericDelay.polyOrder'] = str(inps.poly_order)
     if not inps.outfile:
         fbase = os.path.splitext(inps.timeseries_file)[0]
-        inps.outfile - f'{fbase}_tropHTC.h5'
+        inps.outfile - f'{fbase}_trophtc.h5'
     
     writefile.write(
-        ts_data,
+        ts_htc_data,
         out_file=inps.outfile,
         metadata=meta,
         ref_file=inps.timeseries_file

From 763c6dbad0c39ba74c63382b278d14e2871ebdf5 Mon Sep 17 00:00:00 2001
From: Yohcy <865533709@qq.com>
Date: Thu, 31 Oct 2024 01:58:26 +0800
Subject: [PATCH 03/12] change the interpolate and gaussian filter

---
 src/mintpy/tropo_htc.py | 29 +++++++++++++++--------------
 1 file changed, 15 insertions(+), 14 deletions(-)

diff --git a/src/mintpy/tropo_htc.py b/src/mintpy/tropo_htc.py
index fb96e58ce..2dd00c37a 100644
--- a/src/mintpy/tropo_htc.py
+++ b/src/mintpy/tropo_htc.py
@@ -6,12 +6,11 @@
 
 
 import os
-import cv2
 import numpy as np
-from scipy.interpolate import RectBivariateSpline
+import scipy
 from scipy.interpolate import griddata
+from scipy.interpolate import interp2d
 
-from mintpy.mask import mask_matrix
 from mintpy.objects import timeseries
 from mintpy.utils import readfile, writefile
 
@@ -61,6 +60,8 @@ def estimate_local_slope(dem, ts_data, inps, n_ref, meta):
     rg = 40
     step = 0.0001
     lamda = 0.3284035 # TODO
+    # w2 = 2 * int(truncate * sigma + 0.5) + 1
+    truncate = ((w2 - 1)/2 - 0.5)/w1
     
     # TODO
     #ref_y = int(meta['REF_Y'])
@@ -161,7 +162,7 @@ def estimate_local_slope(dem, ts_data, inps, n_ref, meta):
                 mask_tmp = mask[U-1:D, L-1:R]
                 
                 A = dem[U-1:D, L-1:R]
-                A_LP = cv2.GaussianBlur(A, (w2, w2), sigmaX=w1, borderType=cv2.BORDER_REPLICATE)
+                A_LP = scipy.ndimage.gaussian_filter(A, sigma=w1, truncate=truncate, mode='nearest')
                 A = A - A_LP
                 A_line = A[~np.isnan(mask_tmp)]
                 A_line = A_line / np.linalg.norm(A_line)
@@ -173,7 +174,7 @@ def estimate_local_slope(dem, ts_data, inps, n_ref, meta):
                 phase_ts_Scor_tmp = ts_data[n, :, :] - k_left * dem
                 C = phase_ts_Scor_tmp[U-1:D, L-1:R]
                 C[np.isnan(C)] = 0
-                C_LP = cv2.GaussianBlur(C, (w2, w2), sigmaX=w1, borderType=cv2.BORDER_REPLICATE)
+                C_LP = scipy.ndimage.gaussian_filter(C, sigma=w1, truncate=truncate, mode='nearest')
                 C = C - C_LP
                 C_line = C[~np.isnan(mask_tmp)]
                 C_line = C_line / np.linalg.norm(C_line)
@@ -184,7 +185,7 @@ def estimate_local_slope(dem, ts_data, inps, n_ref, meta):
                 phase_ts_Scor_tmp = ts_data[n, :, :] - k_right * dem
                 C = phase_ts_Scor_tmp[U-1:D, L-1:R]
                 C[np.isnan(C)] = 0
-                C_LP = cv2.GaussianBlur(C, (w2, w2), sigmaX=w1, borderType=cv2.BORDER_REPLICATE)
+                C_LP = scipy.ndimage.gaussian_filter(C, sigma=w1, truncate=truncate, mode='nearest')
                 C = C - C_LP
                 C_line = C[~np.isnan(mask_tmp)]
                 C_line = C_line / np.linalg.norm(C_line)
@@ -202,7 +203,7 @@ def estimate_local_slope(dem, ts_data, inps, n_ref, meta):
                     phase_ts_Scor_tmp = ts_data[n, :, :] - k[i] * dem
                     C = phase_ts_Scor_tmp[U-1:D, L-1:R]
                     C[np.isnan(C)] = 0
-                    C_LP = cv2.GaussianBlur(C, (w2, w2), sigmaX=w1, borderType=cv2.BORDER_REPLICATE)
+                    C_LP = scipy.ndimage.gaussian_filter(C, sigma=w1, truncate=truncate, mode='nearest')
                     C = C - C_LP
                     C_line = C[~np.isnan(mask_tmp)]
                     C_line = C_line / np.linalg.norm(C_line)
@@ -226,7 +227,7 @@ def slope_interpolation(ts_data, inps, k_htc):
     """Filtering parameters for obtaining high-frequency texture correlation"""
     sigma_slope = 7
     w_slope = 7   
-
+    truncate_slope = ((w_slope - 1)/2 - 0.5)/sigma_slope
     N, Na, Nr = ts_data.shape
     #Na, Nr = dem.shape
     W = inps.windowsize
@@ -243,9 +244,9 @@ def slope_interpolation(ts_data, inps, k_htc):
     # K_htc
     k_htc_interp = np.zeros((N, Na, Nr))
     for n in range(0, N):
-        k_htc_filt = cv2.GaussianBlur(k_htc[n, :, :], (w_slope, w_slope), sigmaX=sigma_slope, borderType=cv2.BORDER_REPLICATE)
-        spline = RectBivariateSpline(Na_C, Nr_C, k_htc_filt) #TODO
-        k_htc_interp[n, :, :] = spline.ev(Xq, Yq)
+        k_htc_filt = scipy.ndimage.gaussian_filter(k_htc[n, :, :], sigma=sigma_slope, truncate=truncate_slope, mode='nearest')
+        f = interp2d(Y, X, k_htc_filt, kind='cubic')
+        k_htc_interp[n, :, :] = f(Yq, Xq)
         #f_interp = interp2d(X.flatten(), Y.flatten(), k_htc_filt, kind='spline')
         #k_htc_interp[n, :, :] = f_interp(Xq, Yq).reshape(Xq.shape) #TODO 
     
@@ -255,7 +256,7 @@ def intercept_filtering(dem, ts_data, inps, k_htc_interp, meta):
     """Filtering parameters for obtaining high-frequency texture correlation"""
     sigma_intercept = 251 
     w_intercept = 251 
-
+    truncate_intercept = ((w_intercept - 1)/2 - 0.5)/sigma_intercept
     # TODO
     #ref_y = int(meta['REF_Y'])
     #ref_x = int(meta['REF_X'])
@@ -274,7 +275,7 @@ def intercept_filtering(dem, ts_data, inps, k_htc_interp, meta):
     intercept = np.zeros(phase_ts_htc_low.shape)
     for n in range(0, N):
         tmp = ts_data[n, :, :] - k_htc_interp[n, :, :] * dem
-        tmp_filt = cv2.GaussianBlur(tmp, (w_intercept, w_intercept), sigmaX=sigma_intercept, borderType=cv2.BORDER_REPLICATE)
+        tmp_filt = scipy.ndimage.gaussian_filter(tmp, sigma=w_intercept, truncate=truncate_intercept, mode='nearest')
         tmp = tmp - tmp_filt
         phase_ts_htc_low[n, :, :] = tmp
         intercept[n, :, :] = tmp_filt
@@ -320,7 +321,7 @@ def run_tropo_htc(inps):
     meta = dict(ts_obj.metadata)
     if not inps.outfile:
         fbase = os.path.splitext(inps.timeseries_file)[0]
-        inps.outfile - f'{fbase}_trophtc.h5'
+        inps.outfile = f'{fbase}_trophtc.h5'
     
     writefile.write(
         ts_htc_data,

From c839f0a1195876963b55deb386d8a9f356915723 Mon Sep 17 00:00:00 2001
From: Yohcy <865533709@qq.com>
Date: Mon, 11 Nov 2024 18:57:20 +0800
Subject: [PATCH 04/12] New interpolation function

---
 .../{tropo_htc.py => tropo_local_texture.py}  |  16 +-
 src/mintpy/defaults/smallbaselineApp.cfg      |  10 +-
 src/mintpy/defaults/smallbaselineApp_auto.cfg |   3 +-
 src/mintpy/smallbaselineApp.py                |  35 ++--
 .../{tropo_htc.py => tropo_local_texture.py}  | 179 ++++++++----------
 5 files changed, 111 insertions(+), 132 deletions(-)
 rename src/mintpy/cli/{tropo_htc.py => tropo_local_texture.py} (82%)
 rename src/mintpy/{tropo_htc.py => tropo_local_texture.py} (72%)

diff --git a/src/mintpy/cli/tropo_htc.py b/src/mintpy/cli/tropo_local_texture.py
similarity index 82%
rename from src/mintpy/cli/tropo_htc.py
rename to src/mintpy/cli/tropo_local_texture.py
index 0b97f907e..9f4a83822 100644
--- a/src/mintpy/cli/tropo_htc.py
+++ b/src/mintpy/cli/tropo_local_texture.py
@@ -2,7 +2,7 @@
 ############################################################
 # Program is part of MintPy                                #
 # Copyright (c) 2013, Zhang Yunjun, Heresh Fattahi         #
-# Author: Antonio Valentino, Zhang Yunjun, Aug 2022        #
+# Author: Yang Qingyue, Hu Changyang, 2024                 #
 ############################################################
 
 import argparse
@@ -18,8 +18,7 @@
 """
 #TODO
 EXAMPLE = """example:
-  tropo_htc.py  timeseries_ramp_demErr.h5  -v velocity.h5  -g inputs/geometryRadar.h5  -m maskTempCoh.h5
-  tropo_htc .py  geo_timeseries_demErr.h5  -g geo_geometryRadar.h5     -m geo_maskTempCoh.h5
+  tropo_local_texture.py  timeseries_ramp_demErr.h5  -v velocity.h5  -g inputs/geometryRadar.h5  -m maskTempCoh.h5
 """
 
 def create_parser(subparsers=None):
@@ -35,8 +34,6 @@ def create_parser(subparsers=None):
                         help='DEM file used for correlation calculation.')
     parser.add_argument('-m', '--mask', dest='mask_file', required=True,
                         help='mask file for pixels used for correlation calculation')
-    parser.add_argument('-v', '--velocity', dest='velo_file', required=True,
-                        help='velocity file for generation of deformation mask')
     
     parser.add_argument('-w', '--windowsize', type=int, default=141,
                         help='window size (square window, must be odd number).')
@@ -53,14 +50,13 @@ def cmd_line_parse(iargs=None):
 
     # check: -r / --overlapration option (must be within [0,1])
     if inps.overlapratio and (not 0.0 <= inps.overlapratio <= 1.0):
-        msg = f'overlap ratio {inps.overlapratio} is NOT within [0.1, 1.0]'
+        msg = f'overlap ratio {inps.overlapratio} is NOT within [0.0, 1.0]'
         raise argparse.ArgumentTypeError(msg)
     
-    #check: -w / --windowsize option (must be odd number)
+    # check: -w / --windowsize option (must be odd number)
     if inps.windowsize and (inps.windowsize % 2 == 0):
         msg = f'window size {inps.windowsize} is NOT odd number'
         raise argparse.ArgumentTypeError(msg)
-    #TODO
 
     return inps
 
@@ -71,10 +67,10 @@ def main(iargs=None):
     inps = cmd_line_parse(iargs)
 
     # import
-    from MintPy.src.mintpy.tropo_htc import run_tropo_htc
+    from mintpy.tropo_local_texture import run_tropo_local_texture
 
     # run
-    run_tropo_htc(inps)
+    run_tropo_local_texture(inps)
 
 
 ############################################################################
diff --git a/src/mintpy/defaults/smallbaselineApp.cfg b/src/mintpy/defaults/smallbaselineApp.cfg
index 1899742a2..4cfbbdd57 100644
--- a/src/mintpy/defaults/smallbaselineApp.cfg
+++ b/src/mintpy/defaults/smallbaselineApp.cfg
@@ -224,7 +224,8 @@ mintpy.ionosphericDelay.excludeDate12 = auto  #[20080520_20090817 / no], auto fo
 ##      NARR  - NARR    from NOAA  [need to install PyAPS from Caltech/EarthDef; recommended for N America]
 ## c. gacos - use GACOS with the iterative tropospheric decomposition model (Yu et al., 2018, JGR)
 ##      need to manually download GACOS products at http://www.gacos.net for all acquisitions before running this step
-mintpy.troposphericDelay.method = auto  #[pyaps / height_correlation / gacos / no], auto for pyaps
+## d. local_texture - use high-frequency texture to estimate the local slope to correct the tropospheric delay (Yang et al., 2024, IEEE-TGRS)
+mintpy.troposphericDelay.method = auto  #[pyaps / height_correlation / gacos / local_texture / no], auto for pyaps
 
 ## Notes for pyaps:
 ## a. GAM data latency: with the most recent SAR data, there will be GAM data missing, the correction
@@ -249,11 +250,10 @@ mintpy.troposphericDelay.minCorrelation = auto  #[0.0-1.0], auto for 0
 ## Set the path below to directory that contains the downloaded *.ztd* files
 mintpy.troposphericDelay.gacosDir = auto # [path2directory], auto for "./GACOS"
 
-## Notes for htc:
+## Notes for local_texture:
 ## Set the path below to directory that contains the deformation velocity velocity.h5 files
-mintpy.troposphericDelay.velocityDir = auto # [pah2directory], auto for ".velocity.h5"
-mintpy.troposphericDelay.windowSize = auto
-mintpy.troposphericDelay.overlapRatio = auto
+mintpy.troposphericDelay.windowSize = auto  #[int], auto for 141, window size, must be odd
+mintpy.troposphericDelay.overlapRatio = auto #[float], auto for 0.4, overlap ratio
 
 ########## 9. deramp (optional)
 ## Estimate and remove a phase ramp for each acquisition based on the reliable pixels.
diff --git a/src/mintpy/defaults/smallbaselineApp_auto.cfg b/src/mintpy/defaults/smallbaselineApp_auto.cfg
index 5313e63c2..dc0ec8c33 100644
--- a/src/mintpy/defaults/smallbaselineApp_auto.cfg
+++ b/src/mintpy/defaults/smallbaselineApp_auto.cfg
@@ -107,8 +107,7 @@ mintpy.troposphericDelay.minCorrelation  = 0
 ## gacos
 mintpy.troposphericDelay.gacosDir        = ./GACOS
 
-## htc
-mintpy.troposphericDelay.velocity        = ./velocity.h5
+## local_texture
 mintpy.troposphericDelay.windowSize      = 141
 mintpy.troposphericDelay.overlapRatio    = 0.4
 
diff --git a/src/mintpy/smallbaselineApp.py b/src/mintpy/smallbaselineApp.py
index 509071fdd..4d2ec5cec 100644
--- a/src/mintpy/smallbaselineApp.py
+++ b/src/mintpy/smallbaselineApp.py
@@ -456,11 +456,10 @@ def get_timeseries_filename(template, work_dir='./'):
         Returns:    steps    : dict of dicts, input/output filenames for each step
         """
         work_dir = os.path.abspath(work_dir)
-        #TODO
-        #fname0 = os.path.join(work_dir, 'timeseries.h5')
-        #fname1 = os.path.join(work_dir, 'timeseries.h5')
-        fname0 = os.path.join(work_dir, 'timeseries_ramp_demErr.h5')
-        fname1 = os.path.join(work_dir, 'timeseries_ramp_demErr.h5')
+        fname0 = os.path.join(work_dir, 'timeseries.h5')
+        fname1 = os.path.join(work_dir, 'timeseries.h5')
+        #fname0 = os.path.join(work_dir, 'timeseries_ramp_demErr.h5')
+        #fname1 = os.path.join(work_dir, 'timeseries_ramp_demErr.h5')
         atr = readfile.read_attribute(fname0)
 
         phase_correction_steps = [
@@ -508,8 +507,8 @@ def get_timeseries_filename(template, work_dir='./'):
                     elif method == 'pyaps':
                         fname1 = f'{os.path.splitext(fname0)[0]}_{model}.h5'
 
-                    elif method == 'htc':
-                        fname1 = f'{os.path.splitext(fname0)[0]}_trophtc.h5'
+                    elif method == 'local_texture':
+                        fname1 = f'{os.path.splitext(fname0)[0]}_tropolocaltexture.h5'
                     
                     else:
                         msg = f'un-recognized tropospheric correction method: {method}'
@@ -622,10 +621,10 @@ def run_ionospheric_delay_correction(self, step_name):
     def run_tropospheric_delay_correction(self, step_name):
         """Correct tropospheric delays."""
         #TODO
-        #geom_file = ut.check_loaded_dataset(self.workDir, print_msg=False)[1]
-        geom_file = os.path.join(self.workDir, 'inputs/geometryGeo.h5')
-        #mask_file = os.path.join(self.workDir, 'maskTempCoh.h5')
-        mask_file = os.path.join(self.workDir, 'temporalCoherence.h5')
+        geom_file = ut.check_loaded_dataset(self.workDir, print_msg=False)[1]
+        #geom_file = os.path.join(self.workDir, 'inputs/geometryGeo.h5')
+        mask_file = os.path.join(self.workDir, 'maskTempCoh.h5')
+        #mask_file = os.path.join(self.workDir, 'temporalCoherence.h5')
 
         fnames = self.get_timeseries_filename(self.template, self.workDir)[step_name]
         in_file = fnames['input']
@@ -697,18 +696,16 @@ def get_dataset_size(fname):
 
                         else:
                             raise ValueError(f'un-recognized dataset name: {tropo_model}.')
-            # High-frequency Texture Correlation (Yang et al., 2024)         
-            elif method == 'htc': 
-                #TODO
-                velocity_dir = self.template['mintpy.troposphericDelay.velocityDir']
+            # High-frequency Texture Correction (Yang et al., 2024)         
+            elif method == 'local_texture': 
                 window_size = self.template['mintpy.troposphericDelay.windowSize']
                 overlap_ratio = self.template['mintpy.troposphericDelay.overlapRatio']
-                iargs = [in_file, '-g', geom_file, '-m', mask_file, '-o', out_file, '-v', velocity_dir, '-w' , window_size, '-r', overlap_ratio]
+                iargs = [in_file, '-g', geom_file, '-m', mask_file, '-o', out_file, '-w' , window_size, '-r', overlap_ratio]
                 print('tropospheric delay correction with high frequency texture correlation approach')
-                print('\ntropo_htc.py', ' '.join(iargs))
+                print('\ntropo_local_texture.py', ' '.join(iargs))
                 if ut.run_or_skip(out_file=out_file, in_file=in_file) == 'run':
-                    import mintpy.cli.tropo_htc
-                    mintpy.cli.tropo_htc.main(iargs)
+                    import mintpy.cli.tropo_local_texture
+                    mintpy.cli.tropo_local_texture.main(iargs)
 
         else:
             print('No tropospheric delay correction.')
diff --git a/src/mintpy/tropo_htc.py b/src/mintpy/tropo_local_texture.py
similarity index 72%
rename from src/mintpy/tropo_htc.py
rename to src/mintpy/tropo_local_texture.py
index 2dd00c37a..80a429a8f 100644
--- a/src/mintpy/tropo_htc.py
+++ b/src/mintpy/tropo_local_texture.py
@@ -1,7 +1,7 @@
 ############################################################
 # Program is part of MintPy                                #
 # Copyright (c) 2013, Zhang Yunjun, Heresh Fattahi         #
-# Author: Zhang Yunjun, Heresh Fattahi, 2013               #
+# Author: Yang Qingyue, Hu Changyang, 2024                 #
 ############################################################
 
 
@@ -9,88 +9,70 @@
 import numpy as np
 import scipy
 from scipy.interpolate import griddata
-from scipy.interpolate import interp2d
+from scipy.interpolate import UnivariateSpline
 
 from mintpy.objects import timeseries
 from mintpy.utils import readfile, writefile
 
 ############################################################################
 
-def read_topographic_data(geom_file, meta):
-    print('read height & incidenceAngle from file: '+geom_file)
+def read_topographic_data(geom_file):
+    print('read height from file: '+geom_file)
     dem = readfile.read(geom_file, datasetName='height', print_msg=False)[0]
     Na, Nr = dem.shape
     dataLine = dem.flatten()
-    
-    # 找到非NaN值的索引
+
     index = np.where(~np.isnan(dataLine))[0]
     x = np.ceil(index / Na).astype(int)
     y = np.mod(index, Na).astype(int)
-    
-    # 创建插值网格
+
     xi = np.arange(0, Nr)
     yi = np.arange(0, Na)
-    
-    # 使用三次样条插值
+
     dataInp = griddata((x, y), dataLine[index], (xi[None, :], yi[:, None]), method='cubic')
     dem[np.isnan(dem)] = dataInp[np.isnan(dem)]
     return dem
 
-def read_velocity_data(velo_file, meta):
-    print('read velocity from file: '+velo_file)
-    velocity = readfile.read(velo_file, datasetName='velocity', print_msg=False)[0] #TODO
-
-    return velocity
-
-
 def estimate_local_slope(dem, ts_data, inps, n_ref, meta):
     """Estimate local slope based on texture correlation for each acquisition of timeseries
     Parameters: dem     : 2D array in size of (          length, width)
                 ts_data : 3D array in size of (num_date, length, width)
                 inps    : Namespace
                 n_ref   : Index of reference SAR image in timeseries
-    Returns:    X       : 2D array in size of (poly_num+1, num_date)  
-    """#TODO
-
-    """Filtering parameters for obtaining high-frequency texture correlation"""
-    w1 = 9
-    w2 = 13
-    res_step = 0.5
-    Iteration = 10
-    rg = 40
-    step = 0.0001
-    lamda = 0.3284035 # TODO
+                meta    : Metadata of timeseries
+    Returns:    k_htc   : 3D array in size of (num_date, length, width), length and width depend on the window size and overlap ratio  
+    """
+
+    """Filtering parameters for obtaining high-frequency texture"""
+    w1 = 9  # slope filtering parameters used in gaussian filter
+    w2 = 13 # texture correlation filtering parameters used in gaussian filter
+    res_step = 0.5 # step size for mask updating
+    Iteration = 10 # iteration times for mask updating
+    rg = 40 # range for slope searching
+    step = 0.0001 # step size for slope searching
     # w2 = 2 * int(truncate * sigma + 0.5) + 1
-    truncate = ((w2 - 1)/2 - 0.5)/w1
+    truncate = ((w2 - 1)/2 - 0.5)/w1 # truncation factor for gaussian filter
     
-    # TODO
-    #ref_y = int(meta['REF_Y'])
-    #ref_x = int(meta['REF_X'])
-    ref_y = 282
-    ref_x = 204
-
-    velocity = readfile.read(inps.velo_file, datasetName='velocity')[0]
-
-    print('reading mask from file: '+inps.mask_file)
-    mask_coh = readfile.read(inps.mask_file, datasetName='temporalCoherence')[0]
-    mask_coh = np.where(mask_coh < 0.87, np.nan, 1) #TODO
+    lamda = int(meta['WAVELENGTH'])
+    ref_y = int(meta['REF_Y'])
+    ref_x = int(meta['REF_X'])
     
-    maskdef = 0 #TODO
-    mask_def = create_deformation_mask(dem, velocity, maskdef)
-    mask = mask_def * mask_coh
-
     N, Na, Nr = ts_data.shape
-    #Na, Nr = dem.shape
     W = inps.windowsize
-    w = (W-1)/2
+    w = (W-1)/2 # half window size
     overlap = round(inps.overlapratio*(2*w+1))
+
+    print('reading mask from file: '+inps.mask_file)
+    mask_coh = readfile.read(inps.mask_file, datasetName='temporalCoherence')[0]
+    mask_coh = np.where(mask_coh < 0.87, np.nan, 1)
+    mask = mask_coh
+
     Na_C = np.arange(w + 1, Na + 1, 2 * w - overlap)
     Nr_C = np.arange(w + 1, Nr + 1, 2 * w - overlap)
     
-    k_LLF = np.zeros((N, len(Na_C), len(Nr_C)))
-    d_LLF = np.zeros((N, len(Na_C), len(Nr_C)))
-    k_htc = np.zeros((N, len(Na_C), len(Nr_C)))
-
+    k_LLF = np.zeros((N, len(Na_C), len(Nr_C)))  # LLF: slope values in spatially discrete distributions
+    d_LLF = np.zeros((N, len(Na_C), len(Nr_C)))  # LLF: intercept values in spatially discrete distributions
+    k_htc = np.zeros((N, len(Na_C), len(Nr_C)))  # HTC: slope values in spatially discrete distributions
     
     ts_data = 4 * np.pi / lamda * ts_data[:N, :, :]
     reference_value = ts_data[:, ref_y-1, ref_x-1]
@@ -123,11 +105,10 @@ def estimate_local_slope(dem, ts_data, inps, n_ref, meta):
             mask_process[np.isnan(tmp)] = np.nan
 
             # solve
-            result_compare = np.zeros((N, 4))
+            result_compare = np.zeros((N, 4))  # line0-3: k_LLF, k_ILLF(iteration), k_htc, d_LLF
             for n in range(N):
                 if n == n_ref:
                     continue
-
                 # ----------- initial value ----------- #
                 # mask
                 mask_std = np.ones((Na, Nr)) * mask_process
@@ -212,7 +193,6 @@ def estimate_local_slope(dem, ts_data, inps, n_ref, meta):
                     record[i] = np.abs(conv_AC)
 
                 # slope
-                #index = np.where(record == np.min(record))[0]
                 index = np.argmin(record)
                 result_compare[n, 1] = coe[0]
                 result_compare[n, 2] = k[index]
@@ -221,57 +201,71 @@ def estimate_local_slope(dem, ts_data, inps, n_ref, meta):
             k_LLF[:, na, nr] = result_compare[:, 0]
             d_LLF[:, na, nr] = result_compare[:, 3]
             k_htc[:, na, nr] = result_compare[:, 2]
-    return [k_LLF, d_LLF, k_htc]
+    return k_htc
 
 def slope_interpolation(ts_data, inps, k_htc):
+    """Estimated local slope interpolation to obtain full-scale slope
+    Parameters: ts_data : 3D array in size of (num_date, length, width)
+                inps    : Namespace
+                k_htc   : 3D array in size of (num_date, length, width), length and width depend on the window size and overlap ratio
+    Returns:    k_htc_interp   : 3D array in size of (num_date, length, width)  
+    """
+
     """Filtering parameters for obtaining high-frequency texture correlation"""
-    sigma_slope = 7
-    w_slope = 7   
-    truncate_slope = ((w_slope - 1)/2 - 0.5)/sigma_slope
+    sigma_slope = 7 # standard deviation for gaussian filter
+    w_slope = 7 # window size for gaussian filter
+    truncate_slope = ((w_slope - 1)/2 - 0.5)/sigma_slope # truncation factor for gaussian filter
     N, Na, Nr = ts_data.shape
-    #Na, Nr = dem.shape
     W = inps.windowsize
     w = (W-1)/2
     overlap = round(inps.overlapratio*(2*w+1))
     Na_C = np.arange(w + 1, Na + 1, 2 * w - overlap)
     Nr_C = np.arange(w + 1, Na + 1, 2 * w - overlap)
-    # grid
-    Y = Na_C.reshape(-1, 1) * np.ones((1, len(Nr_C)))
-    X = np.ones((len(Na_C), 1)) * Nr_C
-
-    Xq, Yq = np.meshgrid(np.arange(0, Nr), np.arange(0, Na))
-    
     # K_htc
     k_htc_interp = np.zeros((N, Na, Nr))
     for n in range(0, N):
+        # filtering
         k_htc_filt = scipy.ndimage.gaussian_filter(k_htc[n, :, :], sigma=sigma_slope, truncate=truncate_slope, mode='nearest')
-        f = interp2d(Y, X, k_htc_filt, kind='cubic')
-        k_htc_interp[n, :, :] = f(Yq, Xq)
-        #f_interp = interp2d(X.flatten(), Y.flatten(), k_htc_filt, kind='spline')
-        #k_htc_interp[n, :, :] = f_interp(Xq, Yq).reshape(Xq.shape) #TODO 
-    
+        # interpolation
+        result = np.zeros((Na, Nr))
+        coords = Na_C - 1
+        # Row interpolation
+        for i in range(k_htc_filt.shape[0]):
+            spline_x = UnivariateSpline(coords, k_htc_filt[i, :], k=3, s=0)
+            result[i, :] = spline_x(np.arange(0, Nr))
+        # Column interpolation
+        for j in range(Nr):
+            spline_y = UnivariateSpline(coords, result[:k_htc_filt.shape[0], j], k=3, s=0)
+            result[:, j] = spline_y(np.arange(0, Na))
+        
+        k_htc_interp[n] = result
     return k_htc_interp
 
 def intercept_filtering(dem, ts_data, inps, k_htc_interp, meta):
+    """Estimate and correct tropospheric delay using intercept filtering
+    Parameters: dem     : 2D array in size of (          length, width)
+                ts_data : 3D array in size of (num_date, length, width)
+                inps    : Namespace
+                k_htc_interp : 3D array in size of (num_date, length, width)
+                meta    : Metadata of timeseries
+    Returns:    phase_ts_htc_low   : 3D array in size of (num_date, length, width)  
+    """
+
     """Filtering parameters for obtaining high-frequency texture correlation"""
-    sigma_intercept = 251 
-    w_intercept = 251 
-    truncate_intercept = ((w_intercept - 1)/2 - 0.5)/sigma_intercept
-    # TODO
-    #ref_y = int(meta['REF_Y'])
-    #ref_x = int(meta['REF_X'])
-    ref_y = 282
-    ref_x = 204
+    sigma_intercept = 251 # standard deviation for gaussian filter
+    w_intercept = 251 # window size for gaussian filter
+    truncate_intercept = ((w_intercept - 1)/2 - 0.5)/sigma_intercept # truncation factor for gaussian filter    
+    ref_y = int(meta['REF_Y'])
+    ref_x = int(meta['REF_X'])
+    lamda = int(meta['WAVELENGTH'])
     N, Na, Nr = ts_data.shape
-    lamda = 0.3284035
 
     ts_data = 4 * np.pi / lamda * ts_data[:N, :, :]
     reference_value = ts_data[:, ref_y-1, ref_x-1]
     ts_data = ts_data - reference_value[:, np.newaxis, np.newaxis]
     ts_data[np.isnan(ts_data)] = 0
 
-    phase_ts_htc_low = ts_data
-    
+    phase_ts_htc_low = ts_data.copy()
     intercept = np.zeros(phase_ts_htc_low.shape)
     for n in range(0, N):
         tmp = ts_data[n, :, :] - k_htc_interp[n, :, :] * dem
@@ -284,20 +278,10 @@ def intercept_filtering(dem, ts_data, inps, k_htc_interp, meta):
 
     return phase_ts_htc_low
     
-    
-
-def create_deformation_mask(dem, rate, maskdef):
-    if maskdef == 1:
-        mask_def = rate
-        mask_def[np.abs(rate)>1.5] = np.nan
-        mask_def[~np.isnan(mask_def)] = 1
-    else:
-        mask_def = np.ones(dem.shape) 
-    return mask_def
 
 ############################################################################
 
-def run_tropo_htc(inps):
+def run_tropo_local_texture(inps):
 
     # read time-series data
     ts_obj = timeseries(inps.timeseries_file)
@@ -307,21 +291,24 @@ def run_tropo_htc(inps):
     n_ref = inps.date_list.index(ts_obj.metadata['REF_DATE'])
 
     # read topographic data (DEM)
-    dem = read_topographic_data(inps.geom_file, ts_obj.metadata)
-    #TODO
+    dem = read_topographic_data(inps.geom_file)
+
     # slope estimation
-    k_htc = estimate_local_slope(dem, ts_data, inps, n_ref, ts_obj.metadata)[2]
+    k_htc = estimate_local_slope(dem, ts_data, inps, n_ref, ts_obj.metadata)
+    
     # slope interpolation
     k_htc_interp = slope_interpolation(ts_data, inps, k_htc)
+    
     # intercept filtering
     ts_htc_low = intercept_filtering(dem, ts_data, inps, k_htc_interp, ts_obj.metadata)
-    lamda = 0.3284035 #TODO
+    lamda = ts_obj.metadata['WAVELENGTH']
     ts_htc_data = lamda / 4 /np.pi * ts_htc_low
+    
     # write corrected time-series file
     meta = dict(ts_obj.metadata)
     if not inps.outfile:
         fbase = os.path.splitext(inps.timeseries_file)[0]
-        inps.outfile = f'{fbase}_trophtc.h5'
+        inps.outfile = f'{fbase}_tropolocaltexture.h5'
     
     writefile.write(
         ts_htc_data,

From 4bc9454930c669ec18b9e0558ce7a75638a64f6a Mon Sep 17 00:00:00 2001
From: Yohcy <865533709@qq.com>
Date: Tue, 12 Nov 2024 01:13:26 +0800
Subject: [PATCH 05/12] test

---
 src/mintpy/cli/tropo_local_texture.py | 10 ++---
 src/mintpy/smallbaselineApp.py        | 10 ++---
 src/mintpy/tropo_local_texture.py     | 53 +++++++++++++--------------
 3 files changed, 36 insertions(+), 37 deletions(-)

diff --git a/src/mintpy/cli/tropo_local_texture.py b/src/mintpy/cli/tropo_local_texture.py
index 9f4a83822..3e89f32a2 100644
--- a/src/mintpy/cli/tropo_local_texture.py
+++ b/src/mintpy/cli/tropo_local_texture.py
@@ -12,7 +12,7 @@
 
 ############################################################################
 REFERENCE = """reference:
-    Yang, Q., Yunjun, Z., Wang, R. Y. Heterogeneous InSAR Tropospheric Correction 
+    Yang, Q., Yunjun, Z., Wang, R. Y. Heterogeneous InSAR Tropospheric Correction
     Based on Local Texture Correlation, IEEE Transactions on Geoscience and Remote Sensing, 62,
       doi:10.1109/TGRS.2024.3356749.
 """
@@ -28,18 +28,18 @@ def create_parser(subparsers=None):
     name = __name__.split('.')[-1]
     parser = create_argument_parser(
         name, synopsis=synopsis, description=synopsis, epilog=epilog, subparsers=subparsers)
-    
+
     parser.add_argument('timeseries_file', help='time-series file to be corrected')
     parser.add_argument('-g', '--geometry', dest='geom_file', required=True,
                         help='DEM file used for correlation calculation.')
     parser.add_argument('-m', '--mask', dest='mask_file', required=True,
                         help='mask file for pixels used for correlation calculation')
-    
+
     parser.add_argument('-w', '--windowsize', type=int, default=141,
                         help='window size (square window, must be odd number).')
     parser.add_argument('-r', '--overlapratio', type=float, default=0.4,
                         help='overlap ratio for window filtering')
-    
+
     parser.add_argument('-o', '--outfile', help='output corrected timeseries file name')
     return parser
 
@@ -52,7 +52,7 @@ def cmd_line_parse(iargs=None):
     if inps.overlapratio and (not 0.0 <= inps.overlapratio <= 1.0):
         msg = f'overlap ratio {inps.overlapratio} is NOT within [0.0, 1.0]'
         raise argparse.ArgumentTypeError(msg)
-    
+
     # check: -w / --windowsize option (must be odd number)
     if inps.windowsize and (inps.windowsize % 2 == 0):
         msg = f'window size {inps.windowsize} is NOT odd number'
diff --git a/src/mintpy/smallbaselineApp.py b/src/mintpy/smallbaselineApp.py
index 4d2ec5cec..9dfa8f242 100644
--- a/src/mintpy/smallbaselineApp.py
+++ b/src/mintpy/smallbaselineApp.py
@@ -509,7 +509,7 @@ def get_timeseries_filename(template, work_dir='./'):
 
                     elif method == 'local_texture':
                         fname1 = f'{os.path.splitext(fname0)[0]}_tropolocaltexture.h5'
-                    
+
                     else:
                         msg = f'un-recognized tropospheric correction method: {method}'
                         raise ValueError(msg)
@@ -620,10 +620,10 @@ def run_ionospheric_delay_correction(self, step_name):
 
     def run_tropospheric_delay_correction(self, step_name):
         """Correct tropospheric delays."""
-        #TODO
         geom_file = ut.check_loaded_dataset(self.workDir, print_msg=False)[1]
-        #geom_file = os.path.join(self.workDir, 'inputs/geometryGeo.h5')
         mask_file = os.path.join(self.workDir, 'maskTempCoh.h5')
+        #For test only
+        #geom_file = os.path.join(self.workDir, 'inputs/geometryGeo.h5')
         #mask_file = os.path.join(self.workDir, 'temporalCoherence.h5')
 
         fnames = self.get_timeseries_filename(self.template, self.workDir)[step_name]
@@ -696,8 +696,8 @@ def get_dataset_size(fname):
 
                         else:
                             raise ValueError(f'un-recognized dataset name: {tropo_model}.')
-            # High-frequency Texture Correction (Yang et al., 2024)         
-            elif method == 'local_texture': 
+            # High-frequency Texture Correction (Yang et al., 2024)
+            elif method == 'local_texture':
                 window_size = self.template['mintpy.troposphericDelay.windowSize']
                 overlap_ratio = self.template['mintpy.troposphericDelay.overlapRatio']
                 iargs = [in_file, '-g', geom_file, '-m', mask_file, '-o', out_file, '-w' , window_size, '-r', overlap_ratio]
diff --git a/src/mintpy/tropo_local_texture.py b/src/mintpy/tropo_local_texture.py
index 80a429a8f..c1209272f 100644
--- a/src/mintpy/tropo_local_texture.py
+++ b/src/mintpy/tropo_local_texture.py
@@ -6,10 +6,10 @@
 
 
 import os
+
 import numpy as np
 import scipy
-from scipy.interpolate import griddata
-from scipy.interpolate import UnivariateSpline
+from scipy.interpolate import UnivariateSpline, griddata
 
 from mintpy.objects import timeseries
 from mintpy.utils import readfile, writefile
@@ -40,7 +40,7 @@ def estimate_local_slope(dem, ts_data, inps, n_ref, meta):
                 inps    : Namespace
                 n_ref   : Index of reference SAR image in timeseries
                 meta    : Metadata of timeseries
-    Returns:    k_htc   : 3D array in size of (num_date, length, width), length and width depend on the window size and overlap ratio  
+    Returns:    k_htc   : 3D array in size of (num_date, length, width), length and width depend on the window size and overlap ratio
     """
 
     """Filtering parameters for obtaining high-frequency texture"""
@@ -52,11 +52,11 @@ def estimate_local_slope(dem, ts_data, inps, n_ref, meta):
     step = 0.0001 # step size for slope searching
     # w2 = 2 * int(truncate * sigma + 0.5) + 1
     truncate = ((w2 - 1)/2 - 0.5)/w1 # truncation factor for gaussian filter
-    
-    lamda = int(meta['WAVELENGTH'])
+
+    lamda = float(meta['WAVELENGTH'])
     ref_y = int(meta['REF_Y'])
     ref_x = int(meta['REF_X'])
-    
+
     N, Na, Nr = ts_data.shape
     W = inps.windowsize
     w = (W-1)/2 # half window size
@@ -69,11 +69,11 @@ def estimate_local_slope(dem, ts_data, inps, n_ref, meta):
 
     Na_C = np.arange(w + 1, Na + 1, 2 * w - overlap)
     Nr_C = np.arange(w + 1, Nr + 1, 2 * w - overlap)
-    
+
     k_LLF = np.zeros((N, len(Na_C), len(Nr_C)))  # LLF: slope values in spatially discrete distributions
     d_LLF = np.zeros((N, len(Na_C), len(Nr_C)))  # LLF: intercept values in spatially discrete distributions
     k_htc = np.zeros((N, len(Na_C), len(Nr_C)))  # HTC: slope values in spatially discrete distributions
-    
+
     ts_data = 4 * np.pi / lamda * ts_data[:N, :, :]
     reference_value = ts_data[:, ref_y-1, ref_x-1]
     ts_data = ts_data - reference_value[:, np.newaxis, np.newaxis]
@@ -98,7 +98,7 @@ def estimate_local_slope(dem, ts_data, inps, n_ref, meta):
             U = int(U)
             D = int(D)
             L = int(L)
-            R = int(R)    
+            R = int(R)
             tmp = np.full((Na, Nr), np.nan)
             tmp[U-1:D, L-1:R] = 1
             mask_process = mask.copy()
@@ -130,7 +130,7 @@ def estimate_local_slope(dem, ts_data, inps, n_ref, meta):
                     if iter == 0:
                         result_compare[n, 0] = coe[0]
                         result_compare[n, 3] = coe[1]
-                    
+
                     # mask uploading
                     max_tmp = np.nanmax(np.abs(cor_tmp)) #TODO
                     max_tmp = np.where(max_tmp > res_step, max_tmp - res_step, max_tmp)
@@ -139,15 +139,15 @@ def estimate_local_slope(dem, ts_data, inps, n_ref, meta):
                     if np.nansum(~np.isnan(mask_std)) < np.nansum(~np.isnan(mask_process)) / 10: #TODO
                         break
 
-                # ----------- texture correlation ----------- #        
+                # ----------- texture correlation ----------- #
                 mask_tmp = mask[U-1:D, L-1:R]
-                
+
                 A = dem[U-1:D, L-1:R]
                 A_LP = scipy.ndimage.gaussian_filter(A, sigma=w1, truncate=truncate, mode='nearest')
                 A = A - A_LP
                 A_line = A[~np.isnan(mask_tmp)]
                 A_line = A_line / np.linalg.norm(A_line)
-                
+
                 # range = 40
                 # step = 0.0001
                 # left
@@ -196,7 +196,7 @@ def estimate_local_slope(dem, ts_data, inps, n_ref, meta):
                 index = np.argmin(record)
                 result_compare[n, 1] = coe[0]
                 result_compare[n, 2] = k[index]
-            
+
             # result recording
             k_LLF[:, na, nr] = result_compare[:, 0]
             d_LLF[:, na, nr] = result_compare[:, 3]
@@ -208,7 +208,7 @@ def slope_interpolation(ts_data, inps, k_htc):
     Parameters: ts_data : 3D array in size of (num_date, length, width)
                 inps    : Namespace
                 k_htc   : 3D array in size of (num_date, length, width), length and width depend on the window size and overlap ratio
-    Returns:    k_htc_interp   : 3D array in size of (num_date, length, width)  
+    Returns:    k_htc_interp   : 3D array in size of (num_date, length, width)
     """
 
     """Filtering parameters for obtaining high-frequency texture correlation"""
@@ -237,7 +237,7 @@ def slope_interpolation(ts_data, inps, k_htc):
         for j in range(Nr):
             spline_y = UnivariateSpline(coords, result[:k_htc_filt.shape[0], j], k=3, s=0)
             result[:, j] = spline_y(np.arange(0, Na))
-        
+
         k_htc_interp[n] = result
     return k_htc_interp
 
@@ -248,16 +248,16 @@ def intercept_filtering(dem, ts_data, inps, k_htc_interp, meta):
                 inps    : Namespace
                 k_htc_interp : 3D array in size of (num_date, length, width)
                 meta    : Metadata of timeseries
-    Returns:    phase_ts_htc_low   : 3D array in size of (num_date, length, width)  
+    Returns:    phase_ts_htc_low   : 3D array in size of (num_date, length, width)
     """
 
     """Filtering parameters for obtaining high-frequency texture correlation"""
     sigma_intercept = 251 # standard deviation for gaussian filter
     w_intercept = 251 # window size for gaussian filter
-    truncate_intercept = ((w_intercept - 1)/2 - 0.5)/sigma_intercept # truncation factor for gaussian filter    
+    truncate_intercept = ((w_intercept - 1)/2 - 0.5)/sigma_intercept # truncation factor for gaussian filter
     ref_y = int(meta['REF_Y'])
     ref_x = int(meta['REF_X'])
-    lamda = int(meta['WAVELENGTH'])
+    lamda = float(meta['WAVELENGTH'])
     N, Na, Nr = ts_data.shape
 
     ts_data = 4 * np.pi / lamda * ts_data[:N, :, :]
@@ -273,11 +273,11 @@ def intercept_filtering(dem, ts_data, inps, k_htc_interp, meta):
         tmp = tmp - tmp_filt
         phase_ts_htc_low[n, :, :] = tmp
         intercept[n, :, :] = tmp_filt
-    reference_phase = phase_ts_htc_low[:, ref_y - 1, ref_x - 1]   
+    reference_phase = phase_ts_htc_low[:, ref_y - 1, ref_x - 1]
     phase_ts_htc_low = phase_ts_htc_low - reference_phase[:, np.newaxis, np.newaxis] #TODO
 
     return phase_ts_htc_low
-    
+
 
 ############################################################################
 
@@ -295,21 +295,21 @@ def run_tropo_local_texture(inps):
 
     # slope estimation
     k_htc = estimate_local_slope(dem, ts_data, inps, n_ref, ts_obj.metadata)
-    
+
     # slope interpolation
     k_htc_interp = slope_interpolation(ts_data, inps, k_htc)
-    
+
     # intercept filtering
     ts_htc_low = intercept_filtering(dem, ts_data, inps, k_htc_interp, ts_obj.metadata)
-    lamda = ts_obj.metadata['WAVELENGTH']
+    lamda = float(ts_obj.metadata['WAVELENGTH'])
     ts_htc_data = lamda / 4 /np.pi * ts_htc_low
-    
+
     # write corrected time-series file
     meta = dict(ts_obj.metadata)
     if not inps.outfile:
         fbase = os.path.splitext(inps.timeseries_file)[0]
         inps.outfile = f'{fbase}_tropolocaltexture.h5'
-    
+
     writefile.write(
         ts_htc_data,
         out_file=inps.outfile,
@@ -318,4 +318,3 @@ def run_tropo_local_texture(inps):
     )
 
     return
-

From 9177ecf7967c85bdd7865987a65ebb5891b7758f Mon Sep 17 00:00:00 2001
From: HUCHANGYANG <149142874+Yohcy@users.noreply.github.com>
Date: Tue, 12 Nov 2024 01:27:01 +0800
Subject: [PATCH 06/12] Update tropo_local_texture.py

---
 src/mintpy/tropo_local_texture.py | 53 +++++++++++++++----------------
 1 file changed, 26 insertions(+), 27 deletions(-)

diff --git a/src/mintpy/tropo_local_texture.py b/src/mintpy/tropo_local_texture.py
index 80a429a8f..c1209272f 100644
--- a/src/mintpy/tropo_local_texture.py
+++ b/src/mintpy/tropo_local_texture.py
@@ -6,10 +6,10 @@
 
 
 import os
+
 import numpy as np
 import scipy
-from scipy.interpolate import griddata
-from scipy.interpolate import UnivariateSpline
+from scipy.interpolate import UnivariateSpline, griddata
 
 from mintpy.objects import timeseries
 from mintpy.utils import readfile, writefile
@@ -40,7 +40,7 @@ def estimate_local_slope(dem, ts_data, inps, n_ref, meta):
                 inps    : Namespace
                 n_ref   : Index of reference SAR image in timeseries
                 meta    : Metadata of timeseries
-    Returns:    k_htc   : 3D array in size of (num_date, length, width), length and width depend on the window size and overlap ratio  
+    Returns:    k_htc   : 3D array in size of (num_date, length, width), length and width depend on the window size and overlap ratio
     """
 
     """Filtering parameters for obtaining high-frequency texture"""
@@ -52,11 +52,11 @@ def estimate_local_slope(dem, ts_data, inps, n_ref, meta):
     step = 0.0001 # step size for slope searching
     # w2 = 2 * int(truncate * sigma + 0.5) + 1
     truncate = ((w2 - 1)/2 - 0.5)/w1 # truncation factor for gaussian filter
-    
-    lamda = int(meta['WAVELENGTH'])
+
+    lamda = float(meta['WAVELENGTH'])
     ref_y = int(meta['REF_Y'])
     ref_x = int(meta['REF_X'])
-    
+
     N, Na, Nr = ts_data.shape
     W = inps.windowsize
     w = (W-1)/2 # half window size
@@ -69,11 +69,11 @@ def estimate_local_slope(dem, ts_data, inps, n_ref, meta):
 
     Na_C = np.arange(w + 1, Na + 1, 2 * w - overlap)
     Nr_C = np.arange(w + 1, Nr + 1, 2 * w - overlap)
-    
+
     k_LLF = np.zeros((N, len(Na_C), len(Nr_C)))  # LLF: slope values in spatially discrete distributions
     d_LLF = np.zeros((N, len(Na_C), len(Nr_C)))  # LLF: intercept values in spatially discrete distributions
     k_htc = np.zeros((N, len(Na_C), len(Nr_C)))  # HTC: slope values in spatially discrete distributions
-    
+
     ts_data = 4 * np.pi / lamda * ts_data[:N, :, :]
     reference_value = ts_data[:, ref_y-1, ref_x-1]
     ts_data = ts_data - reference_value[:, np.newaxis, np.newaxis]
@@ -98,7 +98,7 @@ def estimate_local_slope(dem, ts_data, inps, n_ref, meta):
             U = int(U)
             D = int(D)
             L = int(L)
-            R = int(R)    
+            R = int(R)
             tmp = np.full((Na, Nr), np.nan)
             tmp[U-1:D, L-1:R] = 1
             mask_process = mask.copy()
@@ -130,7 +130,7 @@ def estimate_local_slope(dem, ts_data, inps, n_ref, meta):
                     if iter == 0:
                         result_compare[n, 0] = coe[0]
                         result_compare[n, 3] = coe[1]
-                    
+
                     # mask uploading
                     max_tmp = np.nanmax(np.abs(cor_tmp)) #TODO
                     max_tmp = np.where(max_tmp > res_step, max_tmp - res_step, max_tmp)
@@ -139,15 +139,15 @@ def estimate_local_slope(dem, ts_data, inps, n_ref, meta):
                     if np.nansum(~np.isnan(mask_std)) < np.nansum(~np.isnan(mask_process)) / 10: #TODO
                         break
 
-                # ----------- texture correlation ----------- #        
+                # ----------- texture correlation ----------- #
                 mask_tmp = mask[U-1:D, L-1:R]
-                
+
                 A = dem[U-1:D, L-1:R]
                 A_LP = scipy.ndimage.gaussian_filter(A, sigma=w1, truncate=truncate, mode='nearest')
                 A = A - A_LP
                 A_line = A[~np.isnan(mask_tmp)]
                 A_line = A_line / np.linalg.norm(A_line)
-                
+
                 # range = 40
                 # step = 0.0001
                 # left
@@ -196,7 +196,7 @@ def estimate_local_slope(dem, ts_data, inps, n_ref, meta):
                 index = np.argmin(record)
                 result_compare[n, 1] = coe[0]
                 result_compare[n, 2] = k[index]
-            
+
             # result recording
             k_LLF[:, na, nr] = result_compare[:, 0]
             d_LLF[:, na, nr] = result_compare[:, 3]
@@ -208,7 +208,7 @@ def slope_interpolation(ts_data, inps, k_htc):
     Parameters: ts_data : 3D array in size of (num_date, length, width)
                 inps    : Namespace
                 k_htc   : 3D array in size of (num_date, length, width), length and width depend on the window size and overlap ratio
-    Returns:    k_htc_interp   : 3D array in size of (num_date, length, width)  
+    Returns:    k_htc_interp   : 3D array in size of (num_date, length, width)
     """
 
     """Filtering parameters for obtaining high-frequency texture correlation"""
@@ -237,7 +237,7 @@ def slope_interpolation(ts_data, inps, k_htc):
         for j in range(Nr):
             spline_y = UnivariateSpline(coords, result[:k_htc_filt.shape[0], j], k=3, s=0)
             result[:, j] = spline_y(np.arange(0, Na))
-        
+
         k_htc_interp[n] = result
     return k_htc_interp
 
@@ -248,16 +248,16 @@ def intercept_filtering(dem, ts_data, inps, k_htc_interp, meta):
                 inps    : Namespace
                 k_htc_interp : 3D array in size of (num_date, length, width)
                 meta    : Metadata of timeseries
-    Returns:    phase_ts_htc_low   : 3D array in size of (num_date, length, width)  
+    Returns:    phase_ts_htc_low   : 3D array in size of (num_date, length, width)
     """
 
     """Filtering parameters for obtaining high-frequency texture correlation"""
     sigma_intercept = 251 # standard deviation for gaussian filter
     w_intercept = 251 # window size for gaussian filter
-    truncate_intercept = ((w_intercept - 1)/2 - 0.5)/sigma_intercept # truncation factor for gaussian filter    
+    truncate_intercept = ((w_intercept - 1)/2 - 0.5)/sigma_intercept # truncation factor for gaussian filter
     ref_y = int(meta['REF_Y'])
     ref_x = int(meta['REF_X'])
-    lamda = int(meta['WAVELENGTH'])
+    lamda = float(meta['WAVELENGTH'])
     N, Na, Nr = ts_data.shape
 
     ts_data = 4 * np.pi / lamda * ts_data[:N, :, :]
@@ -273,11 +273,11 @@ def intercept_filtering(dem, ts_data, inps, k_htc_interp, meta):
         tmp = tmp - tmp_filt
         phase_ts_htc_low[n, :, :] = tmp
         intercept[n, :, :] = tmp_filt
-    reference_phase = phase_ts_htc_low[:, ref_y - 1, ref_x - 1]   
+    reference_phase = phase_ts_htc_low[:, ref_y - 1, ref_x - 1]
     phase_ts_htc_low = phase_ts_htc_low - reference_phase[:, np.newaxis, np.newaxis] #TODO
 
     return phase_ts_htc_low
-    
+
 
 ############################################################################
 
@@ -295,21 +295,21 @@ def run_tropo_local_texture(inps):
 
     # slope estimation
     k_htc = estimate_local_slope(dem, ts_data, inps, n_ref, ts_obj.metadata)
-    
+
     # slope interpolation
     k_htc_interp = slope_interpolation(ts_data, inps, k_htc)
-    
+
     # intercept filtering
     ts_htc_low = intercept_filtering(dem, ts_data, inps, k_htc_interp, ts_obj.metadata)
-    lamda = ts_obj.metadata['WAVELENGTH']
+    lamda = float(ts_obj.metadata['WAVELENGTH'])
     ts_htc_data = lamda / 4 /np.pi * ts_htc_low
-    
+
     # write corrected time-series file
     meta = dict(ts_obj.metadata)
     if not inps.outfile:
         fbase = os.path.splitext(inps.timeseries_file)[0]
         inps.outfile = f'{fbase}_tropolocaltexture.h5'
-    
+
     writefile.write(
         ts_htc_data,
         out_file=inps.outfile,
@@ -318,4 +318,3 @@ def run_tropo_local_texture(inps):
     )
 
     return
-

From 82f23911751d0edc9725542f1ecbaed9d9ff7cda Mon Sep 17 00:00:00 2001
From: HUCHANGYANG <149142874+Yohcy@users.noreply.github.com>
Date: Tue, 12 Nov 2024 01:29:40 +0800
Subject: [PATCH 07/12] Update smallbaselineApp.py

---
 src/mintpy/smallbaselineApp.py | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/src/mintpy/smallbaselineApp.py b/src/mintpy/smallbaselineApp.py
index 4d2ec5cec..9dfa8f242 100644
--- a/src/mintpy/smallbaselineApp.py
+++ b/src/mintpy/smallbaselineApp.py
@@ -509,7 +509,7 @@ def get_timeseries_filename(template, work_dir='./'):
 
                     elif method == 'local_texture':
                         fname1 = f'{os.path.splitext(fname0)[0]}_tropolocaltexture.h5'
-                    
+
                     else:
                         msg = f'un-recognized tropospheric correction method: {method}'
                         raise ValueError(msg)
@@ -620,10 +620,10 @@ def run_ionospheric_delay_correction(self, step_name):
 
     def run_tropospheric_delay_correction(self, step_name):
         """Correct tropospheric delays."""
-        #TODO
         geom_file = ut.check_loaded_dataset(self.workDir, print_msg=False)[1]
-        #geom_file = os.path.join(self.workDir, 'inputs/geometryGeo.h5')
         mask_file = os.path.join(self.workDir, 'maskTempCoh.h5')
+        #For test only
+        #geom_file = os.path.join(self.workDir, 'inputs/geometryGeo.h5')
         #mask_file = os.path.join(self.workDir, 'temporalCoherence.h5')
 
         fnames = self.get_timeseries_filename(self.template, self.workDir)[step_name]
@@ -696,8 +696,8 @@ def get_dataset_size(fname):
 
                         else:
                             raise ValueError(f'un-recognized dataset name: {tropo_model}.')
-            # High-frequency Texture Correction (Yang et al., 2024)         
-            elif method == 'local_texture': 
+            # High-frequency Texture Correction (Yang et al., 2024)
+            elif method == 'local_texture':
                 window_size = self.template['mintpy.troposphericDelay.windowSize']
                 overlap_ratio = self.template['mintpy.troposphericDelay.overlapRatio']
                 iargs = [in_file, '-g', geom_file, '-m', mask_file, '-o', out_file, '-w' , window_size, '-r', overlap_ratio]

From d458ed27c5dd2b60202766b65cf68ce75b03cb0e Mon Sep 17 00:00:00 2001
From: HUCHANGYANG <149142874+Yohcy@users.noreply.github.com>
Date: Tue, 12 Nov 2024 01:43:31 +0800
Subject: [PATCH 08/12] Update tropo_local_texture.py

---
 src/mintpy/cli/tropo_local_texture.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/mintpy/cli/tropo_local_texture.py b/src/mintpy/cli/tropo_local_texture.py
index 9f4a83822..c8b67df8c 100644
--- a/src/mintpy/cli/tropo_local_texture.py
+++ b/src/mintpy/cli/tropo_local_texture.py
@@ -18,7 +18,7 @@
 """
 #TODO
 EXAMPLE = """example:
-  tropo_local_texture.py  timeseries_ramp_demErr.h5  -v velocity.h5  -g inputs/geometryRadar.h5  -m maskTempCoh.h5
+  tropo_local_texture.py  timeseries_ramp_demErr.h5 -g inputs/geometryRadar.h5  -m maskTempCoh.h5
 """
 
 def create_parser(subparsers=None):

From 4354ab6ff645241faa0b9779201303cc2a497ef1 Mon Sep 17 00:00:00 2001
From: Yohcy <865533709@qq.com>
Date: Tue, 12 Nov 2024 14:56:30 +0800
Subject: [PATCH 09/12] Update tropo_local_texture

---
 src/mintpy/cli/tropo_local_texture.py |  2 +-
 src/mintpy/tropo_local_texture.py     | 19 ++++++++++---------
 2 files changed, 11 insertions(+), 10 deletions(-)

diff --git a/src/mintpy/cli/tropo_local_texture.py b/src/mintpy/cli/tropo_local_texture.py
index 3e89f32a2..38d39e422 100644
--- a/src/mintpy/cli/tropo_local_texture.py
+++ b/src/mintpy/cli/tropo_local_texture.py
@@ -18,7 +18,7 @@
 """
 #TODO
 EXAMPLE = """example:
-  tropo_local_texture.py  timeseries_ramp_demErr.h5  -v velocity.h5  -g inputs/geometryRadar.h5  -m maskTempCoh.h5
+  tropo_local_texture.py  timeseries_ramp_demErr.h5 -g inputs/geometryRadar.h5  -m maskTempCoh.h5
 """
 
 def create_parser(subparsers=None):
diff --git a/src/mintpy/tropo_local_texture.py b/src/mintpy/tropo_local_texture.py
index c1209272f..4da74b153 100644
--- a/src/mintpy/tropo_local_texture.py
+++ b/src/mintpy/tropo_local_texture.py
@@ -43,7 +43,7 @@ def estimate_local_slope(dem, ts_data, inps, n_ref, meta):
     Returns:    k_htc   : 3D array in size of (num_date, length, width), length and width depend on the window size and overlap ratio
     """
 
-    """Filtering parameters for obtaining high-frequency texture"""
+    # Filtering parameters for obtaining high-frequency texture
     w1 = 9  # slope filtering parameters used in gaussian filter
     w2 = 13 # texture correlation filtering parameters used in gaussian filter
     res_step = 0.5 # step size for mask updating
@@ -116,7 +116,7 @@ def estimate_local_slope(dem, ts_data, inps, n_ref, meta):
                 # iterative linear fitting
                 # res+step = 0.5 # step [rad]
                 # Iteration = 10
-                for iter in range(Iteration):
+                for i in range(Iteration):
                     # linear fitting
                     phase_tmp = ts_data[n, :, :] * mask_std
                     dem_tmp = dem * mask_std
@@ -127,7 +127,7 @@ def estimate_local_slope(dem, ts_data, inps, n_ref, meta):
                     cor_tmp = phase_tmp - coe[0]*dem_tmp - coe[1]
 
                     # result recording
-                    if iter == 0:
+                    if i == 0:
                         result_compare[n, 0] = coe[0]
                         result_compare[n, 3] = coe[1]
 
@@ -180,7 +180,7 @@ def estimate_local_slope(dem, ts_data, inps, n_ref, meta):
                     k = np.arange(-rg, 1) * step + coe[0]
                     record = np.zeros(len(k))
 
-                for i in range(0, len(k)):
+                for i in enumerate(k):
                     phase_ts_Scor_tmp = ts_data[n, :, :] - k[i] * dem
                     C = phase_ts_Scor_tmp[U-1:D, L-1:R]
                     C[np.isnan(C)] = 0
@@ -228,14 +228,15 @@ def slope_interpolation(ts_data, inps, k_htc):
         k_htc_filt = scipy.ndimage.gaussian_filter(k_htc[n, :, :], sigma=sigma_slope, truncate=truncate_slope, mode='nearest')
         # interpolation
         result = np.zeros((Na, Nr))
-        coords = Na_C - 1
+        coords_y = Na_C - 1
+        coords_x = Nr_C - 1
         # Row interpolation
         for i in range(k_htc_filt.shape[0]):
-            spline_x = UnivariateSpline(coords, k_htc_filt[i, :], k=3, s=0)
+            spline_x = UnivariateSpline(coords_x, k_htc_filt[i, :], k=3, s=0)
             result[i, :] = spline_x(np.arange(0, Nr))
         # Column interpolation
         for j in range(Nr):
-            spline_y = UnivariateSpline(coords, result[:k_htc_filt.shape[0], j], k=3, s=0)
+            spline_y = UnivariateSpline(coords_y, result[:k_htc_filt.shape[0], j], k=3, s=0)
             result[:, j] = spline_y(np.arange(0, Na))
 
         k_htc_interp[n] = result
@@ -251,14 +252,14 @@ def intercept_filtering(dem, ts_data, inps, k_htc_interp, meta):
     Returns:    phase_ts_htc_low   : 3D array in size of (num_date, length, width)
     """
 
-    """Filtering parameters for obtaining high-frequency texture correlation"""
+    # Filtering parameters for obtaining high-frequency texture correlation
     sigma_intercept = 251 # standard deviation for gaussian filter
     w_intercept = 251 # window size for gaussian filter
     truncate_intercept = ((w_intercept - 1)/2 - 0.5)/sigma_intercept # truncation factor for gaussian filter
     ref_y = int(meta['REF_Y'])
     ref_x = int(meta['REF_X'])
     lamda = float(meta['WAVELENGTH'])
-    N, Na, Nr = ts_data.shape
+    N = ts_data.shape[0]
 
     ts_data = 4 * np.pi / lamda * ts_data[:N, :, :]
     reference_value = ts_data[:, ref_y-1, ref_x-1]

From 51e61be1cd1f254e1a83529adeb7610648994e8f Mon Sep 17 00:00:00 2001
From: HUCHANGYANG <149142874+Yohcy@users.noreply.github.com>
Date: Tue, 12 Nov 2024 15:02:37 +0800
Subject: [PATCH 10/12] Update tropo_local_texture.py

---
 src/mintpy/tropo_local_texture.py | 19 ++++++++++---------
 1 file changed, 10 insertions(+), 9 deletions(-)

diff --git a/src/mintpy/tropo_local_texture.py b/src/mintpy/tropo_local_texture.py
index c1209272f..4da74b153 100644
--- a/src/mintpy/tropo_local_texture.py
+++ b/src/mintpy/tropo_local_texture.py
@@ -43,7 +43,7 @@ def estimate_local_slope(dem, ts_data, inps, n_ref, meta):
     Returns:    k_htc   : 3D array in size of (num_date, length, width), length and width depend on the window size and overlap ratio
     """
 
-    """Filtering parameters for obtaining high-frequency texture"""
+    # Filtering parameters for obtaining high-frequency texture
     w1 = 9  # slope filtering parameters used in gaussian filter
     w2 = 13 # texture correlation filtering parameters used in gaussian filter
     res_step = 0.5 # step size for mask updating
@@ -116,7 +116,7 @@ def estimate_local_slope(dem, ts_data, inps, n_ref, meta):
                 # iterative linear fitting
                 # res+step = 0.5 # step [rad]
                 # Iteration = 10
-                for iter in range(Iteration):
+                for i in range(Iteration):
                     # linear fitting
                     phase_tmp = ts_data[n, :, :] * mask_std
                     dem_tmp = dem * mask_std
@@ -127,7 +127,7 @@ def estimate_local_slope(dem, ts_data, inps, n_ref, meta):
                     cor_tmp = phase_tmp - coe[0]*dem_tmp - coe[1]
 
                     # result recording
-                    if iter == 0:
+                    if i == 0:
                         result_compare[n, 0] = coe[0]
                         result_compare[n, 3] = coe[1]
 
@@ -180,7 +180,7 @@ def estimate_local_slope(dem, ts_data, inps, n_ref, meta):
                     k = np.arange(-rg, 1) * step + coe[0]
                     record = np.zeros(len(k))
 
-                for i in range(0, len(k)):
+                for i in enumerate(k):
                     phase_ts_Scor_tmp = ts_data[n, :, :] - k[i] * dem
                     C = phase_ts_Scor_tmp[U-1:D, L-1:R]
                     C[np.isnan(C)] = 0
@@ -228,14 +228,15 @@ def slope_interpolation(ts_data, inps, k_htc):
         k_htc_filt = scipy.ndimage.gaussian_filter(k_htc[n, :, :], sigma=sigma_slope, truncate=truncate_slope, mode='nearest')
         # interpolation
         result = np.zeros((Na, Nr))
-        coords = Na_C - 1
+        coords_y = Na_C - 1
+        coords_x = Nr_C - 1
         # Row interpolation
         for i in range(k_htc_filt.shape[0]):
-            spline_x = UnivariateSpline(coords, k_htc_filt[i, :], k=3, s=0)
+            spline_x = UnivariateSpline(coords_x, k_htc_filt[i, :], k=3, s=0)
             result[i, :] = spline_x(np.arange(0, Nr))
         # Column interpolation
         for j in range(Nr):
-            spline_y = UnivariateSpline(coords, result[:k_htc_filt.shape[0], j], k=3, s=0)
+            spline_y = UnivariateSpline(coords_y, result[:k_htc_filt.shape[0], j], k=3, s=0)
             result[:, j] = spline_y(np.arange(0, Na))
 
         k_htc_interp[n] = result
@@ -251,14 +252,14 @@ def intercept_filtering(dem, ts_data, inps, k_htc_interp, meta):
     Returns:    phase_ts_htc_low   : 3D array in size of (num_date, length, width)
     """
 
-    """Filtering parameters for obtaining high-frequency texture correlation"""
+    # Filtering parameters for obtaining high-frequency texture correlation
     sigma_intercept = 251 # standard deviation for gaussian filter
     w_intercept = 251 # window size for gaussian filter
     truncate_intercept = ((w_intercept - 1)/2 - 0.5)/sigma_intercept # truncation factor for gaussian filter
     ref_y = int(meta['REF_Y'])
     ref_x = int(meta['REF_X'])
     lamda = float(meta['WAVELENGTH'])
-    N, Na, Nr = ts_data.shape
+    N = ts_data.shape[0]
 
     ts_data = 4 * np.pi / lamda * ts_data[:N, :, :]
     reference_value = ts_data[:, ref_y-1, ref_x-1]

From 7cfdb4d7d7741b1bab20a6f039f1c7a07085c0dd Mon Sep 17 00:00:00 2001
From: HUCHANGYANG <149142874+Yohcy@users.noreply.github.com>
Date: Tue, 12 Nov 2024 15:07:33 +0800
Subject: [PATCH 11/12] Update tropo_local_texture.py

---
 src/mintpy/tropo_local_texture.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/mintpy/tropo_local_texture.py b/src/mintpy/tropo_local_texture.py
index 4da74b153..83678cf60 100644
--- a/src/mintpy/tropo_local_texture.py
+++ b/src/mintpy/tropo_local_texture.py
@@ -211,7 +211,7 @@ def slope_interpolation(ts_data, inps, k_htc):
     Returns:    k_htc_interp   : 3D array in size of (num_date, length, width)
     """
 
-    """Filtering parameters for obtaining high-frequency texture correlation"""
+    # Filtering parameters for obtaining high-frequency texture correlation
     sigma_slope = 7 # standard deviation for gaussian filter
     w_slope = 7 # window size for gaussian filter
     truncate_slope = ((w_slope - 1)/2 - 0.5)/sigma_slope # truncation factor for gaussian filter

From 8b6932b4356d4a5ff8b32c9817a433cd5edc3a8b Mon Sep 17 00:00:00 2001
From: Yohcy <865533709@qq.com>
Date: Tue, 12 Nov 2024 15:10:44 +0800
Subject: [PATCH 12/12] Update tropo_local_texture

---
 src/mintpy/tropo_local_texture.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/mintpy/tropo_local_texture.py b/src/mintpy/tropo_local_texture.py
index 4da74b153..83678cf60 100644
--- a/src/mintpy/tropo_local_texture.py
+++ b/src/mintpy/tropo_local_texture.py
@@ -211,7 +211,7 @@ def slope_interpolation(ts_data, inps, k_htc):
     Returns:    k_htc_interp   : 3D array in size of (num_date, length, width)
     """
 
-    """Filtering parameters for obtaining high-frequency texture correlation"""
+    # Filtering parameters for obtaining high-frequency texture correlation
     sigma_slope = 7 # standard deviation for gaussian filter
     w_slope = 7 # window size for gaussian filter
     truncate_slope = ((w_slope - 1)/2 - 0.5)/sigma_slope # truncation factor for gaussian filter