ARIAtools.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font size=\"1\">Copyright 2021, by the California Institute of Technology. ALL RIGHTS RESERVED. United States Government sponsorship acknowledged. Any commercial use must be negotiated with the Office of Technology Transfer at the California Institute of Technology.</font>\n",
    "    \n",
    "<font size=\"1\">This software may be subject to U.S. export control laws and regulations. By accepting this document, the user agrees to comply with all applicable U.S. export laws and regulations. User has the responsibility to obtain export licenses, or other export authority as may be required, before exporting such information to foreign countries or providing access to foreign persons.<font>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ARIAtools\n",
    "The steps in this example have been extracted from the unit tests found in ARIA-tools config.yml, found [here](https://github.com/aria-tools/ARIA-tools/blob/dev/.circleci/config.yml) and the ARIA-tools documentation, found [here](https://github.com/aria-tools/ARIA-tools-docs)\n",
    "\n",
    "#### Kernel: mintpy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Initial set-up\n",
    "\n",
    "Make sure the ~/.netrc file is updated to include your earthdata login. A line has been inserted in which you substitute your username and password:\n",
    "\n",
    "    machine urs.earthdata.nasa.gov login <earthdata_user> password <earthdata_password>\n",
    "\n",
    "This will enable the downloads from the ASF DAAC.\n",
    "\n",
    "You also need to ensure your AWS credentials are valid."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Step 1 - Make a few imports and define the function that performs the plotting"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "from osgeo import gdal, ogr\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.ticker import FuncFormatter, FormatStrFormatter, StrMethodFormatter\n",
    "from mpl_toolkits.axes_grid1 import make_axes_locatable  \n",
    "import ARIAtools.shapefile_util as shputil\n",
    "\n",
    "# this block makes sure the directory set-up/change is only done once and relative to the notebook's directory\n",
    "try:\n",
    "    start_dir\n",
    "except NameError:\n",
    "    start_dir = os.getcwd()\n",
    "    output_dir = os.path.join(start_dir, 'notebook_output/ARIAtools_dem')\n",
    "    os.makedirs(output_dir, exist_ok=True)\n",
    "    py_dir = os.path.join(start_dir, 'python')\n",
    "\n",
    "os.chdir(output_dir)\n",
    "\n",
    "\n",
    "def plot_layer(path_layer, lay_type=None, cmap=None, **kwargs):\n",
    "    \"\"\" \n",
    "        path_layers is a string to the GDAL compatible dataset to be plotted\n",
    "    \"\"\"\n",
    "    \n",
    "    if not lay_type: \n",
    "        lay_type = os.path.dirname(path_layer)\n",
    "    title = [os.path.basename(lay_type)]\n",
    "    \n",
    "    ## get the lon lat bounds\n",
    "    ds       = gdal.Open(path_layer, gdal.GA_ReadOnly)\n",
    "    trans    = ds.GetGeoTransform()\n",
    "    extent   = [trans[0], trans[0] + ds.RasterXSize * trans[1], trans[3] + ds.RasterYSize*trans[5], trans[3]]\n",
    "    \n",
    "    ## loading the data\n",
    "    n_bands  = ds.RasterCount\n",
    "    lst_arrs = []\n",
    "    \n",
    "    for band in range(n_bands):\n",
    "        raster = ds.GetRasterBand(band+1)\n",
    "        arr    = raster.ReadAsArray()\n",
    "        try:\n",
    "            NoData = raster.GetNoDataValue()\n",
    "            arr = np.ma.masked_where((arr>1e20) |(arr==NoData),arr )\n",
    "        except:\n",
    "            print('Could not find a no-data value...')\n",
    "            arr = np.ma.masked_where(arr>1e20,arr)\n",
    "        \n",
    "        lst_arrs.append(arr)\n",
    "\n",
    "    ds = None\n",
    "    if n_bands < 4:\n",
    "        nrows = 1; ncols = n_bands\n",
    "    else:\n",
    "        raise Exception('Number of bands currently unsupported')\n",
    "        \n",
    "    \n",
    "    ## initializing a figure\n",
    "    fig, axes = plt.subplots(figsize=(12,9), ncols=ncols, nrows=nrows, sharex='col', sharey='row')\n",
    "    axes = axes if isinstance(axes, np.ndarray) else np.array(axes)\n",
    "    axe  = axes.ravel() \n",
    "    cmap = plt.cm.Greys_r\n",
    "    cmap.set_under('black')\n",
    "    \n",
    "    ## definging the plotting options for differnt layer types\n",
    "    # Amplitude:\n",
    "    if lay_type.endswith('amplitude'): \n",
    "        # will fix the maximum amplitude bound\n",
    "        vmin=None\n",
    "        vmax = 2000 \n",
    "    # Coherence:\n",
    "    elif lay_type.endswith('coherence'): \n",
    "        # has fixed range between 0-1\n",
    "        vmin=0\n",
    "        vmax = 1\n",
    "    # Incidence angle:\n",
    "    elif lay_type.endswith('incidenceAngle'): \n",
    "        vmin=None\n",
    "        vmax=None\n",
    "    # water\n",
    "    elif lay_type.startswith('water'):\n",
    "        # no bounds needed will be a 0/1 mask\n",
    "        vmin=0\n",
    "        vmax=1\n",
    "        cmap='Greys'\n",
    "    # deformation or unwrapped phase\n",
    "    elif lay_type.startswith('defo'): \n",
    "        # let the data drive the bounds\n",
    "        vmin=None\n",
    "        vmax=None\n",
    "        # change colormap to a warm type\n",
    "        cmap=plt.cm.coolwarm\n",
    "    elif lay_type.startswith('terr') or lay_type.startswith('topo'): \n",
    "        # let the data drive the bounds\n",
    "        vmin=None\n",
    "        vmax=None\n",
    "        # change colormap to a warm type\n",
    "        cmap=plt.cm.terrain\n",
    "    elif lay_type == 'ENU':\n",
    "        vmin=None\n",
    "        vmax=None\n",
    "        title = ['East', 'North', 'Up']\n",
    "        fig.subplots_adjust(wspace=0.5)\n",
    "\n",
    "    else:\n",
    "        # let the data drive the bounds\n",
    "        vmin=None\n",
    "        vmax=None\n",
    "        # change colormap to a warm type\n",
    "        cmap=plt.cm.coolwarm\n",
    "        \n",
    "    # plotting the data    \n",
    "    for i, ax in enumerate(axe):\n",
    "        im   = ax.imshow(lst_arrs[i], cmap=cmap, vmin=vmin, vmax=vmax, extent=extent,interpolation='nearest')\n",
    "        divider = make_axes_locatable(ax)\n",
    "        cax     = divider.append_axes('right', size='5%', pad=0.25)\n",
    "        if lay_type == 'ENU':\n",
    "            fig.colorbar(im, cax=cax, format=FuncFormatter(lambda x, y: '{:.3f}'.format(x)))\n",
    "        elif lay_type.startswith('water'):\n",
    "            fig.colorbar(im, cax=cax, ticks=[vmin, vmax])\n",
    "        else:\n",
    "            fig.colorbar(im, cax=cax)\n",
    "\n",
    "        ax.set_title(title[i], fontsize=15)\n",
    "        ax.grid(False)\n",
    "\n",
    "    axe[0].set_ylabel('latitude', labelpad=15, fontsize=15)\n",
    "    axe[int(np.floor(n_bands/2))].set_xlabel('longitude', labelpad=15, fontsize=15)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Step 2 : Determine the product(s) from track 42 between 2018/04/02 and 2018/04/14 over the big island.\n",
    "- Output the count of products\n",
    "- Output the product URLs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "test_dir = os.path.join(output_dir, 'test2')\n",
    "prod_dir = os.path.join(test_dir, 'products')\n",
    "os.makedirs(prod_dir, exist_ok=True)\n",
    "os.chdir(test_dir)\n",
    "          \n",
    "!rm -f products/*.txt\n",
    "\n",
    "# there should be 6 products\n",
    "!ariaDownload.py -t 042 -w products --ifg 20180414_20180402 -v -o count\n",
    "\n",
    "# get the urls of the products\n",
    "!ariaDownload.py -t 042 -w products --ifg 20180414_20180402 -v -o url"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Step 3 : Download and process a product from the ASF DAAC"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "os.chdir(prod_dir)\n",
    "!echo Downloading `head -2 download_products_*_0.txt | tail -1` \n",
    "!wget --no-clobber --load-cookies /tmp/cookies.txt --save-cookies /tmp/cookies.txt --keep-session-cookies `head -2 download_products_*_0.txt | tail -1`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Step 4 : Extract the bounding box from the (first) product file\n",
    "\n",
    "Make sure your AWS credentials are valid.\n",
    "\n",
    "#### Notice that we are invoking DEM.py using the isce conda environment! "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%cd $test_dir\n",
    "\n",
    "# get the bounding box of the image\n",
    "!ariaExtract.py -f \"products/*.nc\"\n",
    "\n",
    "!cat user_bbox.json\n",
    "\n",
    "# here we do a bit of black magic to invoke DEM.py (which requires ISCE3) in the isce environment\n",
    "!conda run -n isce python $py_dir/DEM.py -b `python $py_dir/extract_bbox.py -f user_bbox.json`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Step 5 : Perform some extractions of the data; azimuthAngle and amplitude\n",
    "\n",
    "The first invocation of ariaExtract specifies the NISAR dem file to which a cutline will be applied. Subsequent calls will reuse it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%cd $test_dir\n",
    "\n",
    "# extract azimuthAngle layer\n",
    "!echo \"==========  azimuthAngle  ==========\"\n",
    "!ariaExtract.py -f \"products/*.nc\" -l azimuthAngle -d dem.tif\n",
    "\n",
    "# list the result files\n",
    "!echo \"Result files\"\n",
    "!ls -1 azimuthAngle\n",
    "\n",
    "# extract amplitude layer\n",
    "!echo \"==========    amplitude   ==========\"\n",
    "!ariaExtract.py -f \"products/*.nc\" -l amplitude\n",
    "\n",
    "# list the result files\n",
    "!echo \"Result files\"\n",
    "!ls -1 amplitude\n",
    "\n",
    "# extract coherehence layer\n",
    "!echo \"==========    coherence   ==========\"\n",
    "!ariaExtract.py -f \"products/*.nc\" -l coherence\n",
    "\n",
    "# list the result files\n",
    "!echo \"Result files\"\n",
    "!ls -1 coherence\n",
    "\n",
    "# extract unwrappedPhase layer\n",
    "!echo \"========== unwrappedPhase ==========\"\n",
    "!ariaExtract.py -f \"products/*.nc\" -l unwrappedPhase\n",
    "\n",
    "# list the result files\n",
    "!echo \"Result files\"\n",
    "!ls -1 unwrappedPhase"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Step 6 : Show the product bounding box"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!ls -1 productBoundingBox\n",
    "\n",
    "shputil.plot_shapefile(os.path.join('productBoundingBox','productBoundingBox.json'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Step 7 : Plot the azimuth extraction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_layer(os.path.join('azimuthAngle', '20180414_20180402.vrt'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Step 8 : Plot the amplitude extraction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_layer(os.path.join('amplitude', '20180414_20180402.vrt'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Step 9 : Plot the coherence"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_layer(os.path.join('coherence', '20180414_20180402.vrt'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Step 10 : Plot the unwrappedPhase"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_layer(os.path.join('unwrappedPhase', '20180414_20180402.vrt'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Step 11: Time-series prep\n",
    "TBD: Use of watermask switch (i.e. -m download) goes into seemingly non-terminating loop."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!ariaTSsetup.py -d download -f \"products/*.nc\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Step 12 : Run loading of aria products in mintpy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# bit of a kludge to identify where to find prep_aria.py\n",
    "import mintpy\n",
    "mintpy_dir = os.path.split(mintpy.__file__)[0]\n",
    "\n",
    "!python $mintpy_dir/prep_aria.py -s stack/ -d DEM/SRTM_3arcsec.dem -i incidenceAngle/20180414_20180402.vrt -a azimuthAngle/20180414_20180402.vrt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font size=\"1\">This notebook is compatible with NISAR Jupyter Server Stack v1.7.1 and above</font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "mintpy [conda env:.local-mintpy]",
   "language": "python",
   "name": "conda-env-.local-mintpy-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.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}