licsbas_tutorial.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "mighty-uganda",
   "metadata": {},
   "source": [
    "# COMET InSAR Workshop 2024: LiCSBAS Tutorial\n",
    "\n",
    "The aim of this tutorial is to understand and complete the processing steps of open-source LiCSBAS software required to prepare a time series analysis for Campi Flegrei using LiCSAR data from 2016-2021.  \n",
    "\n",
    "\n",
    "This notebook is provided with default parameters using a subset of data to ensure that the practical session can be completed in the allotted time. The practical instructor will talk through the parameters and the operation of the notebook; however, you are not required to make edits for the notebook to run."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e6328ecb-95e6-4501-a26c-bcd731724175",
   "metadata": {},
   "source": [
    "### please run the following cell first to load some necessary modules:  \n",
    "(to run a cell, just click into it and press shift-Enter)  \n",
    "\n",
    "#### Note: in case your notebook gets frozen/kernel disconnected, try clicking Kernel->Restart (rather than Refreshing page). If this runs ok, you may just rerun the following cell, and then skip back to the cell you were in before."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cooperative-saudi",
   "metadata": {},
   "source": [
    "**cell \\[ 00 \\]** - this is a naming of a cell below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e2b51bc9-5d5d-439e-ab6e-b1d714f8281d",
   "metadata": {},
   "outputs": [],
   "source": [
    "from IPython.display import Image\n",
    "import os\n",
    "\n",
    "# below is for the local jupyhub env - should not affect the binder version\n",
    "os.environ['PATH'] = os.environ['PATH']+':/data/notebooks/repositories/LiCSBAS/bin'\n",
    "import sys\n",
    "sys.path.append('/data/notebooks/repositories/LiCSBAS/LiCSBAS_lib')\n",
    "\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "\n",
    "# some extra functions\n",
    "def display(png):\n",
    "    if not os.path.exists(png):\n",
    "        print('This png file does not exist, recheck please')\n",
    "    else:\n",
    "        return Image(png)\n",
    "\n",
    "def printbold(instr='done'):\n",
    "    print('\\033[1m' + instr + '\\033[0m')\n",
    "\n",
    "# selected frame\n",
    "frame = '022D_04826_121209' # for tutorial\n",
    "\n",
    "# Once you identified the frame ID to process and confirmed that \n",
    "# the sufficient interferograms are available on the frame, \n",
    "# make a processing directory and change to the directory\n",
    "framedir = os.path.join(os.environ['HOME'],frame)\n",
    "try:\n",
    "    os.mkdir(framedir)\n",
    "except:\n",
    "    pass\n",
    "\n",
    "os.chdir(framedir)\n",
    "\n",
    "# below check is for the local jupyhub env only, the files should have already been extracted if you run through binder\n",
    "if not os.path.exists('GEOCml2GACOS'):\n",
    "    print('extracting licsbas data')\n",
    "    os.system('7za x /data/notebooks/COMET_InSAR_Workshop_2024/Day2_LiCSBAS/tutorial.7z')\n",
    "    # or this one - a result of time series only:\n",
    "    os.system('7za x /data/notebooks/COMET_InSAR_Workshop_2024/Day2_LiCSBAS/tutorialts.7z')\n",
    "\n",
    "\n",
    "print('prerequisites loaded, thank you. Running in directory:')\n",
    "print(os.getcwd())\n",
    "\n",
    "printbold()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6ec6ebf3-5608-417a-9854-48834d89b0dc",
   "metadata": {},
   "source": [
    "# LiCSBAS Tutorial\n",
    "\n",
    "LiCSBAS is an open-source package in Python and bash to carry out InSAR time series analysis using LiCSAR interfeometric products (primarily unwrapped interferograms and coherence) which are freely available on the [COMET-LiCS web portal](https://comet.nerc.ac.uk/COMET-LiCS-portal/).\n",
    "\n",
    "\n",
    "\n",
    "Users can easily and automatically derive the time series and velocity of the displacement if sufficient LiCSAR products are available in the area of interest. LiCSBAS also contains visualization tools to interactively display the time series of displacement to help investigation and interpretation of the results. Current release offers various improvements including option to perform customized unwrapping procedure using *LiCSBAS02to05_unwrap.py*.  \n",
    "\n",
    "## Documentation\n",
    "\n",
    "The [COMET LiCSBAS repository](https://github.com/comet-licsar/licsbas) is an actively developed fork of the [original LiCSBAS](https://github.com/yumorishita/LiCSBAS) tool developed by Dr. Yu Morishita during his stay at the University of Leeds. For more detailed information, see the [**wiki**](https://github.com/yumorishita/LiCSBAS/wiki) pages and [quick start](https://github.com/yumorishita/LiCSBAS/wiki/2_0_workflow#quick-start).  \n",
    "\n",
    "This tutorial is run through Jupyter Notebook - see short [tutorial](https://www.dataquest.io/blog/jupyter-notebook-tutorial) if you are not familiar with this environment."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c8a7e4e8-230e-4ffc-8f6c-5dc4c0efd1b1",
   "metadata": {},
   "source": [
    "## Standard LiCSBAS Workflow\n",
    "\n",
    "<img src=\"https://raw.githubusercontent.com/wiki/yumorishita/LiCSBAS/images/workflow.png\" width=650>\n",
    "\n",
    "**Note:** Current workflow is extended by several optional processes, such as advanced selection of reference point (LiCSBAS120 by Dr. Qi Ou), pixel-based nullification of unwrapping errors (part of LiCSBAS12, by Dr. Lin Shen), advanced unwrapping procedure (LiCSBAS02to05, by Dr. Milan Lazecky), estimation of non-linear model and offsets (LiCSBAS_cum2vel, by Dr. Yu Morishita and M. Lazecky), ability to load ALOS interferograms (LiCSBAS01_get_geotiff_ALOS included from Dr. Yu Morishita's updated solutions), among various other changes. As we keep using standard LiCSBAS functionality within this tutorial, we did not update flowchart and other figures and leave them for investigation of individual parameters of the LiCSBAS scripts."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6e50647b-76bc-442e-99c3-ee33e12d11bb",
   "metadata": {},
   "source": [
    "# Module 0: Preparing the dataset\n",
    "\n",
    "This tutorial will run on a preprocessed dataset (i.e. modules 0 are complete here) that is clipped in both time and space, and cover covers volcanic area of Campi Flegrei in Napoli region (Italy). The figure below shows GPS network measurements over the Napoli region, published in https://www.mdpi.com/2072-4292/13/14/2725. The Campi Flegrei caldera (in the middle) has been inflating since 2010s, a continuous uplift can be observed from the GPS station time series (middle plot).  \n",
    "\n",
    "<img src=\"tutorial/campi_flegrei_gps.png\" width=\"1000\">\n",
    "\n",
    "You may use the following module 0 steps to process another frame of your interest. **For the purposes of this tutorial, they are marked as skipped**.  \n",
    "\n",
    "\n",
    "Note: you may identify LiCSAR frames covering your own area of interest (AOI) in your free time later. For this, you may just check the [COMET LiCSAR portal](https://comet.nerc.ac.uk/COMET-LiCS-portal/).  \n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6c25fa27-fdb8-49e1-badb-76f151d1858a",
   "metadata": {},
   "source": [
    "## (skipping) Step 0-1: Download GeoTIFF from COMET-LiCS web portal\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7c0c06c4-78d1-4c4e-a35d-21e4770c7f44",
   "metadata": {},
   "source": [
    "### LiCSBAS01_get_geotiff.py\n",
    "This script downloads GeoTIFF files of unw (unwrapped interferogram) and cc (coherence) in the specified frame ID from COMET-LiCS web portal.  \n",
    "The -f option is not necessary when the frame ID can be automatically identified from the name of the working directory.  \n",
    "GACOS data can also be downloaded if available, as well as phase data (and intensity) required for the lately added unwrapping procedure.  \n",
    "Existing GeoTIFF files are not re-downloaded to save time, i.e., only the newly available data will be downloaded.  \n",
    "\n",
    "**Note:** LiCSBAS is capable of running on any network of standard interferograms if these exist in GEOC/YYYYMMDD_YYYYMMDD together with coherence (either in 0-1 or 0-255 range)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "surface-inspector",
   "metadata": {},
   "source": [
    "**cell \\[ 01 \\]**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5d6d32a6-4a56-4627-8153-be06a5409390",
   "metadata": {},
   "outputs": [],
   "source": [
    "! LiCSBAS01_get_geotiff.py -h\n",
    "#e.g. LiCSBAS01_get_geotiff.py -f 167D_04884_131212 --get_gacos --n_para 1\n",
    "#printbold()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "45384849-2ea9-4232-b7c9-68e249f775f7",
   "metadata": {},
   "source": [
    "## (skipping) Step 0-2: Convert and downsample"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5750e7cf-3036-4832-b675-fac77838bd8a",
   "metadata": {},
   "source": [
    "### LiCSBAS02_ml_prep.py\n",
    "This script converts GeoTIFF files of unw and cc to float32 and uint8 format, respectively, for further time series analysis, and also downsamples (multilooks) data if specified.  \n",
    "Existing files are not re-created to save time, i.e., only the newly available data will be processed.  \n",
    "\n",
    "**Note:** Use of LiCSBAS02to05_unwrap.py replaces this and following 0-x optional steps, however they still can be applied to its output."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "modular-species",
   "metadata": {},
   "source": [
    "**cell \\[ 02 \\]**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b245132a-3ae9-40fe-9066-18846478f2ec",
   "metadata": {},
   "outputs": [],
   "source": [
    "! LiCSBAS02_ml_prep.py -h\n",
    "# e.g. LiCSBAS02_ml_prep.py -i GEOC -o GEOCml10 -n 10 --npara 1\n",
    "printbold()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dcaa8ffa-a34c-4362-9e83-e9c4953fc370",
   "metadata": {},
   "source": [
    "# (skipping) Step 0-3: GACOS correction (optional)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "03d2bd57-8475-4954-bff9-5c6c0d29cc77",
   "metadata": {},
   "source": [
    "### LiCSBAS03op_GACOS.py\n",
    "This script applies a tropospheric correction to unw data using GACOS data.  \n",
    "GACOS data may be automatically downloaded from COMET-LiCS web at step01 (if available), or could be externally obtained by requesting on a [GACOS web](http://www.gacos.net/).  \n",
    "(if you request the GACOS data through the GACOS web, the dates and time of interest can be found in baselines and slc.mli.par, respectively. Note that original GACOS data comes as zenith total delay, ZTD. This is automatically converted to the slant direction of the satellite, i.e. into slant total delay, SLTD).  \n",
    "\n",
    " \n",
    "The impact of the correction can be visually checked by showing GACOS_info.png and \\*/\\*.gacos.png:  \n",
    "- unw_org - original unwrapped data (STD = overall standard deviation of all unw values)\n",
    "- unw_cor - unwrapped data after GACOS correction\n",
    "- dsltd - differential (respective to the interferometric pair) correction image from GACOS ZTD data converted to SLTD (see above)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "532cf541-ae12-40a4-b0e6-2241a3301ec6",
   "metadata": {},
   "source": [
    "<img src=\"https://raw.githubusercontent.com/wiki/yumorishita/LiCSBAS/images/20170407_20170501.gacos.png\" width=600>  \n",
    "\n",
    "<img src=\"https://raw.githubusercontent.com/wiki/yumorishita/LiCSBAS/images/GACOS_info.png\" width=600>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "lyric-latest",
   "metadata": {},
   "source": [
    "**cell \\[ 03 \\]**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1b4f6b3d-bc42-41da-bd36-c34f0e51045f",
   "metadata": {},
   "outputs": [],
   "source": [
    "! LiCSBAS03op_GACOS.py -h\n",
    "# e.g. LiCSBAS03op_GACOS.py -i GEOCml10 -o GEOCml10GACOS --fillhole --n_para 1\n",
    "printbold()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "534865b5-903a-47ea-8cfa-98fb1d6fd579",
   "metadata": {
    "tags": []
   },
   "source": [
    "## (skipping) Step 0-4: Mask UNW (optional)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3d9ab168-e1cb-4406-9748-89aca9e2cb4e",
   "metadata": {},
   "source": [
    "### LiCSBAS04op_mask_unw.py\n",
    "This script masks specified areas or low coherence areas in the unw data.  \n",
    "The masking is effective when the unw data include areas which have many unwrapping errors and are not of interest, and can improve the result of Step 1-2 (loop closure).   \n",
    "\n",
    "<img src=\"https://raw.githubusercontent.com/wiki/yumorishita/LiCSBAS/images/20150524_20150617.unw_mask.png\" width=600>\n",
    "\n",
    "**Note:** In practice, an option to mask data based on their individual coherence (parameter *-s*) is valuable for reliable inversion in the later steps."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "partial-monkey",
   "metadata": {},
   "source": [
    "**cell \\[ 04 \\]**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f9527f75-4c3e-4bd8-972a-b13784cdccff",
   "metadata": {},
   "outputs": [],
   "source": [
    "! LiCSBAS04op_mask_unw.py -h\n",
    "printbold()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e3e5e279-9194-4fcc-a5e0-e197de3125da",
   "metadata": {},
   "source": [
    "## (skipping) Step 0-5: Clip UNW (optional)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "63d7da81-c0f9-42cf-a167-9b1b73b2f34e",
   "metadata": {},
   "source": [
    "### LiCSBAS05op_clip_unw.py\n",
    "This script clips a specified rectangular area of interest from unw and cc data.  \n",
    "The clipping can make the data size smaller and processing faster, and improve the quality measures calculated for the whole scene, such as avg loop phase closure result of Step 1-2.  \n",
    "\n",
    "<img src=\"https://raw.githubusercontent.com/wiki/yumorishita/LiCSBAS/images/20150524_20150617.unw_clip.png\" width=400>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "innovative-assistant",
   "metadata": {},
   "source": [
    "**cell \\[ 05 \\]**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0ca628c2-476f-45d9-b8ca-fe8293d363b5",
   "metadata": {},
   "outputs": [],
   "source": [
    "! LiCSBAS05op_clip_unw.py -h\n",
    "# e.g. LiCSBAS05op_clip_unw.py -i GEOCml10GACOS -o GEOCml10GACOSclip -g 14.03/14.22/40.78/40.90 --n_para 1\n",
    "printbold()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ed86fb69-9892-43fc-8db8-d962b5f954b3",
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true
    },
    "tags": []
   },
   "source": [
    "# Running steps 11-16 using existing preprocessed data\n",
    "\n",
    "LiCSBAS processing directory has following folders:\n",
    "- GEOC - location for the original data, (auto-)downloaded from COMET LiCSAR\n",
    "- GEOCmlX - result of Step 0-2 - X-times downsampled converted data\n",
    "- GEOCmlXGACOS - optional - result of GACOS correction (steps 0-4 and 0-5 also generate optional 'final pre-processed dataset'\n",
    "- GEOCmlX\\*clip - optional - spatially clipped dataset\n",
    "- TS_GEOC* - folder containing outputs of time series processing (see further)\n",
    "\n",
    "**Note:** We have preprocessed data for the given frame for you, for the purposes of this tutorial."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "psychological-conservative",
   "metadata": {},
   "source": [
    "**cell \\[ 06 \\]**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "international-irrigation",
   "metadata": {},
   "outputs": [],
   "source": [
    "# check the contents of the folder\n",
    "! pwd\n",
    "! ls\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "practical-example",
   "metadata": {},
   "outputs": [],
   "source": [
    "# check the number of interferograms\n",
    "! echo The folder contains `ls GEOCml2GACOSclip | wc -l` interferograms\n",
    "! echo \"last 30 folders+files:\"\n",
    "! ls GEOCml2GACOSclip | tail -n 30"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "conceptual-fancy",
   "metadata": {},
   "source": [
    "**cell \\[ 07 \\]**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "29650253-cf06-4060-80a2-3b085710915a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# check the interferograms network:\n",
    "! ls GEOCml2GACOSclip | grep '_20' > ifg_list.txt\n",
    "print('generating network plot')\n",
    "! LiCSBAS_plot_network.py -i ifg_list.txt -b GEOCml2GACOSclip/baselines -o network.png >/dev/null 2>/dev/null\n",
    "print('first 10 interferogram pairs:')\n",
    "! head -n 10 ifg_list.txt\n",
    "display('network.png')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fiscal-vienna",
   "metadata": {},
   "source": [
    "**cell \\[ 08 \\]**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f8aadb68-d18a-44b9-b55b-156bff7d2033",
   "metadata": {},
   "outputs": [],
   "source": [
    "# look into some random interferogram:\n",
    "pair = '20161014_20161026'\n",
    "png = os.path.join('GEOCml2GACOSclip',pair,pair+'.unw.png')\n",
    "display(png)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "worse-berry",
   "metadata": {},
   "source": [
    "**cell \\[ 09 \\]**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "10f7fa06-b613-4f0d-bfbd-d68db9dd32d5",
   "metadata": {},
   "outputs": [],
   "source": [
    "# see effect of GACOS of this selected interferogram:\n",
    "# note GACOS corrections are created before clipping, thus they are in the second extracted folder\n",
    "pair='20161008_20161026'\n",
    "#pair='20161002_20161020'  # pair causing higher std after GACOS\n",
    "print('unw_org - original unwrapped data')\n",
    "print('unw_cor - corrected unw data')\n",
    "print('dsltd - differential slant direction total delay estimate from GACOS data')\n",
    "print('colour scale: 6pi/colour cycle')\n",
    "png = os.path.join('GEOCml2GACOS',pair,pair+'.gacos.png')\n",
    "display(png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d46b8057-e1a5-445c-8d21-fc22efa01689",
   "metadata": {},
   "outputs": [],
   "source": [
    "# look into some interferogram - e.g. see effect of GACOS:\n",
    "! head -n 10 GEOCml2GACOSclip/GACOS_info.txt\n",
    "png = 'GEOCml2GACOSclip/GACOS_info.png'\n",
    "display(png)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c5a8c7ad-24b2-4318-a0a6-975dc2dccbdb",
   "metadata": {},
   "source": [
    "# Performing basic time series analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "07435d42-219b-4b67-9ba5-3853e418db32",
   "metadata": {},
   "source": [
    "# Step 1-1: Quality Check"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dff4ce6c-9110-4c7b-97f6-b799fc1e46bd",
   "metadata": {},
   "source": [
    "### LiCSBAS11_check_unw.py\n",
    "This script checks quality of unw data and identifies bad interferograms based on average coherence and coverage of the unw data.  \n",
    "This also prepares a time series working directory (overwrite if already exists).  \n",
    "\n",
    "Examples of low quality (to be removed) interferograms:  \n",
    "<img src=\"https://raw.githubusercontent.com/wiki/yumorishita/LiCSBAS/images/20180108_20180213.unw.geo.png\" width=\"300\">\n",
    "<img src=\"https://raw.githubusercontent.com/wiki/yumorishita/LiCSBAS/images/20150325_20150418.unw.geo.png\" width=\"300\">  \n",
    "\n",
    "**Note:** LiCSAR team works on updated quality checking together to flag and remove interferograms containing coregistration or other errors. An optional parameter **-s** would detect interferograms with azimuth ramp higher than approx. 0.1 mm/km, as often result due to strong ionospheric influence.  \n",
    "<img src=\"tutorial/licsbas11_s.png\" width=800>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "social-fields",
   "metadata": {},
   "source": [
    "**cells \\[ 10 \\]**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "psychological-confidence",
   "metadata": {},
   "outputs": [],
   "source": [
    "! LiCSBAS11_check_unw.py -h\n",
    "# e.g. LiCSBAS11_check_unw.py -d GEOCml2GACOSclip\n",
    "printbold()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "86fb4df1-bd54-4c7a-be7f-8e2be9e39724",
   "metadata": {},
   "outputs": [],
   "source": [
    "print('running the command')\n",
    "! LiCSBAS11_check_unw.py -d GEOCml2GACOSclip\n",
    "printbold()\n",
    "# ETA few seconds"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "olive-arrest",
   "metadata": {},
   "source": [
    "**cells \\[ 11 \\]**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "221ebdcc-c9f7-44f3-8ee5-bc34b5c3b0b2",
   "metadata": {},
   "outputs": [],
   "source": [
    "# check outputs of the step:\n",
    "! ls TS_GEOCml2GACOSclip\n",
    "# outputs of step 1-1 are in folders starting by 11"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6852412b-32f3-432c-91e5-b8c88c214208",
   "metadata": {},
   "outputs": [],
   "source": [
    "# see bad ifgs:\n",
    "! ls TS_GEOCml2GACOSclip/11bad_ifg_ras"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "964ef78d-f1c2-499c-802c-1e7214f75cf0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# check stats - perhaps change thresholds?\n",
    "! head -n 30 TS_GEOCml2GACOSclip/info/11ifg_stats.txt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "63610cd4-c8cf-473b-af82-4d5ee80beed0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# see one of the discarded ifgs:\n",
    "pair='20200414_20200420'\n",
    "png = os.path.join('TS_GEOCml2GACOSclip','11bad_ifg_ras',pair+'.unw.png')\n",
    "#png = os.path.join('GEOCml2GACOSclip',pair,pair+'.unw.png')\n",
    "display(png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "domestic-connection",
   "metadata": {},
   "outputs": [],
   "source": [
    "#check network\n",
    "png='TS_GEOCml2GACOSclip/network/network11.png'\n",
    "display(png)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "93964719-d9c1-4cb0-ae8c-13680582abb0",
   "metadata": {},
   "source": [
    "# Step 1-2: Loop Closure"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "seeing-therapy",
   "metadata": {},
   "source": [
    "### LiCSBAS12_loop_closure.py\n",
    "This script identifies bad unw by checking loop closures of interferometric triplets, i.e. non-zero residuals \n",
    "$$\\phi_{123}=\\phi_{12}+\\phi_{23}-\\phi_{13}$$\n",
    "\n",
    "\n",
    "<img src=\"https://raw.githubusercontent.com/wiki/yumorishita/LiCSBAS/images/20171203_20171215_20171227_loop.png\" width=\"400\">"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "20de2987-7da5-463c-9140-fea9ed0e3d08",
   "metadata": {},
   "source": [
    "This procedure runs 4 iterations where the following is applied, to investigate loop closure: \n",
    "1. iteration - identify loops and mark/not use unw ifgs of higher overall RMSE than given threshold, tolerating $2\\pi$ phase jumps\n",
    "2. iteration - check loop errors per pixel and select stable reference point (all valid unw data and a minimum of loop phase RMSE)\n",
    "3. iteration - use reference area, do not tolerate $2\\pi$ jumps, re-identify unw ifgs of higher overall RMSE than given threshold and remove them as bad\n",
    "4. iteration - calculate statistics (n_loop_err, n_ifg_noloop, and their ratio) and optionally nullify pixels of unws that only appear in non-closed loop -> good pixels of bad ifgs are still used  \n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cathedral-strand",
   "metadata": {},
   "source": [
    "**cells \\[ 12 \\]**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "84ea1a53-3d69-4063-bc1e-6ff7d1662e6b",
   "metadata": {},
   "outputs": [],
   "source": [
    "! LiCSBAS12_loop_closure.py -h\n",
    "# e.g. LiCSBAS12_loop_closure.py -d GEOCml10GACOS -l 3 --n_para 1\n",
    "printbold()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "mechanical-wagner",
   "metadata": {},
   "outputs": [],
   "source": [
    "print('running the command')\n",
    "! LiCSBAS12_loop_closure.py -d GEOCml2GACOSclip -l 1.5 --n_para 1 --multi_prime\n",
    "printbold()\n",
    "# ETA 7 min., rerunning: ETA few seconds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "682a1127-1e42-4548-8a68-5856e59278d3",
   "metadata": {},
   "outputs": [],
   "source": [
    "import LiCSBAS_io_lib as io_lib\n",
    "import LiCSBAS_loop_lib as loop_lib\n",
    "import LiCSBAS_tools_lib as tools_lib\n",
    "import LiCSBAS_inv_lib as inv_lib\n",
    "import LiCSBAS_plot_lib as plot_lib\n",
    "import xarray as xr\n",
    "import SCM\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a78b94fa-fee4-4a76-8e5f-5c51f1eb41d4",
   "metadata": {},
   "outputs": [],
   "source": [
    "tools_lib.get_cmap('SCM.romaO')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "38008f8a-2d15-497d-816a-1a271df0fb7b",
   "metadata": {},
   "outputs": [],
   "source": [
    "! grep import /data/notebooks/repositories/LiCSBAS/bin/LiCSBAS12_loop_closure.py"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eb8f1812-83cd-4ffe-b46d-f58283768d32",
   "metadata": {},
   "source": [
    " "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "responsible-newman",
   "metadata": {},
   "source": [
    "**cells \\[ 13 \\]**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "silver-bristol",
   "metadata": {},
   "outputs": [],
   "source": [
    "! head -n 10 TS_GEOCml2GACOSclip/12loop/loop_info.txt\n",
    "! echo \"...\"\n",
    "! tail -n 30 TS_GEOCml2GACOSclip/12loop/loop_info.txt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "asian-railway",
   "metadata": {},
   "outputs": [],
   "source": [
    "! ls TS_GEOCml2GACOSclip/12loop/bad_loop_png"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "underlying-material",
   "metadata": {},
   "outputs": [],
   "source": [
    "# see one of the triplets of discarded ifgs:\n",
    "triplet = '20171226_20180101_20180113'\n",
    "#triplet = '20171108 20171226 20180101' # example of bad loop\n",
    "#triplet = '20161014 20161020 20161101' # example of candidate to bad loop that is used further\n",
    "triplet = triplet.replace(' ','_')\n",
    "png = os.path.join('TS_GEOCml2GACOSclip','12loop','bad_loop_png',triplet+'_loop.png')\n",
    "if not os.path.exists(png):\n",
    "    png = os.path.join('TS_GEOCml2GACOSclip','12loop','good_loop_png',triplet+'_loop.png')\n",
    "print('Note: the figure below was generated during 1st iteration, i.e. before referencing/without multi-prime')\n",
    "display(png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "treated-conversation",
   "metadata": {},
   "outputs": [],
   "source": [
    "# check also no loop ifgs:\n",
    "! echo \"There are \"`ls TS_GEOCml2GACOSclip/12no_loop_ifg_ras | wc -l`\" pairs identified as not forming a complete loop triplet\"\n",
    "! ls TS_GEOCml2GACOSclip/12no_loop_ifg_ras"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "understanding-suite",
   "metadata": {},
   "outputs": [],
   "source": [
    "#check network\n",
    "png='TS_GEOCml2GACOSclip/network/network12.png'\n",
    "display(png)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "41cf1de6-c375-4e9d-b1b5-c71bc41209a1",
   "metadata": {},
   "source": [
    "# Step 1-3: Small Baseline Inversion"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "electric-jewelry",
   "metadata": {},
   "source": [
    "### LiCSBAS13_sb_inv.py\n",
    "This script inverts the SB network of unw to obtain the time series and velocity using NSBAS ([López-Quiroz et al., 2009](https://www.sciencedirect.com/science/article/pii/S0926985109000251); [Doin et al., 2011](https://earth.esa.int/documents/10174/1573056/Presentation_small_baseline_NSBAS_Etna_deformation_Envisat.pdf)) approach.  \n",
    "A stable reference point is determined after the inversion.  \n",
    "RMS of the time series wrt median among all points is calculated for each point.  \n",
    "Then the point with minimum RMS and minimum n_gap is selected as new stable reference point."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "spiritual-corpus",
   "metadata": {},
   "source": [
    "**cells \\[ 14 \\]**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b3212738-2cad-43ca-a182-1102acf55bb0",
   "metadata": {},
   "outputs": [],
   "source": [
    "! LiCSBAS13_sb_inv.py -h \n",
    "# e.g. LiCSBAS13_sb_inv.py -d GEOCml10GACOS --n_para 1 --mem_size 1000\n",
    "printbold()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "incorporated-benjamin",
   "metadata": {},
   "outputs": [],
   "source": [
    "print('running the command')\n",
    "! LiCSBAS13_sb_inv.py -d GEOCml2GACOSclip --n_para 1 --mem_size 1000 --singular --nopngs #--only_sb\n",
    "#without --singular:  (original constrained SBAS method aka NSBAS)\n",
    "#ETA 25 min. ... although, on your local computer this would take approx. 5 minutes\n",
    "#\n",
    "#with --singular:\n",
    "#ETA 7 minutes, but separate subsets can be connected in lower precision\n",
    "#\n",
    "#with --only_sb\n",
    "#ETA 1 minute, but would skip all points that contain some NaN (only pixels with full SBAS network are inverted)\n",
    "\n",
    "printbold()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "increasing-phenomenon",
   "metadata": {},
   "source": [
    "LiCSBAS procedure will first invert matrix for points having all unwrapped values (no NaNs) using the full NSBAS set of equations (same result as for $\\bf d = \\bf G \\bf m$ where $\\bf d$ is set of unwrapped data related to $\\bf m$ that is a vector of incremental displacements, through design matrix $\\bf G$) and then runs for incomplete points to find $\\bf m$, filling empty values under assumption that the displacement is linear in time, i.e. solving points with missing values using set of equations:\n",
    "$$\\bf d = \\bf G \\bf m$$ $$0 = 0.0001(B m - vt - c)$$\n",
    "\n",
    "for details (e.g. definition of matrix B) see https://doi.org/10.3390/RS12030424  \n",
    "\n",
    "In our case, using **--singular** approach, we simplify the calculation of the gaps (missing values) by using only the partial range of SBAS-extracted cumulated increments, and filling the gaps assuming constant linear velocity connecting epochs between each gap. Although this leads to similar outputs, the computation is significantly simpler and keeps non-linear information better as connected clusters are not aligned along the estimated linear trend. Yet it is considered more prone to mistakes due to errors at the gap connecting epochs. In future release this can be replaced by extrapolation using Gaussian window in time.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "derived-durham",
   "metadata": {},
   "source": [
    "**cells \\[ 15 \\]**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "institutional-baking",
   "metadata": {},
   "outputs": [],
   "source": [
    "#check network\n",
    "png='TS_GEOCml2GACOSclip/network/network13.png'\n",
    "display(png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "caring-bidding",
   "metadata": {},
   "outputs": [],
   "source": [
    "! ls TS*/results/*.png\n",
    "png='TS_GEOCml2GACOSclip/results/vel.png'\n",
    "display(png)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "163ebdcf-5239-4148-9258-ed81e551861a",
   "metadata": {},
   "source": [
    "# Step 1-4: Calculate STD of Velocity by Bootstrap"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4969127f-5363-426e-a632-b6aa94f2d14b",
   "metadata": {},
   "source": [
    "### LiCSBAS14_vel_std.py\n",
    "This script primarily calculates the standard deviation of the velocity by the bootstrap method and STC (spatio-temporal consistency; [Hanssen et al., 2008](https://raw.githubusercontent.com/wiki/yumorishita/LiCSBAS/documents/Hanssen_2008.pdf)) but can also be used to optionally re-estimate velocity with outlier detection/removal."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "tropical-death",
   "metadata": {},
   "source": [
    "**cells \\[ 16 \\]**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "88d9a7d5-f230-424d-b88e-15b2a91246c3",
   "metadata": {},
   "outputs": [],
   "source": [
    "! LiCSBAS14_vel_std.py -h\n",
    "# e.g. LiCSBAS14_vel_std.py -t TS_GEOCml10GACOS\n",
    "printbold()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "present-smart",
   "metadata": {},
   "outputs": [],
   "source": [
    "print('running the command')\n",
    "! LiCSBAS14_vel_std.py -t TS_GEOCml2GACOSclip --mem_size 1000\n",
    "# ETA 5 sec\n",
    "printbold()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cleared-differential",
   "metadata": {},
   "source": [
    "**cells \\[ 17 \\]**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "extended-difference",
   "metadata": {},
   "outputs": [],
   "source": [
    "! ls TS*/results/*.png -alh\n",
    "png='TS_GEOCml2GACOSclip/results/vstd.png'\n",
    "display(png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "blessed-maryland",
   "metadata": {},
   "outputs": [],
   "source": [
    "png='TS_GEOCml2GACOSclip/results/stc.png'\n",
    "display(png)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7393457a-49dc-493f-bb10-335457822326",
   "metadata": {},
   "source": [
    "# Step 1-5: Mask Time Series"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cfa1dbd4-4413-46c7-8469-4a8e7d6f9c98",
   "metadata": {},
   "source": [
    "### LiCSBAS15_mask_ts.py\n",
    "\n",
    "This script makes a mask for time series using several noise indices.  \n",
    "The pixel is masked if any of the values of the noise indices for a pixel is worse (larger or smaller) than a specified threshold.  \n",
    "\n",
    "<img src=\"https://raw.githubusercontent.com/wiki/yumorishita/LiCSBAS/images/mask_ts.png\" width=\"900\">"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2e7de9fd-d4de-48e2-b3f2-38319fd5139a",
   "metadata": {},
   "source": [
    "| Noise index | Meaning |\n",
    "| --- | ----------- |\n",
    "| coh_avg | Average value of the interferometric coherence across the stack (0–1) |\n",
    "| n_unw | Number of unwrapped data used in the time series calculation (as a ratio to the number of images; i.e., if 1.5, n_unw=1.5*n_im) |\n",
    "| vstd | Standard deviation of the velocity (mm/yr) estimated in Step 1-4 |\n",
    "| maxTlen | Maximum time length of the connected network (years). The larger the value, the better the precision in the velocity estimate. |\n",
    "| n_gap | Number of gaps in the interferogram network and time series breaks. |\n",
    "| **n_gap_merged** | (only if set as parameter *--n_gap_use_merged*) Modified n_gap counting neighbouring gaps as one. |\n",
    "| stc | Spatio-temporal consistency (mm) (Hanssen et al., 2008), estimated in Step 1-4. |\n",
    "| n_ifg_noloop | Number of interferograms with no loops that cannot be checked by the loop closure and possibly have unidentified unwrapping errors. |\n",
    "| n_loop_err | Number of the unclosed loops after the network refinement in Step 1-2. |\n",
    "| **n_loop_err_rat** | Ratio of n_loop_err per number of checked loops. Will replace n_loop_err in the next major release. |\n",
    "| **loop_ph_avg_abs** | Average of the absolute wrapped phase loop closure (aka phase bias) (rad). |\n",
    "| resid_rms | RMS of residuals in the small baseline (SB) inversion (mm). |\n",
    "\n",
    "(bold: newly added signatures since the last release)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "proof-legislature",
   "metadata": {},
   "source": [
    "**cells \\[ 18 \\]**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "67ff8e49-e164-4ac8-b055-766d05318ca0",
   "metadata": {},
   "outputs": [],
   "source": [
    "! LiCSBAS15_mask_ts.py -h\n",
    "# e.g. LiCSBAS15_mask_ts.py -t TS_GEOCml10GACOS -c 0.05 -u 1.5 -v 100 -T 1 -g 10 -s 5  -i 50 -l 5 -r 2\n",
    "printbold()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "collect-california",
   "metadata": {},
   "outputs": [],
   "source": [
    "print('running the command')\n",
    "! LiCSBAS15_mask_ts.py -t TS_GEOCml2GACOSclip \n",
    "printbold()\n",
    "png = 'TS_GEOCml2GACOSclip/results/vel.mskd.png'\n",
    "display(png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "common-cruise",
   "metadata": {},
   "outputs": [],
   "source": [
    "print('Remind Step 1-2 - there were 64 interferograms not completing a loop! thus n_ifg_noloop=64 ')\n",
    "png = 'TS_GEOCml2GACOSclip/mask_ts.png'\n",
    "display(png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "martial-speech",
   "metadata": {},
   "outputs": [],
   "source": [
    "print('running the command')\n",
    "! LiCSBAS15_mask_ts.py -t TS_GEOCml2GACOSclip -i 70 -l 10 -r 2\n",
    "# ETA 2 sec\n",
    "printbold()\n",
    "png = 'TS_GEOCml2GACOSclip/results/vel.mskd.png'\n",
    "display(png)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "advanced-samoa",
   "metadata": {},
   "source": [
    "**cell \\[ 19 \\]**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "handled-segment",
   "metadata": {},
   "outputs": [],
   "source": [
    "png = 'TS_GEOCml2GACOSclip/mask_ts.png'\n",
    "display(png)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bb2d396b-bc38-4501-82d5-aeda0d388d8f",
   "metadata": {},
   "source": [
    "# Step 1-6: Filter (& Deramp) Time Series"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ed22b202-176e-4f6f-a5ac-a224b0c7549e",
   "metadata": {},
   "source": [
    "### LiCSBAS16_filt_ts.py\n",
    "This script applies spatio-temporal filter (HP in time and LP in space with gaussian kernel, same as StaMPS) to the time series of displacement.  \n",
    "Deramping (1D, bilinear, or 2D polynomial) can also be applied if -r option is used.  \n",
    "Topography-correlated components (linear with elevation) can also be subtracted with --hgt_linear option simultaneously with deramping before spatio-temporal filtering.  \n",
    "Further on, we incorporated DEM error estimation implemented in [LiCSBAS2](https://github.com/yumorishita/LiCSBAS2).\n",
    "The impact of filtering (deramp and linear elevation as well) can be visually checked by showing 16filt\\*/\\*png.  \n",
    "A stable reference point is determined after the filtering in the same manner as in Step 1-3.  \n",
    "\n",
    "### Previews:\n",
    "Deramping:  \n",
    "<img src=\"https://raw.githubusercontent.com/wiki/yumorishita/LiCSBAS/images/20180102_deramp.png\" width=\"600\">  \n",
    "\n",
    "Removing phase correlated to topography:  \n",
    "<img src=\"https://raw.githubusercontent.com/wiki/yumorishita/LiCSBAS/images/20180102_hgt_linear.png\" width=\"600\">  \n",
    "<img src=\"https://raw.githubusercontent.com/wiki/yumorishita/LiCSBAS/images/20180102_hgt_corr.png\" width=\"400\">  \n",
    "\n",
    "Spatio-temporal filtering:  \n",
    "<img src=\"https://raw.githubusercontent.com/wiki/yumorishita/LiCSBAS/images/20180102_filt.png\" width=\"600\">"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "anonymous-tucson",
   "metadata": {},
   "source": [
    "**cells \\[ 20 \\]**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bc8cb384-24ad-452f-b8f3-a227d08e5e73",
   "metadata": {},
   "outputs": [],
   "source": [
    "! LiCSBAS16_filt_ts.py -h\n",
    "# e.g. LiCSBAS16_filt_ts.py -t TS_GEOCml2GACOSclip -r 2 --n_para 1\n",
    "printbold()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "detected-litigation",
   "metadata": {},
   "outputs": [],
   "source": [
    "print('running the command')\n",
    "! LiCSBAS16_filt_ts.py -t TS_GEOCml2GACOSclip -r 1 --n_para 1 --hgt_linear --nofilter\n",
    "# ETA (on typical local computer): 10 min. with hgt_linear, 7 minutes without\n",
    "# Binder: >20 minutes if using spatio-temporal filter -> here we only deramp and use hgt_linear\n",
    "printbold()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "asian-advancement",
   "metadata": {},
   "source": [
    "**cells \\[ 21 \\]**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "confident-sugar",
   "metadata": {},
   "outputs": [],
   "source": [
    "print('first 30 files of filtered dataset')\n",
    "! ls TS_GEOCml2GACOSclip/16filt_cum | head -n 30"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "upset-blind",
   "metadata": {},
   "outputs": [],
   "source": [
    "#epoch = '20161026'\n",
    "epoch = '20161008'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "chief-paper",
   "metadata": {},
   "outputs": [],
   "source": [
    "print('effect of deramping:')\n",
    "png='TS_GEOCml2GACOSclip/16filt_cum/{}_deramp.png'.format(epoch)\n",
    "display(png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "hungry-limitation",
   "metadata": {},
   "outputs": [],
   "source": [
    "print('effect of removing linear correlation to height:')\n",
    "print('colour scale: 6pi/cycle')\n",
    "png='TS_GEOCml2GACOSclip/16filt_cum/{}_hgt_linear.png'.format(epoch)\n",
    "display(png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "greatest-brother",
   "metadata": {},
   "outputs": [],
   "source": [
    "print('correlation to height')\n",
    "png='TS_GEOCml2GACOSclip/16filt_cum/{}_hgt_corr.png'.format(epoch)\n",
    "display(png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "forced-browser",
   "metadata": {},
   "outputs": [],
   "source": [
    "print('Note: we have avoided spatio-temporal filter using --nofilter option. If run, you would see results here:')\n",
    "print('effect of spatio-temporal filter:')\n",
    "png='TS_GEOCml2GACOSclip/16filt_cum/{}_filt.png'.format(epoch)\n",
    "display(png)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a958550e-7ae6-4b79-a56c-d05aa75d9a29",
   "metadata": {},
   "source": [
    "# Viewing the final time series"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "wrapped-sunglasses",
   "metadata": {},
   "source": [
    "**cells \\[ 22 \\]**"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e87e6501-7b75-4565-bb84-82582410c644",
   "metadata": {},
   "source": [
    "The original LiCSBAS script offers advanced visualisation of various output parameters, allows for re-referencing the data etc."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5f12efed-6957-487a-a70f-22c197a59c73",
   "metadata": {},
   "outputs": [],
   "source": [
    "! LiCSBAS_plot_ts.py -h\n",
    "# e.g. LiCSBAS_plot_ts.py -i TS_GEOCml2GACOSclip/cum_filt.h5"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8e3aae0d-1184-4fca-8349-4cc9c541a367",
   "metadata": {},
   "source": [
    "<img src=\"https://raw.githubusercontent.com/wiki/yumorishita/LiCSBAS/images/sample_vel.png\" width=\"550\"> <img src=\"https://raw.githubusercontent.com/wiki/yumorishita/LiCSBAS/images/sample_ts.png\" width=\"500\">"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b85673f2-7373-45f8-88b7-c97f68bc0b66",
   "metadata": {},
   "outputs": [],
   "source": [
    "print('see the static example of the viewer')\n",
    "%run ../bin/LiCSBAS_plot_ts.py -i TS_GEOCml2GACOSclip/cum.h5 -p 53/39\n",
    "# %run ../bin/LiCSBAS_plot_ts.py -i TS_GEOCml2GACOSclip/cum_filt.h5 -p 53/39"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5cd9a702-4603-4d92-914b-bf08a897dc73",
   "metadata": {},
   "outputs": [],
   "source": [
    "print('use the viewer to see filtered (deramped) dataset - note the time series lines')\n",
    "%run ../bin/LiCSBAS_plot_ts.py -i TS_GEOCml2GACOSclip/cum_filt.h5 -p 53/39"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9084ec5c-cbb1-4932-9b58-5fc1a93fc891",
   "metadata": {},
   "outputs": [],
   "source": [
    "print('generating png file of time series for given coordinates')\n",
    "#! LiCSBAS_plot_ts.py -i TS_GEOCml2GACOSclip/cum.h5 --ts_png test_ts.png -p 53/39 >/dev/null\n",
    "! LiCSBAS_plot_ts.py -i TS_GEOCml2GACOSclip/cum_filt.h5 --ts_png test_ts.png -p 53/39 >/dev/null\n",
    "display('test_ts.png')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a7ac74be-9200-480a-b276-39014e2686bd",
   "metadata": {},
   "source": [
    "It is, however, convenient to transform the output data to a standard NetCDF in geographic coordinates, so we can use any standard program for further processing. Following optional scripts will export the results to a standard GeoTIFF file or a NetCDF datacube that you may use both in any GIS environment, or e.g. to load to python via xarray etc.\n",
    "\n",
    "For more information (and example on QGIS Time Series viewer) see [LiCSBAS wiki](https://github.com/yumorishita/LiCSBAS/wiki/4_other_tools) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0dbc3a5a-ddbc-473e-a01b-b94fac666da1",
   "metadata": {},
   "outputs": [],
   "source": [
    "# to export velocities to a geotiff\n",
    "! LiCSBAS_flt2geotiff.py -h\n",
    "# e.g. LiCSBAS_flt2geotiff.py -i TS_GEOCml2*/results/vel -p GEOCml2GACOSclip/EQA.dem_par -o 022D_04826_121209.vel.geo.tif\n",
    "\n",
    "# to export to Google Earth kmz tiles:\n",
    "#! LiCSBAS_color_geotiff2tiles.py -h"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4a25e7e9-7c87-45fb-8e51-ca250761c2bb",
   "metadata": {},
   "outputs": [],
   "source": [
    "# to export all relevant data to a NetCDF datacube\n",
    "! LiCSBAS_out2nc.py -h\n",
    "# e.g. LiCSBAS_out2nc.py -i TS_GEOCml2*/cum.h5 -o 022D_04826_121209.nc\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ef8e344e-0b0d-4eff-8274-e79011f584b0",
   "metadata": {},
   "source": [
    "Let's export LiCSBAS output to the standard NetCDF datacube, and then view it through jupyter notebook. You may also download the NetCDF file and load into QGIS, ncview etc."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "98a701d4-1614-42e4-8b38-1d9d080bc44a",
   "metadata": {},
   "outputs": [],
   "source": [
    "! LiCSBAS_out2nc.py -i TS_GEOCml2*/cum_filt.h5 -o 022D_04826_121209.nc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ddbf3e60-3526-415e-9b13-d79fe7e33769",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# re-define in case we run this at the end without finalising the processing\n",
    "import os\n",
    "frame = '022D_04826_121209'\n",
    "framedir = os.path.join(os.environ['HOME'],frame)\n",
    "os.chdir(framedir)\n",
    "\n",
    "# loading the netcdf as xarray dataset\n",
    "import xarray as xr\n",
    "ncfile='022D_04826_121209.nc'\n",
    "if not os.path.exists(ncfile):\n",
    "    os.system('LiCSBAS_out2nc.py -i TS_GEOCml2GACOSclip/cum_filt.h5 -o 022D_04826_121209.nc')\n",
    "\n",
    "cube=xr.open_dataset(ncfile)\n",
    "\n",
    "# explore the contents\n",
    "cube"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bda1c4a5-02b4-42fa-aca6-7476fa375180",
   "metadata": {},
   "source": [
    "Finally we can visualise the results using an original approach of making the great [pyGMT](https://www.generic-mapping-tools.org/egu22pygmt/intro.html) interactive - with great appreciation to Miss Rochelle Pun, an outstanding student that was selected for the [COMET Summer MSc Internship 2024](https://comet.nerc.ac.uk/internships) during which she developed her [IntPyGmt tool](https://github.com/chelle0425/IntPyGMT) as a side activity."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "166f0416-d33e-4188-aeab-196ae786eba6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# visualise time series\n",
    "# (c) 2024, Rochelle Pun during her short COMET Summer MSc Internship 2024\n",
    "# originally published in https://github.com/chelle0425/IntPyGMT\n",
    "# adapted to LiCSBAS\n",
    "\n",
    "from LiCSBAS_extras import pygmt_plot_interactive\n",
    "\n",
    "title = 'LiCSBAS Tutorial: time series'\n",
    "lims = [-15,60]\n",
    "label='velocity [mm/year]'\n",
    "cmap='polar'\n",
    "\n",
    "# note this will work only with ipympl installed:\n",
    "%matplotlib widget\n",
    "\n",
    "pygmt_plot_interactive(cube, title, label, lims, cmap,\n",
    "                           photobg=False, plotvec=None)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1ad553d7-1341-48ae-a8ce-4e0ab6dd75d4",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "bc503173-9368-4abe-bbae-4ae945359d5b",
   "metadata": {},
   "source": [
    "## Citations\n",
    "\n",
    "Morishita, Y.; Lazecky, M.; Wright, T.J.; Weiss, J.R.; Elliott, J.R.; Hooper, A. LiCSBAS: An Open-Source InSAR Time Series Analysis Package Integrated with the LiCSAR Automated Sentinel-1 InSAR Processor. *Remote Sens.* **2020**, *12*, 424, https://doi.org/10.3390/RS12030424.\n",
    "\n",
    "Morishita, Y.: Nationwide urban ground deformation monitoring in Japan using Sentinel-1 LiCSAR products and LiCSBAS. *Prog. Earth Planet. Sci.* **2021**, *8*, 6,  https://doi.org/10.1186/s40645-020-00402-7.\n",
    "\n",
    "Lazecký, M.; Spaans, K.; González, P.J.; Maghsoudi, Y.; Morishita, Y.; Albino, F.; Elliott, J.; Greenall, N.; Hatton, E.; Hooper, A.; Juncu, D.; McDougall, A.; Walters, R.J.; Watson, C.S.; Weiss, J.R.; Wright, T.J. LiCSAR: An Automatic InSAR Tool for Measuring and Monitoring Tectonic and Volcanic Activity. *Remote Sens.* **2020**, *12*, 2430, https://doi.org/10.3390/rs12152430.\n",
    "\n",
    "Lazecký, M.; Ou, Q.; McGrath, J.; Payne, J.; Espin, P.; Hooper, A.; Wright, T. Strategies for improving and correcting unwrapped interferograms implemented in LiCSBAS. *Procedia Computer Science* **2024**, *239*, 2408-2412, https://doi.org/10.1016/j.procs.2024.06.435.\n",
    "\n",
    "## Acknowledgements\n",
    "\n",
    "The core and main functionality of LiCSBAS has been accomplished by Dr. Yu Morishita (Geospatial Information Authority of Japan) during his visit at the University of Leeds, funded by JSPS Overseas Research Fellowship (June 2018-March 2020). LiCSBAS is further developed by the team at the University of Leeds, under the supervision of Prof. Andy Hooper and Prof. Tim Wright.  \n",
    "\n",
    "COMET is the UK Natural Environment Research Council's Centre for the Observation and Modelling of Earthquakes, Volcanoes and Tectonics. LiCSAR is developed as part of the NERC large grant, \"Looking inside the continents from Space\" (NE/K010867/1). LiCSAR contains modified Copernicus Sentinel-1 data analysed by the COMET. LiCSAR uses [JASMIN](http://jasmin.ac.uk), the UK’s collaborative data analysis environment.  \n",
    "\n",
    "The [Scientific Colour Maps](http://www.fabiocrameri.ch/colourmaps.php) ([Crameri, 2018](https://doi.org/10.5194/gmd-11-2541-2018)) is used in LiCSBAS.\n",
    "\n",
    "\n",
    "[<img src=\"https://raw.githubusercontent.com/wiki/yumorishita/LiCSBAS/images/COMET_logo.png\"  width=\"100\">](https://comet.nerc.ac.uk/)   [<img src=\"https://raw.githubusercontent.com/wiki/yumorishita/LiCSBAS/images/logo-leeds.png\"  width=\"250\">](https://environment.leeds.ac.uk/see/)  [<img src=\"https://raw.githubusercontent.com/wiki/yumorishita/LiCSBAS/images/LiCS_logo.jpg\"  width=\"150\">](https://comet.nerc.ac.uk/COMET-LiCS-portal/)  [<img src=\"https://mybinder.org/static/logo.svg\" width=\"200\">](https://mybinder.org/)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.12.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}