EarthScope2024/1.4_Stripmap_Data_Processing_Interferometry/stripmapApp.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Processing your Data with stripmapApp.py\n",
    "\n",
    "**Author**: Heresh Fattahi, Piyush Agram, updated for 2022, 2023, and 2024 Eric Fielding\n",
    "\n",
    "In this notebook, we will walk through the various steps of processing with the ISCE2 stripmapApp.py. \n",
    "\n",
    "stripmapApp.py is a pair-by-pair interferometric processor that takes as input two SAR acquisitions acquired in stripmap mode. stripmapApp.py will not work for other acquisition mode such as TOPS or ScanSAR. At this time, ISCE2's topsApp supports TOPS mode acquisitions in SLC data format from Sentinel-1 A and B. Processing of stripmap data in ISCE2 is supported from RAW or SLC level data inputs for a number of satellites.\n",
    "\n",
    "To illustrate the usage of stripApp.py, we will use a JAXA ALOS dataset capturing the surface deformation as result of the March 2011 Kamoamoa dike intrusion and fissure eruption that occurred on Kilauea in Hawaii. JAXA calls the \"raw\" data that we are processing here Level 1.0. Other satellite systems call this Level 0 data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# stripmap acquisition mode"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "jp-MarkdownHeadingCollapsed": true,
    "tags": []
   },
   "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 ISCE2 **stripmapApp.py**. \n",
    "\n",
    "The stripmap mode has been used by several 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 acquiring 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",
    "The NISAR mission will use a SAR acquisition mode called SweepSAR, which enables a wide 240 km swath with a digital beam forming method. The processing of NISAR raw (Level 0) data is complicated, but after focusing it is close to stripmap data for Level 1 SLC and Level 2 products because the squint angle does not change. NISAR data processing requires the ISCE3 package.\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 literature:\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",
    " \n",
    "Rosen, P. A., E. Gurrola, G. F. Sacco, and H. Zebker (2012), The InSAR Scientific Computing Environment, paper presented at 9th European Conference on Synthetic Aperture Radar, Nuremberg, Germany, 23-26 April.\n",
    "\n",
    "Rosen, P. A., E. M. Gurrola, P. Agram, J. Cohen, M. Lavalle, B. V. Riel, H. Fattahi, M. A. G. Aivazis, M. Simons, and S. M. Buckley (2018), The InSAR Scientific Computing Environment 3.0: A Flexible Framework for NISAR Operational and User-Led Science Processing, paper presented at IGARSS 2018 - 2018 IEEE International Geoscience and Remote Sensing Symposium, 22-27 July 2018. (doi:10.1109/IGARSS.2018.8517504)\n",
    "\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 ISCE2 application, designed for interferometric processing of SAR data acquired with stripmap mode onboard platforms with precise orbits. This ISCE2 application is equivalent to insarApp (an older ISCE2 application which was 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": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "import os\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from osgeo import gdal\n",
    "from tqdm import tqdm\n",
    "import urllib.request\n",
    "import boto3                      # For talking to s3 bucket\n",
    "from botocore import UNSIGNED # NEW Change\n",
    "from botocore.client import Config # NEW Change\n",
    "import netrc\n",
    "\n",
    "\n",
    "# create .netrc if it does not exist--adapted from ariaExtract notebook\n",
    "if not os.path.exists(os.path.expanduser('~/.netrc')):\n",
    "    print('NEEDED To Download ALOS and SRTM data: \\n Link to create account : https://urs.earthdata.nasa.gov/')\n",
    "    earthdata_user = input('Please type your Earthdata username:')\n",
    "    earthdata_user = str(earthdata_user)\n",
    "    earthdata_password = input('Please type your Earthdata password:')\n",
    "    earthdata_password = str(earthdata_password)\n",
    "    os.system('echo machine urs.earthdata.nasa.gov login \"{usern}\" password \"{passwd}\" > ~/.netrc; chmod 600 ~/.netrc'.format( \\\n",
    "              usern = earthdata_user, passwd = earthdata_password))\n",
    "#\n",
    "\n",
    "# use the netrc module to read the file\n",
    "machine_name = \"urs.earthdata.nasa.gov\"\n",
    "credentials = netrc.netrc()\n",
    "login, _, password = credentials.authenticators(machine_name)\n",
    "\n",
    "ASF_USER = login\n",
    "ASF_PASS = password\n",
    "\n",
    "# the working directory:\n",
    "home_dir = os.path.join(os.getenv(\"HOME\"), \"work\")\n",
    "PROCESS_DIR = os.path.join(home_dir, \"Hawaii_ALOS1\")\n",
    "DATA_DIR =  os.path.join(PROCESS_DIR, \"data\")\n",
    "\n",
    "if len(ASF_PASS)==0 | len(ASF_USER)==0:\n",
    "    print(\"ERROR: The ASF USER pass needs to be set either by running this cell or should be included in ~/.netrc file\") \n",
    "\n",
    "# defining backup dirs in case of download issues on the local server\n",
    "# boto3 is package to access AWS S3 bucket from Python\n",
    "# not presently used\n",
    "s3 = boto3.resource(\"s3\", config=Config(signature_version=UNSIGNED))\n",
    "data_backup_bucket = s3.Bucket(\"asf-jupyter-data-west\")\n",
    "data_backup_dir = \"stripmap_Hawaii\"\n",
    "\n",
    "# Utility to copy data from \n",
    "def copy_from_bucket(file_in_bucket, dest_file,\n",
    "                    bucket=data_backup_bucket):\n",
    "    if os.path.exists(dest_file):\n",
    "        print(\"Destination file {0} already exists. Skipping download...\".format(dest_file))\n",
    "    else:\n",
    "        bucket.download_file(file_in_bucket, dest_file)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def configure_inputs(outDir): \n",
    "\n",
    "    \"\"\"Write Configuration files for ISCE2 stripmapApp to process ALOS Level 1.0 products\"\"\"\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",
    "    \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>stripmap_Hawaii/demLat_N18_N21_Lon_W157_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": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "\n",
    "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",
    "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": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "tags": []
   },
   "outputs": [],
   "source": [
    "cmd = \"wget  https://datapool.asf.alaska.edu/L1.0/A3/ALPSRP265743230-L1.0.zip --user={0} --password='{1}'\".format(ASF_USER, ASF_PASS)\n",
    "if not os.path.exists(os.path.join(DATA_DIR, \"ALPSRP265743230-L1.0.zip\")):\n",
    "    os.system(cmd)\n",
    "else:\n",
    "    print(\"ALPSRP265743230-L1.0.zip already exists\")\n",
    "    \n",
    "cmd = \"wget  https://datapool.asf.alaska.edu/L1.0/A3/ALPSRP272453230-L1.0.zip --user={0} --password='{1}'\".format(ASF_USER, ASF_PASS)\n",
    "if not os.path.exists(os.path.join(DATA_DIR, \"ALPSRP272453230-L1.0.zip\")):\n",
    "    os.system(cmd)\n",
    "else:\n",
    "    print(\"ALPSRP272453230-L1.0.zip already exists\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now the data should be downloading from ASF. If the download did not start, then check your `.netrc` file has your EarthData login. You can also download with the SSARA command below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "ls "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "unzip the downloaded files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "if not os.path.exists(os.path.join(DATA_DIR, \"ALPSRP265743230-L1.0\")):\n",
    "    !unzip ALPSRP265743230-L1.0.zip\n",
    "if not os.path.exists(os.path.join(DATA_DIR, \"ALPSRP272453230-L1.0\")):\n",
    "    !unzip ALPSRP272453230-L1.0.zip"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " looking at the unzipped directories there are multiple files:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "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 (or two in this example of dual-polarization data) 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 (SLC) or 1.5 (ground-range detected) 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": null,
   "metadata": {
    "scrolled": true,
    "tags": []
   },
   "outputs": [],
   "source": [
    "!cat ALPSRP265743230-L1.0/ALPSRP265743230.l0.workreport"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "!grep Img_SceneCenterDateTime ALPSRP265743230-L1.0/ALPSRP265743230.l0.workreport\n",
    "!grep Img_SceneCenterDateTime ALPSRP272453230-L1.0/ALPSRP272453230.l0.workreport"
   ]
  },
  {
   "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": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "!mkdir 20110119\n",
    "!mkdir 20110306\n",
    "!mv ALPSRP265743230-L1.0 20110119\n",
    "!mv ALPSRP272453230-L1.0 20110306"
   ]
  },
  {
   "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": null,
   "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": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "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": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "configure_inputs(PROCESS_DIR)"
   ]
  },
  {
   "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>reference.xml</catalog>\n",
    "    </component>\n",
    "    <component name=\"secondary\">\n",
    "        <catalog>secondary.xml</catalog>\n",
    "    </component>\n",
    "\n",
    "    <!--  \n",
    "    <property name=\"demFilename\">\n",
    "        <value>stripmap_Hawaii/demLat_N18_N21_Lon_W157_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": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!stripmapApp.py --help --steps"
   ]
  },
  {
   "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": {
    "tags": []
   },
   "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": [
    "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. The data is still the raw data (which JAXA calls level 1.0 and other providers call level 0) but it is converted with all the required metadatato the ISCE format."
   ]
  },
  {
   "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": {
    "tags": []
   },
   "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 single-look complex (SLC) images. JAXA calls the SLC products level 1.1 and other providers call it level 1. Many satellite providers only provide SLC products, not raw data, so for many satellites, the stripmap processing starts with the SLC data."
   ]
  },
  {
   "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": [
    "!gdalinfo 20110119_slc/20110119.slc"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Display a subset of SLC's amplitude and phase. In the amplitude image, we can see the coast of the big island of Hawaii that is covered by this scene."
   ]
  },
  {
   "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": [
    "Checks if the DEM was given in the input xml file. If the DEM is not given, then the app downloads SRTM DEM from the NASA Land Processes DAAC."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "!stripmapApp.py stripmapApp.xml  --start=verifyDEM --end=verifyDEM"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This step will check the DEM file specified in the stripmapApp.xml. If no DEM file has been specified, stripmapApp.py will download the DEM on the fly, mosaic the individual DEM tiles into one file, remove the geoid, and track the filename for subsequent processing. \n",
    "\n",
    "<div class=\"alert alert-danger\">\n",
    "   \n",
    "<b>POTENTIAL ISSUE:</b> \n",
    "You did not set your earthdata credentials as detailed in the ISCE installation. Find instructions here: [Step 2 here](https://wiki.earthdata.nasa.gov/display/EL/How+To+Access+Data+With+cURL+And+Wget)\n",
    "</div>\n",
    "\n",
    "<div class=\"alert alert-danger\">\n",
    "<b>POTENTIAL ISSUE:</b> \n",
    "You did set up your earthdata credentials as detailed in the ISCE installation but you have special characters in your password. \"Escape\" these characters by adding a backslash in front of them. \n",
    "</div>\n",
    "\n",
    "<div class=\"alert alert-danger\">\n",
    "<b>POTENTIAL ISSUE:</b> \n",
    "The DEM download site is down and returns no or partial data tiles.\n",
    "</div></div>  \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### ONLY if download fails"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If your DEM download failed and you are unable to resolve the issue, you can use the DEM provided with this tutorial. This requires three steps to complete, before you can rerun the **verifyDEM** step:\n",
    "\n",
    "1) Copy over the DEM from the backup folder to your processing location"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fileZip = ['stripmap_Hawaii.zip']\n",
    "\n",
    "# zip file with the DEM files\n",
    "#'demLat_N18_N21_Lon_W157_W154.dem.wgs84',       #Binary raster\n",
    "#         'demLat_N18_N21_Lon_W157_W154.dem.wgs84.xml',   #XML needed by ISCE\n",
    "#         'demLat_N18_N21_Lon_W157_W154.dem.wgs84.vrt']   #VRT needed by ISCE\n",
    "\n",
    "for file in fileZip:\n",
    "    if not os.path.exists(os.path.join(PROCESS_DIR,file)):\n",
    "        copy_from_bucket(os.path.join(\"unavco2022\",file),\n",
    "                 os.path.join(PROCESS_DIR,file))\n",
    "        print(file + \" done\")\n",
    "    else:\n",
    "        print(file + \" already exists\")\n",
    "        \n",
    "if not os.path.exists(os.path.join(PROCESS_DIR, \"stripmap_Hawaii\")):\n",
    "    !unzip stripmap_Hawaii.zip\n",
    "else:\n",
    "    print (\"stripmap_Hawaii DEM directory already exists\")\n",
    "    !ls stripmap_Hawaii\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "2) To use this DEM, you will need to edit the **stripmapApp.xml** file and uncomment the property for the **demFilename**. You can use the terminal editor `vim` or the built-in Jupyter Lab editor.\n",
    "\n",
    "``` vim\n",
    "    vim work/Hawaii_ALOS1/stripmapApp.xml\n",
    "```\n",
    "\n",
    "Remember the *xml* guidelines:\n",
    "```xml\n",
    "<!--<property> ..COMMENTED.. </property>--> \n",
    "<property> ..UNCOMMENTED.. </property>\n",
    "```\n",
    "If the staged file name is exactly the same as what `strimapApp.py` was trying to download internally, then the edit may not be necessary, but the unpacked staged files are in a different location in this case."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Now rerun the **verifyDEM** step to ensure the new DEM information is correctly loaded into the processing:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "!stripmapApp.py --dostep=verifyDEM"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### topo (mapping radar coordinates to geo coordinates)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "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), platform's trajectory, and an existing DEM providing topography, 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": [
    "!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": {},
   "outputs": [],
   "source": [
    "!stripmapApp.py stripmapApp.xml --start=geo2rdr --end=geo2rdr"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "At this step, given the geo-coordinates 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 of the secondary image from the reference, 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": [
    "!ls offsets"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "azimuth.off: contains the geometric offsets between reference and secondary images in azimuth direction\n",
    "\n",
    "range.off:   contains the geometric offsets between reference and secondary images in range direction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "from osgeo 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": {
    "tags": []
   },
   "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 approximately 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": [
    "from osgeo 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": {
    "tags": []
   },
   "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. Coregistration can also be affected by changes in the ionosphere, especially for L-band data like this ALOS data 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 by the ionosphere 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 and the average offset over the whole scene is estimated. \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": [
    "!ls misreg/"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order to extract the estimated misregistration offsets we can use the ISCE internal objects to load the information saved in the XML files:"
   ]
  },
  {
   "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": [
    "For this pair, there is a range misregistration of 0.99 pixels or about 9 meters."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### refine_resample (resampling using geometrical offsets + misregistration)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "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": [
    "!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": {},
   "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": "markdown",
   "metadata": {},
   "source": [
    "Now we are finally ready to calculate the interferogram between the two precisely coregistered SLC images."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!stripmapApp.py stripmapApp.xml --start=interferogram --end=interferogram"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "At this step the reference SLC 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": [
    "!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": {},
   "outputs": [],
   "source": [
    "!stripmapApp.py stripmapApp.xml --start=filter --end=filter"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A power spectrum filter (Goldstein and Werner, 1998) is applied to the multi-looked interferogram to reduce noise.\n",
    "\n",
    "Now we can display the interferogram before and after the filtering and see the ground displacement signals in this interferogram."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "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. There is also some subsidence at the summit of Kīlauea due to magma withdrawal.\n",
    "<br>\n",
    "The overall ramp of fringes across the scene is due to the ionospheric variations that we will discuss later.\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": {},
   "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 in radar coordinates. The unwrapped interferogram is a two band data with magnitude and phase components."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "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": "markdown",
   "metadata": {},
   "source": [
    "We can look at profiles of the phase in the wrapped and unwrapped interferograms. The first profile is a small distance where the phase does not vary more than 2-pi. The second profile is a longer profile that shows the effects of the phase unwrapping required when the phase varies by more than 2-pi."
   ]
  },
  {
   "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": "markdown",
   "metadata": {},
   "source": [
    "Now we can geocode the interferogram and put it into the geographic coordinates that can be combined with other datasets."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "!stripmapApp.py stripmapApp.xml --start=geocode  --end=geocode"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Final Geocoded Unwrapped phase\n",
    "Now we can plot the geocoded unwrapped phase that is sometimes called a GUNW product. We will work with the GRFN ARIA-processed Sentinel-1 GUNW products later in the course."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# 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": [
    "The dike injection signal is partly hidden by the large phase ramp across the scene that is caused by ionospheric variations in this L-band pair."
   ]
  },
  {
   "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  reference.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 (Raw and SLC)\n",
    "-   Radarsat-1\n",
    "-  Radarsat-2\n",
    "-  TerraSARX\n",
    "-  TanDEMX (or PAZ)\n",
    "-  Sentinel1 (stripmap mode)\n",
    "    \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",
    "The above example is for Envisat SLC (ASA_IMS_1) images. StripmapApp also supports Envisat Raw (ASA_IM__0) images.\n",
    "\n",
    "```\n",
    "\n",
    "<br>\n",
    "<div class=\"alert alert-info\">\n",
    "<b>Note :</b> \n",
    "Note that a directory that contains the orbits and another with the Instrument calibration files is required for processing ENVISAT data. \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",
    "\n",
    "ALOS2 data is only distributed as SLC products.\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",
    "TerraSAR-X data is only distributed as SLC products\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",
    "# go to the processing directory in case not already there\n",
    "os.chdir(PROCESS_DIR)\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 in orbit = {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": [
    "# range pixel size in slant range can be converted to ground range pixel size\n",
    "# here we calculate it at the middle of the swath\n",
    "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": [
    "#### Example 3: compute ground azimuth pixels size"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# the azimuth pixel size calculated in Example 1 was at the orbit altitude\n",
    "# here we do an approximate calculation for the pixel spacing on the ground\n",
    "az_pixel_ground = az_pixel_size*Re/(Re+H)\n",
    "print(\"Azimuth Pixel Size on ground = {0} m\".format(az_pixel_ground))"
   ]
  },
  {
   "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 [conda env:.local-earthscope_insar]",
   "language": "python",
   "name": "conda-env-.local-earthscope_insar-py"
  },
  "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": 4
}