Iterative outlier rejection for timeseries fits now generalized (#76)

* Iterative outlier rejection for timeseries fits now generalized for GNSS and InSAR data. Applied only to GNSS TS fits for now.

* Clear outputs

---------

Co-authored-by: rzinke <robert.zinke@jpl.nasa.gov>

Author: rzinke

methods/coseismic/Coseismic_Requirement_Validation.ipynb

@@ -113,7 +113,7 @@
     "from mintpy.objects import gnss, timeseries\n",
     "from mintpy.smallbaselineApp import TimeSeriesAnalysis\n",
     "\n",
-    "from solid_utils.gnss_utils import remove_gnss_outliers, model_gnss_timeseries\n",
+    "from solid_utils.fitting import IterativeOutlierFit\n",
     "from solid_utils.sampling import load_geo, samp_pair, profile_samples, haversine_distance\n",
     "from solid_utils.plotting import display_coseismic_validation, display_validation_table"
    ]
@@ -904,35 +904,26 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "# Get InSAR geometry file\n",
-    "geom_file = ut.get_geometry_file(['incidenceAngle','azimuthAngle'],\n",
-    "                                 work_dir=mintpy_dir, coord='geo')\n",
+    "# Outlier ID and removal parameters\n",
+    "outlier_thresh = 3  # standard deviations\n",
+    "outlier_iters = 1  # iterations\n",
     "\n",
     "# Loop through valid GNSS stations\n",
     "for i, site_name in enumerate(site_names):\n",
     "    # Retrieve station information\n",
     "    gnss_stn = gnss_stns[site_name]\n",
-    "    gnss_stn.get_los_displacement(geom_file,\n",
+    "    gnss_stn.get_los_displacement(insar_metadata,\n",
     "                                  start_date=start_date,\n",
     "                                  end_date=end_date)\n",
     "\n",
-    "    # Prepare dictionary of displacement measurements\n",
-    "    gnss_stn.dis_dict = {}  # empty for now\n",
-    "    \n",
-    "    # Save unfiltered data for plotting\n",
-    "    unfilt_dates = gnss_stn.dates\n",
-    "    unfilt_dis = gnss_stn.dis_los\n",
-    "    \n",
     "    # Outlier detection and removal\n",
-    "    remove_gnss_outliers(gnss_stn, 'LOS', model=ts_functions,\n",
-    "                         threshold=3, max_iter=2, verbose=False)\n",
+    "    gnss_fit = IterativeOutlierFit(gnss_stn.dates, gnss_stn.dis_los,\n",
+    "                                   model=ts_functions, threshold=3,\n",
+    "                                   max_iter=1)\n",
     "    \n",
-    "    # Model outlier-removed time-series (estimate method 2)\n",
-    "    dis_los_hat, m_hat, mhat_se = model_gnss_timeseries(gnss_stn, 'LOS', ts_functions)\n",
-    "    \n",
-    "    # Record station velocity\n",
-    "    gnss_stn.dis_dict['dis_mm'] = m_hat[2] * 1000\n",
-    "    gnss_stn.dis_dict['dis_error_mm'] = mhat_se[2] * 1000\n",
+    "    # Record station displacement\n",
+    "    gnss_stn.dis_dict = {'dis_mm': gnss_fit.m_hat[2] * 1000,\n",
+    "                         'dis_error_mm': gnss_fit.mhat_se[2] * 1000}\n",
     "\n",
     "    # Report\n",
     "    if i == 0 :\n",
@@ -950,20 +941,18 @@
     "    if ax_nb == 0:\n",
     "        fig, axes = plt.subplots(figsize=(10, 4), ncols=2)\n",
     "\n",
+    "    # Plot outliers\n",
+    "    if gnss_fit.n_outliers > 0:\n",
+    "        axes[ax_nb].scatter(gnss_fit.outlier_dates,\n",
+    "                            1000*gnss_fit.outlier_dis,\n",
+    "                            2**2, 'firebrick', label='outliers')\n",
+    "\n",
     "    # Plot filtered data and model fit\n",
-    "    axes[ax_nb].scatter(gnss_stn.dates, 1000*gnss_stn.dis_los,\n",
-    "                        3**2, 'k', label='filt data')\n",
-    "    axes[ax_nb].plot(gnss_stn.dates, 1000*dis_los_hat,\n",
-    "                     'c', linewidth=3, label='model fit')\n",
+    "    axes[ax_nb].scatter(gnss_fit.dates, 1000*gnss_fit.dis, 3**2,\n",
+    "                        'dimgrey', zorder=1)\n",
+    "    axes[ax_nb].plot(gnss_fit.dates, 1000*gnss_fit.dis_hat,\n",
+    "                     'c', linewidth=3, label='model fit', zorder=2)\n",
     "\n",
