3_Access_HyP3_Data.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "ca471365-111f-4b70-8911-8d282003db90",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "# Access & Subset HyP3 SBAS Stack (InSAR or Burst-InSAR)\n",
    "\n",
    "*Author: Alex Lewandowski; Alaska Satellite Facility*\n",
    "\n",
    "This notebook assumes that you have already ordered processing for an InSAR Short Baseline Subset (SBAS) stack using [ASF HyP3](https://hyp3-docs.asf.alaska.edu/), available in a web browser at [ASF Vertex](https://search.asf.alaska.edu/) or programmatically with [hyp3-sdk](https://github.com/ASFHyP3/hyp3-sdk).\n",
    "\n",
    "---\n",
    "\n",
    "<div class=\"alert alert-info\" style=\"display: flex; align-items: center; font-family: 'Times New Roman', Times, serif; background-color: #d1ecf1;\">\n",
    "  <div style=\"display: flex; align-items: center; width: 10%;\">\n",
    "    <a href=\"https://github.com/ASFOpenSARlab/opensarlab_MintPy_Recipe_Book/issues\">\n",
    "      <img src=\"github_issues.png\" alt=\"GitHub logo over the word Issues\" style=\"width: 100%;\">\n",
    "    </a>\n",
    "  </div>\n",
    "  <div style=\"width: 95%;\">\n",
    "    <b>Did you find a bug? Do you have a feature request?</b>\n",
    "    <br/>\n",
    "    Explore GitHub Issues on this Jupyter Book's GitHub repository. Find solutions, add to the discussion, or start a new bug report or feature request: <a href=\"https://github.com/ASFOpenSARlab/opensarlab_MintPy_Recipe_Book/issues\">opensarlab_MintPy_Recipe_Book Issues</a>\n",
    "  </div>\n",
    "</div>\n",
    "\n",
    "<div class=\"alert alert-info\" style=\"display: flex; align-items: center; justify-content: space-between; font-family: 'Times New Roman', Times, serif; background-color: #d1ecf1;\">\n",
    "  <div style=\"display: flex; align-items: center; width: 10%; margin-right: 10px;\">\n",
    "    <a href=\"mailto:uso@asf.alaska.edu\">\n",
    "      <img src=\"ASF_support_logo.png\" alt=\"ASF logo\" style=\"width: 100%\">\n",
    "    </a>\n",
    "  </div>\n",
    "  <div style=\"width: 95%;\">\n",
    "    <b>Have a question related to SAR, ASF data access, or performing SBAS time series analyses with MintPy?</b>\n",
    "    <br/>\n",
    "    Contact ASF User Support: <a href=\"mailto:uso@asf.alaska.edu\">uso@asf.alaska.edu</a>\n",
    "  </div>\n",
    "</div>\n",
    "\n",
    "---\n",
    "\n",
    "## 0. Import Required Software "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "98210e5a-6f29-4a48-9021-18080b1a5761",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "from datetime import datetime\n",
    "from pathlib import Path\n",
    "import shutil\n",
    "import sys\n",
    "from tqdm.auto import tqdm\n",
    "\n",
    "import contextily as ctx\n",
    "import geopandas as gpd\n",
    "from hyp3_sdk import Batch, HyP3\n",
    "from ipyfilechooser import FileChooser\n",
    "from IPython.display import Markdown, display, HTML\n",
    "import matplotlib.pyplot as plt \n",
    "from matplotlib.lines import Line2D\n",
    "import numpy as np\n",
    "import opensarlab_lib as osl\n",
    "from osgeo import gdal, ogr\n",
    "gdal.UseExceptions()\n",
    "from shapely.geometry import box\n",
    "from rasterio.warp import transform_bounds\n",
    "\n",
    "current = Path(\"..\").resolve()\n",
    "sys.path.append(str(current))\n",
    "import util.util as util\n",
    "\n",
    "%matplotlib widget"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4453ed21-0573-4764-b53d-fe6c656be756",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "---\n",
    "\n",
    "## 1. Select or create a working directory for the analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d09e7a9e-3073-4c43-b2f1-9d5c75d407c7",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "age = osl.select_parameter(\n",
    "    [\n",
    "        \"Access a new SBAS stack\",\n",
    "        \"Add to existing SBAS stack\"\n",
    "    ]\n",
    ")\n",
    "display(age)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cc504661-8853-4904-8f86-0c87d950bfd7",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "new = 'new' in age.value\n",
    "\n",
    "if new:\n",
    "    print(f'Current working directory: {Path.cwd()}')\n",
    "    print('Create a new directory to hold your data:')\n",
    "    data_path = input(f'Enter an unused path for a new data directory:  {Path.home()}/')\n",
    "    try:\n",
    "        data_path = Path.home() / data_path.strip()\n",
    "        data_path.mkdir()\n",
    "    except:\n",
    "        raise\n",
    "else:\n",
    "    path = Path.home()\n",
    "    fc = FileChooser(path)\n",
    "    display(fc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c6168373-7bc7-4f26-91cd-e27d77672ffd",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "if not new:\n",
    "    data_path = Path.home()/fc.selected_path"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3ec94add-ab46-42a2-b784-a0a90851b779",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "---\n",
    "\n",
    "## 2. Migrate SBAS Stack from HyP3\n",
    "\n",
    "**Create a HyP3 object and authenticate**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c012ca68-6c09-46c5-94b8-32c7a633229f",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "hyp3 = HyP3(prompt=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ad54c61e-fe22-4b7e-8138-8950e1b7a29c",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "**You may search for InSAR projects in your own account or migrate data from any user's account**\n",
    "\n",
    "- Retrieving data from another user's account only requires their username and the project name.\n",
    "- It does **not** require the other user's password. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ef02d298-9f0b-4406-b960-ae53a3d6a5f2",
   "metadata": {},
   "outputs": [],
   "source": [
    "hyp3_project = osl.select_parameter(\n",
    "    [\n",
    "        'Access InSAR data with any valid HyP3 username and HyP3 Project Name',\n",
    "        'Search your Projects for available InSAR data'\n",
    "    ]\n",
    ")\n",
    "display(hyp3_project)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "666b8617-4195-4d48-bad4-3dd8e4272212",
   "metadata": {},
   "source": [
    "**Select your SBAS stack's HyP3 product type**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e0dec948-2ea6-40a2-bebe-77434774beea",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "product_select = osl.select_parameter(\n",
    "    [\n",
    "        'INSAR_GAMMA',\n",
    "        'INSAR_ISCE_BURST'\n",
    "    ]\n",
    ")\n",
    "print(\"Select your SBAS stack's HyP3 product type\")\n",
    "display(product_select)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6cd265a1-634c-4e7d-b60b-381f874a8dae",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "product_type = product_select.value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1b2dad03-6ade-4bd2-8103-79641edcc757",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "search = \"Search\" in hyp3_project.value\n",
    "if search:\n",
    "    my_hyp3_info = hyp3.my_info()\n",
    "    active_projects = dict()\n",
    "    \n",
    "    print(\"Checking all HyP3 projects for current INSAR_GAMMA jobs\")\n",
    "    for project in tqdm(my_hyp3_info['job_names']):\n",
    "            batch = Batch()\n",
    "            batch = hyp3.find_jobs(\n",
    "                name=project, \n",
    "                job_type=product_type\n",
    "            ).filter_jobs(running=False, include_expired=False)\n",
    "            if len(batch) > 0:\n",
    "                active_projects.update({batch.jobs[0].name: batch})\n",
    "    \n",
    "    if len(active_projects) > 0:\n",
    "        display(Markdown(\"<text style='color:darkred;'>Note: After selecting a project, you must select the next cell before hitting the 'Run' button or typing Shift/Enter.</text>\"))\n",
    "        display(Markdown(\"<text style='color:darkred;'>Otherwise, you will rerun this code cell.</text>\"))\n",
    "        print('\\nSelect a Project:')\n",
    "        project_select = osl.select_parameter(active_projects.keys())\n",
    "        display(project_select)\n",
    "    else:\n",
    "        print(\"Found no active projects containing InSAR products\")\n",
    "else:\n",
    "    username = input(\"enter the HyP3 username on the account containing an SBAS stack to migrate\")\n",
    "    project_name = input(\"Enter the HyP3 project name\")\n",
    "    batch = Batch()\n",
    "    batch = hyp3.find_jobs(\n",
    "        name=project_name, \n",
    "        job_type=product_type, \n",
    "        user_id=username\n",
    "    ).filter_jobs(running=False, include_expired=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a18ba167-933d-4ffc-bca4-ed4334c0a318",
   "metadata": {},
   "source": [
    "**Select a date range of products to migrate:**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5b0413d4-266d-422a-af7d-8f0be64ce439",
   "metadata": {},
   "outputs": [],
   "source": [
    "if search:\n",
    "    jobs = active_projects[project_select.value]\n",
    "else:\n",
    "    jobs = batch\n",
    "\n",
    "display(Markdown(\"<text style='color:darkred;'>Note: After selecting a date range, you should select the next cell before hitting the 'Run' button or typing Shift/Enter.</text>\"))\n",
    "display(Markdown(\"<text style='color:darkred;'>Otherwise, you may simply rerun this code cell.</text>\"))\n",
    "print('\\nSelect a Date Range:')\n",
    "dates = osl.get_job_dates(jobs)\n",
    "date_picker = osl.gui_date_picker(dates)\n",
    "display(date_picker)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f7ec3cd7-27fd-46b4-bf20-b2adbe14027e",
   "metadata": {},
   "source": [
    "**Save the selected date range and remove products falling outside of it:**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d054f377-520b-4015-9c43-6c7641d09ebd",
   "metadata": {},
   "outputs": [],
   "source": [
    "date_range = osl.get_slider_vals(date_picker)\n",
    "date_range[0] = date_range[0].date()\n",
    "date_range[1] = date_range[1].date()\n",
    "print(f\"Date Range: {str(date_range[0])} to {str(date_range[1])}\")\n",
    "jobs = osl.filter_jobs_by_date(jobs, date_range)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1b2facfa-175a-49ef-a9bb-efd7ec72f791",
   "metadata": {},
   "source": [
    "**Gather the available paths and orbit directions for the remaining products:**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e119636a-93a5-479f-8a33-594073104077",
   "metadata": {},
   "outputs": [],
   "source": [
    "display(Markdown(\"<text style='color:darkred;'><text style='font-size:150%;'>This may take some time for projects containing many jobs...</text></text>\"))\n",
    "osl.set_paths_orbits(jobs)\n",
    "paths = set()\n",
    "for p in jobs:\n",
    "    paths.add(p.path)\n",
    "display(Markdown(f\"<text style=color:blue><text style='font-size:175%;'>Done.</text></text>\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "836f377b-cb38-418c-b362-821574276a53",
   "metadata": {},
   "source": [
    "**Select a path:**\n",
    "\n",
    "- Sentinel-1 has a 12-day repeat cycle so it is not appropriate to merge interferograms across multiple orbital paths. \n",
    "- If multiple paths are represented in the SBAS stack, select one."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c84cfbac-454d-41c6-bde0-4f6332247674",
   "metadata": {},
   "outputs": [],
   "source": [
    "display(Markdown(\"<text style='color:darkred;'>Note: After selecting a path, you must select the next cell before hitting the 'Run' button or typing Shift/Enter.</text>\"))\n",
    "display(Markdown(\"<text style='color:darkred;'>Otherwise, you will simply rerun this code cell.</text>\"))\n",
    "print('\\nSelect a Path:')\n",
    "path_choice = osl.select_parameter(paths)\n",
    "display(path_choice)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fe8443a8-f5c0-49ea-8ff5-289892339d80",
   "metadata": {},
   "source": [
    "**Save the selected flight path/s:**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8faddbeb-dd8e-4019-883a-f4074e18b00e",
   "metadata": {},
   "outputs": [],
   "source": [
    "flight_path = path_choice.value\n",
    "if flight_path:\n",
    "    if flight_path:\n",
    "        print(f\"Flight Path: {flight_path}\")\n",
    "    else:\n",
    "        print('Flight Path: All Paths')\n",
    "else:\n",
    "    print(\"WARNING: You must select a flight path in the previous cell, then rerun this cell.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a2674277-2205-42ad-bfca-31864bdf8b2d",
   "metadata": {},
   "source": [
    "**Filter jobs by path:**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9c502067-f99d-494b-8a76-61c79a97695e",
   "metadata": {},
   "outputs": [],
   "source": [
    "jobs = osl.filter_jobs_by_path(jobs, [flight_path])\n",
    "print(f\"There are {len(jobs)} products to migrate.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "860a58f0-4c08-4ff9-8cf8-6b83dc0d757a",
   "metadata": {},
   "source": [
    "**Migrate the products, unzip them into a directory named after the product type, and delete the zip files:**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "494360ca-c07e-46db-8a0b-d383a19df815",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "print(f\"\\nProject: {jobs.jobs[0].name}\")\n",
    "project_zips = jobs.download_files(data_path)\n",
    "for z in project_zips:\n",
    "    osl.asf_unzip(str(data_path), str(z))\n",
    "    z.unlink()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3346ad90-3426-480a-ae36-29366788e45d",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "## 3. Confirm Presence of a DEM, Azimuth Angle Map, and Incidence Angle Map\n",
    "\n",
    "- These are optional addon products for HyP3, which are necessary for MintPy\n",
    "    - Incidence angle maps are included with HyP3 jobs when the `Include Look Vectors` option is selected.\n",
    "    - DEMs are included with HyP3 jobs when the `Include DEM` option is selected\n",
    "- This is an optional addon product for HyP3, which is necessary for MintPy if running the correct_SET (Solid Earth Tides) step\n",
    "    - Azimuth angle maps are included with HyP3 jobs when the `Include Look Vectors` option is selected\n",
    "\n",
    "**All of the above mentioned files will be included in an InSAR project if Set MintPy Options is selected when adding InSAR jobs to a project in ASF-Search (Vertex)**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e9799fd1-b252-4515-8b23-40780376b58f",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "dem = sorted(list(data_path.glob('*/*dem*.tif')))\n",
    "lv_phi = sorted(list(data_path.glob('*/*lv_phi*.tif')))\n",
    "lv_theta = sorted(list(data_path.glob('*/*lv_theta*.tif')))\n",
    "water_mask = sorted(list(data_path.glob('*/*_water_mask*.tif')))\n",
    "unw = sorted(list(data_path.glob('*/*_unw_phase*.tif')))\n",
    "corr = sorted(list(data_path.glob('*/*_corr*.tif')))\n",
    "conn_comp = sorted(list(data_path.glob('*/*_conncomp*.tif')))\n",
    "tiff_path = dem + lv_phi + lv_theta + water_mask + unw + corr + conn_comp\n",
    "\n",
    "if len(dem) > 0:\n",
    "    print(\"Success: Found at least 1 DEM.\")\n",
    "else:\n",
    "    raise FileNotFoundError(\"Failed to find at least 1 DEM.\\n\"\n",
    "                            \"You will not be able to successfully run a MintPy time-series unless you\"\n",
    "                            \"reorder your HyP3 project with DEMS or provide one from another source.\")\n",
    "                            \n",
    "if len(lv_phi) > 0:\n",
    "    print(\"Success: Found at least 1 lv_phi look vector file.\")\n",
    "else:\n",
    "    raise FileNotFoundError(\"Failed to find at least 1 lv_phi look vector file.\\n\"\n",
    "                            \"You will not be able to successfully run a MintPy time-series unless your\"\n",
    "                            \"reorder your HyP3 project with 'Include Look Vectors' option selected.\")\n",
    "    \n",
    "if len(lv_theta) > 0:\n",
    "    print(\"Success: Found at least 1 lv_theta look vector file.\")\n",
    "else:\n",
    "    raise FileNotFoundError(\"Failed to find at least 1 lv_theta look vector file.\\n\"\n",
    "                            \"You will not be able to successfully run a MintPy time-series unless your\"\n",
    "                            \"reorder your HyP3 project with 'Include Look Vectors' option selected.\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "69bae960-a5d5-42a2-8b90-5dd57f1f878c",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "gdf = gpd.GeoDataFrame(\n",
    "    {\n",
    "    'tiff_path': tiff_path,\n",
    "    'EPSG': [util.get_epsg(p) for p in tiff_path],\n",
    "    'geometry': [util.get_geotiff_bbox(p) for p in tiff_path],\n",
    "    }\n",
    ")\n",
    "\n",
    "# check for multiple projections and project to the predominant EPSG \n",
    "if gdf['EPSG'].nunique() > 1:\n",
    "    proj_count = gdf['EPSG'].value_counts()\n",
    "    predominant_epsg = proj_count.idxmax()\n",
    "    print(f'reprojecting to predominant EPSG: {predominant_epsg}')\n",
    "    for _, row in gdf.loc[gdf['EPSG'] != predominant_epsg].iterrows():\n",
    "        pth = row['tiff_path']\n",
    "        no_data_val = util.get_no_data_val(pth)\n",
    "        res = util.get_res(pth)\n",
    "        \n",
    "        temp = pth.parent/f\"temp_{pth.stem}.tif\"\n",
    "        pth.rename(temp)\n",
    "        src_epsg = row['EPSG']\n",
    "\n",
    "        warp_options = {\n",
    "            \"dstSRS\":f\"EPSG:{predominant_epsg}\", \"srcSRS\":f\"EPSG:{src_epsg}\",\n",
    "            \"targetAlignedPixels\":True,\n",
    "            \"xRes\":res, \"yRes\":res,\n",
    "            \"dstNodata\": no_data_val\n",
    "        }\n",
    "        gdal.Warp(str(pth), str(temp), **warp_options)\n",
    "        temp.unlink()\n",
    "\n",
    "    gdf = gpd.GeoDataFrame(\n",
    "    {\n",
    "    'tiff_path': tiff_path,\n",
    "    'EPSG': [util.get_epsg(p) for p in tiff_path],\n",
    "    'geometry': [util.get_geotiff_bbox(p) for p in tiff_path],\n",
    "    }\n",
    "    )\n",
    "        \n",
    "display(gdf)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "48e233a7-32fe-4404-a536-a4b2827460b8",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "---\n",
    "## 4. Subset the Stack\n",
    "\n",
    "**Select how to define your AOI for subsetting**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b0ff5364-8f2c-4107-afde-83d42c8db86c",
   "metadata": {
    "editable": true,
    "scrolled": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "subset_option = osl.select_parameter([\n",
    "    'Draw a bounding box on a map',\n",
    "    'Subset to area of common coverage, shared by all data in the stack',\n",
    "    'Provide a polygon in Well-Known Text (WKT)',\n",
    "    'Provide a shapefile'\n",
    "])\n",
    "\n",
    "display(subset_option)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cc8c848d-3ae3-4b04-aef5-51d346cd88d9",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "**Determine the maximum and common extents of the stack and plot an Area-of_Interest Selector:**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cde455e0-be45-4dd5-b48d-1fec72086c7c",
   "metadata": {
    "editable": true,
    "scrolled": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "from shapely.geometry import Polygon\n",
    "import shapely.wkt\n",
    "\n",
    "draw = 'Draw' in subset_option.value\n",
    "wkt_poly = 'WKT' in subset_option.value\n",
    "common_coverage = 'common coverage' in subset_option.value\n",
    "shapefile = 'shapefile' in subset_option.value\n",
    "\n",
    "max_extents = osl.get_max_extents(unw)\n",
    "xmin, ymin, xmax, ymax = transform_bounds(int(osl.get_projection(str(unw[0]))), 3857, *max_extents)\n",
    "max_extents_3857 = [xmin, ymin, xmax, ymax]\n",
    "\n",
    "common_extents = osl.get_common_coverage_extents(unw)\n",
    "xmin, ymin, xmax, ymax = transform_bounds(int(osl.get_projection(str(unw[0]))), 3857, *common_extents)\n",
    "common_extents_3857 = [xmin, ymin, xmax, ymax]\n",
    "\n",
    "if draw:   \n",
    "    print('Maximum Extents: Pixels in this area are guaranteed to be included in at least one interferogram in the stack.\\n\\n')\n",
    "    \n",
    "    print('Common Extents: Pixels in this area are guaranteed to included in every interferogram in the stack.\\n\\n')\n",
    "\n",
    "    img_src = '<img src=\"https://opensarlab-docs.asf.alaska.edu/opensarlab-notebook-assets/notebook_images/common_max_stack_coverage.png\" alt=\"\" width=\"500\"/>'\n",
    "    \n",
    "    display(HTML(img_src))\n",
    "    \n",
    "    print(f\"Select an AOI inside the common area covered by the stack.\")\n",
    "    \n",
    "    aoi = osl.AOI_Selector(max_extents_3857, common_extents_3857, figsize=(10, 8)) \n",
    "elif shapefile:\n",
    "    print('Select a shapefile (*.shp)')\n",
    "    shp_fc = FileChooser(Path.home())\n",
    "    display(shp_fc)\n",
    "else:\n",
    "    correct_wkt_input = False\n",
    "    while not correct_wkt_input:\n",
    "        if common_coverage:\n",
    "            wkt = (f'POLYGON(({common_extents[0]} {common_extents[1]}, {common_extents[2]} {common_extents[1]}, {common_extents[2]} '\n",
    "                   f'{common_extents[3]}, {common_extents[0]} {common_extents[3]}, {common_extents[0]} {common_extents[1]}))')\n",
    "            wkt_shapely_geom = shapely.wkt.loads(wkt)\n",
    "        else:\n",
    "            wkt, wkt_shapely_geom = util.get_valid_wkt()\n",
    "            \n",
    "        wkt_ogr_geom = ogr.CreateGeometryFromWkt(wkt)\n",
    "        \n",
    "        if not util.check_within_bounds(wkt_shapely_geom, gdf):\n",
    "            print('WKT exceeds bounds of at least one dataset')\n",
    "            if common_coverage:\n",
    "                raise Exception('Error determining area of common coverage')\n",
    "            else:\n",
    "                continue\n",
    "    \n",
    "        correct_wkt_input = True\n",
    "        \n",
    "    shp_path = data_path / f'shape_{datetime.strftime(datetime.now(), \"%Y%m%dT%H%M%S\")}.shp'\n",
    "    epsg = int(gdf.iloc[0]['EPSG'])\n",
    "    util.save_shapefile(wkt_ogr_geom, epsg, shp_path)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bc3bcbb4-ec34-4298-a685-52572e5a648d",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "**If providing an AOI with WKT or a shapefile, or if subsetting to a common extent, confirm the AOI location on a plot before subsetting**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d835cb9c-829b-4ae9-b1e8-ef5d1557ac55",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "plt.close()\n",
    "if shapefile:\n",
    "    shp_path = Path(shp_fc.selected)\n",
    "    if shp_path.suffix != '.shp':\n",
    "        raise Exception(f'Selected file suffix not \".shp\"')\n",
    "if not draw:   \n",
    "    shp_gdf = gpd.read_file(shp_path)\n",
    "    shp_gdf = shp_gdf.to_crs(crs='EPSG:4326')\n",
    "\n",
    "    box1 = gpd.GeoDataFrame({\"geometry\": [box(*max_extents)]}, crs=f'EPSG:{epsg}')\n",
    "    box2 = gpd.GeoDataFrame({\"geometry\": [box(*common_extents)]}, crs=f'EPSG:{epsg}')\n",
    "    box1 = box1.to_crs(crs='EPSG:4326')\n",
    "    box2 = box2.to_crs(crs='EPSG:4326')\n",
    "    \n",
    "    fig, ax = plt.subplots(1, 1, figsize=(10, 10))\n",
    "    shp_gdf.plot(ax=ax, color='blue', alpha=0.25, edgecolor='k')\n",
    "    box1.plot(ax=ax, color='none', edgecolor='red', linewidth=4)\n",
    "    box2.plot(ax=ax, color='none', edgecolor='green', linewidth=4)\n",
    "\n",
    "    ctx.add_basemap(ax, crs=shp_gdf.crs.to_string(), source=ctx.providers.Esri.WorldImagery)\n",
    "    ctx.add_basemap(ax, crs=shp_gdf.crs.to_string(), source=ctx.providers.OpenStreetMap.Mapnik, alpha=0.5)\n",
    "    \n",
    "    minx, miny, maxx, maxy = box1.total_bounds\n",
    "    buffer_percentage = 0.05\n",
    "    buffer_x = (maxx - minx) * buffer_percentage\n",
    "    buffer_y = (maxy - miny) * buffer_percentage\n",
    "    buffered_minx = minx - buffer_x\n",
    "    buffered_maxx = maxx + buffer_x\n",
    "    buffered_miny = miny - buffer_y\n",
    "    buffered_maxy = maxy + buffer_y\n",
    "    \n",
    "    ax.set_xlim(buffered_minx, buffered_maxx)\n",
    "    ax.set_ylim(buffered_miny, buffered_maxy)\n",
    "    ax.set_title('Loaded Shape for Subsetting')\n",
    "    \n",
    "    legend_elements = [\n",
    "        Line2D([0], [0], color='red', lw=4, label='Max area covered by data stack'),\n",
    "        Line2D([0], [0], color='green', lw=4, label='Common area covered by data stack')\n",
    "    ]\n",
    "    ax.legend(handles=legend_elements)\n",
    "    plt.show()\n",
    "    print('Max area covered by data stack: Every pixel guaranteed to be present in at least 1 image')\n",
    "    print('Common area covered by data stack: Every pixel guaranteed to be present in all images')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a13ad4d8-6550-45bc-8a09-31e644ee6a8e",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "**Subset the data** "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4c6cf820-6a21-4b94-8d97-ebb1845be0a3",
   "metadata": {
    "editable": true,
    "scrolled": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "for pth in tqdm(gdf['tiff_path']):\n",
    "    print(f'Subsetting: {pth}')\n",
    "\n",
    "    if draw:\n",
    "        try:\n",
    "            xmin, ymin, xmax, ymax = transform_bounds(\n",
    "                3857,\n",
    "                int(gdf['EPSG'].iloc[0]), \n",
    "                *[aoi.x1, aoi.y1, aoi.x2, aoi.y2]\n",
    "            )\n",
    "            ul = [xmin, ymax]\n",
    "            lr = [xmax, ymin]\n",
    "        except TypeError:\n",
    "            print('TypeError')\n",
    "            display(Markdown(f'<text style=color:red>This error may occur if an AOI was not selected.</text>'))\n",
    "            display(Markdown(f'<text style=color:red>Note that the square tool icon in the AOI selector menu is <b>NOT</b> the selection tool. It is the zoom tool.</text>'))\n",
    "        temp_pth = pth.parent/f'subset_{pth.name}'\n",
    "        gdal.Translate(destName=str(temp_pth), srcDS=str(pth), projWin=[ul[0], ul[1], lr[0], lr[1]])\n",
    "        pth.unlink()\n",
    "        temp_pth.rename(pth)\n",
    "    elif common_coverage:\n",
    "        temp_pth = pth.parent/f'subset_{pth.name}'\n",
    "        gdal.Translate(destName=str(temp_pth), srcDS=str(pth), projWin=[common_extents[0], common_extents[3], common_extents[2], common_extents[1]])\n",
    "        pth.unlink()\n",
    "        temp_pth.rename(pth)\n",
    "    else:\n",
    "        warp_args = gdal.WarpOptions(cutlineDSName=shp_path, cropToCutline=True) \n",
    "        gdal.Warp(str(pth), str(pth), options=warp_args)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "11d2af15-039e-4194-bb69-a606c8b5945d",
   "metadata": {
    "editable": true,
    "scrolled": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "**Select whether to convert to WGS84 (lat/lon) prior to loading into MintPy**\n",
    "\n",
    "- HyP3 data are delivered in UTM.\n",
    "- This notebook intentionally only converts to WGS84, not back to UTM.\n",
    "- Repeatedly converting data back and forth between coordinate systems will lead to a loss in precision.\n",
    "- If you have already converted the data to WGS84 and wish to go back to UTM, delete the data and rerun this notebook to access fresh copies from HyP3."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a12faccc-8b8d-4ff2-9485-37176bbf4dc3",
   "metadata": {},
   "outputs": [],
   "source": [
    "coord_sys_option = osl.select_parameter(['UTM', 'WGS84 (lat/lon)'], 'Select a coordinate system for your data')\n",
    "display(coord_sys_option)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "72e4bb1f-e835-4bcb-bb5f-7c1d6645b97c",
   "metadata": {},
   "source": [
    "**If desired, convert to WGS84**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7b10b607-8d3e-4b49-a029-ebb0a770fa1f",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "wgs84 = 'WGS84' in coord_sys_option.value\n",
    "if wgs84:\n",
    "    for pth in tqdm(gdf['tiff_path']):\n",
    "        print(f'Converting {pth} to WGS84')\n",
    "        gdal.Warp(str(pth), str(pth), dstSRS='EPSG:4326')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "da7df9d9-04bf-494e-b162-8cfa7416ba4e",
   "metadata": {},
   "source": [
    "**Remove any subset scenes containing no data:**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "20888840-21ea-4496-b4f4-ad183f9ef7c7",
   "metadata": {},
   "outputs": [],
   "source": [
    "removed = []\n",
    "# check unw_phase for no-data filled rasters and remove parent directories accordingly\n",
    "for pth in tqdm(gdf.loc[gdf['tiff_path'].astype(str).str.contains('unw_phase', na=False)]['tiff_path']):\n",
    "    try:\n",
    "        raster = gdal.Open(str(pth))\n",
    "        band = raster.ReadAsArray()\n",
    "    except RuntimeError:\n",
    "        raise\n",
    "        \n",
    "    if np.count_nonzero(band) < 1:\n",
    "        shutil.rmtree(pth.parent)\n",
    "        removed.append(pth)\n",
    "\n",
    "if len(removed) == 0:\n",
    "    print(\"No interferograms were removed\")\n",
    "else:\n",
    "    print(f\"{len(removed)} interferograms removed:\")\n",
    "    for pth in removed:\n",
    "        print(pth.parent)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:.local-opensarlab_mintpy_recipe_book]",
   "language": "python",
   "name": "conda-env-.local-opensarlab_mintpy_recipe_book-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.9.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}