smallbaselineApp_aria.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p align=\"left\">\n",
    "  <img width=\"100%\" src=\"docs/Banner.png\">\n",
    "</p>    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "# InSAR Time Series Analysis: MintPy + ARIA GUNW products\n",
    "\n",
    "**Author:** Zhang Yunjun, Heresh Fattahi, August 12-16, 2019 @UNAVCO\n",
    "\n",
    "The Miami Insar Timeseries software in PYthon (MintPy) is an open-source package for InSAR time-series analysis. Mintpy currently starts from stacks of unwrapped interferograms (in either geo or radar coordinates) and estimates ground displacement time-series. Mintpy is primerily consistent with stacks of interferograms processed with ISCE. However, the software also supports interfarograms processed with other InSAR processors such as GAMMA, SNAP and ROI_PAC. \n",
    "\n",
    "MintPy is available on Github from the following page: https://github.com/insarlab/MintPy\n",
    "\n",
    "References: The detailed algorithms implemented in MintPy can be found in the following manuscript: \n",
    "\n",
    "+ Yunjun, Z., H. Fattahi, F. Amelung (2019), Small baseline InSAR time series analysis: Unwrapping error correction and noise reduction, Computers & Geosciences, 133, 104331, doi:[10.1016/j.cageo.2019.104331](https://doi.org/10.1016/j.cageo.2019.104331), [ArXiv](https://eartharxiv.org/9sz6m/)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "# 0. Initial setup of the notebook"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "The cell below performs the intial setup of the notebook and must be run every time the notebook (re)starts. It defines the processing location and check the example dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import os\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from mintpy import view, tsview, plot_network, plot_transection, plot_coherence_matrix\n",
    "\n",
    "# define the work directory\n",
    "work_dir = os.path.abspath(os.getcwd())\n",
    "\n",
    "#The Path to the folder where ARIA data after running ariaTSsetup.py is stored\n",
    "aria_dir = os.path.join(os.getenv(\"HOME\"),'work/04_ariaTSprep')\n",
    "\n",
    "print(\"Work directory: \", work_dir)\n",
    "print(\"The path to the directory with ARIA data after running ariaTSsetup:\", aria_dir)\n",
    "\n",
    "if not os.path.isdir(work_dir):\n",
    "    os.makedirs(work_dir)\n",
    "    print('Create directory: {}'.format(work_dir))\n",
    "\n",
    "print('Go to work directory: {}'.format(work_dir))\n",
    "os.chdir(work_dir)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "# 1. smallbaselineApp.py (general overview)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "This application provides a workflow which includes several steps to invert a stack of unwrapped interferograms and apply different corrections to obtain ground displacement timeseries.  \n",
    "The workflow consists of two main blocks:\n",
    "\n",
    "* correcting unwrapping errors and inverting for the raw phase time-series (blue ovals),\n",
    "* correcting for noise from different sources to obtain the displacement time-series (green ovals).\n",
    "\n",
    "Some steps are optional, which are switched off by default (marked by dashed boundaries). Configuration parameters for each step are initiated with default values in a customizable text file: [smallbaselineApp.cfg](https://github.com/insarlab/MintPy/blob/master/mintpy/defaults/smallbaselineApp.cfg). In this notebook, we will walk through the various steps."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "## 1.1 Processing steps of smallbaselineApp.py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "The smallbaselineApp.py workflow can be called with a single command-line call; by default it will run all the required processing steps with options pulled from the template files. However, in this notebook, we will use the \"step\" processing, which allows to re-start the processing from a given step. More detailed usage can be found in help."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "!smallbaselineApp.py --help"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "The app includes the following processing steps:\n",
    "\n",
    "'load_data', 'modify_network', 'reference_point', 'correct_unwrap_error', 'stack_interferograms',    \n",
    "'invert_network', 'correct_LOD', 'correct_troposphere', 'deramp', 'correct_topography',    \n",
    "'residual_RMS', 'reference_date', 'velocity', 'geocode', 'google_earth', 'hdfeos5'\n",
    "\n",
    "#### Input\n",
    "\n",
    "* `load_data:` loads the stack unwrapped interferograms, coherence files, connected componenents, baseline, incidence angle and other metadata in to HDF5 files with multiple datasets and attributes.  \n",
    "\n",
    "#### Network inversion\n",
    "\n",
    "* `modify_network:` this step allows to modify the network of interferograms if needed. e.g., the network can be modified based on temporal and spatial baselines threshold, or can be modified by removing specific pairs.\n",
    "* `reference_point:` The unwrapped interferograms may be relative to differnt reference pixels. This step introduces a common refernce pixel to all interferograms. For intuitive interpretation, one may choose a stable coherent non-deforming pixel. However, since the estimated InSAR displacement time-series is relative in both time and space, then choosing a deforming pixel does not change the results.\n",
    "* `correct_unwrap_error:` Input unwrapped interferograms may be affected by phase unwrapping errors (wrong $2\\pi$ phase jumps added during phase unwrapping). This step offers three methods to possibly identify and correct phase unwrapping errors.\n",
    "* `stack_interferograms:` This step allos to simply average the unwrapped interferograms. This can be used to have a quick look of expected rate maps even before inversion. \n",
    "* `invert_network:` inverts the stack of unwrapped interferograms to form InSAR phase time-series. This is equivalent to transforming the network of small-baseline interferograms to a single-master network of interferogram (i.e., the unwrapped phase timeseries). \n",
    "\n",
    "#### Noise reduction of displacement time-series\n",
    "\n",
    "* `correct_LOD:` this step is specific to Envisat data and applies an empirucal correction to account for possible local oscilator drift of the radar.\n",
    "* `correct_troposphere:` corrects tropospheric delay using atmospheric models and with empirical phase elevation approach estimated from InSAR data.\n",
    "* `deramp:` This step (if requested) removes a ramp from each acquisition. Note that deramping removes residual long-wavelength interferometric phase components which may be dute to noise (geometrical residual, atmospheric delay) or signal (due to tectonic deformation).  \n",
    "* `correct_topography:` estimates residual topographic effects which are correlated with temporal variation of perpendicular baseline.  \n",
    "* `residual_RMS:` estimates the average noise level for each acquisition by calculating the RMS of the residual phase\n",
    "* `reference_date:` change reference date\n",
    "* `velocity:` estimates the average velocity\n",
    "\n",
    "#### Output\n",
    "\n",
    "* `geocode:` if the original stack in radar-coordinates, convert it to geo-coordinates in lat/lon\n",
    "* `google_earth:` output the average velocity into an Google Earth KMZ file.\n",
    "* `hdfeos5:` output the displacement time-series with geometry info into one file in [HDF-EOS5](http://hdfeos.org) format."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "## 1.2 Configuring the processing parameters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "The processing parameters for the smallbaselineApp.py are configured through two configuration files. At least one configuration is required to run smallbaselineApp.py.\n",
    "\n",
    "* `default configuration`: [smallbaselineApp.cfg](https://github.com/insarlab/MintPy/blob/master/mintpy/defaults/smallbaselineApp.cfg). It contains all configuration parameters, grouped by steps, with default auto values (which are defined in [smallbaselineApp_auto.cfg](https://github.com/insarlab/MintPy/blob/master/mintpy/defaults/smallbaselineApp_auto.cfg)). This file is copied over to the current working directory and read every time smallbaselineApp.py runs.\n",
    "* `custom configuration` (optional but recommended): `SanFranSenDT42.txt` in the example dataset. It constains selective, manually modified configuration parameters. The custom template file name is arbitrary. Custom template has higher priority than the default template; if custom template is input, smallbaselineApp.py will update the default smallbaselineApp.cfg file accordingly."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "### custom configuration for the dataset in this notebook"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Create a text file named SanFranSenDT42.txt and copy the following few lines in that file: \n",
    "\n",
    "```cfg\n",
    "mintpy.reference.yx                          = 490,630\n",
    "mintpy.network.coherenceBased                = no\n",
    "mintpy.networkInversion.weightFunc           = no\n",
    "mintpy.topographicResidual.pixelwiseGeometry = no\n",
    "mintpy.troposphericDelay.method              = no\n",
    "mintpy.deramp                                = no\n",
    "mintpy.topographicResidual                   = no\n",
    "```\n",
    "    \n",
    "Check **more examples for custom template** on [here](https://github.com/insarlab/MintPy/tree/master/docs/examples/input_files).      \n",
    "Run **load_data.py -H** for example input file paths for ISCE direct outputs as well as outputs from GAMMA, SNAP and ROI_PAC.   "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "# 2. Time Series Analysis (small baseline approach)  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "## Load the ARIA data into MintPy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "MintPy is most consistent with the ISCE direct outputs. However, it supports interferograms processed with other InSAR software including Gamma and SNAP. In this tutorial we are not using the direct ISCE outputs in their flat binary formats, but rather we use the ISCE outputs packaged in an HDF5 file, which have been produced by ARIA. Therefore we skip the conventional load data step in MintPy and instead manually run a script called prep_aria.py to load the ARIA products into MintPy. In other words prep_aria.py creates HDF5 files which contains the stack of unwrapped interferograms, coherence, connected components and required metdata which are used as inputs to MintPy.\n",
    "\n",
    "If the user has already prepared the ARIA data, then the data preparation can be skipped. Please see the [example data directory](https://mintpy.readthedocs.io/en/latest/dir_structure/) about setting up the stack for MintPy processing. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "#ariaDownload.py -b '37.25 38.1 -122.6 -121.75' --track 42\n",
    "#ariaTSsetup.py -f 'products/*.nc' -b '37.25 38.1 -122.6 -121.75' --mask Download\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "<br>\n",
    "<div class=\"alert alert-info\">\n",
    "<b>Note :</b> \n",
    "smallbaslineApp.py currently does not load ARIA products. Therefore users need to run prep_aria.py script before running the smallbaselineApp.py inorder to load the products in a single HDF5 file that MintPy expects.  \n",
    "</div>\n",
    "To load the ARIA interferograms into mintpy run the following command:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "!prep_aria.py -s ../04_ariaTSprep/stack/ -i ../04_ariaTSprep/incidenceAngle/20150605_20150512.vrt -d ../04_ariaTSprep/DEM/SRTM_3arcsec.dem -a ../04_ariaTSprep/azimuthAngle/20150605_20150512.vrt --water-mask ../04_ariaTSprep/mask/watermask.msk"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "By running this command, the \"inputs\" directory inside the working directory is created and two HDF5 files are produced as"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "ls inputs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "ifgramStack.h5: This file contains 6 dataset cubes and multiple metadata. The datasets include Perpendicular baselines (average value), coherence data cube, connected commponents cube, unwrapped interferometric phases, \"dropIfgram\" a boolean dataset which indicates weather an interferogram is used for inversion or ignored and the \"date\" dataset with a shape of \"number of interferograms\" by 2 where the first column is the reference date and the 2nd column is the secondary dates of the interferogram. \n",
    "\n",
    "geometryGeo.h5: This file contains geometrical datasets including incidence angle, azimuth angle, shadow layover mask, slant range distance and height. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "<div class=\"alert alert-info\">\n",
    "<b>info.py :</b> \n",
    "To get general infomation about a MintPy product, run info.py on the file.   \n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "!info.py inputs/ifgramStack.h5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "!info.py inputs/geometryGeo.h5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "## Plot the interferogram network"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Before inversion, its useful to take a look at the network of interferograms. Running plot_network.py gives an overview of the network and the average coherence of the stack. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "plot_network.main(['inputs/ifgramStack.h5'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Note that with the --nodisplay argument, the plots won't be displayed but saved as pdf files in the current directory. Running this command creates multiple files as follows:\n",
    "\n",
    "ifgramStack_coherence_spatialAvg.txt: A simple text file that provides an overview to the stack and contains the interferogram dates, average coherence temporal and spatial baseline separation.\n",
    "\n",
    "Network.pdf: Displays the network of interferograms on time-baseline coordinates, colorcoded by avergae coherence of the interferograms. Circles represent the acquisition dates and lines represent the interferograms. Solid lines are the interferograms used for time-series analysis and dashed line are the interferograms ignored in the time-series analysis. \n",
    "\n",
    "CoherenceMatrix.pdf shows the avergae coherence pairs between all available pairs in the stack.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "## Generate masks before runing time-series"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Before running the time-series inversion, one may want to looks at average coherence in the stack. To create a map of average spatial coherence use \"temporal_average.py\":"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "!temporal_average.py ./inputs/ifgramStack.h5 -d coherence -o avgSpatialCoh.h5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "view.main(['avgSpatialCoh.h5'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Also one may optionally extract the water mask from the geometry file to be used later for masking the time-series results:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "!generate_mask.py inputs/geometryGeo.h5 waterMask --nonzero -o waterMask.h5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "view.main(['waterMask.h5'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "!generate_mask.py  inputs/ifgramStack.h5  --nonzero  -o maskConnComp.h5  --update\n",
    "view.main(['maskConnComp.h5'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "## reference_point"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "The interferometric phase is relative observation by nature. The phases of each unwrapped interferogram are relative with respect to an arbitrary pixel. Therfore we need to reference all interferograms to a common reference pixel.\n",
    "The step \"reference_point\" selects a common reference pixel for the stack of interferograms. The default approach of mintpy is to choose a pixel with highest spatial coherence in the stack. Other options include specifying the longitude and latitude of the desired reference pixel or the line and column number of the refence pixel.    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "!smallbaselineApp.py SanFranSenDT42.txt --dostep reference_point"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Running the \"reference_step\" adds additional attributes \"REF_X, REF_Y\" and \"REF_LON, REF_LAT\" to the ifgramStack.h5 file. To see the attributes of the file run info.py:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "!info.py inputs/ifgramStack.h5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Note that reference_point does not change the actual values of the unwrapped phase dataset. However, MintPy takes into account the phase at the reference point while performing the time-series inversion. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "## invert_network"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "In the next step we invert the network of differential unwrapped interferograms to estimate the time-series of unwrapped phase with respect to a reference acquisition date, which by default is the first acquisition. The estimated time-series is converted to distance change from radar to target and is provided in meters.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "!smallbaselineApp.py SanFranSenDT42.txt --dostep invert_network"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "The main product generated after running inver_network step, is timeseries.h5. To see the general content of the file run info.py"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "!info.py timeseries.h5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "The timeseries file contains three datasets, the time-series which is the interferometric range change for each acquisition relative to the reference acquisition, the \"date\" dataset which contains the acquisition date for each acquisition and the bperp dataset which contains the timeseries of the perpendicular baseline.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "# for complex options\n",
    "scp_args = 'timeseries.h5 -v -5 5 --noaxis '\n",
    "view.main(scp_args.split())\n",
    "# for all-default or no options in one-line : view.main(['timeseries.h5'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "<div class=\"alert alert-info\">\n",
    "<b>Question :</b> \n",
    "Why we invert a network of small baseline interferograms to estimate a single master time-series, instead of forming the single-master network of interferograms at first place? \n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "## velocity"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "The ground deformation caused by many geophysical or anthropogenic processes are linear at first order approximation. Therefore it is common to estimate the rate of the ground deformation which is the slope of linear fit to the time-series. The step \"velocity\" estimates the rate of the displacement. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "!smallbaselineApp.py SanFranSenDT42.txt --dostep velocity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "scp_args = 'velocity.h5 velocity -v -1 1'\n",
    "view.main(scp_args.split())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "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": {
    "hidden": true
   },
   "source": [
    "Obvious features in the estimated velocity map:\n",
    "\n",
    "1) The general pattern of displacement is consistent with tectonic setting of the region. Pacific plate is moving north-west with respect to North american plate. The satellites is on a descending track and the estimated displacement shows the blue region is moving away from the red region. \n",
    "\n",
    "2) The magnitude of the relative movement (blue region relative to red region) is about ~15 mm/yr in radar Line-Of-Sight (LOS) direction which is consistent with ~40 mm/yr horizontal movements of pacific relative to north america.\n",
    "\n",
    "3) The estimated velocity shows a linear feature almost aligned with the south-east north-west diagonal of the map. This linear feature shows the aseismic fault creep on Hayward fault. \n",
    "\n",
    "4) Further north and around latitude 38N another linear feature represents Concord fault parallel to Hayward fault.\n",
    "\n",
    "5) Around latitude 37.9 N, a east-west linear discontinuity is evident. This is most likely caused by missing bursts in some interferograms. \n",
    "\n",
    "6) The red region at the south-east corenr (around 37.35 N, 121.9W) seems to be a hydrological signal showing ground uplift caused by acquifer recharge.\n",
    "\n",
    "7) The block box at 37.7N, 122.3W is the reference pixel for this map. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "<div class=\"alert alert-info\">\n",
    "<b>Note:</b> \n",
    "For an interactive plot of tinme-series for each pixel, use the tsview.ipynb jupyter notebook. From a terminal one may directly call the tsview.py to display InSAR time-series. \n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "The estimated velocity also comes with an expression of unecrtainty which is simply based on the goodness of fit while fitting a linear model to the time-series. This quantity is saved in \"velocity.h5\" under the velocityStd dataset. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "scp_args = 'velocity.h5 velocityStd -v 0 0.2'\n",
    "view.main(scp_args.split())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "The estimated standard deviation only represents the goodness of fit and can be biased or maybe under-estimating the actual uncertainty of the product. However, the spatial pattern of the estimated standard deviation is interesting and clearly shows the spatial correltion of noise in the time-series. The uncertainty is distance dependent and increases with increasing distance between pixels. This map shows the uncertainty for each pixle relative to the reference pixel. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "<div class=\"alert alert-info\">\n",
    "<b>Question :</b> \n",
    "What are the sources of errors that can potentially increase the uncertainty or bias the estimated velocity at this stage?  \n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "# 3. Error analysis (what is signal, what is noise!)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Uncertainty of the ground displacement products derived from InSAR time-series, depends on the quality of the inversion of the stack of interferograms and the accuracy in separating the ground displacement from other components of the InSAR data. Therefore the definition of signal vs noise is different at the two main steps in mintpy:  \n",
    "\n",
    "1) During the inversion: \n",
    "    At this step all systematic components of the interferometric phase (e.g., ground displacement, propagation delay, geometrical residuals caused by DEM or platform's orbit inaccuracy) are considered signal, while the interferometric phase decorrelation, phase unwrapping error and phase inconsistency are considered noise. \n",
    "    \n",
    "2) After inversion: the ground displacement component of the time-serieses is signal, and everything else (including the propagation delay and geometrical residuals) are considered noise\n",
    "\n",
    "Therefore we first discuss the possible sources of error during the inversion and the existing ways in MintPy to evaluate the quality of inversion and to improve the uncertainty of the inversion. Afterwards we explain the different components of the time-series and the different processing steps in MintPy to separate them from ground displacement signal.  \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "## 3.1 Quality of the inversion"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "The main sources of noise during the time-series inversion includes decorrelation, phase unwrapping error and the inconsistency of triplets of interferofgrams. Here we mainly focus on the decorrelation and unwrapping errors. We first show the existing quantities in MintPy to evaluate decorrelation and unwrapping errors and then discuss the existing ways in MintPy to reduce the decorrelation and unwrapping errors on the time-series inversion."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "### 3.1.1 Average spatial coherence"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Mintpy computes temporal average of spatial coherence of the entire stack as a potential ancillary measure to choose reliable pixels after time-series inversion. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "view.main(['avgSpatialCoh.h5'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "### 3.1.2 Temporal coherence"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "In addition to timeseries.h5 which contains the time-series dataset, invert_network produces other quantities, which contain metrics to evaluate the quality of the inversion including temporalCoherence.h5. Temporal coherence represents the consistency of the timeseries with the network of interferograms. \n",
    "\n",
    "Temporal coherence varies from 0 to 1. Pixels with values closer to 1 are considered reliable and pixels with values closer to zero are considered unreliable. For a dense network of interferograms, a threshold of 0.7 may be used (Yunjun et al, 2019)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "view.main(['temporalCoherence.h5'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "<div class=\"alert alert-info\">\n",
    "<b>Question :</b> \n",
    "For a network of interferograms with only one interferogram between each acquisition and subsequent acquisition (sequential network of interferograms), how the temporal coherence looks like? Can the temporal coherence be trusted in this case?\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "<div class=\"alert alert-info\">\n",
    "<b>Question :</b> \n",
    "Generally, the temporal coherence may resamble the spatial pattern of the average spatial coherence. Under what condition, an area with very high average spatial coherence may show very low temporal coherence?  \n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "### 3.1.3 Modify the Network of interferograms"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "This step enables modifying the network of interferograms before the network inversion. Several options exist to modify a network. The network modification can be based on choosing interferogram based on an average coherence value over the entire interferogram (optionally masked area), or over an area of interest. The network can be also modified based on temporal and perpendicular baseline thresholds, maximum number of connections between each acquisition and following acquisitions, or simply by excluding dates or interferogram pairs or their indices. All different network modification options can be configured through the configuration file:\n",
    "\n",
    "```cfg\n",
    "## 1) Coherence-based network modification = Threshold + MST, by default\n",
    "## It calculates a average coherence value for each interferogram using spatial coherence and input mask (with AOI)\n",
    "## Then it finds a minimum spanning tree (MST) network with inverse of average coherence as weight (keepMinSpanTree)\n",
    "## For all interferograms except for MST's, exclude those with average coherence < minCoherence.\n",
    "mintpy.network.coherenceBased  = auto  #[yes / no], auto for no, exclude interferograms with coherence < minCoherence\n",
    "mintpy.network.keepMinSpanTree = auto  #[yes / no], auto for yes, keep interferograms in Min Span Tree network\n",
    "mintpy.network.minCoherence    = auto  #[0.0-1.0], auto for 0.7\n",
    "mintpy.network.maskFile        = auto  #[file name, no], auto for waterMask.h5 or no (use all pixels)\n",
    "mintpy.network.aoiYX           = auto  #[y0:y1,x0:x1 / no], auto for no, area of interest for coherence calculation\n",
    "mintpy.network.aoiLALO         = auto  #[lat0:lat1,lon0:lon1 / no], auto for no - use the whole area\n",
    "## 2) Network modification based on temporal/perpendicular baselines, date etc.\n",
    "mintpy.network.tempBaseMax     = auto  #[1-inf, no], auto for no, maximum temporal baseline in days\n",
    "mintpy.network.perpBaseMax     = auto  #[1-inf, no], auto for no, maximum perpendicular spatial baseline in meter\n",
    "mintpy.network.connNumMax      = auto  #[1-inf, no], auto for no, maximum number of neighbors for each acquisition\n",
    "mintpy.network.referenceFile   = auto  #[date12_list.txt / ifgramStack.h5 / no], auto for no\n",
    "mintpy.network.excludeDate     = 20160107  #[20080520,20090817 / no], auto for no\n",
    "mintpy.network.excludeIfgIndex = auto  #[1:5,25 / no], auto for no, list of ifg index (start from 0)\n",
    "mintpy.network.startDate       = auto  #[20090101 / no], auto for no\n",
    "mintpy.network.endDate         = auto  #[20110101 / no], auto for no\n",
    "```\n",
    "\n",
    "<div class=\"alert alert-info\">\n",
    "<b>TIP:</b> \n",
    "For coherence-based network modification, a customized area of interest (AOI) that includes the low coherence pixels surrounding is recommended. \n",
    "</div>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "<div class=\"alert alert-info\">\n",
    "<b>TIP:</b> \n",
    "For best results, we recommend (if possible) users process a redundant network of interferograms to have more freedom to modify the network before the inversion. \n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "By looking at the temporal coherence we suspect that some of the interferograms may have missing frames or may have wrong stitching. Since we have many interefrograms, we have the flexibility to remove few if possible. Let's first plot all interferograms and visually investigate any problem with some pairs. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "scp_args = 'inputs/ifgramStack.h5 -v -10 10'\n",
    "view.main(scp_args.split())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "By investigating all the interferograms it's obvious that the interferograms with the following indices have clear jumps at a burst boundary most likely caused by missing bursts in some acquisitions or due to problems in stitching:\n",
    "\n",
    "25, 30, 34, 40, 46, 52, 56, 57, 61, 62, 67, 70, 74, 80\n",
    "\n",
    "Therfore we decide to remove these interferograms from the stack and repeat the inversion. To remove the interferograms with specific indices we need to exclude the following option in the configuration file"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "\n",
    "\n",
    "```cfg\n",
    "mintpy.network.excludeIfgIndex = 25, 30, 34, 40, 46, 52, 56, 57, 61, 62, 67, 70, 74, 80\n",
    "mintpy.network.keepMinSpanTree = yes\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "<div class=\"alert alert-info\">\n",
    "<b>NOTE:</b> \n",
    "By modifying the network, the interferograms are not physically removed. However the value for that interferogram in the dropIfgram dataset in ifgramStack.h5 turns to False.  \n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "<div class=\"alert alert-info\">\n",
    "<b>NOTE:</b> \n",
    "By forcing the mintpy.network.keepMinSpanTree = yes option, the connectivity of the network is always gauranteed.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "!smallbaselineApp.py SanFranSenDT42.txt --dostep modify_network"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Plotting the network after network modification shows how the network looks like"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "plot_network.main(['inputs/ifgramStack.h5'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "After modifying the network, we need to manually remove the timeseries and re-run the inversion to see the impact of the network modification"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "!rm timeseries.h5\n",
    "!smallbaselineApp.py SanFranSenDT42.txt --dostep invert_network"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Let's check how the temporal coherence looks like after removing those interferograms. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "view.main(['temporalCoherence.h5'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Note that the temporal coherence at the top has increased and the discontinuity has disapeared."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "<div class=\"alert alert-info\">\n",
    "<b>Question:</b> \n",
    "If someone is interested only in the timeseries at one of the islands arround latitude 37.9N, where the reference pixel should be put in?</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "### 3.1.4 Phase unwrapping error correction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "The interferometric phases are wrapped (modulo 2pi) and integration of the wrapped phase,\n",
    "commonly called phase unwrapping, is required to obtain a field of relative phase with\n",
    "respect to a given pixel. The phase unwrapping algorithms add integer number of 2PI phase jumps to recover the unwrapped phase. Interferometric phase noise and discontinuities among different coherent regions may lead to wrong 2PI jumps added to the phase field known as unwrapping error. Unwrapping errors can bias the estimated time-series. \n",
    "\n",
    "MintPy provides three methods to possibly detect and correct the phase unwrapping errors. The first method is based on the phase closure of the triplets of the interferograms. The second methods is automating the traditional manual bridging method in which coherent components with the smallest distance from each other are assumed connected and therefore the a smooth phase variation across them are enforced. The third approach is a hybrid approach and simply uses the both approached mentioned before.\n",
    "\n",
    "Note that to use the phase closure approach a dense network of interferograms should be available.\n",
    "\n",
    "To use the phase unwrapping error correction methods usually a common mask is generated, which shows pixels with valid unwrapped phase in all interferograms. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "The phase unwrapping error correction can be configured using the following options:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "```cfg\n",
    "## reference: Yunjun et al., 2019\n",
    "## supported methods:\n",
    "## a. phase_closure           - suitable for highly redundant network\n",
    "## b. bridging                - suitable for islands or areas with steep topography\n",
    "## c. bridging+phase_closure  - recommended\n",
    "## Options for the bridging method:\n",
    "## 1. ramp - a phase ramp could be removed before estimating the phase difference between reliable regions\n",
    "##    and added back afterwards.\n",
    "## 2. bridgePtsRadius - half size of the window used to calculate the median value of phase difference\n",
    "mintpy.unwrapError.method          = auto  #[bridging / phase_closure / bridging+phase_closure / no], auto for no\n",
    "mintpy.unwrapError.waterMaskFile   = auto  #[waterMask.h5 / no], auto for waterMask.h5 or no [if no waterMask.h5 found]\n",
    "mintpy.unwrapError.ramp            = auto  #[linear / quadratic], auto for no; recommend linear for L-band data\n",
    "mintpy.unwrapError.bridgePtsRadius = auto  #[1-inf], auto for 50, half size of the window around end points\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "### 3.1.5 Weighted least squares"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "The following **weighted least squares (WLS) inversions** methods are supported:\n",
    "\n",
    "+ Inverse of covariance\n",
    "+ Finisher Information Matrix (FIM)\n",
    "+ Spatial coherence\n",
    "+ Uniform / no weighting\n",
    "\n",
    "The corresponding template options are: \n",
    "\n",
    "```cfg\n",
    "## Invert network of interferograms into time-series using weighted least sqaure (WLS) estimator.\n",
    "## weighting options for least square inversion [fast option available but not best]:\n",
    "## a. var - use inverse of covariance as weight (Tough et al., 1995; Guarnieri & Tebaldini, 2008) [recommended]\n",
    "## b. fim - use Fisher Information Matrix as weight (Seymour & Cumming, 1994; Samiei-Esfahany et al., 2016).\n",
    "## c. coh - use coherence as weight (Perissin & Wang, 2012)\n",
    "## d. no  - uniform weight (Berardino et al., 2002)\n",
    "\n",
    "mintpy.networkInversion.weightFunc      = auto #[var / fim / coh / no], auto for var\n",
    "\n",
    "```\n",
    "\n",
    "By default MintPy uses the same network of interferograms for all the pixels. However, it is possible to use variable networks for different pixels by specifying a mask dataset, a threshold value for the mask and a minimum redundance value: \n",
    "  \n",
    "```cfg\n",
    "\n",
    "## mask options for unwrapPhase of each interferogram before inversion (recommed if weightFunct=no):\n",
    "## a. coherence        - mask out pixels with spatial coherence < maskThreshold\n",
    "## b. connectComponent - mask out pixels with False/0 value\n",
    "## c. no               - no masking [recommended].\n",
    "## SBAS (Berardino et al., 2002) = minNormVelocity (yes) + weightFunc (no)\n",
    "\n",
    "mintpy.networkInversion.maskDataset     = auto #[coherence / connectComponent / no], auto for no\n",
    "mintpy.networkInversion.maskThreshold   = auto #[0-1], auto for 0.4\n",
    "mintpy.networkInversion.minRedundancy   = auto #[1-inf], auto for 1.0, min num_ifgram for every SAR acquisition\n",
    "```\n",
    "\n",
    "<div class=\"alert alert-info\">\n",
    "<b>NOTE:</b> \n",
    "If you choose to use variable network for different pixels, then interpret the temporal coherence with caution as the temporal coherence for a network with only obe connection between each acquisition and next acquisition is 1. A by product named numInvIfgram.h5 shows the number of interferograms used in the inversion for each pixel.\n",
    "</div>\n",
    "\n",
    "MintPy by default masks the estimated time-series using waterMask.h5. However, the mask can be specified or can be turned off:\n",
    "\n",
    "```cfg\n",
    "mintpy.networkInversion.waterMaskFile   = auto #[filename / no], auto for waterMask.h5 or no [if no waterMask.h5 found]\n",
    "```\n",
    "\n",
    "The least squares inversion can be performed using phase (A matrix in Berardino et al, 2002) or using phase velocity (B matrix in Berardino at al, 2002). The latter allows to invert a disconnected network. For a connected network, the design matrix is full-rank and the inversion using either methods are the same.\n",
    "The minimization is primarily based on L2 norm minimization. The L1 norm minimization has not been extensively tested. \n",
    "\n",
    "```cfg\n",
    "\n",
    "mintpy.networkInversion.minNormVelocity = auto #[yes / no], auto for yes, min-norm deformation velocity or phase\n",
    "mintpy.networkInversion.residualNorm    = auto #[L2 ], auto for L2, norm minimization solution\n",
    "```\n",
    "\n",
    "MintPy offers some options to use dask for parallel computing on High Performance Computing\n",
    "\n",
    "```cfg\n",
    "## Parallel processing with Dask for HPC\n",
    "mintpy.networkInversion.parallel        = auto #[yes / no], auto for no, parallel processing using dask\n",
    "mintpy.networkInversion.numWorker       = auto #[int > 0], auto for 40, number of works for dask cluster to use\n",
    "mintpy.networkInversion.walltime        = auto #[HH:MM], auto for 00:40, walltime for dask workers\n",
    "## Temporal coherence is calculated and used to generate final mask (Pepe & Lanari, 2006, IEEE-TGRS)\n",
    "mintpy.networkInversion.minTempCoh      = auto #[0.0-1.0], auto for 0.7, min temporal coherence for mask\n",
    "mintpy.networkInversion.minNumPixel     = auto #[int > 0], auto for 100, min number of pixels in mask above\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "## 3.2 Noise reduction of InSAR displacement time-series"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "After inversion of the network of interferograms, the estimated time-series contains different components including tropospheric delay, topographic residuals ground displacement and other possible geophysical components (e.g., tide, ionosphere if they have not been corrected) or instrumental effects (e.g., the local oscilator drift of Envisat). Given ground displacement as our signal of interest, the following processing steps attempt to separate signal from noise and provide a ground displacement time-series for selected coherent pixels."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "### 3.2.1 Correct local oscillator drift (for Envisat) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Data from Envisat’s Advanced Synthetic Aperture Radar (ASAR) instrument include a phase ramp in range direction due to timing errors. This step corrects the local oscillator drift using the empirical model given by Marinkovic and Larsen (2013). It's automatically turned ON for Envisat data and OFF for all the other dataset."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "### 3.2.2 Tropospheric delay correction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "This step corrects the tropospheric phase delay. Two methods are supported:\n",
    "\n",
    "+ Global Atmospheric Models (GAMs) (Jolivet et al., 2011; 2014; PyAPS needs to be installed).\n",
    "+ Empirical relationship between stratified tropospheric delay and topography (Doin et al., 2009).\n",
    "\n",
    "The corresponding template options are:\n",
    "\n",
    "```cfg\n",
    "## correct tropospheric delay using the following methods:\n",
    "## a. height_correlation - correct stratified tropospheric delay (Doin et al., 2009, J Applied Geop)\n",
    "## b. pyaps - use Global Atmospheric Models (GAMs) data (Jolivet et al., 2011; 2014)\n",
    "##      ERA5  - ERA-5 from ECMWF [default; need to install pyaps3 on GitHub]\n",
    "##      ECMWF - ERA-Interim from ECMWF [need to install pyaps on Caltech/EarthDef]\n",
    "##      MERRA - MERRA-2 from NASA Goddard [need to install pyaps on Caltech/EarthDef]\n",
    "##      NARR  - NARR from NOAA [recommended for areas in North America; need to install pyaps on Caltech/EarthDef]\n",
    "mintpy.troposphericDelay.method = auto  #[pyaps / height_correlation / no], auto for pyaps\n",
    "\n",
    "## Notes for pyaps: \n",
    "## a. GAM data latency: with the most recent SAR data, there will be GAM data missing, the correction\n",
    "## will be applied to dates with GAM data available and skipped for the others.\n",
    "## b. WEATHER_DIR: if you define an environmental variable named WEATHER_DIR to contain the path to a \n",
    "## directory, then MintPy applications will download the GAM files into the indicated directory. Also MintPy\n",
    "## application will look for the GAM files in the directory before downloading a new one to prevent downloading\n",
    "## multiple copies if you work with different dataset that cover the same date/time.\n",
    "mintpy.troposphericDelay.weatherModel = auto  #[ERA5 / ECMWF / MERRA / NARR], auto for ERA5, for pyaps method\n",
    "mintpy.troposphericDelay.weatherDir   = auto  #[path2directory], auto for WEATHER_DIR or \"./\"\n",
    "\n",
    "## Notes for height_correlation:\n",
    "## Extra multilooking is applied to estimate the empirical phase/elevation ratio ONLY.\n",
    "## For an dataset with 5by15 looks, looks=8 will generate phase with (5*8)by(15*8) looks\n",
    "## to estimate the empirical parameter; then apply the correction to original phase (with 5by15 looks),\n",
    "## if the phase/elevation correlation is larger than minCorrelation.\n",
    "mintpy.troposphericDelay.polyOrder      = auto  #[1 / 2 / 3], auto for 1, for height_correlation method\n",
    "mintpy.troposphericDelay.looks          = auto  #[1-inf], auto for 8, for height_correlation, extra multilooking number\n",
    "mintpy.troposphericDelay.minCorrelation = auto  #[0.0-1.0], auto for 0, for height_correlation\n",
    "```\n",
    "\n",
    "It outputs:\n",
    "+ timeseries_ECMWF.h5: the troposphere-corrected time-series (ECMWF data is selected in this case).\n",
    "+ ECMWF.h5: the estimated tropospheric phase delay time-series (for pyaps method only)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "### 3.2.3 deramping (optional)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "This step estimate and removes a linear or quadratic ramp from each acquisition based on the phase of the reliable pixels. It's recommended for localized deformation signals, such as volcanic deformation, landslides and city subsidence; but not recommeded for long spatial wavelength deformation signals, such as interseismic deformation.\n",
    "\n",
    "The cooresponding template options are:\n",
    "\n",
    "```cfg\n",
    "## estimate and remove a phase ramp for each acquisition based on the reliable pixels.\n",
    "## recommended for localized deformation signals, i.e. volcanic deformation, landslide and subsidence, etc.\n",
    "mintpy.deramp          = auto  #[no / linear / quadratic], auto for no - no ramp will be removed\n",
    "mintpy.deramp.maskFile = auto  #[filename / no], auto for maskTempCoh.h5, mask file for ramp estimation\n",
    "```\n",
    "\n",
    "It outputs a new time-series HDF5 file with suffix _ramp_: timeseries_ECMWF_ramp.h5 in this example."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "### 3.2.4 Topographic residual correction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "This step corrects the phase residual caused by the inaccuracy of DEM (DEM error) using its relationship with the perpendicular baseline time-series (Fattahi and Amelung, 2013, IEEE-TGRS). The corresponding template options are:\n",
    "\n",
    "```cfg\n",
    "## reference: Fattahi and Amelung, 2013, IEEE-TGRS\n",
    "## Notes on options:\n",
    "## stepFuncDate      - Specify stepFuncDate option if you know there are sudden displacement jump in your area,\n",
    "##    i.e. volcanic eruption, or earthquake, and check timeseriesStepModel.h5 afterward for their estimation.\n",
    "## excludeDate       - Dates excluded for error estimation only\n",
    "## pixelwiseGeometry - Use pixel-wise geometry info, such as incidence angle and slant range distance for error estimation\n",
    "##    yes - use pixel-wise geometry when they are available [slow; used by default]\n",
    "##    no  - use mean geometry [fast]\n",
    "mintpy.topographicResidual                    = auto  #[yes / no], auto for yes\n",
    "mintpy.topographicResidual.polyOrder          = auto  #[1-inf], auto for 2, poly order of temporal deformation model\n",
    "mintpy.topographicResidual.phaseVelocity      = auto  #[yes / no], auto for no - phase, use phase velocity for error estimation\n",
    "mintpy.topographicResidual.stepFuncDate       = auto  #[20080529,20100611 / no], auto for no, date of step jump\n",
    "mintpy.topographicResidual.excludeDate        = auto  #[20070321 / txtFile / no], auto for exclude_date.txt\n",
    "mintpy.topographicResidual.pixelwiseGeometry  = auto  #[yes / no], auto for yes, use pixel-wise geometry info\n",
    "```\n",
    "\n",
    "It outputs:\n",
    "+ timeseries_ECMWF_ramp_demErr.h5: the topographic residual corrected time-series.\n",
    "+ timeseriesResidual.h5: residual phase time-series of this least square inversion.\n",
    "+ demErr.h5: the estimated DEM error."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "!smallbaselineApp.py smallbaselineApp.cfg --dostep correct_topography"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "To view the estimated DEM error:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "view.main(['demErr.h5'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "### 3.2.5 Residual RMS for noise evaluation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "This step calculates the Root Mean Square (RMS) of the residual phase time-series for each acquisition; then it:\n",
    "1. selects the date with the minimum RMS value as the optimal reference date.\n",
    "2. detects the noisy acquisitions with RMS beyond the outlier detection threshold.\n",
    "\n",
    "The corresponding template options are:\n",
    "```cfg\n",
    "## 1) Residual Phase Root Mean Square\n",
    "## calculate the Root Mean Square (RMS) of residual phase time-series for each acquisition\n",
    "## To get rid of long wavelength component in space, a ramp is removed for each acquisition\n",
    "## Set optimal reference date to date with min RMS\n",
    "## Set exclude dates (outliers) to dates with RMS > cutoff * median RMS (Median Absolute Deviation)\n",
    "mintpy.residualRMS.maskFile = auto  #[file name / no], auto for maskTempCoh.h5, mask for ramp estimation\n",
    "mintpy.residualRMS.deramp   = auto  #[quadratic / linear / no], auto for quadratic\n",
    "mintpy.residualRMS.cutoff   = auto  #[0.0-inf], auto for 3\n",
    "```\n",
    "\n",
    "It outputs:\n",
    "+ rms_timeseriesResidual_ramp.txt: for RMS value of each acquisition\n",
    "+ rms_timeseriesResidual_ramp.pdf: plot of the rms_timeseriesResidual_ramp.txt\n",
    "+ reference_date.txt: date in YYYYMMDD format for the optional reference date\n",
    "+ exclude_date.txt: date(s) in YYYYMMDD format for the noisy acquisitions (if at least one is detected)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "!smallbaselineApp.py smallbaselineApp.cfg --dostep residual_RMS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "cat rms_timeseriesResidual_ramp.txt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "cat reference_date.txt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "### 3.2.6 Changing the reference date"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "This step changes the reference date of all phase time-series files, based on the input template option:\n",
    "\n",
    "```cfg\n",
    "## reference all time-series to one date in time\n",
    "## no     - do not change the default reference date (1st date)\n",
    "mintpy.reference.date  = auto   #[reference_date.txt / 20090214 / no], auto for reference_date.txt\n",
    "```\n",
    "\n",
    "This step operates on the existing time-series files and does not output new files.\n",
    "\n",
    "<div class=\"alert alert-warning\">\n",
    "<b>Note:</b> \n",
    "The optimal reference date (default option) gives the time-series plot a \"clean\" looks only. Changing the reference is equivalent to adding a constant to the displacement time series, which does not change the velocity or any other information derived from the displacement time series.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "!smallbaselineApp.py smallbaselineApp.cfg --dostep reference_date"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "### 3.2.7 Re-estimating Velocity after noise reduction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "This step estimates the slope of the best fitting line to the displacement time-series and its standard deviation. Noisy acquisitions (identified in \"residual_RMS\" step) from exclude_date.txt file are excluded by default during the estimation. \n",
    "\n",
    "The corresponding template options are:\n",
    "\n",
    "```cfg\n",
    "## estimate linear velocity from time-series, and from tropospheric delay file if exists.\n",
    "mintpy.velocity.excludeDate = auto   #[exclude_date.txt / 20080520,20090817 / no], auto for exclude_date.txt\n",
    "mintpy.velocity.startDate   = auto   #[20070101 / no], auto for no\n",
    "mintpy.velocity.endDate     = auto   #[20101230 / no], auto for no\n",
    "```\n",
    "\n",
    "It outputs:\n",
    "+ velocity.h5: the estimated average velocity and its standard deviation from displacement time-series.\n",
    "+ velocityECMWF.h5: same as above but from the tropospheric delay time-series, to see the potential bias introduced by troposphere if it was not corrected."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "!smallbaselineApp.py smallbaselineApp.cfg --dostep velocity\n",
    "!smallbaselineApp.py smallbaselineApp.cfg --dostep google_earth"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "scp_args = 'velocity.h5 velocity -v -1 1'\n",
    "view.main(scp_args.split())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "scp_args = 'velocity.h5 --start-lalo 37.5307 -122.3631 --end-lalo 37.9828 -121.7707 '\n",
    "scp_args = 'velocity.h5 --start-lalo 37.6523 -122.1163 --end-lalo 37.7047 -122.0497 '\n",
    "plot_transection.main(scp_args.split())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "The transect across Hayward fault shows ~2.5 mm/yr fault creep in LOS direction."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "# 4. Validation (comparing InSAR with GPS)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "MintPy's analysis is independent of GPS observations. This allows validating InSAR products with GPS data when they are available. For this purpose MintPy automatically downloads GPS data over the region of interest from Nevada Geodetic Laboratory at University of Nevada, Reno. To download the GPS data, project the components to InSAR LOS direction and plot them on the InSAR velocity map one can use the following command:    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "scp_args = 'velocity.h5 velocity --show-gps --ref-gps P225 --gps-comp enu2los'\n",
    "view.main(scp_args.split())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "In order to display the GPS station names on the plot add --gps-label to the plot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "scp_args = 'velocity.h5 velocity --show-gps --ref-gps P225 --gps-comp enu2los --gps-label --fontsize 8'\n",
    "view.main(scp_args.split())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hideCode": true,
    "hideOutput": false
   },
   "source": [
    "# 3.3 Other MintPy functionalities\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "### 3.3.1 Stand alone scripts for every processing step"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Almost all processing steps are callable through running individual command line scripts. As an example one can invert the interferograms using `ifgram_inversion.py  inputs/ifgramStack.h5 -t smallbaselineApp.cfg`.\n",
    "\n",
    "Or perform the atmospheric delay correction using `tropo_pyaps3.py` or `tropo_phase_elevation.py`.\n",
    "\n",
    "Or do the topographic residual correction using `dem_error.py`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "### 3.3.2 Projecting line-of-sight to vertical and horizontal components"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "If geocoded displacement products from ascending and descending paths are available, `asc_desc2horz_vert.py` script can be called to project the LOS components to Horizontal and vertical components"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "### 2.3.3 Generate a mask"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "generate_mask.py can be used to generate a mask based on existing datasets with differnt options to choose a threshold. For example the following command creates a boolean mask which is True for pixels with temporal coherence greater or equal 0.7 and False everywhere else."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "generate_mask.py  temporalCoherence.h5 -m 0.7 -o maskTempCoh.h5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "### 2.3.4 Geocoding products which are in radar coordinates"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "If the stack of interferograms are in radar coordinate system, the time-series inversion and produced products will be generated in radar coordinate system. To geocode the products, one may use geocode.py.\n",
    "\n",
    "Note that the same command with an additional \"--geo2radar\" option can be used to transform a map in geo-coordinate to the radar coordinate.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "### 2.3.5 Mosaiking multiple geocoded files"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "match.py Merges two or more geocoded datasets sharing common area into one. It finds the overlap area and calculates the average offset between the datasets. For example:\n",
    "\n",
    "match.py  vel_AlosAT422.h5  vel_AlosAT423.h5  vel_AlosAT424.h5  vel_AlosAT425.h5\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "<div class=\"alert alert-info\">\n",
    "<b>Question :</b> \n",
    "Let's assume one has generated Velocity fields for two adjacent tracks. After mosaiking would you expect the two adjacent track be always consistent at the overlap? What may cause inconsistency at the overalp of two adjacent tracks?\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "### 2.3.6 Plot transect"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "plots transect along a line specified with lat-lon coordinates. Check the link below for interactive plot:\n",
    "https://nbviewer.jupyter.org/github/insarlab/MintPy/blob/master/docs/tutorials/plot_transection.ipynb"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "scp_args = 'velocity.h5 --start-lalo 37.7629 -122.4929 --end-lalo 37.9504 -121.9296 '\n",
    "plot_transection.main(scp_args.split())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "### Plot coherence matrix (pixelwise)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Interactive notebook:\n",
    "https://nbviewer.jupyter.org/github/insarlab/MintPy/blob/master/docs/tutorials/plot_transection.ipynb"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "scp_args = 'inputs/ifgramStack.h5 --yx 858 512'\n",
    "plot_coherence_matrix.main(scp_args.split())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "### 2.3.7 Filters (spatial and temporal)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Standard time-series analysis with smallbaselineApp.py avoids any filtering. However, individual scripts are available for users to apply spatial and temporal filtering on the time-series products."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "### 2.3.8 Export to other formats"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true,
    "scrolled": true
   },
   "source": [
    "MintPy provides tools to export the products in HDF5 format to other formats including ...\n",
    "\n",
    "+ ROI_APC (save_roipac.py)\n",
    "+ GBIS (save_gbis.py)\n",
    "+ GMT (save_gmt.py)\n",
    "+ HDF-EOS5 (save_hdfeos5.py)\n",
    "+ [Google Earth raster (save_kmz.py)](https://mintpy.readthedocs.io/en/latest/google_earth/#2_raster_image)\n",
    "+ [Google Earth points for displacement time-series (save_kmz_timeseries.py)](https://mintpy.readthedocs.io/en/latest/google_earth/#1_displacement_time-series)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "### 2.3.9 view timeseries and other products"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "view.py : to view and plot the stack of interferogram, coherence, etc\n",
    "          velocity, temporal coherence, etc\n",
    "        DEM, products overlaid on DEM and many more"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "tsview.py : a time-series viewer"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "timeseries in kmz:\n",
    "\n",
    "save_kmz_timeseries.py timeseries_demErr.h5 --vel velocity.h5  --tcoh temporalCoherence.h5 --mask maskTempCoh.h5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Check out University of Miami online time-series viewer:\n",
    "    https://insarmaps.miami.edu/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "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.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}