methods/secular/Secular_Requirement_Validation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Workflow to Validate NISAR L2 Secular Displacement Requirement\n",
    "\n",
    "**Original code authored by:** David Bekaert, Heresh Fattahi, Eric Fielding, and Zhang Yunjun \n",
    "\n",
    "Extensive modifications by Adrian Borsa and Amy Whetter 2022 <br>\n",
    "Subsequent updates and reformatting by Rob Zinke, Katia Tymofyeveva and Adrian Borsa 2024\n",
    "\n",
    "<div class=\"alert alert-warning\">\n",
    "Step 1 of the validation workflow is executed in the <b>ARIA_prep</b> notebook. All processing steps in this notebook MUST be run sequentially, although several are optional and can be skipped.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<hr/>\n",
    "\n",
    "## Table of Contents: <a id='TOC'></a>\n",
    "\n",
    "[**Environment Setup**](#setup)\n",
    "- [Load Python Packages](#load_packages)\n",
    "- [Define CalVal Site and Parameters](#set_calval_params)\n",
    "- [Set Directories and Files](#set_directories)\n",
    "\n",
    "[**1. Download and Prepare Interferograms**](#prep_ifg)\n",
    "[Executed in ARIA_prep]\n",
    "\n",
    "[**2. Generate Time Series from Interferograms**](#gen_ts)\n",
    "- [2.1. Validate/Modify Interferogram Network](#validate_network)\n",
    "- [2.2. Generate Quality Control Mask](#generate_mask)\n",
    "- [2.3. Reference Interferograms To Common Lat/Lon](#common_latlon)\n",
    "- [2.4. Invert for SBAS Line-of-Sight Timeseries](#invert_SBAS)\n",
    "\n",
    "[**3. Optional Corrections**](#opt_correction)\n",
    "- [3.1. Solid Earth Tide Correction](#solid_earth)\n",
    "- [3.2. Tropospheric Delay Correction](#tropo_corr)\n",
    "- [3.3. Phase Deramping ](#phase_deramp)\n",
    "- [3.4. Topographic Residual Correction ](#topo_corr) \n",
    "\n",
    "[**4. Decompose InSAR and GNSS Time Series**](#decomp_ts)\n",
    "- [4.1. Estimate InSAR LOS Velocities](#insar_vel)\n",
    "- [4.2. Not Used](#empty)\n",
    "- [4.3. Find Collocated GNSS Stations](#find_gps)  \n",
    "- [4.4. Get GNSS Position Time Series](#gps_ts) \n",
    "- [4.5. Estimate GNSS LOS Velocities](#gps_vel)\n",
    "- [4.6. Re-reference GNSS and InSAR](#reference)\n",
    "\n",
    "[**5. Validation Method 1: GNSS-InSAR Direct Comparison**](#validation1)\n",
    "- [5.1. Make GNSS-InSAR Velocity Residuals at GNSS Station Locations](#make_resid)\n",
    "- [5.2. Make Double-differenced Velocity Residuals](#make_ddiff)\n",
    "- [5.3. Secular Requirement Validation: Method 1](#secular_validation_method1)\n",
    "\n",
    "[**6. Validation Method 2: InSAR-only Structure Function**](#validation2)\n",
    "- [6.1. Read InSAR Array and Mask Pixels with No Data](#array_mask)\n",
    "- [6.2. Randomly Sample and Pair Pixels](#array_sample)\n",
    "- [6.3. Secular Requirement Validation: Method 2](#secular_validation_method2)\n",
    "\n",
    "[**Appendix: Supplementary Comparisons and Plots**](#appendix)\n",
    "- [A.1. Compare Raw Velocities](#plot_vel)\n",
    "- [A.2. Plot Velocity Residuals](#plot_residual)\n",
    "- [A.3. Plot Double-differenced Residuals](#plot_ddiff)\n",
    "- [A.4. GNSS Timeseries Plots](#plot_gps_tseries)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "<br>\n",
    "<hr>\n",
    "\n",
    "<a id='#setup'></a>\n",
    "## Environment Setup"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "### Load Python Packages <a id='#load_packages'></a>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import glob\n",
    "import json\n",
    "import numpy as np\n",
    "import os\n",
    "import pandas as pd\n",
    "import subprocess\n",
    "\n",
    "from datetime import datetime as dt\n",
    "from matplotlib import pyplot as plt\n",
    "from mintpy.cli import view, plot_network\n",
    "from mintpy.objects import gnss, timeseries\n",
    "from mintpy.smallbaselineApp import TimeSeriesAnalysis\n",
    "from mintpy.utils import ptime, readfile, time_func, utils as ut\n",
    "from pathlib import Path\n",
    "from scipy import signal\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})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "### Define Calval Site and Parameters <a id='set_calval_params'></a>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Specify a calval location ID from my_sites.txt\n",
    "site = 'test' \n",
    "\n",
    "# Choose the requirement to validate\n",
    "# Options: 'Secular' 'Coseismic' 'Transient'\n",
    "requirement = 'Secular' \n",
    "\n",
    "# What dataset are you processing?\n",
    "# Options: \n",
    "#'ARIA_S1' (old directory structure for Sentinel-1 testing with aria-tools)\n",
    "#'ARIA_S1_new' (new directory structure for Sentinel-1 testing with aria-tools)\n",
    "dataset = 'ARIA_S1_new'\n",
    "aria_gunw_version = '3_0_1'\n",
    "\n",
    "# The date and version of this Cal/Val run\n",
    "rundate = '20240909'\n",
    "version = '1'\n",
    "\n",
    "# Provide the file where you keep your customized list of sites.\n",
    "custom_sites = '/home/jovyan/my_sites.txt'\n",
    "\n",
    "# Enter a username for storing your outputs\n",
    "if os.path.exists('/home/jovyan/me.txt'):\n",
    "    with open('/home/jovyan/me.txt') as m:\n",
    "        you = m.readline().strip()\n",
    "else:\n",
    "    you = input('Please type a username for your calval outputs:')\n",
    "    with open ('/home/jovyan/me.txt', 'w') as m: \n",
    "        m.write(you)\n",
    "\n",
    "# Load metadata for calval locations\n",
    "with open('/home/jovyan/my_sites.txt','r') as fid:\n",
    "    sitedata = json.load(fid)\n",
    "#sitedata['sites'][site]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Set Directories and Files <a id='set_directories'></a>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Directory location for Cal/Val data (do not change)\n",
    "start_directory = '/scratch/nisar-st-calval-solidearth' \n",
    "\n",
    "# Site directory\n",
    "site_dir = os.path.join(start_directory, dataset, site)\n",
    "\n",
    "# Working directory for calval processing\n",
    "work_dir = os.path.join(site_dir, requirement, you, rundate, 'v' + version)\n",
    "print(\"  Work directory:\", work_dir)\n",
    "\n",
    "# Directory for storing GUNW interferograms\n",
    "gunw_dir = os.path.join(site_dir,'products')\n",
    "print(\"  GUNW directory:\", gunw_dir) \n",
    "\n",
    "# Directory for storing MintPy outputs\n",
    "mintpy_dir = os.path.join(work_dir,'MintPy')\n",
    "print(\"MintPy directory:\", mintpy_dir)\n",
    "\n",
    "# Flag missing MintPy directory\n",
    "if not os.path.exists(mintpy_dir):\n",
    "    print()\n",
    "    print('ERROR: Stop! Your MintPy processing directory does not exist for this requirement, site, version, or date of your ATBD run.')\n",
    "    print('You may need to run the prep notebook first!')\n",
    "    print()\n",
    "else:\n",
    "    os.chdir(mintpy_dir)\n",
    "\n",
    "# Set MintPy filenames\n",
    "vel_file = os.path.join(mintpy_dir, 'velocity.h5')\n",
    "msk_file = os.path.join(mintpy_dir, 'maskConnComp.h5')  # maskTempCoh.h5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br>\n",
    "<hr>\n",
    "\n",
    "<a id='#prep_ifg'></a>\n",
    "## 1. Download and Prepare Interferograms \n",
    "Executed in *ARIA_prep* notebook"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "<br>\n",
    "<hr>\n",
    "\n",
    "<a id='#gen_ts'></a>\n",
    "## 2. Generation of Time Series from Interferograms"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "### 2.1. Validate/Modify Interferogram Network <a id='validate_network'></a>\n",
    "\n",
    "Add additional parameters to config_file in order to remove selected interferograms, change minimum coherence, etc."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "config_file = os.path.join(mintpy_dir, sitedata['sites'][site]['calval_location'] + '.cfg')\n",
    "command = 'smallbaselineApp.py ' + str(config_file) + ' --dostep modify_network'\n",
    "process = subprocess.run(command, shell=True)\n",
    "plot_network.main(['inputs/ifgramStack.h5'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "### 2.2. Generate Quality Control Mask <a id='generate_mask'></a>\n",
    "\n",
    "Mask files can be can be used to mask pixels in the time-series processing. Below we generate an initial mask file `maskConnComp.h5` based on the connected components for all the interferograms, which is a metric for unwrapping quality. After time-series analysis is complete, we will calculate a mask from the temporal coherence or variation of phase or displacement with time to make `maskTempCoh.h5`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "tags": []
   },
   "outputs": [],
   "source": [
    "command='generate_mask.py inputs/ifgramStack.h5  --nonzero  -o maskConnComp.h5  --update'\n",
    "process = subprocess.run(command, shell=True)\n",
    "view.main(['maskConnComp.h5', 'mask'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "### 2.3. Reference Interferograms To Common Lat/Lon <a id='common_latlon'></a>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "tags": []
   },
   "outputs": [],
   "source": [
    "command = 'smallbaselineApp.py ' + str(config_file) + ' --dostep reference_point'\n",
    "process = subprocess.run(command, shell=True)\n",
    "os.system('info.py inputs/ifgramStack.h5 | egrep \"REF_\"');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "### 2.4. Invert for SBAS Line-of-Sight Timeseries <a id='invert_SBAS'></a>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "tags": []
   },
   "outputs": [],
   "source": [
    "command = 'smallbaselineApp.py ' + str(config_file) + ' --dostep invert_network'\n",
    "process = subprocess.run(command, shell=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "<br>\n",
    "<hr>\n",
    "\n",
    "<a id='opt_correction'></a>\n",
    "## 3. Optional Corrections\n",
    "\n",
    "Phase distortions related to solid earth and ocean tidal effects, as well as those due to temporal variations in the vertical stratification of the atmosphere, can be mitigated using the approaches described below. At this point, it is expected that these corrections will not be needed to validate the mission requirements, but they may be used to produce the highest quality data products. Typically, these are applied to the estimated time series product rather than to the individual interferograms, since they are a function of the time of each radar acquisition."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.1. Solid Earth Tides Correction <a id='solid_earth'></a>\n",
    "\n",
    "[MintPy provides functionality for this correction, but it is not part of this notebook]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "### 3.2. Tropospheric Delay Correction <a id='tropo_corr'></a>\n",
    "\n",
    "Optional atmospheric correction utilizes the PyAPS (Jolivet et al., 2011, Jolivet and Agram, 2012) module within GIAnT (or eventually a merged replacement for GIAnT and MintPy). PyAPS is well documented, maintained and can be freely downloaded. PyAPS is included in GIAnT distribution). PyAPS currently includes support for ECMWF’s ERA-Interim, NOAA’s NARR and NASA’s MERRA weather models. A final selection of atmospheric models to be used for operational NISAR processing will be done during Phase C.\n",
    "\n",
    "[T]ropospheric delay maps are produced from atmospheric data provided by Global Atmospheric Models. This method aims to correct differential atmospheric delay correlated with the topography in interferometric phase measurements. Global Atmospheric Models (hereafter GAMs)... provide estimates of the air temperature, the atmospheric pressure and the humidity as a function of elevation on a coarse resolution latitude/longitude grid. In PyAPS, we use this 3D distribution of atmospheric variables to determine the atmospheric phase delay on each pixel of each interferogram.\n",
    "\n",
    "The absolute atmospheric delay is computed at each SAR acquisition date. For a pixel a_i at an elevation z at acquisition date i, the four surrounding grid points are selected and the delays for their respective elevations are computed. The resulting delay at the pixel a_i is then the bilinear interpolation between the delays at the four grid points. Finally, we combine the absolute delay maps of the InSAR partner images to produce the differential delay maps used to correct the interferograms.\n",
    "\n",
    "REFERENCE : https://github.com/insarlab/pyaps#2-account-setup-for-era5\n",
    "<br> Read Section 2 for ERA5 [link above] to create an account on the CDS website."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "if 'do_tropo_correction' in sitedata['sites']:\n",
    "    do_tropo_correction = sitedata['sites'][site]['do_tropo_correction']\n",
    "else:\n",
    "    do_tropo_correction = 'False'\n",
    "\n",
    "if do_tropo_correction == 'True':\n",
    "    if not Use_Staged_Data and not os.path.exists(Path.home()/'.cdsapirc'):\n",
    "        print('NEEDED to download ERA5, link: https://cds.climate.copernicus.eu/user/register')\n",
    "        UID = input('Please type your CDS_UID:')\n",
    "        CDS_API = input('Please type your CDS_API:')\n",
    "        \n",
    "        cds_tmp = '''url: https://cds.climate.copernicus.eu/api/v2\n",
    "        key: {UID}:{CDS_API}'''.format(UID=UID, CDS_API=CDS_API)\n",
    "        os.system('echo \"{cds_tmp}\" > ~/.cdsapirc; chmod 600 ~/.cdsapirc'.format(cds_tmp = str(cds_tmp)))\n",
    "    \n",
    "    command = 'smallbaselineApp.py ' + str(config_file) + ' --dostep correct_troposphere'\n",
    "    process = subprocess.run(command, shell=True)\n",
    "    \n",
    "    view.main(['inputs/ERA5.h5'])\n",
    "    timeseries_filename = 'timeseries_ERA5.h5'\n",
    "else:\n",
    "    timeseries_filename = 'timeseries.h5'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "### 3.3. Phase Deramping <a id='phase_deramp'></a>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "command = 'smallbaselineApp.py ' + str(config_file) + ' --dostep deramp'\n",
    "process = subprocess.run(command, shell=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "### 3.4. Topographic Residual Correction <a id='topo_corr'></a>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "command = 'smallbaselineApp.py ' + str(config_file) + ' --dostep correct_topography'\n",
    "process = subprocess.run(command, shell=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "<br>\n",
    "<hr>\n",
    "\n",
    "<a id='decomp_ts'></a>\n",
    "## 4. Estimate InSAR and GNSS Velocities\n",
    "The approach that will be used for the generation of NISAR L3 products for Requirements 660 and 663 allows for an explicit inclusion of key basis functions (e.g., Heaviside functions, secular rate, etc.) in the InSAR inversion. Modifications to this algorithm may be identified and implemented in response to NISAR Phase C activities. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "### 4.1. Estimate InSAR LOS Velocities <a id='insar_vel'></a>\n",
    "\n",
    "Given a time series of InSAR LOS displacements, the observations for a given pixel, $U(t)$, can be parameterized as:\n",
    "\n",
    "$$U(t) = a \\;+\\; vt \\;+\\; c_1 cos (\\omega_1t - \\phi_{1,}) \\;+\\; c_2 cos (\\omega_2t - \\phi_2) \\;+\\; \\sum_{j=1}^{N_{eq}} \\left( h_j+f_j F_j (t-t_j) \\right)H(t - t_j) \\;+\\; \\frac{B_\\perp (t)}{R sin \\theta}\\delta z \\;+\\; residual$$ \n",
    "\n",
    "which includes a constant offset $(a)$, velocity $(v)$, and amplitudes $(c_j)$ and phases $(\\phi_j)$ of annual $(\\omega_1)$ and semiannual $(\\omega_2)$ sinusoidal terms.  Where needed we can include additional complexity, such as coseismic and postseismic processes parameterized by Heaviside (step) functions $H$ and postseismic functions $F$ (the latter typically exponential and/or logarithmic).   $B_\\perp(t)$, $R$, $\\theta$, and $\\delta z$ are, respectively, the perpendicular component of the interferometric baseline relative to the first date, slant range distance, incidence angle and topography error correction for the given pixel. \n",
    "\n",
    "Thus, given either an ensemble of interferograms or the output of SBAS (displacement vs. time), we can write the LSQ problem as \n",
    "\n",
    "$$ \\textbf{G}\\textbf{m} = \\textbf{d}$$\n",
    "\n",
    "where $\\textbf{G}$ is the design matrix (constructed out of the different functional terms in Equation 2 evaluated either at the SAR image dates for SBAS output, or between the dates spanned by each pair for interferograms), $\\textbf{m}$ is the vector of model parameters (the coefficients in Equation 2) and $\\textbf{d}$ is the vector of observations.  (In the case of GNSS displacement timeseries, $\\textbf{G}, \\textbf{d}, \\textbf{m}$ are constructed using values evaluated at single epochs corresponding to the GNSS solution times, otherwise the estimate is the same.) \n",
    "\n",
    "With this formulation, we can obtain InSAR and GNSS velocity estimates and their formal uncertainties (including in areas where zero deformation is expected). The default InSAR (and GNSS) velocity fit in MintPy is to estimate a mean linear velocity $(v)$ in in the equation. To estimate a secular deformation rate, one can also estimate additional parameters for annual and semi-annual terms to account for seasonal periodic fluctuations in position. Below, we define a set of basis functions that MintPy will invert for, for both InSAR and GNSS data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Specify basis functions\n",
    "# Changing these might change the meaning of this exercise! Use caution if altering.\n",
    "ts_functions = {\n",
    "                'polynomial': 1,\n",
    "                'periodic': [1.0, 0.5],  # Periodic terms are in units of years\n",
    "                'stepDate': [],\n",
    "                'polyline': [],\n",
    "                'exp': {},\n",
    "                'log': {}\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Parse the TS model into argument string for command line execution\n",
    "function_str = ''\n",
    "\n",
    "# Loop through basis functions\n",
    "for fcn in ts_functions.keys():\n",
    "    arg = ts_functions[fcn]\n",
    "    \n",
    "    if fcn == 'polynomial':\n",
    "        function_str += f' --polynomial {arg:d}'\n",
    "    \n",
    "    else:\n",
    "        # Check whether values are passed\n",
    "        if len(arg) > 0:\n",
    "            # Append basis function name\n",
    "            function_str += f' --{fcn:s}'\n",
    "\n",
    "            # Append function arguments (e.g, period, decay, etc.)\n",
    "            if type(arg) == list:\n",
    "                function_str += ' {:s}'.format(' '.join([str(a) for a in arg]))\n",
    "            elif type(arg) == dict:\n",
    "                for event_date in arg.keys():\n",
    "                    function_str += f' {event_date:s}'\n",
    "                    for val in arg[event_date]:\n",
    "                        function_str += f' {val:f}'\n",
    "\n",
    "# Run command\n",
    "command = 'timeseries2velocity.py ' + timeseries_filename + function_str\n",
    "process = subprocess.run(command, shell=True)\n",
    "\n",
    "# load velocity file\n",
    "insar_velocities,_ = readfile.read(vel_file, datasetName = 'velocity')  # read velocity file\n",
    "insar_velocities = insar_velocities * 1000.  # convert InSAR velocities from m/yr to mm/yr\n",
    "\n",
    "# set masked pixels to NaN\n",
    "msk,_ = readfile.read(msk_file)\n",
    "insar_velocities[msk == 0] = np.nan\n",
    "insar_velocities[insar_velocities == 0] = np.nan"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we plot the mean linear velocity fit. The MintPy `view` module automatically reads the temporal coherence mask `maskTempCoh.h5` and applies that to mask out pixels with unreliable velocities (white)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "scp_args = 'velocity.h5 velocity -v -25 25 --colormap RdBu_r --figtitle LOS_Velocity --unit mm/yr -m ' + msk_file;\n",
    "view.main(scp_args.split());"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4.2. Not Used <a id='empty'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "### 4.3. Find Collocated GNSS Stations <a id='find_gps'></a>\n",
    "\n",
    "The project will have access to L2 position data for continuous GNSS stations in third-party networks (e.g., NSF’s Plate Boundary Observatory, the HVO network for Hawaii, GEONET-Japan, and GEONET-New Zealand), located in target regions for NISAR solid earth calval. Station data will be post-processed by one or more analysis centers, will be freely available, and will have latencies of several days to weeks, as is the case with positions currently produced by the NSF’s GAGE Facility and separately by the University of Nevada Reno. Networks will contain one or more areas of high-density station coverage (2~20 km nominal station spacing over 100 x 100 km or more) to support validation of L2 NISAR requirements at a wide range of length scales.\n",
    "\n",
    "MintPy supports several options for downloading and handling GNSS data. Currently supported data sets are:\n",
    "* UNR National Geodetic Laboratory\n",
    "* SIO/JPL MEaSUREs ESESES\n",
    "* JPL SIDESHOW"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# GNSS processing source\n",
    "if 'gnss_source' in sitedata['sites']:\n",
    "    gnss_source = sitedata['sites'][site]['gnss_source']\n",
    "else:\n",
    "    gnss_source = 'UNR'\n",
    "print(f\"GNSS processing source: {gnss_source:s}\")\n",
    "\n",
    "# Get analysis metadata from InSAR velocity file\n",
    "insar_metadata = readfile.read_attribute(vel_file)\n",
    "lat_step = float(insar_metadata['Y_STEP'])\n",
    "lon_step = float(insar_metadata['X_STEP'])\n",
    "(S,N,W,E) = ut.four_corners(insar_metadata)\n",
    "start_date = insar_metadata.get('START_DATE', None).split('T')[0]\n",
    "end_date = insar_metadata.get('END_DATE', None).split('T')[0]\n",
    "start_date_gnss = dt.strptime(start_date, \"%Y%m%d\")\n",
    "end_date_gnss = dt.strptime(end_date, \"%Y%m%d\")\n",
    "\n",
    "# Get center lat/lon of analysis region\n",
    "bbox = sitedata['sites'][site]['analysis_region'].replace(\"'\",\"\").split()\n",
    "lat0 = (float(bbox[0]) + float(bbox[1]))/2\n",
    "lon0 = (float(bbox[2]) + float(bbox[3]))/2\n",
    "# Identify geometry file x/y location associated with center lat/lon\n",
    "geom_obj = mintpy_dir + '/inputs/geometryGeo.h5'\n",
    "atr = readfile.read_attribute(geom_obj)\n",
    "coord = ut.coordinate(atr, lookup_file=geom_obj)\n",
    "y, x = coord.geo2radar(lat0, lon0)[0:2]\n",
    "y = max(0, y);  y = min(int(atr['LENGTH'])-1, y)\n",
    "x = max(0, x);  x = min(int(atr['WIDTH'])-1, x)\n",
    "# Write out inclination/azimuth at specified x/y location\n",
    "kwargs = dict(box=(x,y,x+1,y+1))\n",
    "inc_angle = readfile.read(geom_obj, datasetName='incidenceAngle', **kwargs)[0][0,0]\n",
    "az_angle  = readfile.read(geom_obj, datasetName='azimuthAngle',   **kwargs)[0][0,0]\n",
    "\n",
    "# Set GNSS Parameters\n",
    "gnss_completeness_threshold = 0.9    #0.9  #percent of data timespan with valid GNSS epochs\n",
    "gnss_residual_stdev_threshold = 10.  #0.03  #max threshold standard deviation of residuals to linear GNSS fit\n",
    "\n",
    "# Search for collocated GNSS stations\n",
    "site_names, site_lats, site_lons = gnss.search_gnss(SNWE=(S,N,W,E),\n",
    "                                                    start_date=start_date,\n",
    "                                                    end_date=end_date,\n",
    "                                                    source=gnss_source)\n",
    "site_names = [str(stn) for stn in site_names]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "### 4.4. Get GNSS Position Time Series <a id='gps_ts'></a>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "use_stn = []  #stations to keep\n",
    "bad_stn = []  #stations to toss\n",
    "use_lats = [] \n",
    "use_lons = []\n",
    "\n",
    "for counter, stn in enumerate(site_names):\n",
    "    gnss_obj = gnss.get_gnss_class(gnss_source)(site = stn, data_dir = os.path.join(mintpy_dir,'GNSS'))\n",
    "    gnss_obj.open(print_msg=False)\n",
    "    \n",
    "    # count number of dates in time range\n",
    "    dates = gnss_obj.dates\n",
    "    range_days = (end_date_gnss - start_date_gnss).days\n",
    "    gnss_count = np.histogram(dates, bins=[start_date_gnss, end_date_gnss])\n",
    "    gnss_count = int(gnss_count[0][0])\n",
    "\n",
    "    # for this quick screening check of data quality, we use constant incidence and azimuth angles \n",
    "    # get standard deviation of residuals to linear fit\n",
    "    disp_los = ut.enu2los(gnss_obj.dis_e, gnss_obj.dis_n, gnss_obj.dis_u, inc_angle, az_angle)\n",
    "    disp_detrended = signal.detrend(disp_los)\n",
    "    stn_stdv = np.std(disp_detrended)\n",
    "   \n",
    "    # select GNSS stations based on data completeness and scatter of residuals\n",
    "    disp_detrended = signal.detrend(disp_los)\n",
    "    if range_days*gnss_completeness_threshold <= gnss_count:\n",
    "        if stn_stdv > gnss_residual_stdev_threshold:\n",
    "            bad_stn.append(stn)\n",
    "        else:\n",
    "            use_stn.append(stn)\n",
    "            use_lats.append(site_lats[counter])\n",
    "            use_lons.append(site_lons[counter])\n",
    "    else:\n",
    "        bad_stn.append(stn)\n",
    "\n",
    "site_names = use_stn\n",
    "site_lats = use_lats\n",
    "site_lons = use_lons\n",
    "\n",
    "# [optional] manually remove additional stations\n",
    "gnss_to_remove=[]\n",
    "\n",
    "for i, gnss_site in enumerate(gnss_to_remove):\n",
    "    if gnss_site in site_names:\n",
    "        site_names.remove(gnss_site)\n",
    "    if gnss_site not in bad_stn:\n",
    "        bad_stn.append(gnss_site)\n",
    "\n",
    "print(\"\\nFinal list of {} stations used in analysis:\".format(len(site_names)))\n",
    "print(site_names)\n",
    "print(\"List of {} stations removed from analysis\".format(len(bad_stn)))\n",
    "print(bad_stn)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "### 4.5. Estimate GNSS LOS Velocities <a id='gps_vel'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We estimate the GNSS velocities in InSAR line of sight, by projecting the 3D displacements into LOS and fitting the resulting time-series using the same basis functions (i.e., linear and seasonal terms) that were used to fit the InSAR time-series. Outliers are removed from the GNSS position data and the functional fits are recomputed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Empty dicts and lists to store GNSS data\n",
    "gnss_stns = {}\n",
    "gnss_velocities = []\n",
    "gnss_velocity_errors = []\n",
    "\n",
    "# Loop through valid GNSS stations\n",
    "for i, site_name in enumerate(site_names):\n",
    "    # Retrieve station information\n",
    "    gnss_stn = gnss.get_gnss_class(gnss_source)(site = site_name)\n",
    "    gnss_stn.get_los_displacement(insar_metadata,\n",
    "                                  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",
    "    # Record filtered station to dictionary\n",
    "    gnss_stns[site_name] = gnss_stn\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",
    "    # Record station velocity\n",
    "    gnss_velocities.append(m_hat[1])\n",
    "    gnss_velocity_errors.append(mhat_se[1])\n",
    "\n",
    "    # Plotting parameters\n",
    "    n_dates = len(gnss_stn.dates)\n",
    "    label_skips = n_dates//6\n",
    "    ax_nb = i % 2\n",
    "\n",
    "    # Spawn figure for even numbers\n",
    "    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",
    "    \n",
    "    # Format plot\n",
    "    axes[ax_nb].legend()\n",
    "    axes[ax_nb].set_xticks(gnss_stn.dates[::label_skips])\n",
    "    axes[ax_nb].set_xticklabels([date.strftime('%Y-%m-%d') \\\n",
    "                                 for date in gnss_stn.dates[::label_skips]],\n",
    "                                rotation=80)\n",
    "    axes[ax_nb].set_ylabel('LOS displacement (mm)')\n",
    "    axes[ax_nb].set_title(f'{site_name:s} ({1000*gnss_velocities[-1]:.1f} mm/yr)')\n",
    "    fig.tight_layout()\n",
    "\n",
    "# Convert list of GNSS velocities to numpy array\n",
    "gnss_velocities = np.array(gnss_velocities)\n",
    "gnss_velocity_errors = np.array(gnss_velocity_errors)\n",
    "\n",
    "# Convert GNSS velocity m/yr to mm/yr\n",
    "gnss_velocities *= 1000\n",
    "gnss_velocity_errors *= 1000\n",
    "\n",
    "# Report GNSS velocities\n",
    "print('site velocity(mm/yr)')\n",
    "for site_name, vel, se in zip(site_names, gnss_velocities, gnss_velocity_errors):\n",
    "    print(f'{site_name:s} {vel:.1f} +- {se:.2f}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "### 4.6. Re-reference GNSS and InSAR <a id='reference'></a>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# reference GNSS stations to GNSS reference site\n",
    "print(\"Using GNSS reference station: \", sitedata['sites'][site]['gps_ref_site_name'])\n",
    "ref_site_ind = site_names.index(sitedata['sites'][site]['gps_ref_site_name'])\n",
    "gnss_velocities = gnss_velocities - gnss_velocities[ref_site_ind]\n",
    "\n",
    "# reference InSAR to GNSS reference site\n",
    "ref_site_lat = float(site_lats[ref_site_ind])\n",
    "ref_site_lon = float(site_lons[ref_site_ind])\n",
    "ref_y, ref_x = ut.coordinate(insar_metadata).geo2radar(ref_site_lat, ref_site_lon)[:2]\n",
    "insar_velocities = insar_velocities - insar_velocities[ref_y, ref_x]\n",
    "\n",
    "# plot GNSS stations on InSAR velocity field\n",
    "vmin, vmax = -25, 25\n",
    "cmap = plt.get_cmap('RdBu_r')\n",
    "\n",
    "fig, ax = plt.subplots(figsize=[18, 5.5])\n",
    "cax = ax.imshow(insar_velocities, cmap=cmap, vmin=vmin, vmax=vmax, interpolation='nearest', extent=(W, E, S, N))\n",
    "cbar = fig.colorbar(cax, ax=ax)\n",
    "cbar.set_label('LOS velocity [mm/year]')\n",
    "\n",
    "for lat, lon, obs in zip(site_lats, site_lons, gnss_velocities):\n",
    "    color = cmap((obs - vmin)/(vmax - vmin))\n",
    "    ax.scatter(lon, lat, color=color, s=8**2, edgecolors='k')\n",
    "for i, label in enumerate(site_names):\n",
    "     plt.annotate(label, (site_lons[i], site_lats[i]), color='black')\n",
    "\n",
    "out_fig = os.path.abspath('vel_insar_vs_gnss.png')\n",
    "fig.savefig(out_fig, bbox_inches='tight', transparent=True, dpi=300)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-info\">\n",
    "<b>Note :</b> \n",
    "Negative values indicates that target is moving away from the radar (i.e., subsidence in case of vertical deformation).\n",
    "Positive values indicates that target is moving towards the radar (i.e., uplift in case of vertical deformation). \n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "<br>\n",
    "<hr>\n",
    "\n",
    "<a id='validation1'></a>\n",
    "## 5. Validation Method 1: GNSS-InSAR Direct Comparison "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "### 5.1. Make GNSS-InSAR Velocity Residuals at GNSS Station Locations <a id='make_resid'></a>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "#Set Parameters\n",
    "pixel_radius = 5   #number of InSAR pixels to average for comparison with GNSS\n",
    "\n",
    "#Create dictionary with the stations as the key and all their info as an array \n",
    "stn_dict = {}\n",
    "\n",
    "#Loop over GNSS station locations\n",
    "for i in range(len(site_names)): \n",
    "    \n",
    "    # convert GNSS station lat/lon information to InSAR x/y grid\n",
    "    stn_lat = site_lats[i]\n",
    "    stn_lon = site_lons[i]\n",
    "    x_value = round((stn_lon - W)/lon_step)\n",
    "    y_value = round((stn_lat - N)/lat_step)\n",
    "    \n",
    "    # get velocities and residuals\n",
    "    gnss_site_vel = gnss_velocities[i]\n",
    "    #Caution: If you expand the radius parameter farther than the bounding grid it will break. \n",
    "    #To fix, remove the station in section 4 when the site_names list is filtered\n",
    "    vel_px_rad = insar_velocities[y_value-pixel_radius:y_value+1+pixel_radius, \n",
    "                     x_value-pixel_radius:x_value+1+pixel_radius]\n",
    "    insar_site_vel = np.median(vel_px_rad)\n",
    "    residual = gnss_site_vel - insar_site_vel\n",
    "\n",
    "    # populate data structure\n",
    "    values = [x_value, y_value, insar_site_vel, gnss_site_vel, residual, stn_lat, stn_lon]\n",
    "    stn = site_names[i]\n",
    "    stn_dict[stn] = values\n",
    "\n",
    "# extract data from structure\n",
    "res_list = []\n",
    "insar_site_vels = []\n",
    "gnss_site_vels = []\n",
    "lat_list = []\n",
    "lon_list = []\n",
    "for i in range(len(site_names)): \n",
    "    stn = site_names[i]\n",
    "    insar_site_vels.append(stn_dict[stn][2])\n",
    "    gnss_site_vels.append(stn_dict[stn][3])\n",
    "    res_list.append(stn_dict[stn][4])\n",
    "    lat_list.append(stn_dict[stn][5])\n",
    "    lon_list.append(stn_dict[stn][6])\n",
    "num_stn = len(site_names) \n",
    "print('Finish creating InSAR residuals at GNSS sites')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "### 5.2. Make Double-Differenced Velocity Residuals <a id='make_ddiff'></a>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "diff_res_list = []\n",
    "stn_dist_list = []\n",
    "n_gnss_sites = len(site_names)\n",
    "\n",
    "# loop over stations\n",
    "for i in range(n_gnss_sites - 1):\n",
    "    stn1 = site_names[i]\n",
    "    for j in range(i + 1, n_gnss_sites):\n",
    "        stn2 = site_names[j]\n",
    "\n",
    "        # calculate GNSS and InSAR velocity differences between stations\n",
    "        gnss_vel_diff = stn_dict[stn1][3] - stn_dict[stn2][3]\n",
    "        insar_vel_diff = stn_dict[stn1][2] - stn_dict[stn2][2]\n",
    "\n",
    "        # calculate GNSS vs InSAR differences (double differences) between stations\n",
    "        diff_res = gnss_vel_diff - insar_vel_diff\n",
    "        diff_res_list.append(diff_res)\n",
    "\n",
    "        # get distance (km) between stations using Haversine formula\n",
    "        # index 5 is lat, 6 is lon\n",
    "        stn_dist = haversine_distance(stn_dict[stn1][6], stn_dict[stn1][5], stn_dict[stn2][6], stn_dict[stn2][5])\n",
    "        stn_dist_list.append(stn_dist)\n",
    "\n",
    "# Write data for statistical tests\n",
    "gnss_site_dist = np.array(stn_dist_list)\n",
    "double_diff_rel_measure = np.array(np.abs(diff_res_list))\n",
    "ndx = np.argsort(gnss_site_dist)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "### 5.3. Secular Requirement Validation: Method 1 <a id='secular_validation_method1'></a>\n",
    "\n",
    "We assume that the distribution of residuals is Gaussian and that the requirement success threshold represents a 1-sigma limit within which we expect 68.3% of residuals to lie."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define requirement\n",
    "secular_threshold_rqmt = 2  # mm/yr\n",
    "secular_distance_rqmt = [0.1, 50.0]  # km\n",
    "\n",
    "# Validation parameters\n",
    "n_bins = 10\n",
    "threshold = 0.683  \n",
    "\n",
    "# Validation figure and assessment\n",
    "validation_table, fig = display_validation(gnss_site_dist,                      # binned distance for point\n",
    "                                           double_diff_rel_measure,             # binned double-difference velocities mm/yr\n",
    "                                           site,                                # cal/val site name\n",
    "                                           start_date,                          # start date of InSAR dataset\n",
    "                                           end_date,                            # end date of InSAR dataset \n",
    "                                           requirement=secular_threshold_rqmt,  # measurement requirement to meet, e.g 2 mm/yr for 3 years of data over 0.1-50km\n",
    "                                           distance_rqmt=secular_distance_rqmt, # distance over requirement is to meet, e.g. over length scales of 0.1-50 km [0.1, 50] \n",
    "                                           n_bins=n_bins,                       # number of bins, to collect statistics \n",
    "                                           threshold=threshold,                 # quantile threshold for point-pairs that pass requirement, e.g. 0.683 - we expect 68.3% of residuals to lie. \n",
    "                                           sensor='Sentinel-1',                 # sensor that is validated, Sentinel-1 or NISAR\n",
    "                                           validation_type='secular',           # validation for: secular, transient, coseismic requirement\n",
    "                                           validation_data='GNSS')              # validation method: GNSS - Method 1, InSAR - Method 2\n",
    "\n",
    "out_fig = os.path.abspath('secular_insar-gnss_velocity_vs_distance.png')\n",
    "fig.savefig(out_fig, bbox_inches='tight', transparent=True, dpi=30)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "display_validation_table(validation_table)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "<div class=\"alert alert-warning\">\n",
    "Validation Method 1: success for a baseline distance bin occurs when the percentage of residuals within the requirement success threshold is >0.683\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br>\n",
    "<hr>\n",
    "\n",
    "<a id='validation2'></a>\n",
    "## 6. Validation Method 2: InSAR-only Structure Function"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In Validation Method 2, we use a time interval and area where we assume no deformation.  As with Method 1, we assume that the distribution of residuals is Gaussian and that the requirement success threshold represents a 1-sigma limit within which we expect 68.3% of residuals to lie."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 6.1. Read InSAR Array and Mask Pixels with no Data <a id='array_mask'></a>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# use the assumed non-earthquake displacement as the insar_displacment for statistics and convert to mm\n",
    "insar_velocities,_ = readfile.read(vel_file, datasetName = 'velocity')  #read velocity\n",
    "velStart = sitedata['sites'][site]['download_start_date']\n",
    "insar_velocities = insar_velocities * 1000.  # convert velocities from m to mm\n",
    "\n",
    "# set masked pixels to NaN\n",
    "msk,_ = readfile.read(msk_file)\n",
    "insar_velocities[msk == 0] = np.nan\n",
    "insar_velocities[insar_velocities == 0] = np.nan\n",
    "\n",
    "# display map of velocity data after masking\n",
    "cmap = plt.get_cmap('RdBu_r')\n",
    "fig, ax = plt.subplots(figsize=[18, 5.5])\n",
    "img1 = ax.imshow(insar_velocities, cmap=cmap, vmin=-20, vmax=20, interpolation='nearest', extent=(W, E, S, N))\n",
    "ax.set_title(\"Secular \\n Date \"+velStart)\n",
    "cbar1 = fig.colorbar(img1, ax=ax)\n",
    "cbar1.set_label('LOS velocity [mm/year]')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 6.2. Randomly Sample Pixels and Pair Them Up with Option to Remove Trend <a id='remove_trend'></a>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define requirement\n",
    "secular_threshold_rqmt = 2  # mm/yr\n",
    "secular_distance_rqmt = [0.1, 50.0]  # km\n",
    "\n",
    "sample_mode = 'points'  # 'points' or 'profile'\n",
    "# note that the 'profile' method may take significantly longer\n",
    "\n",
    "# Collect samples using the specified method\n",
    "if sample_mode in ['points']:\n",
    "    X0,Y0 = load_geo(insar_metadata)\n",
    "    X0_2d,Y0_2d = np.meshgrid(X0,Y0)\n",
    "\n",
    "    insar_sample_dist, insar_rel_measure = samp_pair(X0_2d, Y0_2d, insar_velocities, num_samples=1000000)\n",
    "\n",
    "elif sample_mode in ['profile']:\n",
    "    # Sample grid setup\n",
    "    length, width = int(insar_metadata['LENGTH']), int(insar_metadata['WIDTH'])\n",
    "    X = np.linspace(W+lon_step, E-lon_step, width)  # longitudes\n",
    "    Y = np.linspace(N+lat_step, S-lat_step, length)  # latitudes\n",
    "    X_coords, Y_coords = np.meshgrid(X, Y)\n",
    "\n",
    "    # Draw random samples from map (without replacement)\n",
    "    num_samples = 20000\n",
    "    \n",
    "    # Retrieve profile samples\n",
    "    insar_sample_dist, insar_rel_measure = profile_samples(\\\n",
    "                    x=X_coords.reshape(-1,1),\n",
    "                    y=Y_coords.reshape(-1,1),\n",
    "                    data=insar_velocities,\n",
    "                    metadata=insar_metadata,\n",
    "                    len_rqmt=secular_distance_rqmt,\n",
    "                    num_samples=num_samples)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 6.3. Secular Requirement Validation: Method 2  <a id='secular_validation_method2'></a>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Validation parameters\n",
    "n_bins = 10\n",
    "threshold = 0.683  \n",
    "\n",
    "# Validation figure and assessment\n",
    "validation_table, fig = display_validation(insar_sample_dist,              # binned distance for point\n",
    "                                           insar_rel_measure,              # binned relative velocities mm/yr\n",
    "                                           site,                           # cal/val site name\n",
    "                                           start_date,                     # start date of InSAR dataset\n",
    "                                           end_date,                       # end date of InSAR dataset \n",
    "                                           requirement=secular_threshold_rqmt,  # measurement requirement to meet, e.g 2 mm/yr for 3 years of data over 0.1-50km\n",
    "                                           distance_rqmt=secular_distance_rqmt,  # distance over requirement is to meet, e.g. over length scales of 0.1-50 km [0.1, 50] \n",
    "                                           n_bins=n_bins,                  # number of bins, to collect statistics \n",
    "                                           threshold=threshold,            # quantile threshold for point-pairs that pass requirement, e.g. 0.683 - we expect 68.3% of residuals to lie. \n",
    "                                           sensor='Sentinel-1',            # sensor that is validated, Sentinel-1 or NISAR\n",
    "                                           validation_type='secular',      # validation for: secular, transient, coseismic requirement\n",
    "                                           validation_data='InSAR')        # validation method: GNSS - Method 1, InSAR - Method 2\n",
    "\n",
    "out_fig = os.path.abspath('secular_insar-only_vs_distance_'+site+'_date'+velStart+'.png')\n",
    "fig.savefig(out_fig, bbox_inches='tight', transparent=True, dpi=300)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "display_validation_table(validation_table)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-warning\">\n",
    "Final result Method 2—\n",
    "    68% of points below the requirements line is success\n",
    "</div>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "<br>\n",
    "<hr>\n",
    "\n",
    "<a id='appendix'></a>\n",
    "## Appendix: Supplementary Comparisons and Plots"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "### A.1. Compare Raw Velocities <a id='compare_raw'></a>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "vmin, vmax = -25, 25\n",
    "plt.figure(figsize=(10,3))\n",
    "plt.hist(insar_site_vels, range=[vmin, vmax], bins=50, color=\"green\", edgecolor='grey', label='V_InSAR')\n",
    "plt.hist(gnss_site_vels, range=[vmin, vmax], bins=50, color=\"orange\", edgecolor='grey', label='V_gnss', alpha=0.5)\n",
    "plt.legend(loc='upper right')\n",
    "plt.title(f\"Velocities \\n Date range {start_date}-{end_date} \\n Reference stn: {sitedata['sites'][site]['gps_ref_site_name']} \\n Number of stations used: {num_stn}\")\n",
    "plt.xlabel('LOS Velocity (mm/year)')\n",
    "plt.ylabel('N Stations')\n",
    "plt.ylim(0,20)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "### A.2. Plot Velocity Residuals <a id='plot_residuals'></a>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "vmin, vmax = -10, 10\n",
    "plt.figure(figsize=(10,3))\n",
    "plt.hist(res_list, bins = 40, range=[vmin,vmax], edgecolor='grey', color=\"darkblue\", linewidth=1, label='V_gnss - V_InSAR (area average)')\n",
    "plt.legend(loc='upper right')\n",
    "plt.title(f\"Residuals \\n Date range {start_date}-{end_date} \\n Reference stn: {sitedata['sites'][site]['gps_ref_site_name']} \\n Number of stations used: {num_stn}\")\n",
    "plt.xlabel('Velocity Residual (mm/year)')\n",
    "plt.ylabel('N Stations')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "### A.3. Plot Double Difference Residuals <a id='plot_ddiff'></a>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "plt.figure(figsize=(10,3))\n",
    "plt.hist(diff_res_list, range = [vmin, vmax],bins = 40, color = \"darkblue\",edgecolor='grey',label='V_gnss_(s1-s2) - V_InSAR_(s1-s2)')\n",
    "plt.legend(loc='upper right')\n",
    "plt.title(f\"Difference Residualts \\n Date range {start_date}-{end_date} \\n Reference stn: {sitedata['sites'][site]['gps_ref_site_name']} \\n Number of stations used: {num_stn}\")\n",
    "plt.xlabel('Double Differenced Velocity Residual (mm/year)')\n",
    "plt.ylabel('N Stations')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "### A.4. GNSS Timeseries Plots <a id='plot_gps_tseries'></a>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# grab the time-series file used for time function estimation given the template setup\n",
    "template = readfile.read_template(os.path.join(mintpy_dir, 'smallbaselineApp.cfg'))\n",
    "template = ut.check_template_auto_value(template)\n",
    "insar_ts_file = TimeSeriesAnalysis.get_timeseries_filename(template, mintpy_dir)['velocity']['input']\n",
    "\n",
    "# read the time-series file\n",
    "insar_ts, atr = readfile.read(insar_ts_file, datasetName='timeseries')\n",
    "mask = readfile.read(os.path.join(mintpy_dir, 'maskTempCoh.h5'))[0]\n",
    "print(f'reading timeseries from file: {insar_ts_file}')\n",
    "\n",
    "# Get date list\n",
    "date_list = timeseries(insar_ts_file).get_date_list()\n",
    "num_date = len(date_list)\n",
    "date0, date1 = date_list[0], date_list[-1]\n",
    "insar_dates = ptime.date_list2vector(date_list)[0]\n",
    "\n",
    "# spatial reference\n",
    "coord = ut.coordinate(atr)\n",
    "ref_site = sitedata['sites'][site]['gps_ref_site_name']\n",
    "ref_gnss_obj = gnss_stns[ref_site]\n",
    "ref_lat, ref_lon = ref_gnss_obj.get_site_lat_lon()\n",
    "ref_y, ref_x = coord.geo2radar(ref_lat, ref_lon)[:2]\n",
    "if not mask[ref_y, ref_x]:\n",
    "    raise ValueError(f'Given reference GNSS site ({ref_site}) is in mask-out unrelible region in InSAR! Change to a different site.')\n",
    "ref_insar_dis = insar_ts[:, ref_y, ref_x]\n",
    "\n",
    "# Plot displacements and velocity timeseries at GNSS station locations\n",
    "num_site = len(site_names)\n",
    "prog_bar = ptime.progressBar(maxValue=num_site)\n",
    "for i, site_name in enumerate(site_names):\n",
    "    prog_bar.update(i+1, suffix=f'{site_name} {i+1}/{num_site}')\n",
    "\n",
    "    ## read data\n",
    "    # recall gnss station displacements with outliers removed\n",
    "    gnss_obj = gnss_stns[site_name]\n",
    "    gnss_lalo = (gnss_obj.site_lat, gnss_obj.site_lon)\n",
    "\n",
    "    # get relative LOS displacement on common dates\n",
    "    gnss_dates = np.array(sorted(list(set(gnss_obj.dates) & set(ref_gnss_obj.dates))))\n",
    "    gnss_dis = np.zeros(gnss_dates.shape, dtype=np.float32)\n",
    "    for i, date_i in enumerate(gnss_dates):\n",
    "        idx1 = np.where(gnss_obj.dates == date_i)[0][0]\n",
    "        idx2 = np.where(ref_gnss_obj.dates == date_i)[0][0]\n",
    "        gnss_dis[i] = gnss_obj.dis_los[idx1] - ref_gnss_obj.dis_los[idx2]\n",
    "    \n",
    "    # shift GNSS to zero-mean in time [for plotting purpose]\n",
    "    gnss_dis -= np.nanmedian(gnss_dis)\n",
    "\n",
    "    # read InSAR\n",
    "    y, x = coord.geo2radar(gnss_lalo[0], gnss_lalo[1])[:2]\n",
    "    insar_dis = insar_ts[:, y, x] - ref_insar_dis\n",
    "    # apply a constant shift in time to fit InSAR to GNSS\n",
    "    comm_dates = sorted(list(set(gnss_dates) & set(insar_dates)))\n",
    "    if comm_dates:\n",
    "        insar_flag = [x in comm_dates for x in insar_dates]\n",
    "        gnss_flag = [x in comm_dates for x in gnss_dates]\n",
    "        insar_dis -= np.nanmedian(insar_dis[insar_flag] - gnss_dis[gnss_flag])\n",
    "\n",
    "    ## plot figure\n",
    "    if gnss_dis.size > 0 and np.any(~np.isnan(insar_dis)):\n",
    "        fig, ax = plt.subplots(figsize=(10, 3))\n",
    "        ax.axhline(color='grey',linestyle='dashed', linewidth=2)\n",
    "        ax.scatter(gnss_dates, gnss_dis*1000, s=2**2, label=\"GNSS Daily Positions\")\n",
    "        ax.scatter(insar_dates, insar_dis*1000, label=\"InSAR Positions\")\n",
    "        # axis format\n",
    "        ax.set_title(f\"Station Name: {site_name}\") \n",
    "        ax.set_ylabel('LOS displacement [mm]')\n",
    "        ax.legend()\n",
    "prog_bar.close()\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "hide_input": false,
  "kernelspec": {
   "display_name": "solid_earth_atbd",
   "language": "python",
   "name": "solid_earth_atbd"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.5"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": false,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  },
  "toc-autonumbering": false,
  "toc-showmarkdowntxt": false
 },
 "nbformat": 4,
 "nbformat_minor": 4
}