stripmapApp_Alos.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"../MathBasics/Figs/Banner.png\" width=\"100%\" />\n",
    "<font face=\"Calibri\">\n",
    "<br>\n",
    "<font size=\"6\"> <b> InSAR pair processing with stripmapApp  <b> </font>\n",
    "\n",
    "<font size=\"4\"> <b> <font color='rgba(200,0,0,0.2)'> </font> </b> </font>\n",
    "\n",
    "<br> <img src=\"../MathBasics/Figs/NASALogo.png\" width=\"250\" align=\"right\" /> <br> \n",
    "<font size=\"4\"> <b> Heresh Fattahi</b> \n",
    "<font size=\"3\">  <br>\n",
    "<font> <b>Date: </b> Aug 10, 2020 </font>\n",
    "</font>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In conventional stripmap Synthetic Aperture Radar(SAR) imaging mode, the radar antenna\n",
    "is fixed to a specific direction, illuminating a single swath of the scene with a fixed squint angle (i.e., the angle between the radar beam and the cross-track direction). The imaging swath width can be increased using the scanning SAR (ScanSAR) or Terrain Observation by Progressive Scan(TOPS). In this notebook we focus on interferometric processing of stripmap data using stripmapApp.py. \n",
    "\n",
    "The stripmap mode has been used by sevreal SAR missions, such as Envisat, ERS, RadarSAT-1, Radarsat-2, ALOS-1, Cosmo Sky-Med and TerraSAR-X. Although Sentinel-1 A/B and ALOS-2 are capable of acuqiring SAR data with stripmap mode, their operational imaging modes are TOPS and ScanSAR respectively. Both missions have been acquiring stripmap data over certain regions.\n",
    "\n",
    "For processing TOPS data using topsApp, please see the topsApp notebook. However, we recommend that new InSAR users may start with the stripmapApp notebook first, and then try the topsApp notebook. \n",
    "\n",
    "The detailed algorithms for stripmap processing and TOPS processing implemented in ISCE software can be found in the following literatures:\n",
    "\n",
    "### stripmapApp:\n",
    "\n",
    "H. Fattahi, M. Simons, and P. Agram,  \"InSAR Time-Series Estimation of the Ionospheric Phase Delay: An Extension of the Split Range-Spectrum Technique\", IEEE Trans. Geosci. Remote Sens., vol. 55, no. 10, 5984-5996, 2017.\n",
    "(https://ieeexplore.ieee.org/abstract/document/7987747/)\n",
    "\n",
    "### topsApp:\n",
    "\n",
    "H. Fattahi, P. Agram, and M. Simons, “A network-based enhanced spectral diversity approach for TOPS time-series analysis,” IEEE Trans. Geosci. Remote Sens., vol. 55, no. 2, pp. 777–786, Feb. 2017. (https://ieeexplore.ieee.org/abstract/document/7637021/)\n",
    "\n",
    "### ISCE framework:\n",
    "Rosen et al, IGARSS 2018 [Complete reference here] \n",
    "\n",
    "![title](notebook_docs/Stripmap_Tops.png)\n",
    "\n",
    "(Figure from Fattahi et. al., 2017)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "# stripmapApp (general overview)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "stripmapApp.py is an ISCE application, designed for interferometric processing of SAR data acquired with stripmap mode onboard platforms with precise orbits. This ISCE application is equivalent to insarApp (an older ISCE application which has been widely used by ISCE users), in this sense that both apps process stripmap data. Although the naming convention of the products and interface of the App is similiar to insarApp, the core modules used in stripmapApp are fundamentally different. The main features of stripmapApp includes the following:\n",
    "\n",
    "#### a) Focusing RAW data to native Doppler:\n",
    "\n",
    "If processing starts from RAW data, JPL's ROI software is used for focusing the raw data to SLC(Single Look Complex) SAR images in native Doppler geometry. stripmapApp does not use the motion compensation algorithm which was used in insarApp. \n",
    "\n",
    "#### b) Interferometric processing of SLCs in native or zero Doppler geometry \n",
    "If the input data are SLC images, then focusing is not required and will be skipped. stripmapApp can process SLCs focused to zero or native Dopplers\n",
    "\n",
    "#### c) Coregistration using SAR acquisition geometry (Orbit + DEM) \n",
    "The geometry module of ISCE is used for coregistration of SAR images, i.e., range and azimuth offsets are computed for each pixel using SAR acquisition geometry, orbit information and an existing Digital Elevation Model(DEM). The geometrical offsets are refined with a small constant shift in range and azimuth directions. The constant shifts are estimated using incoherent cross-correlation of the two SAR images already coregistered using pure geometrical information. \n",
    "\n",
    "#### d) More optional precise coregistration step\n",
    "An optional step called \"rubbersheeting\" is available for more precise coregistration. If \"rubbersheeting\" is requested, a dense azimuth offsets is computed using incoherent cross-correlation between the two SAR images, and is added to the geometrical offsets for more precise coregistration. Rubbersheeting may be required if SAR images are affected by ionospheric scintillation. \n",
    "\n",
    "#### e) Ionospheric phase estimation\n",
    "Split Range-Spectrum technique and ionospheric phase estimation are available as optional processing steps.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Prepare directories, download raw data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Importing some python modules and setting up some variables: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'/scratch/05861/tg851601/Alos'"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# If you want to use the staged data on S3 bucket, change this flag to True. \n",
    "# when the falg is False, the notebook downloads the ALOS raw data from ASF and \n",
    "# ISCE automatically downloads an SRTM DEM for processing\n",
    "Use_Staged_Data = False\n",
    "\n",
    "import os\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "# import gdal\n",
    "from osgeo import gdal\n",
    "from tqdm import tqdm\n",
    "import urllib.request\n",
    "\n",
    "\n",
    "\n",
    "#\n",
    "ASF_USER = \"\"\n",
    "ASF_PASS = \"\"\n",
    "\n",
    "# the working directory:\n",
    "# home_dir = os.path.join(os.getenv(\"HOME\"), \"work\")\n",
    "home_dir = os.path.join(os.getenv(\"SCRATCHDIR\"), \"Alos\")\n",
    "PROCESS_DIR = os.path.join(home_dir, \"Hawaii_ALOS1\")\n",
    "DATA_DIR =  os.path.join(PROCESS_DIR, \"data\")\n",
    "\n",
    "\n",
    "# if the ASF user/pass is not provided above, try to read it from ~/.netrc file\n",
    "if (len(ASF_PASS)==0 | len(ASF_USER)==0) & (os.path.exists(os.path.join(os.getenv(\"HOME\"), \".netrc\"))):\n",
    "    netrc_path = os.path.join(os.getenv(\"HOME\"), \".netrc\")\n",
    "    count = len(open(netrc_path).readlines(  ))\n",
    "    if count == 1:\n",
    "        file = open(os.path.join(os.getenv(\"HOME\"), \".netrc\"), \"r\")\n",
    "        contents = file.read().split(\" \")\n",
    "        ASF_USER = contents[3]\n",
    "        ASF_PASS = contents[5]\n",
    "        file.close()\n",
    "    else:\n",
    "        ASF_USER = np.loadtxt(os.path.join(os.getenv(\"HOME\"), \".netrc\"), skiprows=1, usecols=1, dtype=str)[0]\n",
    "        ASF_PASS = np.loadtxt(os.path.join(os.getenv(\"HOME\"), \".netrc\"), skiprows=1, usecols=1, dtype=str)[1]\n",
    "\n",
    "if (len(ASF_PASS)==0 | len(ASF_USER)==0) | (not os.path.exists(os.path.join(os.getenv(\"HOME\"), \".netrc\"))) :\n",
    "    print(\"WARNING: The ASF USER pass needs to be included in ~/.netrc file.\")\n",
    "    print(\" The ~/.netrc file does not exixt or is not setup properly.\")\n",
    "    print(\"Follow this link for instructions on setting up your ~/.netrc file\")\n",
    "    print(\"\"\"If you wish to download the data from ASF please make sure: \n",
    "             1) you have valid earthdata login and password stored in your ~/.netrc file\n",
    "             2) make sure that you have logged into ASF vertex page and have accepted the EULA agreement.\"\"\")\n",
    "    print(\"Without further actions you will still be able to run this notebook using already available data on S3 bucket.\")\n",
    "    print(\"Using data on S3 bucket ...\")\n",
    "    Use_Staged_Data = True\n",
    "\n",
    "if Use_Staged_Data:\n",
    "    print(\"Using the staged data for this notebook has been turned on.\")\n",
    "\n",
    "import shutil\n",
    "# shutil.copytree(os.getenv('WORKDIR') + '/TESTDATA_STRIPMAPAPP_ALOS/',home_dir)\n",
    "shutil.copytree('/work2/05861/tg851601/stampede2/insarlab/TESTDATA_STRIPMAPAPP_ALOS',home_dir,dirs_exist_ok=True)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def configure_inputs(outDir, Use_Staged_Data): \n",
    "\n",
    "    \"\"\"Wraite Configuration files for ISCE2 stripmapApp to process NISAR sample products\"\"\"\n",
    "    # FA 5/2024: It seems hardwired \n",
    "    cmd_reference_config = '''<component>\n",
    "    <property name=\"IMAGEFILE\">\n",
    "       <value>[data/20110119/ALPSRP265743230-L1.0/IMG-HH-ALPSRP265743230-H1.0__D]</value>\n",
    "    </property>\n",
    "    <property name=\"LEADERFILE\">\n",
    "        <value>[data/20110119/ALPSRP265743230-L1.0/LED-ALPSRP265743230-H1.0__D]</value>\n",
    "    </property>\n",
    "    <property name=\"OUTPUT\">\n",
    "        <value>20110119</value>\n",
    "    </property>\n",
    "</component>'''\n",
    "\n",
    "    print(\"writing reference.xml\")\n",
    "    with open(os.path.join(outDir,\"reference.xml\"), \"w\") as fid:\n",
    "        fid.write(cmd_reference_config)\n",
    "    \n",
    "    cmd_secondary_config = '''<component>\n",
    "    <property name=\"IMAGEFILE\">\n",
    "        <value>[data/20110306/ALPSRP272453230-L1.0/IMG-HH-ALPSRP272453230-H1.0__D]</value>\n",
    "    </property>\n",
    "    <property name=\"LEADERFILE\">\n",
    "        <value>[data/20110306/ALPSRP272453230-L1.0/LED-ALPSRP272453230-H1.0__D]</value>\n",
    "    </property>\n",
    "    <property name=\"OUTPUT\">\n",
    "        <value>20110306</value>\n",
    "    </property>\n",
    "\n",
    "</component>'''\n",
    "    \n",
    "    print(\"writing secondary.xml\")\n",
    "    with open(os.path.join(outDir,\"secondary.xml\"), \"w\") as fid:\n",
    "        fid.write(cmd_secondary_config)\n",
    "\n",
    "    if Use_Staged_Data:\n",
    "        cmd_stripmap_config = '''<?xml version=\"1.0\" encoding=\"UTF-8\"?>\n",
    "<stripmapApp>\n",
    "  <component name=\"insar\">\n",
    "    <property name=\"sensor name\">ALOS</property>\n",
    "    <component name=\"reference\">\n",
    "        <catalog>reference.xml</catalog>\n",
    "    </component>\n",
    "    <component name=\"secondary\">\n",
    "        <catalog>secondary.xml</catalog>\n",
    "    </component>\n",
    "\n",
    "    <property name=\"demFilename\">\n",
    "        <value>demLat_N18_N21_Lon_W156_W154.dem.wgs84</value>\n",
    "    </property>\n",
    "\n",
    "    <property name=\"unwrapper name\">icu</property>\n",
    "\n",
    "    <property name=\"do split spectrum\">False</property>\n",
    "\n",
    "    <property name=\"do dispersive\">False</property>\n",
    "\n",
    "</component>\n",
    "</stripmapApp>'''\n",
    "    else:\n",
    "            \n",
    "        cmd_stripmap_config = '''<?xml version=\"1.0\" encoding=\"UTF-8\"?>\n",
    "<stripmapApp>\n",
    "  <component name=\"insar\">\n",
    "    <property name=\"sensor name\">ALOS</property>\n",
    "    <component name=\"reference\">\n",
    "        <catalog>reference.xml</catalog>\n",
    "    </component>\n",
    "    <component name=\"secondary\">\n",
    "        <catalog>secondary.xml</catalog>\n",
    "    </component>\n",
    "\n",
    "    <!--\n",
    "    <property name=\"demFilename\">\n",
    "        <value>demLat_N18_N21_Lon_W156_W154.dem.wgs84</value>\n",
    "    </property>\n",
    "    -->\n",
    "\n",
    "    <property name=\"unwrapper name\">icu</property>\n",
    "\n",
    "    <property name=\"do split spectrum\">False</property>\n",
    "\n",
    "    <property name=\"do dispersive\">False</property>\n",
    "\n",
    "</component>\n",
    "</stripmapApp>'''\n",
    "\n",
    "    print(\"writing stripmapApp.xml\")\n",
    "    with open(os.path.join(outDir,\"stripmapApp.xml\"), \"w\") as fid:\n",
    "        fid.write(cmd_stripmap_config)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Check if the PROCESS_DIR and DATA_DIR already exist. If they don't exist, we create them:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/scratch/05861/tg851601/Alos/Hawaii_ALOS1  already exists!\n",
      "/scratch/05861/tg851601/Alos/Hawaii_ALOS1/data  already exists!\n"
     ]
    }
   ],
   "source": [
    "if not os.path.exists(PROCESS_DIR):\n",
    "    print(\"create \", PROCESS_DIR)\n",
    "    os.makedirs(PROCESS_DIR)\n",
    "else:\n",
    "    print(PROCESS_DIR, \" already exists!\")\n",
    "\n",
    "if not os.path.exists(DATA_DIR):\n",
    "    print(\"create \", DATA_DIR)\n",
    "    os.makedirs(DATA_DIR)\n",
    "else:\n",
    "    print(DATA_DIR, \" already exists!\")\n",
    "\n",
    "\n",
    "os.chdir(DATA_DIR)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this tutorial we will process two ALOS1 PALSAR acquistions over Hawaii. The two acquisitions cover a dike opening even in March 2011.\n",
    "![title](notebook_docs/Hawaii.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Download two ALOS-1 acquistions from ASF using the following command:\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now the data should be downloading from ASF. Check the terminal that you strated your jupyter notebook from and you should see the progress bar for the download "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#Alternative ssara command\n",
    "#!ssara_federated_query.py --platform=ALOS --intersectsWith='POLYGON((-155.3 19.5, -155.3 19.8,-155.0 19.8,-155.0 19.5, -155.3 19.5 ))' --print --kml --flightDirection=D --beamMode=FBS,FBD --relativeOrbit=598   -s 2011-01-17 -e 2011-03-07 --download"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'/scratch/05861/tg851601/Alos'"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "os.chdir(home_dir)\n",
    "os.getcwd()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "unzip the downloaded files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Archive:  ALPSRP265743230-L1.0.zip\n",
      "  inflating: ALPSRP265743230-L1.0/IMG-HH-ALPSRP265743230-H1.0__D  \n",
      "  inflating: ALPSRP265743230-L1.0/IMG-HV-ALPSRP265743230-H1.0__D  \n",
      "  inflating: ALPSRP265743230-L1.0/VOL-ALPSRP265743230-H1.0__D  \n",
      "  inflating: ALPSRP265743230-L1.0/LED-ALPSRP265743230-H1.0__D  \n",
      "  inflating: ALPSRP265743230-L1.0/TRL-ALPSRP265743230-H1.0__D  \n",
      "  inflating: ALPSRP265743230-L1.0/workreport  \n",
      "  inflating: ALPSRP265743230-L1.0/ALPSRP265743230.l0.workreport  \n",
      "Archive:  ALPSRP272453230-L1.0.zip\n",
      "  inflating: ALPSRP272453230-L1.0/IMG-HH-ALPSRP272453230-H1.0__D  \n",
      "  inflating: ALPSRP272453230-L1.0/IMG-HV-ALPSRP272453230-H1.0__D  \n",
      "  inflating: ALPSRP272453230-L1.0/VOL-ALPSRP272453230-H1.0__D  \n",
      "  inflating: ALPSRP272453230-L1.0/LED-ALPSRP272453230-H1.0__D  \n",
      "  inflating: ALPSRP272453230-L1.0/TRL-ALPSRP272453230-H1.0__D  \n",
      "  inflating: ALPSRP272453230-L1.0/workreport  \n",
      "  inflating: ALPSRP272453230-L1.0/ALPSRP272453230.l0.workreport  \n"
     ]
    }
   ],
   "source": [
    "!unzip -o ALPSRP265743230-L1.0.zip\n",
    "!unzip -o ALPSRP272453230-L1.0.zip"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " looking at the unzipped directories there are multiple files:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ALPSRP265743230.l0.workreport\tLED-ALPSRP265743230-H1.0__D  workreport\n",
      "IMG-HH-ALPSRP265743230-H1.0__D\tTRL-ALPSRP265743230-H1.0__D\n",
      "IMG-HV-ALPSRP265743230-H1.0__D\tVOL-ALPSRP265743230-H1.0__D\n"
     ]
    }
   ],
   "source": [
    "!ls ALPSRP265743230-L1.0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When you download PALSAR data from a data provider, each frame comprises an image data file and an image leader file, as well as possibly some other ancillary files that are not used by ISCE.  \n",
    "\n",
    "Files with IMG as prefix are images. \n",
    "Files with LED as prefix are leaders. \n",
    "\n",
    "The leader file contains parameters of the sensor that are relevant to the imaging mode, all the information necessary to process the data.  The data file contains the raw data samples if Level 1.0 raw data (this is just a different name from what other satellites call Level 0) and processed imagery if Level 1.1 or 1.5 image data.  The naming convention for these files is standardized across data archives, and has the following taxonomy:\n",
    "\n",
    "\n",
    "![title](notebook_docs/ALOS1_PALSAR.png)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To see the acquisition date of this PALSAR acquisition we can look at the following file:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Odi_ProductManagementNo            = \"X0200065\"\n",
      "Odi_ProductManagementBranchNo      = \"815\"\n",
      "Scs_SceneID                        = \"ALPSRP265743230\"\n",
      "Pds_ProductID                      = \"H1.0__D\"\n",
      "Pds_OrbitDataPrecision             = \"Precision\"\n",
      "Pds_AttitudeDataPrecision          = \"Standard\"\n",
      "Img_SceneCenterDateTime            = \"20110119 20:33:44.268\"\n",
      "Img_SceneStartDateTime             = \"20110119 20:33:36.068\"\n",
      "Img_SceneEndDateTime               = \"20110119 20:33:52.468\"\n",
      "Img_NoOfPSR10LossLines             = \"0\"\n",
      "Pdi_ProductDataSize                = \"756.4\"\n",
      "Pdi_CntOfL10ProductFileName        = \"5\"\n",
      "Pdi_L10ProductFileName01           = \"VOL-ALPSRP265743230-H1.0__D\"\n",
      "Pdi_L10ProductFileName02           = \"LED-ALPSRP265743230-H1.0__D\"\n",
      "Pdi_L10ProductFileName03           = \"IMG-HH-ALPSRP265743230-H1.0__D\"\n",
      "Pdi_L10ProductFileName04           = \"IMG-HV-ALPSRP265743230-H1.0__D\"\n",
      "Pdi_L10ProductFileName05           = \"TRL-ALPSRP265743230-H1.0__D\"\n",
      "Pdi_NoOfPixels                     = \"5968\"\n",
      "Pdi_NoOfLines                      = \"31478\"\n",
      "Pdi_ProductFormat                  = \"CEOS\"\n",
      "Ach_TimeCheck                      = \"OK\"\n",
      "Ach_AttitudeCheck                  = \"OK\"\n",
      "Ach_AbsoluteNavigationStatus       = \"OK\"\n",
      "Ach_HouseKeepingDataCheck          = \"OK\"\n",
      "Ach_OrbitCheck                     = \"OK\"\n",
      "Ach_OnBoardAttitudeCheck           = \"OK\"\n",
      "Ach_LossLines                      = \"OK\"\n",
      "Ach_AbsoluteNavigationTime         = \"OK\"\n",
      "Ach_CalFlagCheck                   = \"OK\"\n",
      "Rst_PSR10PrcStatus                 = \"0000\"\n",
      "Rst_PSR10PrcStartTime              = \"20150112 05:49:19\"\n",
      "Rst_PSR10PrcEndTime                = \"20150112 05:51:48\"\n",
      "Rst_PSR_AbsoluteNavigationTime     = \"0000\"\n",
      "Ver_PSR10Prc                       = \"09.03\"\n",
      "Ver_PSR10AutoCheck                 = \"08.00\"\n",
      "Ver_PSR10ConstantForProduct        = \"08.00\"\n",
      "Ver_PSR10Product                   = \"08.00\"\n",
      "Ver_OS_VersionInDataProcessingUnit = \"Red Hat Enterprise Linux Server release 5.3\"\n"
     ]
    }
   ],
   "source": [
    "!cat ALPSRP265743230-L1.0/ALPSRP265743230.l0.workreport"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Img_SceneCenterDateTime            = \"20110119 20:33:44.268\"\n",
      "Img_SceneCenterDateTime            = \"20110306 20:32:42.519\"\n"
     ]
    }
   ],
   "source": [
    "!grep Img_SceneCenterDateTime ALPSRP265743230-L1.0/ALPSRP265743230.l0.workreport\n",
    "!grep Img_SceneCenterDateTime ALPSRP272453230-L1.0/ALPSRP272453230.l0.workreport"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "for clarity let's create two directories for the two acquisition dates and move the unziped folders there:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'/scratch/05861/tg851601/Alos/Hawaii_ALOS1/data/20110306'"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "!mkdir 20110119\n",
    "!mkdir 20110306\n",
    "!mv ALPSRP265743230-L1.0 20110119\n",
    "!mv ALPSRP272453230-L1.0 20110306\n",
    "shutil.move('20110119', DATA_DIR)\n",
    "shutil.move('20110306', DATA_DIR)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have the data ready let's cd to the main PROCESS directroy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "os.chdir(PROCESS_DIR)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To make sure where we are, run pwd:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/scratch/05861/tg851601/Alos/Hawaii_ALOS1\n"
     ]
    }
   ],
   "source": [
    "!pwd"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Setting up input xml files for processing with stripmapApp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create the input configuration files (refernce.xml, secondary.xml, stripmapApp.xml) to configure the inputs and the processing parameters.\n",
    "The configurations files can be created using your favorit editor or by calling the \"configure\" funstion which is defined at the top of this notebook:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Will not be using S3 pre-staged data, The DEM will be automatically downloaded by ISCE\n"
     ]
    }
   ],
   "source": [
    "if Use_Staged_Data:\n",
    "    # Check if a stage file from S3 already exist, if not try and download it\n",
    "    if not os.path.isfile('demLat_N18_N21_Lon_W156_W154.dem.wgs84'):\n",
    "        !aws s3 cp s3://asf-jupyter-data/fattahi-backup/demLat_N18_N21_Lon_W156_W154.dem.wgs84 demLat_N18_N21_Lon_W156_W154.dem.wgs84\n",
    "        !aws s3 cp s3://asf-jupyter-data/fattahi-backup/demLat_N18_N21_Lon_W156_W154.dem.wgs84.vrt demLat_N18_N21_Lon_W156_W154.dem.wgs84.vrt\n",
    "        !aws s3 cp s3://asf-jupyter-data/fattahi-backup/demLat_N18_N21_Lon_W156_W154.dem.wgs84.xml demLat_N18_N21_Lon_W156_W154.dem.wgs84.xml\n",
    "    else:\n",
    "        print(\"The DEM already exists\")\n",
    "else:\n",
    "    print(\"Will not be using S3 pre-staged data, The DEM will be automatically downloaded by ISCE\")\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "writing reference.xml\n",
      "writing secondary.xml\n",
      "writing stripmapApp.xml\n"
     ]
    }
   ],
   "source": [
    "configure_inputs(PROCESS_DIR, Use_Staged_Data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here is an example refernce.xml file for this tutorial:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### reference.xml"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```xml\n",
    "<component>\n",
    "    <property name=\"IMAGEFILE\">\n",
    "       <value>[data/20110119/ALPSRP265743230-L1.0/IMG-HH-ALPSRP265743230-H1.0__D]</value>\n",
    "    </property>\n",
    "    <property name=\"LEADERFILE\">\n",
    "        <value>[data/20110119/ALPSRP265743230-L1.0/LED-ALPSRP265743230-H1.0__D]</value>\n",
    "    </property>\n",
    "    <property name=\"OUTPUT\">\n",
    "        <value>20110119</value>\n",
    "    </property>\n",
    "</component>\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### secondary.xml"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```xml\n",
    "<component>\n",
    "    <property name=\"IMAGEFILE\">\n",
    "        <value>[data/20110306/ALPSRP272453230-L1.0/IMG-HH-ALPSRP272453230-H1.0__D]</value>\n",
    "    </property>\n",
    "    <property name=\"LEADERFILE\">\n",
    "        <value>[data/20110306/ALPSRP272453230-L1.0/LED-ALPSRP272453230-H1.0__D]</value>\n",
    "    </property>\n",
    "    <property name=\"OUTPUT\">\n",
    "        <value>20110306</value>\n",
    "    </property>\n",
    "\n",
    "</component>\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### stripmapApp.xml\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```xml\n",
    "<?xml version=\"1.0\" encoding=\"UTF-8\"?>\n",
    "<stripmapApp>\n",
    "  <component name=\"insar\">\n",
    "    <property name=\"sensor name\">ALOS</property>\n",
    "    <component name=\"reference\">\n",
    "        <catalog>refernce.xml</catalog>\n",
    "    </component>\n",
    "    <component name=\"secondary\">\n",
    "        <catalog>secondary.xml</catalog>\n",
    "    </component>\n",
    "\n",
    "    <!--  \n",
    "    <property name=\"demFilename\">\n",
    "        <value>demLat_N18_N21_Lon_W156_W154.dem.wgs84</value>\n",
    "    </property>\n",
    "    -->\n",
    "      \n",
    "    <property name=\"unwrapper name\">icu</property>\n",
    "     \n",
    "</component>\n",
    "</stripmapApp>\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br>\n",
    "<div class=\"alert alert-info\">\n",
    "<b>Note :</b> \n",
    "\n",
    "demFilename is commented out in the stripmapApp.xml. This means that user has not specified the DEM. Therefore, isce looks online and download the SRTM dem.\n",
    "\n",
    "</div>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After downloading the data to process, and setting up the input xml files, we are ready to start processing with stripmapApp. To see a full list of the processing steps run the following command:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2024-05-10 16:01:14,760 - isce.insar - INFO - ISCE VERSION = 2.6.3, RELEASE_SVN_REVISION = ,RELEASE_DATE = 20230418, CURRENT_SVN_REVISION = \n",
      "ISCE VERSION = 2.6.3, RELEASE_SVN_REVISION = ,RELEASE_DATE = 20230418, CURRENT_SVN_REVISION = \n",
      "\n",
      "    Insar Application:\n",
      "    Implements InSAR processing flow for a pair of scenes from\n",
      "    sensor raw data to geocoded, flattened interferograms.\n",
      "    \n",
      "The currently supported sensors are:  ['ALOS', 'ALOS2', 'ALOS_SLC', 'COSMO_SKYMED', 'COSMO_SKYMED_SLC', 'ENVISAT', 'ENVISAT_SLC', 'ERS', 'ERS_ENVISAT', 'ERS_ENVISAT_SLC', 'ERS_SLC', 'ICEYE_SLC', 'KOMPSAT5', 'RADARSAT1', 'RADARSAT2', 'RISAT1', 'RISAT1_SLC', 'ROI_PAC', 'SAOCOM_SLC', 'SENTINEL1', 'SICD_RGZERO', 'TERRASARX', 'UAVSAR_HDF5_SLC', 'UAVSAR_POLSAR', 'UAVSAR_RPI', 'UAVSAR_STACK']\n",
      "\n",
      "Usages: \n",
      "stripmapApp.py <input-file.xml>\n",
      "stripmapApp.py --steps\n",
      "stripmapApp.py --help\n",
      "stripmapApp.py --help --steps\n",
      "\n",
      "\n",
      "    Insar Application:\n",
      "    Implements InSAR processing flow for a pair of scenes from\n",
      "    sensor raw data to geocoded, flattened interferograms.\n",
      "    \n",
      "A description of the individual steps can be found in the README file\n",
      "and also in the ISCE.pdf document\n",
      "\n",
      "Command line options for steps processing are formed by\n",
      "combining the following three options as required:\n",
      "\n",
      "'--start=<step>', '--end=<step>', '--dostep=<step>'\n",
      "\n",
      "The step names are chosen from the following list:\n",
      "\n",
      "['startup', 'preprocess', 'cropraw', 'formslc', 'cropslc']\n",
      "['verifyDEM', 'topo', 'geo2rdr', 'coarse_resample', 'misregistration']\n",
      "['refined_resample', 'dense_offsets', 'rubber_sheet_range', 'rubber_sheet_azimuth', 'fine_resample']\n",
      "['split_range_spectrum', 'sub_band_resample', 'interferogram', 'sub_band_interferogram', 'filter']\n",
      "['filter_low_band', 'filter_high_band', 'unwrap', 'unwrap_low_band', 'unwrap_high_band']\n",
      "['ionosphere', 'geocode', 'geocodeoffsets', 'endup']\n",
      "\n",
      "If --start is missing, then processing starts at the first step.\n",
      "If --end is missing, then processing ends at the final step.\n",
      "If --dostep is used, then only the named step is processed.\n",
      "\n",
      "In order to use either --start or --dostep, it is necessary that a\n",
      "previous run was done using one of the steps options to process at least\n",
      "through the step immediately preceding the starting step of the current run.\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# !export PATH=$ISCE_STACK/stripmapStack:$PATH\n",
    "# add /stripmapStack to ISCE_STACK and prepend it to PATH\n",
    "current_path = os.environ.get('PATH')\n",
    "isce_stack = os.environ.get('ISCE_STACK')\n",
    "os.environ['PATH'] = f\"{isce_stack}/stripmapStack:{current_path}\"\n",
    "current_path = os.environ.get('PATH')\n",
    "os.environ['PATH'] = f\"{current_path}:/work2/05861/tg851601/stampede2/code/rsmas_insar/tools/miniconda3/lib/python3.11/site-packages/isce/applications\"\n",
    "!stripmapApp.py --help --steps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/work2/05861/tg851601/stampede2/code/rsmas_insar/tools/miniconda3/share/isce2/stripmapStack\n",
      "/work2/05861/tg851601/stampede2/code/rsmas_insar/tools/miniconda3/bin\n",
      "/work2/05861/tg851601/stampede2/code/rsmas_insar/tools/miniconda3/condabin\n",
      "/home1/05861/tg851601/.vscode-server/cli/servers/Stable-b58957e67ee1e712cebf466b995adf4c5307b2bd/server/bin/remote-cli\n",
      "/usr/local/bin\n",
      "/usr/bin\n",
      "/usr/local/sbin\n",
      "/usr/sbin\n",
      "/work2/05861/tg851601/stampede2/code/rsmas_insar/tools/miniconda3/lib/python3.11/site-packages/isce/applications\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "'/scratch/05861/tg851601/Alos/Hawaii_ALOS1'"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# display Path\n",
    "path = os.environ.get('PATH')\n",
    "paths = path.split(os.pathsep)\n",
    "for p in paths:\n",
    "    print(p)\n",
    "os.getcwd()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# stripmapApp processing steps"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The default setting of this App includes the following steps to generate a geocoded interferogram from raw data or SLC images:\n",
    "\n",
    "'startup', 'preprocess', \n",
    "'formslc',\n",
    "'verifyDEM', \n",
    "'topo', \n",
    "'geo2rdr', \n",
    "'coarse_resample', \n",
    "'misregistration', \n",
    "'refined_resample', \n",
    "'interferogram', \n",
    "'filter', \n",
    "'unwrap', \n",
    "'geocode'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br>\n",
    "<div class=\"alert alert-info\">\n",
    "<b>Note (to process the interferogram with one command):</b> \n",
    "\n",
    "stripmapApp.py stripmapApp.xml --start=startup --end=endup\n",
    "\n",
    "</div>\n",
    "\n",
    "However in this tutorial we process the interferogram step by step.\n",
    "\n",
    "<br>\n",
    "<div class=\"alert alert-info\">\n",
    "<b>At the end of each step, you will see a mesage showing the remaining steps:</b> \n",
    "\n",
    "The remaining steps are (in order):  [.....]\n",
    "</div>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### pre-processing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "!stripmapApp.py stripmapApp.xml --start=startup --end=preprocess"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By the end of \"preprocess\", the following folders are created:\n",
    "\n",
    "20110119_raw\n",
    "\n",
    "20110306_raw\n",
    "\n",
    "If you look into one of these folders:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "os.chdir(PROCESS_DIR)\n",
    "!ls 20110119_raw"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "20110119.raw contains the raw data (I/Q real and imaginary parts of each pulse, sampled along track (azimuth direction) with Pulse Repitition Frequency (PRF) and across track(range direction) with Range Sampling Frequency. stripmapApp currently only handles data acquired (or resampled) to a constant PRF. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### crop raw data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "!stripmapApp.py stripmapApp.xml --start=cropraw --end=cropraw"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The \"cropraw\" step would crop the raw data based on the region of interest if it was requested in the stripmapApp.xml. The region of interest can be added to stripmapApp.xml as:\n",
    "```xml\n",
    "<property name=\"regionOfInterest\">[19.0, 19.9, -155.4, -154.7]</property>\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since we have not sopecified the region of interest, then \"cropraw\" will be ignored and the whole frame will be processed."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### focusing "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "!stripmapApp.py stripmapApp.xml --start=formslc --end=formslc"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By the end of \"formslc\", the raw data for both reference and secondary images are focused to SLC images."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "!ls 20110119_slc"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "20110119.slc: Single Look Comlex image for 20110119 acquisition. \n",
    "\n",
    "20110119.slc.vrt: A gdal VRT file which contains the size, data type, etc.\n",
    "\n",
    "20110119.slc.xml: ISCE xml metadat file"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order to see the number of lines and pixels for an SLC image (or any data readable by GDAL):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "os.chdir(PROCESS_DIR)\n",
    "!gdalinfo 20110119_slc/20110119.slc"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Display a subset of SLC's amplitude and phase"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "\n",
    "ds = gdal.Open(\"20110119_slc/20110119.slc\", gdal.GA_ReadOnly)\n",
    "# extract a part of the SLC to display\n",
    "x0 = 0\n",
    "y0 = 10000\n",
    "x_offset = 5000\n",
    "y_offset = 10000\n",
    "slc = ds.GetRasterBand(1).ReadAsArray(x0, y0, x_offset, y_offset)\n",
    "ds = None\n",
    "\n",
    "fig = plt.figure(figsize=(14, 12))\n",
    "\n",
    "# display amplitude of the slc\n",
    "ax = fig.add_subplot(1,2,1)\n",
    "ax.imshow(np.abs(slc), vmin = -2, vmax=2, cmap='gray')\n",
    "ax.set_title(\"amplitude\")\n",
    "\n",
    "#display phase of the slc\n",
    "ax = fig.add_subplot(1,2,2)\n",
    "ax.imshow(np.angle(slc))\n",
    "ax.set_title(\"phase\")\n",
    "\n",
    "plt.show()\n",
    "\n",
    "slc = None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### crop SLC"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "!stripmapApp.py stripmapApp.xml  --start=cropslc --end=cropslc"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Similar to crop raw data but for SLC. Since region of interest has not been specified, the whole frame is processed."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### verifyDEM"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Cheks if the DEM was given in the input xml file. If the DEM is not given, then the app downloads SRTM DEM."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "!stripmapApp.py stripmapApp.xml  --start=verifyDEM --end=verifyDEM"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### topo (mapping radar coordinates to geo coordinates)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "!stripmapApp.py stripmapApp.xml  --start=topo --end=topo"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "At this step, based on the SAR acquisition geometry of the reference Image (including Doppler information), platforms trajectory and an existing DEM, each pixel of the reference image is geolocated. The geolocated coordinates will be at the same coordinate system of the platforms state vectors, which are usually given in WGS84 coordinate system. Moreover the incidence angle and heading angles will be computed for each pixel. \n",
    "\n",
    "![title](notebook_docs/Topo.png)\n",
    "\n",
    "Outputs of the step \"topo\" are written to \"geometry\" directory:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "\n",
    "os.chdir(PROCESS_DIR)\n",
    "!ls geometry"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "lat.rdr.full: latitude of each pixel on the ground. \"full\" stands for full SAR image resolution grid (before multi-looking)\n",
    "\n",
    "lon.rdr.full: longitude\n",
    "\n",
    "z.rdr.full: height\n",
    "\n",
    "los.rdr.full: incidence angle and heading angle"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Read a bounding box of latitude\n",
    "ds = gdal.Open('geometry/lat.rdr.full', gdal.GA_ReadOnly)\n",
    "lat = ds.GetRasterBand(1).ReadAsArray(0,10000,5000, 5000)\n",
    "ds = None\n",
    "\n",
    "# Read a bounding box of longitude\n",
    "ds = gdal.Open('geometry/lon.rdr.full', gdal.GA_ReadOnly)\n",
    "lon = ds.GetRasterBand(1).ReadAsArray(0,10000,5000, 5000)\n",
    "ds = None\n",
    "\n",
    "# Read a bounding box of height\n",
    "ds = gdal.Open('geometry/z.rdr.full', gdal.GA_ReadOnly)\n",
    "hgt = ds.GetRasterBand(1).ReadAsArray(0,10000,5000, 5000)\n",
    "ds = None\n",
    "\n",
    "fig = plt.figure(figsize=(18, 16))\n",
    "\n",
    "ax = fig.add_subplot(1,3,1)\n",
    "cax=ax.imshow(lat)\n",
    "ax.set_title(\"latitude\")\n",
    "ax.set_axis_off()\n",
    "cbar = fig.colorbar(cax, orientation='horizontal')\n",
    "\n",
    "\n",
    "ax = fig.add_subplot(1,3,2)\n",
    "cax=ax.imshow(lon)\n",
    "ax.set_title(\"longitude\")\n",
    "ax.set_axis_off()\n",
    "cbar = fig.colorbar(cax, orientation='horizontal')\n",
    "\n",
    "\n",
    "ax = fig.add_subplot(1,3,3)\n",
    "cax=ax.imshow(hgt, vmin = -100, vmax=1000)\n",
    "ax.set_title(\"height\")\n",
    "ax.set_axis_off()\n",
    "cbar = fig.colorbar(cax, orientation='horizontal')\n",
    "\n",
    "plt.show()\n",
    "\n",
    "lat = None\n",
    "lon = None\n",
    "hgt = None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### geo2rdr (mapping from geo coordinates to radar coordinates)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "!stripmapApp.py stripmapApp.xml --start=geo2rdr --end=geo2rdr"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "At this step, given the geo-ccordinates of each pixel in the reference image (outputs of topo), the range and azimuth time (radar coordinates) is computed given the acquisition geometry and orbit information of the secondary image.  \n",
    "\n",
    "![title](notebook_docs/Geo2rdr.png)\n",
    "\n",
    "The computed range and azimuth time for the secondary image, gives the pure geometrical offset, required for resampling the secondary image to the reference image in the next step.\n",
    "\n",
    "![title](notebook_docs/deltaR.png)\n",
    "\n",
    "After running this step, the geometrical offsets are available in \"offsets\" folder:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "os.chdir(PROCESS_DIR)\n",
    "!ls offsets"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "azimuth.off: contains the offsets betwen reference ans secondary images in azimuth direction\n",
    "\n",
    "range.off:   contains the offsets betwen reference ans secondary images in range direction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "from osgeo import gdal\n",
    "#import gdal\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "ds = gdal.Open('offsets/azimuth.off', gdal.GA_ReadOnly)\n",
    "# extract only part of the data to display\n",
    "az_offsets = ds.GetRasterBand(1).ReadAsArray(100,100,5000,5000)\n",
    "ds = None\n",
    "\n",
    "ds = gdal.Open('offsets/range.off', gdal.GA_ReadOnly)\n",
    "# extract only part of the data to display\n",
    "rng_offsets = ds.GetRasterBand(1).ReadAsArray(100,100,5000,5000)\n",
    "ds = None\n",
    "\n",
    "fig = plt.figure(figsize=(14, 12))\n",
    "\n",
    "ax = fig.add_subplot(1,2,1)\n",
    "cax=ax.imshow(az_offsets)\n",
    "ax.set_title(\"azimuth offsets\")\n",
    "ax.set_axis_off()\n",
    "cbar = fig.colorbar(cax, orientation='horizontal')\n",
    "\n",
    "ax = fig.add_subplot(1,2,2)\n",
    "cax = ax.imshow(rng_offsets)\n",
    "ax.set_title(\"range offsets\")\n",
    "ax.set_axis_off()\n",
    "cbar = fig.colorbar(cax, orientation='horizontal')\n",
    "\n",
    "plt.show()\n",
    "\n",
    "az_offsets = None\n",
    "rng_offsets = None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### resampling (using only geometrical offsets)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "!stripmapApp.py stripmapApp.xml --start=coarse_resample --end=coarse_resample"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "At this step, the gemetrical offsets are used to resample the secondary image to the same grid as the reference image, i.e., the secondary image is co-registered to the reference image. The output of this step is written to \"coregisteredSlc\" folder."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "!ls coregisteredSlc/"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "coarse_coreg.slc: is the secondary SLC coregistered to the reference image"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "os.chdir(PROCESS_DIR)\n",
    "from osgeo import gdal\n",
    "#import gdal\n",
    "\n",
    "ds = gdal.Open(\"coregisteredSlc/coarse_coreg.slc\", gdal.GA_ReadOnly)\n",
    "slc = ds.GetRasterBand(1).ReadAsArray(0, 10000, 5000, 10000)\n",
    "ds = None\n",
    "\n",
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(1,2,1)\n",
    "ax.imshow(np.abs(slc), vmin = -2, vmax=2, cmap='gray')\n",
    "ax.set_title(\"amplitude\")\n",
    "\n",
    "slc = None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### misregistration (estimating constant offsets in range and azimuth directions)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "!stripmapApp.py stripmapApp.xml --start=misregistration --end=misregistration"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The range and azimuth offsets derived from pure geometry can be potentially affected by inaccuracy of orbit information or inaccurate DEMs, or inaccurate SAR metadata. The current available DEMs (e.g., SRTM DEMs) are accurate enough to estimate offsets with accuracies of 1/100 of a pixel. The Orbit information of most modern SAR sensors are also precise enough to obtain the same order of accuracy. However, inaccurate metadata (such as timing error, constant range bias), or range bulk delay may affect the estimated offsets. To account for such sources of errors the misregistration step is performed to estimate possible constant offsets between coarse coregistered SLC and reference SLC. For this purpose an incoherent cross correlation is performed. \n",
    "\n",
    "The results of the \"misregistration\" step is written to the \"misreg\" folder.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "os.chdir(PROCESS_DIR)\n",
    "!ls misreg/"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order to extract the estimated misregistration offsets:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "import isce\n",
    "import isceobj.StripmapProc.StripmapProc as St\n",
    "\n",
    "stObj=St()\n",
    "stObj.configure()\n",
    "\n",
    "az = stObj.loadProduct(\"misreg/misreg_az.xml\")\n",
    "rng = stObj.loadProduct(\"misreg/misreg_rg.xml\")\n",
    "\n",
    "print(\"azimuth misregistration: \", az._coeffs)\n",
    "print(\"range misregistration: \", rng._coeffs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### refine_resample (resampling using geometrical offsets + misregistration)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "!stripmapApp.py stripmapApp.xml --start=refined_resample --end=refined_resample"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "At this step resampling is re-run to account for the misregistration estimated at the previous step. The new coregisterd SLC (named refined_coreg.slc) is written to the \"coregisteredSlc\" folder."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "os.chdir(PROCESS_DIR)\n",
    "!ls coregisteredSlc/"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### optional steps ('dense_offsets', 'rubber_sheet', 'fine_resample', 'split_range_spectrum' , 'sub_band_resample')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "!stripmapApp.py stripmapApp.xml --start=dense_offsets --end=sub_band_resample"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These steps are optional and will be skipped if user does not request them in the input xml file. We will get back to these steps in a different session where we estimate ionospheric phase."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### interferogram"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "!stripmapApp.py stripmapApp.xml --start=interferogram --end=interferogram"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "At this step the reference image and refined_coreg.slc is used to generate the interferogram. The generated interferogram is multi-looked based on the user inputs in the input xml file. If user does not specify the number of looks in range and azimuth directions, then they will be estimated based on posting. The default posting is 30 m which can be also specified in the input xml file.\n",
    "\n",
    "The results of the interferogram step is written to the \"interferogram\" folder:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "os.chdir(PROCESS_DIR)\n",
    "!ls interferogram/"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "topophase.flat: flattened (geometrical phase removed) and multi-looked interferogram.(one band complex64 data).\n",
    "\n",
    "topophase.cor: coherence and magnitude for the flattened multi-looked interferogram. (two bands float32 data).\n",
    "\n",
    "topophase.cor.full: similar to topophase.cor but at full SAR resolution.\n",
    "\n",
    "topophase.amp: amplitudes of reference amd secondary images. (two bands float32) \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### sub-band interferogram"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "!stripmapApp.py stripmapApp.xml --start=sub_band_interferogram --end=sub_band_interferogram"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This step will be skipped as we have not asked for ionospheric phase estimation. We will get back to this step in the ionospheric phase estimation notebook."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### filter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "!stripmapApp.py stripmapApp.xml --start=filter --end=filter"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A power spectral filter is applied to the multi-looked interferogram to reduce noise.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "from osgeo import gdal\n",
    "# import gdal\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "# reading the multi-looked wrapped interferogram\n",
    "ds = gdal.Open(\"interferogram/topophase.flat\", gdal.GA_ReadOnly)\n",
    "igram = ds.GetRasterBand(1).ReadAsArray()\n",
    "ds = None\n",
    "\n",
    "# reading the multi-looked un-wrapped interferogram\n",
    "ds = gdal.Open(\"interferogram/filt_topophase.flat\", gdal.GA_ReadOnly)\n",
    "filt_igram = ds.GetRasterBand(1).ReadAsArray()\n",
    "ds = None\n",
    "\n",
    "fig = plt.figure(figsize=(18, 16))\n",
    "\n",
    "ax = fig.add_subplot(1,3,1)\n",
    "ax.imshow(np.abs(igram), vmin = 0 , vmax = 60.0, cmap = 'gray')\n",
    "ax.set_title(\"magnitude\")\n",
    "ax.set_axis_off()\n",
    "\n",
    "ax = fig.add_subplot(1,3,2)\n",
    "ax.imshow(np.angle(igram), cmap='jet')\n",
    "ax.plot([1000,2800,2800,1000,1000],[3000,3000,2000,2000,3000],'-k')\n",
    "ax.set_title(\"multi-looked interferometric phase\")\n",
    "ax.set_axis_off()\n",
    "\n",
    "ax = fig.add_subplot(1,3,3)\n",
    "ax.imshow(np.angle(filt_igram), cmap='jet')\n",
    "ax.plot([1000,2800,2800,1000,1000],[3000,3000,2000,2000,3000],'-k')\n",
    "ax.set_title(\"multi-looked & filtered phase\")\n",
    "#ax.set_axis_off()\n",
    "\n",
    "fig = plt.figure(figsize=(18, 16))\n",
    "\n",
    "ax = fig.add_subplot(1,3,1)\n",
    "ax.imshow(np.abs(igram[2000:3000, 1000:2800]), vmin = 0 , vmax = 60.0, cmap = 'gray')\n",
    "ax.set_title(\"magnitude\")\n",
    "ax.set_axis_off()\n",
    "\n",
    "ax = fig.add_subplot(1,3,2)\n",
    "ax.imshow(np.angle(igram[2000:3000, 1000:2800]), cmap='jet')\n",
    "ax.plot([600,1200,1200,600,600],[600,600,100,100,600],'--k')\n",
    "ax.set_title(\"multi-looked interferometric phase\")\n",
    "ax.set_axis_off()\n",
    "\n",
    "ax = fig.add_subplot(1,3,3)\n",
    "ax.imshow(np.angle(filt_igram[2000:3000, 1000:2800]), cmap='jet')\n",
    "ax.plot([600,1200,1200,600,600],[600,600,100,100,600],'--k')\n",
    "ax.set_title(\"multi-looked & filtered phase\")\n",
    "ax.set_axis_off()\n",
    "\n",
    "fig = plt.figure(figsize=(18, 16))\n",
    "\n",
    "ax = fig.add_subplot(1,3,1)\n",
    "ax.imshow(np.abs(igram[2100:2600, 1600:2200]), vmin = 0 , vmax = 60.0, cmap = 'gray')\n",
    "ax.set_title(\"magnitude\")\n",
    "ax.set_axis_off()\n",
    "\n",
    "ax = fig.add_subplot(1,3,2)\n",
    "ax.imshow(np.angle(igram[2100:2600, 1600:2200]), cmap='jet')\n",
    "ax.set_title(\"multi-looked interferometric phase\")\n",
    "ax.set_axis_off()\n",
    "\n",
    "ax = fig.add_subplot(1,3,3)\n",
    "ax.imshow(np.angle(filt_igram[2100:2600, 1600:2200]), cmap='jet')\n",
    "ax.set_title(\"multi-looked & filtered phase\")\n",
    "ax.set_axis_off()\n",
    "\n",
    "filt_igram = None\n",
    "igram = None\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br>\n",
    "<div class=\"alert alert-info\">\n",
    "<b>Note :</b> \n",
    "The interferometric phase, shows ground displacement caused by dike opening event in March 2011, along the east rift zone of Kīlauea Volcano, Hawaii.\n",
    "\n",
    "</div>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### optional steps ('filter_low_band', 'filter_high_band')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "!stripmapApp.py stripmapApp.xml --start=filter_low_band  --end=filter_high_band"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These steps will be skipped since we have not asked for ionospheric phase estimation in the input xml file."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### unwrap"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "!stripmapApp.py stripmapApp.xml --start=unwrap  --end=unwrap"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "At this step the wrapped phase of the filtered and multi-looked interferogram is unwrapped. The unwrapped interferogram is a two band data with magnitude and phase components."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "from osgeo import gdal\n",
    "# import gdal\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# reading the multi-looked wrapped interferogram\n",
    "ds = gdal.Open(\"interferogram/filt_topophase.flat\", gdal.GA_ReadOnly)\n",
    "igram = ds.GetRasterBand(1).ReadAsArray()\n",
    "ds = None\n",
    "\n",
    "# reading the multi-looked unwrapped interferogram\n",
    "ds = gdal.Open(\"interferogram/filt_topophase.unw\", gdal.GA_ReadOnly)\n",
    "igram_unw = ds.GetRasterBand(2).ReadAsArray()\n",
    "ds = None\n",
    "\n",
    "# reading the connected component file\n",
    "ds = gdal.Open(\"interferogram/filt_topophase.conncomp\", gdal.GA_ReadOnly)\n",
    "connected_components = ds.GetRasterBand(1).ReadAsArray()\n",
    "ds = None\n",
    "\n",
    "\n",
    "fig = plt.figure(figsize=(18, 16))\n",
    "\n",
    "ax = fig.add_subplot(1,3,1)\n",
    "cax=ax.imshow(np.angle(igram), cmap='jet')\n",
    "ax.set_title(\"wrapped\")\n",
    "#ax.set_axis_off()\n",
    "cbar = fig.colorbar(cax, ticks=[-3.14,0,3.14],orientation='horizontal')\n",
    "cbar.ax.set_xticklabels([\"$-\\pi$\",0,\"$\\pi$\"])\n",
    "\n",
    "ax = fig.add_subplot(1,3,2)\n",
    "cax = ax.imshow(igram_unw, vmin = -15 , vmax = 15.0, cmap = 'jet')\n",
    "ax.set_title(\"unwrapped\")\n",
    "ax.set_axis_off()\n",
    "cbar = fig.colorbar(cax, ticks=[-15,0, 15], orientation='horizontal')\n",
    "\n",
    "\n",
    "ax = fig.add_subplot(1,3,3)\n",
    "cax = ax.imshow(connected_components, cmap = 'jet')\n",
    "ax.set_title(\"components\")\n",
    "ax.set_axis_off()\n",
    "cbar = fig.colorbar(cax, ticks=[0, 1] , orientation='horizontal')\n",
    "cbar.ax.set_xticklabels([0,1])\n",
    "\n",
    "\n",
    "connected_components = None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br>\n",
    "<div class=\"alert alert-info\">\n",
    "<b>Note (wrapped vs unwrapped) :</b> \n",
    "Note the colorscale for the wrapped and unwrapped interferograms. The wrapped interferometric phase varies from $-\\pi$ to $\\pi$, while the unwrapped interferogram varies from -15 to 15 radians.\n",
    "</div>\n",
    "\n",
    "<br>\n",
    "<div class=\"alert alert-info\">\n",
    "<b>Note :</b> \n",
    "The connected components file is a product of the phase unwrapping. Each interferogram may have several connected compoenets. The unwrapped phase within each component is expected to be correctly unwrapped. However, there might be $2\\pi$ phase jumps between the components. Advanced ISCE users may use the 2-stage unwrapping to adjust ambiguities among different components. stripmapApp currently does not support 2-stage unwrapping. Look for this option in future releases.  \n",
    "</div>\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "\n",
    "profile_wrapped_1 = np.angle(igram[2400,1500:1650])\n",
    "profile_unwrapped_1 = igram_unw[2400,1500:1650]\n",
    "profile_wrapped_2 = np.angle(igram[2400,1500:2000])\n",
    "profile_unwrapped_2 = igram_unw[2400,1500:2000]\n",
    "\n",
    "fig = plt.figure(figsize=(20,8))\n",
    "\n",
    "ax = fig.add_subplot(2,3,1)\n",
    "cax=ax.plot(profile_wrapped_1)\n",
    "ax.set_title(\"wrapped\")\n",
    "\n",
    "ax = fig.add_subplot(2,3,2)\n",
    "cax=ax.plot(profile_unwrapped_1)\n",
    "ax.set_title(\"unwrapped\")\n",
    "\n",
    "ax = fig.add_subplot(2,3,3)\n",
    "cax=ax.plot(np.round((profile_unwrapped_1-profile_wrapped_1)/2.0/np.pi))\n",
    "ax.set_title(\"(unwrapped - wrapped)/(2$\\pi$)\")\n",
    "\n",
    "ax = fig.add_subplot(2,3,4)\n",
    "cax=ax.plot(profile_wrapped_2)\n",
    "ax.set_title(\"wrapped\")\n",
    "\n",
    "ax = fig.add_subplot(2,3,5)\n",
    "cax=ax.plot(profile_unwrapped_2)\n",
    "ax.set_title(\"unwrapped\")\n",
    "\n",
    "ax = fig.add_subplot(2,3,6)\n",
    "cax=ax.plot((profile_unwrapped_2-profile_wrapped_2)/2.0/np.pi)\n",
    "ax.set_title(\"(unwrapped - wrapped)/(2$\\pi$)\")\n",
    "\n",
    "\n",
    "igram = None\n",
    "igram_unw = None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### optional steps ('unwrap_low_band', 'unwrap_high_band', 'ionosphere')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "!stripmapApp.py stripmapApp.xml --start=unwrap_low_band  --end=ionosphere"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since we have not asked for ionospheric phase estimation, all these steps will be skipped."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### geocoding"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "!stripmapApp.py stripmapApp.xml --start=geocode  --end=geocode"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "from osgeo import gdal\n",
    "# import gdal\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# reading the multi-looked wrapped interferogram\n",
    "ds = gdal.Open(\"interferogram/filt_topophase.unw.geo\", gdal.GA_ReadOnly)\n",
    "unw_geocoded = ds.GetRasterBand(2).ReadAsArray()\n",
    "ds = None\n",
    "\n",
    "\n",
    "fig = plt.figure(figsize=(14,12))\n",
    "\n",
    "ax = fig.add_subplot(1,1,1)\n",
    "cax = ax.imshow(unw_geocoded, vmin = -15 , vmax = 15.0, cmap = 'jet')\n",
    "ax.set_title(\"geocoded unwrapped\")\n",
    "ax.set_axis_off()\n",
    "cbar = fig.colorbar(cax, ticks=[-15,0, 15], orientation='horizontal')\n",
    "\n",
    "plt.show()\n",
    "unw_geocoded = None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Supplementary Information"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "### understanding xml files"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "The format of this type of file may seem unfamiliar or strange to you, but with the following description of the basics of the format, it will hopefully become more familiar.   The first thing to point out is that the indentations and line breaks seen above are not required and are simply used to make the structure more clear and the file more readable to humans.  The xml file provides structure to data for consumption by a computer.  As far as the computer is concerned the data structure is equally readable if all of the information were contained on a single very long line, but human readers would have a hard time reading it in that format. \n",
    "\n",
    "The next thing to point out is the method by which the data are structured through the use of tags and attributes.  An item enclosed in the < (less-than) and > (greater-than) symbols is referred to as a tag.  The name enclosed in the < and > symbols is the name of the tag.  Every tag in an xml file must have an associated closing tag that contains the same name but starts with the symbol </ and ends with the symbol >.  This is the basic unit of structure given to the data.  Data are enclosed inside of opening and closing tags that have names identifying the enclosed data.  This structure is nested to any order of nesting necessary to represent the data.  The Python language (in which the ISCE user interface is written) provides powerful tools to parse the xml structure into a data structure object and to very easily “walk” through the structure of that object.  \n",
    "\n",
    "In the above xml file the first and last tags in the file are a tag pair: <stripmapApp> and </stripmapApp> (note again, tags must come in pairs like this).  The first of these two tags, or the opening tag,  marks the beginning of the contents of the tag and the second of these two tags, or the closing tag, marks the end of the contents of the tag.  ISCE expects a “file tag” of this nature to bracket all inputs contained in the file.  The actual name of the file tag, as far as ISCE is concerned, is user selectable.  In this example it is used, as a convenience to the user, to document the ISCE application, named insarApp.py, for which it is meant to provide inputs; it could have been named <foo> and insarApp.py would have been equally happy provided that the closing tag were </foo>.  \n",
    "\n",
    "The next tag  is <component name=\"insar\">.  Its closing tag </component> is located at the penultimate line of the file (one line above the </insar> tag). The name of this tag is component and it has an attribute called name with value “insarApp”. The component tags bound a collection of information that is used by a computational element within ISCE that has the name specified by the name attribute.   The name “insarApp”  in the first component tag tells ISCE that the enclosed information correspond to a functional component in ISCE named “insarApp”, which in this case is actually the application that is run at the command line.  \n",
    "\n",
    "In general, component tags contain information in the form of other component tags or property tags, all of which can be nested to any required level. In this example the  insarApp component contains a property tag and two other component tags.\n",
    "\n",
    "The first tag we see in the insarApp component tag is the property tag with attribute name=“sensor name”. The property tag contains a value tag that contains the name of the sensor, ALOS in this case.  The next tag is a component tag with attribute name=”reference”. This tag contains a catalog tag containing  refernce.xml.  The catalog tag in general informs ISCE to look in the named file (reference.xml in this case) for the contents of the current tag.  The next component tag has the same structure with the catalog  tag containing a different file named secondary.xml."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "### Extra configuration parameters "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "The input configuration file in this tutorial only included mandatory parameters including the reference and secondary images, which are enough to run the application. This means that the application is configured with default parameters hardwired in the code or computed during processing. \n",
    "For custom processing, user may want to set parameters in the input configuration file. In the following a few more parameters are shown that can be added to stripmapApp.xml. \n",
    "\n",
    "### regionOfInterest\n",
    "\n",
    "To specify a region of interest to process:\n",
    "\n",
    "```xml\n",
    "<property name=\"regionOfInterest\">[South, North, West, East]</property>\n",
    "```\n",
    "\n",
    "Example: \n",
    "\n",
    "```xml\n",
    "<property name=\"regionOfInterest\">[19.0, 19.9, -155.4, -154.7]</property>\n",
    "```\n",
    "\n",
    "Default: Full frame is processed.\n",
    "\n",
    "### range looks\n",
    "number of looks in range direction \n",
    "\n",
    "```xml\n",
    "<property name=\"range looks\">USER_INPUT</property>\n",
    "```\n",
    "\n",
    "Deafult: is computed based on the posting parameter.\n",
    "\n",
    "### azimuth looks\n",
    "number of looks in azimuth direction \n",
    "\n",
    "Deafult: is computed based on the posting parameter.\n",
    "\n",
    "\n",
    "### posting\n",
    "Interferogram posting in meters.\n",
    "\n",
    "```xml\n",
    "<property name=\"posting\">USER_INPUT</property>\n",
    "```\n",
    "\n",
    "Default: 30\n",
    "\n",
    "<br>\n",
    "<div class=\"alert alert-info\">\n",
    "<b>Note :</b> \n",
    "If \"range looks\" and \"azimuth looks\" have not been specified, then posting is used to compute them such that the interferogram is generated with a roughly square pixels size with each dimension close to the \"posting\" parameter.\n",
    "\n",
    "</div>\n",
    "\n",
    "### filter strength\n",
    "\n",
    "strength of the adaptive filter used for filtering the wrapped interferogram\n",
    "\n",
    "```xml\n",
    "<property name=\"filter strength\">USER_INPUT</property>\n",
    "```\n",
    "\n",
    "Default: 0.5\n",
    "\n",
    "\n",
    "### useHighResolutionDemOnly\n",
    "```xml\n",
    "<property name=\"useHighResolutionDemOnly\">True</property>\n",
    "```\n",
    "\n",
    "If True and a dem is not specified in input, it will only\n",
    "    download the SRTM highest resolution dem if it is available\n",
    "    and fill the missing portion with null values (typically -32767)\n",
    "\n",
    "Default: False\n",
    "\n",
    "### do unwrap\n",
    "\n",
    "To turn phase unwrapping off\n",
    "```xml\n",
    "<property name=\"do unwrap\">False</property>\n",
    "```\n",
    "\n",
    "Default: True\n",
    "\n",
    "\n",
    "### unwrapper name\n",
    "To choose the name of the phase unwrapping method. e.g., to choose \"snaphu\" for phase unwrapping\n",
    "```xml\n",
    "<property name=\"unwrapper name\">snaphu</property>\n",
    "```\n",
    "\n",
    "Default: \"icu\".\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "### do rubbersheeting\n",
    "\n",
    "To turn on rubbersheeting (estimating azimuth offsets caused by strong ionospheric scentilation)\n",
    "\n",
    "```xml\n",
    "<property name=\"do rubbersheeting\">True</property>\n",
    "```\n",
    "Default : False\n",
    "\n",
    "### rubber sheet SNR Threshold\n",
    "\n",
    "```xml\n",
    "<property name=\"rubber sheet SNR Threshold\">USER_INPUT</property>\n",
    "```\n",
    "If \"do rubbersheeting\" is turned on, then this values is used to mask out azimuth offsets with SNR less that the input threshold. \n",
    "\n",
    "Default: 5\n",
    "\n",
    "### rubber sheet filter size\n",
    "the size of the median filter used for filtering the azimuth offsets\n",
    "\n",
    "```xml\n",
    "<property name=\"rubber sheet filter size\">USER_INPUT</property>\n",
    "```\n",
    "\n",
    "Default: 8\n",
    "\n",
    "### do denseoffsets\n",
    "turn on the dense offsets computation from cross correlation\n",
    "\n",
    "```xml\n",
    "<property name=\"do denseoffsets\">True</property>\n",
    "```\n",
    "\n",
    "Default: False\n",
    "\n",
    "<br>\n",
    "<div class=\"alert alert-info\">\n",
    "<b>Note :</b> \n",
    "\n",
    "If \"do rubbersheeting\" is turned on, then dense offsets computation is turned on regardless of the user input for \"do denseoffsets\"\n",
    "\n",
    "</div>\n",
    "\n",
    "### setting the dense offsets parameters \n",
    "\n",
    "```xml\n",
    "<property name=\"dense window width\">USER_INPUT</property>\n",
    "<property name=\"dense window height\">USER_INPUT</property>\n",
    "<property name=\"dense search width\">USER_INPUT</property>\n",
    "<property name=\"dense search height\">USER_INPUT</property>\n",
    "<property name=\"dense skip width\">USER_INPUT</property>\n",
    "<property name=\"dense skip height\">USER_INPUT</property>\n",
    "```\n",
    "\n",
    "Default values:\n",
    "<br>\n",
    "    dense window width  = 64 \n",
    "<br>\n",
    "    dense window height = 64\n",
    "<br>\n",
    "    dense search width  = 20\n",
    "<br>\n",
    "    dense search height = 20\n",
    "<br>\n",
    "    dense skip width    = 32\n",
    "<br>\n",
    "    dense skip height   = 32\n",
    "\n",
    "\n",
    "### geocode list\n",
    "\n",
    "List of products to be geocoded.\n",
    "```xml\n",
    "<property name=\"geocode list\">\"a list of files to geocode\">\n",
    "```\n",
    "Default: multilooked, filtered wrapped and unwrapped interferograms, coherence, ionospehric phase\n",
    "\n",
    "### offset geocode list\n",
    "List of offset-specific files to geocode\n",
    "```xml\n",
    "<property name=\"offset geocode list\">\"a list of offset files to geocode\">\n",
    "```\n",
    "\n",
    "\n",
    "### do split spectrum\n",
    "\n",
    "turn on split spectrum \n",
    "\n",
    "```xml\n",
    "<property name=\"do split spectrum\">True</property>\n",
    "```\n",
    "\n",
    "Default: False\n",
    "\n",
    "### do dispersive\n",
    "turn on disperive phase estimation\n",
    "\n",
    "```xml\n",
    "<property name=\"do dispersive\">True</property>\n",
    "```\n",
    "\n",
    "Default: False\n",
    "\n",
    "<br>\n",
    "<div class=\"alert alert-info\">\n",
    "<b>Note :</b> \n",
    "\n",
    "By turning on \"do dispersive\", the user input for \"do split spectrum\" is ignored and the split spectrum will be turned on as it is needed for dispersive phase estimation. \n",
    "\n",
    "</div>\n",
    "\n",
    "\n",
    "### control the filter kernel for filtering the dispersive phase  \n",
    "```xml\n",
    "<property name=\"dispersive filter kernel x-size\">800</property>\n",
    "<property name=\"dispersive filter kernel y-size\">800</property>\n",
    "<property name=\"dispersive filter kernel sigma_x\">100</property>\n",
    "<property name=\"dispersive filter kernel sigma_y\">100</property>\n",
    "<property name=\"dispersive filter kernel rotation\">0</property>\n",
    "<property name=\"dispersive filter number of iterations\">5</property>\n",
    "<property name=\"dispersive filter mask type\">coherence</property>\n",
    "<property name=\"dispersive filter coherence threshold\">0.6</property>\n",
    "    \n",
    "```    \n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "### processing data from other stripmap sensors"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "stripmapApp.py is able to process the stripmap data from the following sensors. So far it has been sucessfully tested on the following sensors: \n",
    "    <br>\n",
    "    ALOS1 (Raw and SLC)\n",
    "    ALOS2 (SLC, one frame)\n",
    "    COSMO_SkyMed (Raw and SLC)\n",
    "    ERS\n",
    "    ENVISAT ()\n",
    "    Radarsat-1\n",
    "    Radarsat-2\n",
    "    TerraSARX\n",
    "    TanDEMX\n",
    "    Sentinel1\n",
    "    \n",
    "    \n",
    "    envisat_slc\n",
    "    \n",
    "    \n",
    "### Sample input data xml for different sensors:\n",
    "\n",
    "#### Envisat: \n",
    "```xml\n",
    "\n",
    "<component name=\"reference\">\n",
    "    <property name=\"IMAGEFILE\">data/ASA_IMS_1PNESA20050519_140259_000000172037_00239_16826_0000.N1</property>\n",
    "    <property name=\"INSTRUMENT_DIRECTORY\">/u/k-data/agram/sat_metadata/ENV/INS_DIR</property>\n",
    "    <property name=\"ORBIT_DIRECTORY\">/u/k-data/agram/sat_metadata/ENV/Doris/VOR</property>\n",
    "    <property name=\"OUTPUT\">\n",
    "        20050519\n",
    "    </property>\n",
    "</component>\n",
    "\n",
    "```\n",
    "\n",
    "<br>\n",
    "<div class=\"alert alert-info\">\n",
    "<b>Note :</b> \n",
    "Note that for processing the ENVISAT data a directory that contains the orbits is required. \n",
    "</div>\n",
    "\n",
    "\n",
    "### Sentinel-1 stripmap:\n",
    "```xml\n",
    "    <component name=\"reference\">\n",
    "      <property name=\"orbit directory\">/u/data/sat_metadata/S1/aux_poeorb/</property>\n",
    "      <property name=\"output\">20151024</property>\n",
    "      <property name=\"safe\">/u/data/S1A_S1_SLC__1SSV_20151024T234201_20151024T234230_008301_00BB43_068C.zip</property>\n",
    "    </component>\n",
    "    <component name=\"secondary\">\n",
    "      <property name=\"orbit directory\">/u/k-raw/sat_metadata/S1/aux_poeorb/</property>\n",
    "      <property name=\"output\">20150930</property>\n",
    "      <property name=\"safe\">/u/data/S1A_S1_SLC__1SSV_20150930T234200_20150930T234230_007951_00B1CC_121C.zip</property>\n",
    "    </component>\n",
    "```\n",
    "\n",
    "<br>\n",
    "<div class=\"alert alert-info\">\n",
    "<b>Note :</b> \n",
    "Note that for processing the Sentinel-1 data a directory that contains the orbits is required. \n",
    "</div>\n",
    "\n",
    "### ALOS2 SLC\n",
    "```xml\n",
    "<component>\n",
    "    <property name=\"IMAGEFILE\">\n",
    "        <value>data/20141114/ALOS2025732920-141114/IMG-HH-ALOS2025732920-141114-UBSL1.1__D</value>\n",
    "    </property>\n",
    "    <property name=\"LEADERFILE\">\n",
    "        <value>data/20141114/ALOS2025732920-141114/LED-ALOS2025732920-141114-UBSL1.1__D</value>\n",
    "    </property>\n",
    "    <property name=\"OUTPUT\">\n",
    "        <value>20141114</value>\n",
    "    </property>\n",
    "</component>\n",
    "```\n",
    "\n",
    "### ALOS1 raw data\n",
    "``` xml\n",
    "<component>\n",
    "    <property name=\"IMAGEFILE\">\n",
    "       <value>[data/20110119/ALPSRP265743230-L1.0/IMG-HH-ALPSRP265743230-H1.0__D]</value>\n",
    "    </property>\n",
    "    <property name=\"LEADERFILE\">\n",
    "        <value>[data/20110119/ALPSRP265743230-L1.0/LED-ALPSRP265743230-H1.0__D]</value>\n",
    "    </property>\n",
    "    <property name=\"OUTPUT\">\n",
    "        <value>20110119</value>\n",
    "    </property>\n",
    "</component>\n",
    "```\n",
    "\n",
    "### CosmoSkyMed raw or SLC data\n",
    "\n",
    "```xml\n",
    "<component name=\"reference\">\n",
    "    <property name=\"HDF5\">data/CSKS3_RAW_B_HI_03_HH_RD_SF_20111007021527_20111007021534.h5</property>\n",
    "    <property name=\"OUTPUT\">\n",
    "        20111007\n",
    "    </property>\n",
    "</component>\n",
    "\n",
    "```\n",
    "\n",
    "### TerraSAR-X and TanDEM-X\n",
    "\n",
    "```xml\n",
    "<component name=\"reference\">\n",
    "    <property name=\"xml\">PATH_TO_TSX_DATA_XML</property>\n",
    "    <property name=\"OUTPUT\">OUTPUT_NAME</property> \n",
    "</component>    \n",
    "```\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Using ISCE as a python library"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "ISCE can be used a python library. Users can develop their own workflows within ISCE framework. Here are few simple examples where we try to call isce modules:\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Example 1: (extract metadata, range and azimuth pixel size)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "import isce\n",
    "import isceobj\n",
    "import isceobj.StripmapProc.StripmapProc as St\n",
    "from isceobj.Planet.Planet import Planet\n",
    "\n",
    "stObj = St()\n",
    "stObj.configure()\n",
    "frame = stObj.loadProduct(\"20110119_slc.xml\")\n",
    "print(\"Wavelength = {0} m\".format(frame.radarWavelegth))\n",
    "print(\"Slant Range Pixel Size = {0} m\".format(frame.instrument.rangePixelSize))\n",
    "\n",
    "#For azimuth pixel size we need to multiply azimuth time interval by the platform velocity along the track\n",
    "\n",
    "# the acquisition time at the middle of the scene\n",
    "t_mid = frame.sensingMid\n",
    "\n",
    "#get the orbit for t_mid\n",
    "st_mid=frame.orbit.interpolateOrbit(t_mid)\n",
    "\n",
    "# platform velocity\n",
    "Vs = st_mid.getScalarVelocity()\n",
    "\n",
    "# pulse repitition frequency\n",
    "prf = frame.instrument.PRF\n",
    "\n",
    "#Azimuth time interval \n",
    "ATI = 1.0/prf\n",
    "\n",
    "#Azimuth Pixel size\n",
    "az_pixel_size = ATI*Vs\n",
    "print(\"Azimuth Pixel Size = {0} m\".format(az_pixel_size))\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Example 2: compute ground range pixels size"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "r0 = frame.startingRange\n",
    "rmax = frame.getFarRange()\n",
    "rng =(r0+rmax)/2\n",
    "\n",
    "elp = Planet(pname='Earth').ellipsoid\n",
    "tmid = frame.sensingMid\n",
    "\n",
    "sv = frame.orbit.interpolateOrbit( tmid, method='hermite') #.getPosition()\n",
    "llh = elp.xyz_to_llh(sv.getPosition())\n",
    "\n",
    "\n",
    "hdg = frame.orbit.getENUHeading(tmid)\n",
    "elp.setSCH(llh[0], llh[1], hdg)\n",
    "sch, vsch = elp.xyzdot_to_schdot(sv.getPosition(), sv.getVelocity())\n",
    "\n",
    "Re = elp.pegRadCur\n",
    "H = sch[2]\n",
    "cos_beta_e = (Re**2 + (Re + H)**2 -rng**2)/(2*Re*(Re+H))\n",
    "sin_bet_e = np.sqrt(1 - cos_beta_e**2)\n",
    "sin_theta_i = sin_bet_e*(Re + H)/rng\n",
    "print(\"incidence angle at the middle of the swath: \", np.arcsin(sin_theta_i)*180.0/np.pi)\n",
    "groundRangeRes = frame.instrument.rangePixelSize/sin_theta_i\n",
    "print(\"Ground range pixel size: {0} m \".format(groundRangeRes))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br>\n",
    "<div class=\"alert alert-info\">\n",
    "<b>Note :</b> \n",
    "One can easily get the incidence angle from the los.rdr file inside geometry folder. Even without opening the file, here is a way to get the statistics and the average value of the incidence angle:  gdalinfo geometry/los.rdr -stats\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": 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.11.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}