-    "    # Plot outliers\n",
-    "    outlier_ndxs = [i for i, date in enumerate(unfilt_dates) if date not in gnss_stn.dates]\n",
-    "    if len(outlier_ndxs) > 0:\n",
-    "        outlier_dates = np.array([unfilt_dates[i] for i in outlier_ndxs])\n",
-    "        outlier_dis = np.array([unfilt_dis[i] for i in outlier_ndxs])\n",
-    "        axes[ax_nb].scatter(outlier_dates, 1000*outlier_dis,\n",
-    "                            3**2, 'firebrick', label='outliers')\n",
-    "    \n",
     "    # Format plot\n",
     "    axes[ax_nb].legend()\n",
     "    axes[ax_nb].set_xticks(gnss_stn.dates[::label_skips])\n",

methods/secular/Secular_Requirement_Validation.ipynb

@@ -115,9 +115,10 @@
     "from mintpy.utils import ptime, readfile, time_func, utils as ut\n",
     "from pathlib import Path\n",
     "from scipy import signal\n",
+    "\n",
+    "from solid_utils.fitting import IterativeOutlierFit\n",
     "from solid_utils.plotting import display_validation, display_validation_table\n",
     "from solid_utils.sampling import load_geo, samp_pair, profile_samples, haversine_distance\n",
-    "from solid_utils.gnss_utils import remove_gnss_outliers, model_gnss_timeseries\n",
     "\n",
     "#Set Global Plot Parameters\n",
     "plt.rcParams.update({'font.size': 12})"
