diff --git a/src/mintpy/cli/tropo_local_texture.py b/src/mintpy/cli/tropo_local_texture.py
new file mode 100644
index 000000000..38d39e422
--- /dev/null
+++ b/src/mintpy/cli/tropo_local_texture.py
@@ -0,0 +1,78 @@
+#!/usr/bin/env python3
+############################################################
+# Program is part of MintPy                                #
+# Copyright (c) 2013, Zhang Yunjun, Heresh Fattahi         #
+# Author: Yang Qingyue, Hu Changyang, 2024                 #
+############################################################
+
+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_local_texture.py  timeseries_ramp_demErr.h5 -g inputs/geometryRadar.h5  -m 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('-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.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'
+        raise argparse.ArgumentTypeError(msg)
+
+    return inps
+
+
+############################################################################
+def main(iargs=None):
+    # parse
+    inps = cmd_line_parse(iargs)
+
+    # import
+    from mintpy.tropo_local_texture import run_tropo_local_texture
+
+    # run
+    run_tropo_local_texture(inps)
+
+
+############################################################################
+if __name__ == '__main__':
+    main(sys.argv[1:])
diff --git a/src/mintpy/defaults/smallbaselineApp.cfg b/src/mintpy/defaults/smallbaselineApp.cfg
index 82f56464a..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,6 +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 local_texture:
+## Set the path below to directory that contains the deformation velocity velocity.h5 files
+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 8c5bf907f..dc0ec8c33 100644
--- a/src/mintpy/defaults/smallbaselineApp_auto.cfg
+++ b/src/mintpy/defaults/smallbaselineApp_auto.cfg
@@ -107,6 +107,10 @@ mintpy.troposphericDelay.minCorrelation  = 0
 ## gacos
 mintpy.troposphericDelay.gacosDir        = ./GACOS
 
+## local_texture
+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 9e607093a..9dfa8f242 100644
--- a/src/mintpy/smallbaselineApp.py
+++ b/src/mintpy/smallbaselineApp.py
@@ -458,6 +458,8 @@ def get_timeseries_filename(template, work_dir='./'):
         work_dir = os.path.abspath(work_dir)
         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,6 +507,9 @@ def get_timeseries_filename(template, work_dir='./'):
                     elif method == 'pyaps':
                         fname1 = f'{os.path.splitext(fname0)[0]}_{model}.h5'
 
+                    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)
@@ -617,6 +622,9 @@ 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')
+        #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]
         in_file = fnames['input']
@@ -688,6 +696,16 @@ 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':
+                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]
+                print('tropospheric delay correction with high frequency texture correlation approach')
+                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_local_texture
+                    mintpy.cli.tropo_local_texture.main(iargs)
 
         else:
             print('No tropospheric delay correction.')
diff --git a/src/mintpy/tropo_local_texture.py b/src/mintpy/tropo_local_texture.py
new file mode 100644
index 000000000..83678cf60
--- /dev/null
+++ b/src/mintpy/tropo_local_texture.py
@@ -0,0 +1,321 @@
+############################################################
+# Program is part of MintPy                                #
+# Copyright (c) 2013, Zhang Yunjun, Heresh Fattahi         #
+# Author: Yang Qingyue, Hu Changyang, 2024                 #
+############################################################
+
+
+import os
+
+import numpy as np
+import scipy
+from scipy.interpolate import UnivariateSpline, griddata
+
+from mintpy.objects import timeseries
+from mintpy.utils import readfile, writefile
+
+############################################################################
+
+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()
+
+    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 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
+                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 # truncation factor for gaussian filter
+
+    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
+    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)))  # 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]
+    ts_data[np.isnan(ts_data)] = 0
+
+    for na, nac in enumerate(Na_C):
+        for nr, nrc in enumerate(Nr_C):
+            # patch
+            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) - 1:
+                D = Na
+            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-1:D, L-1:R] = 1
+            mask_process = mask.copy()
+            mask_process[np.isnan(tmp)] = np.nan
+
+            # solve
+            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
+
+                # iterative linear fitting
+                # res+step = 0.5 # step [rad]
+                # Iteration = 10
+                for i in range(Iteration):
+                    # linear fitting
+                    phase_tmp = ts_data[n, :, :] * mask_std
+                    dem_tmp = dem * mask_std
+                    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 i == 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)
+                    mask_std[np.abs(cor_tmp) > max_tmp] = np.nan
+
+                    if np.nansum(~np.isnan(mask_std)) < np.nansum(~np.isnan(mask_process)) / 10: #TODO
+                        break
+
+                # ----------- 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
+                k_left = -rg * step + coe[0]
+                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 = 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)
+                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-1:D, L-1:R]
+                C[np.isnan(C)] = 0
+                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)
+                conv_AC = np.nansum(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 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
+                    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)
+
+                    conv_AC = np.nansum(A_line * C_line)
+                    record[i] = np.abs(conv_AC)
+
+                # slope
+                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_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 # 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
+    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)
+    # 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')
+        # interpolation
+        result = np.zeros((Na, Nr))
+        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_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_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
+    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 # 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 = ts_data.shape[0]
+
+    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.copy()
+    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 = 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
+    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
+
+
+############################################################################
+
+def run_tropo_local_texture(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)
+
+    # 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 = 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,
+        metadata=meta,
+        ref_file=inps.timeseries_file
+    )
+
+    return