@@ -686,7 +687,7 @@
    "outputs": [],
    "source": [
     "# GNSS processing source\n",
-    "if 'gnss_source' in sitedata['sites']:\n",
+    "if 'gnss_source' in sitedata['sites'][site]:\n",
     "    gnss_source = sitedata['sites'][site]['gnss_source']\n",
     "else:\n",
     "    gnss_source = 'UNR'\n",
@@ -861,6 +862,10 @@
    },
    "outputs": [],
    "source": [
+    "# Outlier ID and removal parameters\n",
+    "outlier_thresh = 3  # standard deviations\n",
+    "outlier_iters = 1  # iterations\n",
+    "\n",
     "# Loop through valid GNSS stations\n",
     "for i, site_name in enumerate(site_names):\n",
     "    # Retrieve station information\n",
@@ -869,27 +874,21 @@
     "                                  start_date=start_date,\n",
     "                                  end_date=end_date)\n",
     "\n",
-    "    # Save unfiltered data for plotting\n",
-    "    unfilt_dates = gnss_stn.dates\n",
-    "    unfilt_dis = gnss_stn.dis_los\n",
-    "    \n",
     "    # Outlier detection and removal\n",
-    "    remove_gnss_outliers(gnss_stn, 'LOS', model=ts_functions,\n",
-    "                         threshold=3, max_iter=2, verbose=False)\n",
-    "    \n",
-    "    # Model outlier-removed time-series\n",
-    "    dis_los_hat, m_hat, mhat_se = model_gnss_timeseries(gnss_stn, 'LOS', ts_functions)\n",
-    "    \n",
+    "    gnss_fit = IterativeOutlierFit(gnss_stn.dates, gnss_stn.dis_los,\n",
+    "                                   model=ts_functions, threshold=outlier_thresh,\n",
+    "                                   max_iter=outlier_iters)\n",
+    "\n",
     "    # Record station velocity\n",
-    "    gnss_stn.vel_dict = {'vel_mm_yr': m_hat[1] * 1000,\n",
-    "                         'vel_error_mm_yr': mhat_se[1] * 1000}\n",
+    "    gnss_stn.vel_dict = {'vel_mm_yr': gnss_fit.m_hat[1] * 1000,\n",
+    "                         'vel_error_mm_yr': gnss_fit.mhat_se[1] * 1000}\n",
     "\n",
     "    # Report\n",
     "    if i == 0 :\n",
-    "        print('site velocity(mm/yr)')\n",
-    "    print('{:s} {:.1f} +- {:.2f}'.format(gnss_stn.site,\n",
-    "                                         gnss_stn.vel_dict['vel_mm_yr'],\n",
-    "                                         gnss_stn.vel_dict['vel_error_mm_yr']))\n",
+    "        print('site velocity (mm/yr)')\n",
+    "    print(f\"{gnss_stn.site:s} \"\n",
+    "          f\"{gnss_stn.vel_dict['vel_mm_yr']:.1f} \"\n",
+    "          f\"+- {gnss_stn.vel_dict['vel_error_mm_yr']:.2f}\")\n",
     "\n",
     "    # Plotting parameters\n",
     "    n_dates = len(gnss_stn.dates)\n",
@@ -900,19 +899,17 @@
     "    if ax_nb == 0:\n",
     "        fig, axes = plt.subplots(figsize=(10, 4), ncols=2)\n",
     "\n",
-    "    # Plot filtered data and model fit\n",
-    "    axes[ax_nb].scatter(gnss_stn.dates, 1000*gnss_stn.dis_los,\n",
-    "                        3**2, 'k', label='filt data')\n",
-    "    axes[ax_nb].plot(gnss_stn.dates, 1000*dis_los_hat,\n",
-    "                     'c', linewidth=3, label='model fit')\n",
-    "\n",
     "    # Plot outliers\n",
-    "    outlier_ndxs = [i for i, date in enumerate(unfilt_dates) if date not in gnss_stn.dates]\n",
-    "    if len(outlier_ndxs) > 0:\n",
-    "        outlier_dates = np.array([unfilt_dates[i] for i in outlier_ndxs])\n",
-    "        outlier_dis = np.array([unfilt_dis[i] for i in outlier_ndxs])\n",
-    "        axes[ax_nb].scatter(outlier_dates, 1000*outlier_dis,\n",
-    "                            3**2, 'firebrick', label='outliers')\n",
+    "    if gnss_fit.n_outliers > 0:\n",
+    "        axes[ax_nb].scatter(gnss_fit.outlier_dates,\n",
+    "                            1000*gnss_fit.outlier_dis,\n",
+    "                            2**2, 'firebrick', label='outliers')\n",
+    "\n",
+    "    # Plot filtered data and model fit\n",
+    "    axes[ax_nb].scatter(gnss_fit.dates, 1000*gnss_fit.dis, 3**2,\n",
+    "                        'dimgrey', zorder=1)\n",
+    "    axes[ax_nb].plot(gnss_fit.dates, 1000*gnss_fit.dis_hat,\n",
+    "                     'c', linewidth=3, label='model fit', zorder=2)\n",
     "    \n",
     "    # Format plot\n",
     "    axes[ax_nb].legend()\n",
@@ -922,8 +919,7 @@
     "                                rotation=80)\n",
     "    axes[ax_nb].set_ylabel('LOS displacement (mm)')\n",
     "    axes[ax_nb].set_title(f'{gnss_stn.site:s} ({gnss_stn.vel_dict['vel_mm_yr']:.1f} mm/yr)')\n",
-    "    fig.tight_layout()\n",
-    "    "
+    "    fig.tight_layout()\n"
    ]
   },
   {

solid_utils/fitting.py

@@ -0,0 +1,149 @@
+# Author: R Zinke
+# March, 2025
+
+import numpy as np
+from scipy import stats
+
+from mintpy.utils import time_func, utils as ut
+
+
+class TSModelFit:
+    def __init__(self, dates:np.ndarray, dis:np.ndarray, model:dict, conf=95.45):
+        """Model a displacement time-series.
+    
+        Parameters: dates     - np.ndarray, dates as Python datetime objects
+                    dis       - np.ndarray, displacements
+                    model     - dict, time function model
+        Attributes: dis_hat   - np.ndarray, array of predicted displacement
+                                values
+                    m_hat     - np.ndarray, model fit parameters
+                    mhat_se   - np.ndarray, standard errors for model fit params
+                    err_envelope - list of np.ndarray, upper and lower confidence
+                                bounds
+        """
+        # Record parameters
+        self.model = model
+        if len(dates) != len(dis):
+            raise ValueError('Must have same number of dates and displacements')
+        self.dates = dates
+        self.dis = dis
+        self.n = len(dates)
+
+        # Construct design matrix from dates and model
+        date_list = [date.strftime('%Y%m%d') for date in dates]
+        G = time_func.get_design_matrix4time_func(date_list, self.model)
+    
+        # Invert for model parameters
+        self.m_hat = np.linalg.pinv(G).dot(self.dis)
+    
+        # Predict displacements
+        self.dis_hat = np.dot(G, self.m_hat)
+    
+        # Quantify error on model parameters
+        resids = self.dis - self.dis_hat
+        sse = np.sum(resids**2)
+        dof = len(self.m_hat)
+        pcov = sse/(self.n - dof) * np.linalg.inv(np.dot(G.T, G))
+        self.mhat_se = np.sqrt(np.diag(pcov))
+    
+        # Propagate uncertainty
+        dcov = G.dot(pcov).dot(G.T)
+        derr = np.sqrt(np.diag(dcov))
+    
+        # Error envelope
+        err_scale = stats.t.interval(conf/100, dof)
+        err_lower = self.dis_hat + err_scale[0] * derr
+        err_upper = self.dis_hat + err_scale[1] * derr
+        self.err_envelope = [err_lower, err_upper]
+
+
+class IterativeOutlierFit:
+    @staticmethod
+    def outliers_zscore(dis:np.ndarray, dis_hat:np.ndarray, threshold:float):
+        """Identify outliers using the z-score metric.
+    
+        Compute the number of standard deviations the data are from the mean
+        and return the indices of values greater than the specified threshold.
+    
+        Parameters: dis          - np.ndarray, array of displacement values
+                    dis_hat      - np.ndarray, array of predicted displacement
+                                   values
+                    threshold    - float, z-score value (standard deviation)
+                                   beyond which to exclude data
+        Returns:    outlier_ndxs - np.ndarray, boolean array where True
+                                   indicates an outlier
+                    n_outliers   - int, number of outliers
+        """
+        zscores = (dis - dis_hat) / np.std(dis - dis_hat)
+        outlier_ndxs = np.abs(zscores) > threshold
+        n_outliers = np.sum(outlier_ndxs)
+    
+        return outlier_ndxs, n_outliers
+    
+    def __init__(self, dates, dis, model, threshold=3, max_iter=2):
+        """Determine which data points are outliers based on the z-score
+        metric and remove those points.
+
+        Parameters: dates     - np.ndarray, dates as Python datetime objects
+                    dis       - np.ndarray, displacements
+                    model     - dict, time function model
+                    threshold - float, standard deviations beyond which values
+                                are considered outliers
+                    max_iter  - int, maximutm number of iterations before
+                                stopping
+        Attributes: dis_hat   - np.ndarray, array of predicted displacement
+                                values
+                    mhat      - np.ndarray, model fit parameters
+                    mhat_se   - np.ndarray, standard errors for model fit
+                                params
+                    err_envelope - list of np.ndarray, upper and lower confidence
+                                bounds
+        """
+        # Record parameters
+        self.dates = dates
+        self.dis = dis
+        self.model = model
+        self.outlier_threshold = threshold
+
+        # Initialize outlier removal
+        self.iters = 0
+        self.outlier_dates = np.array([])
+        self.outlier_dis = np.array([])
+
+        # Initial fit to data
+        ts_model = TSModelFit(self.dates, self.dis, self.model)
+
+        # Determine outliers based on z-score
+        (outlier_ndxs,
+         n_outliers) = self.outliers_zscore(ts_model.dis, ts_model.dis_hat,
+                                            self.outlier_threshold)
+        self.n_outliers = n_outliers
+
+        # Remove outliers from data set
+        while (n_outliers > 0) and (self.iters < max_iter):
+            # Update time series
+            self.outlier_dates = np.append(self.outlier_dates,
+                                           self.dates[outlier_ndxs])
+            self.outlier_dis = np.append(self.outlier_dis,
+                                         self.dis[outlier_ndxs])
+
+            self.dates = self.dates[~outlier_ndxs]
+            self.dis = self.dis[~outlier_ndxs]
+
+            # Update timeseries model
+            ts_model = TSModelFit(self.dates, self.dis, self.model)
+
+            # Determine outliers based on z-score
+            (outlier_ndxs,
+             n_outliers) = self.outliers_zscore(ts_model.dis, ts_model.dis_hat,
+                                                self.outlier_threshold)
+            self.n_outliers += n_outliers
+
+            # Update iteration counter
+            self.iters += 1
+
+        # Record final parameters
+        self.m_hat = ts_model.m_hat
+        self.mhat_se = ts_model.mhat_se
+        self.dis_hat = ts_model.dis_hat
+        self.err_envelope = ts_model.err_envelope