Output_GeoTiff.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "842af782-0b88-4fca-921c-9f79e2c759dd",
   "metadata": {},
   "source": [
    "# Output Results to GeoTiff\n",
    "\n",
    "*Author: Alex Lewandowski; Alaska Satellite Facility*\n",
    "\n",
    "Output the results of your MintPy Time Series to GeoTiff\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": "8cde9dcb-a7ce-459d-ae9a-2ba6e61adb69",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "from pathlib import Path\n",
    "import re\n",
    "from typing import Union\n",
    "\n",
    "from affine import Affine\n",
    "import h5py\n",
    "from ipyfilechooser import FileChooser\n",
    "import mintpy.utils.readfile\n",
    "import mintpy.view\n",
    "import numpy as np\n",
    "import opensarlab_lib as osl\n",
    "from osgeo import gdal\n",
    "gdal.UseExceptions()\n",
    "import rasterio\n",
    "from rasterio.crs import CRS\n",
    "import rioxarray as rxr\n",
    "import shutil\n",
    "from tqdm.notebook import tqdm\n",
    "from util import util\n",
    "import xarray as xr"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5ebb8d04-d0a8-4c6f-97d8-ebd58aabab6e",
   "metadata": {},
   "source": [
    "---\n",
    "## 1. Select Your Project's Custom Config File and Define Paths\n",
    "\n",
    "- The custom config is located in your project's `MintPy` directory\n",
    "- It is a text file named after your project\n",
    "  - `path/to/MinPy/my_project.txt`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2942f23e-60f9-4d99-8357-7b2ab5eb07e8",
   "metadata": {},
   "outputs": [],
   "source": [
    "file_chooser_path = util.get_recent_mintpy_config_path()\n",
    "if file_chooser_path:\n",
    "    fc = FileChooser(path=file_chooser_path.parent, filename=file_chooser_path.name, select_default=True)\n",
    "else:\n",
    "    file_chooser_path = Path.home()\n",
    "    fc = FileChooser(file_chooser_path, select_default=False)\n",
    "    \n",
    "print(\"Select your custom MintPy config file (MintPy/my_project_name.txt):\")\n",
    "display(fc)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e9de6112-e508-4ca4-ada5-69776a3ec941",
   "metadata": {},
   "source": [
    "**Define paths**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5ab2b462-44e2-4791-8feb-7faa16b4bf89",
   "metadata": {},
   "outputs": [],
   "source": [
    "if Path(fc.selected) != file_chooser_path:\n",
    "    util.write_recent_mintpy_config_path(Path(fc.selected))\n",
    "\n",
    "config_path = Path(fc.selected)\n",
    "mint_path = config_path.parent\n",
    "inputs_path = mint_path/ 'inputs'\n",
    "ifgramstack = inputs_path / 'ifgramStack.h5'\n",
    "geotiff_path = mint_path/'GeoTiffs'\n",
    "ts_demErr_path = list(mint_path.glob('timeseries*_demErr.h5'))[0]\n",
    "disp_path = geotiff_path / 'displacement_maps'\n",
    "unwrapped_path = disp_path / 'unwrapped'\n",
    "wrapped_path = disp_path / 'wrapped'"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b53da7d8-258e-4656-a6af-b900fe856311",
   "metadata": {},
   "source": [
    "---\n",
    "## 2. Create a List of Dates for All Timesteps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8cd63c9c-fcaa-40dd-92f2-8f7078bfab20",
   "metadata": {},
   "outputs": [],
   "source": [
    "ifgramstack = inputs_path/\"ifgramStack.h5\"\n",
    "\n",
    "with h5py.File(ifgramstack, \"r\") as f:\n",
    "    dates = f[\"date\"][()]\n",
    "    dates = list(set([d.decode(\"utf-8\") for insar in dates for d in insar]))\n",
    "    dates.sort()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e4825e58-9a7b-41c4-8449-ab6960995e42",
   "metadata": {},
   "source": [
    "---\n",
    "## 3. Save the Cumulative Displacement Timeseries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "df413ca9-62b4-4ffa-96b6-f229c284799f",
   "metadata": {},
   "outputs": [],
   "source": [
    "ds = f'{dates[0]}_{dates[-1]}'\n",
    "!save_gdal.py $ts_demErr_path -d $ds --of GTIFF -o $geotiff_path/\"save_gdal_ts_demErr.tif\""
   ]
  },
  {
   "cell_type": "markdown",
   "id": "15340154-9517-49e3-965b-130fba259582",
   "metadata": {},
   "source": [
    "---\n",
    "## 4. Set the Scale for Unwrapped Interferograms"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ce365391-1374-42c9-b0b7-a1b8dde4d6bf",
   "metadata": {},
   "outputs": [],
   "source": [
    "scale_choice = osl.select_parameter(['Meters', 'Centimeters'], description='Select a scale for your unwrapped interferograms:')\n",
    "display(scale_choice)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1d09d01f-ec7a-4f89-90b7-4f1db095364e",
   "metadata": {},
   "source": [
    "---\n",
    "## 5. Save the Unwrapped Displacement GeoTiffs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cbfa2c59-51ef-4912-a125-559e7c30e1e6",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "meters = scale_choice.value == 'Meters'\n",
    "scaled_unwrapped_path = unwrapped_path / scale_choice.value\n",
    "scaled_unwrapped_path.mkdir(parents=True, exist_ok=True)\n",
    "\n",
    "_, unw_info = mintpy.utils.readfile.read(ifgramstack)\n",
    "\n",
    "geotrans = (\n",
    "    float(unw_info['X_FIRST']),\n",
    "    float(unw_info['X_STEP']),\n",
    "    0.0, \n",
    "    float(unw_info['Y_FIRST']), \n",
    "    0.0, \n",
    "    float(unw_info['Y_STEP'])\n",
    ")\n",
    "geotrans = Affine.from_gdal(*geotrans)\n",
    "\n",
    "crs = CRS.from_epsg(unw_info['EPSG'])\n",
    "\n",
    "for i, d in enumerate(tqdm(dates)):\n",
    "    date_range = f'{dates[0]}_{dates[i]}'\n",
    "    cmd = f'view.py {ts_demErr_path} {date_range} --notitle --notick --noaxis'\n",
    "    data, _, _ = mintpy.view.prep_slice(cmd)\n",
    "\n",
    "    if meters:\n",
    "        data = data / 100 # cm -> meters\n",
    "        \n",
    "\n",
    "    with rasterio.open(f'{scaled_unwrapped_path}/{date_range}_{ts_demErr_path.stem}_{scale_choice.value}_unwrapped.tif', 'w', driver='GTiff',\n",
    "                  height = data.shape[0], width = data.shape[1],\n",
    "                  count=1, dtype=str(data.dtype),\n",
    "                  crs=crs,\n",
    "                  transform=geotrans,\n",
    "                  nodata=np.nan) as ds:\n",
    "        \n",
    "        ds.write(data.astype(rasterio.float32), 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e549867b-f10a-45df-b570-00f5d7a5d786",
   "metadata": {},
   "source": [
    "---\n",
    "## 6. Generate and Save Wrapped Interferograms as GeoTiff with a matplotlib 'jet'-like Colormap Symbology\n",
    "\n",
    "**Write a function to add a color ramp to single band GeoTiff**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "caba04f1-c3c7-4e84-a67f-fb47b08ca1a0",
   "metadata": {},
   "outputs": [],
   "source": [
    "def colorize_wrapped_insar(tif_path: Union[str, os.PathLike]):\n",
    "    \"\"\"\n",
    "    Blue: 0 and 2π\n",
    "    Red: π/2\n",
    "    Yellow: π\n",
    "    Green 3/2π\n",
    "    \"\"\"\n",
    "    ds = gdal.Open(str(tif_path), 1)\n",
    "    band = ds.GetRasterBand(1)\n",
    "\n",
    "    # create color table\n",
    "    colors = gdal.ColorTable()\n",
    "    \n",
    "    colors.CreateColorRamp(0, (0, 0, 255),  64, (255, 0, 0)) \n",
    "    colors.CreateColorRamp(64, (255, 0, 0),   128, (255, 255, 0))\n",
    "    colors.CreateColorRamp(128, (255, 255, 0), 192, (0, 255, 0))\n",
    "    colors.CreateColorRamp(192, (0, 255, 0),   255, (0, 0, 255))\n",
    "\n",
    "    # set color table and color interpretation\n",
    "    band.SetRasterColorTable(colors)\n",
    "    band.SetRasterColorInterpretation(gdal.GCI_PaletteIndex)\n",
    "\n",
    "    # close and save file\n",
    "    del band, ds"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b0f7116f-f34e-434e-bc09-f04543b67172",
   "metadata": {},
   "source": [
    "**Collect paths to unwrapped displacement maps**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1099018e-98e9-46f6-a1ef-e312c510e2b5",
   "metadata": {},
   "outputs": [],
   "source": [
    "unwrapped_paths = list((unwrapped_path / scale_choice.value).rglob('*_unwrapped.tif'))\n",
    "unwrapped_paths"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e97f624a-351b-40c6-b07b-5d2c78e25d63",
   "metadata": {},
   "source": [
    "**Generate the wrapped interferogram GeoTiffs**\n",
    "\n",
    "- Please note that the wrapped range used below is currently under review and may not yet correctly correspond to the Sentinel-1 wavelength "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "105eccc1-5c95-4ad9-8012-6e92a3408521",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "meters = scale_choice.value == 'Meters'\n",
    "wrapped_path.mkdir(parents=True, exist_ok=True)\n",
    "\n",
    "sentinel_c_band_lambda = 0.055465763\n",
    "if not meters:\n",
    "    sentinel_c_band_lambda *= 100\n",
    "\n",
    "for unw_path in tqdm(unwrapped_paths):\n",
    "    date_range_regex = f'(?<=/unwrapped/{scale_choice.value}/)\\\\d{{8}}_\\\\d{{8}}'\n",
    "    date_range = re.search(date_range_regex, str(unw_path)).group(0)\n",
    "    \n",
    "    with rxr.open_rasterio(unw_path, masked=True).squeeze() as ds:\n",
    "        # convert unwrapped raster to radians\n",
    "        with xr.set_options(keep_attrs=True):\n",
    "            unw_rad = (ds * 4 * np.pi) / sentinel_c_band_lambda\n",
    "          \n",
    "    wrap_range = [-np.pi, np.pi]\n",
    "       \n",
    "    # wrap the interferogram\n",
    "    with xr.set_options(keep_attrs=True):\n",
    "        wrap = mintpy.utils.utils0.wrap(unw_rad, wrap_range=wrap_range)\n",
    "\n",
    "    # collect crs and transform\n",
    "    with rasterio.open(unw_path, 'r', driver='GTiff') as ds:\n",
    "        unw_crs = ds.read_crs()\n",
    "        unw_transform = ds.transform\n",
    "    \n",
    "    # Save wrapped interferogram as a GeoTiff\n",
    "    wrp_path = wrapped_path/f'{date_range}_{ts_demErr_path.stem}_wrapped_unscaled.tif'\n",
    "    with rasterio.open(wrp_path, 'w', driver='GTiff',\n",
    "                      height = wrap.shape[0], width = wrap.shape[1],\n",
    "                      count=1, dtype=str(wrap.dtype),\n",
    "                      crs=crs,\n",
    "                      transform=geotrans,\n",
    "                      nodata=np.nan) as ds:\n",
    "        ds.write(wrap.astype(rasterio.float32), 1)\n",
    "\n",
    "    # scale wrapped interferogram (0 to 255)\n",
    "    scaled_path = wrapped_path/f'{wrp_path.stem.split(\"_unscaled\")[0]}_scaled.tif'\n",
    "    !gdal_translate -of GTiff -scale -ot BYTE $wrp_path $scaled_path\n",
    "    wrp_path.unlink()\n",
    "    \n",
    "    # add color ramp\n",
    "    colorize_wrapped_insar(scaled_path)\n",
    "    \n",
    "    # convert to 3-band rgb\n",
    "    three_band_path = wrapped_path/f'{scaled_path.stem.split(\"_scaled\")[0]}.tif'\n",
    "    !gdal_translate -of GTiff -expand rgb $scaled_path $three_band_path\n",
    "    scaled_path.unlink()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f49174f3-0515-48db-a620-6969c90e1958",
   "metadata": {},
   "source": [
    "---\n",
    "## 6. Save the Temporal Coherence GeoTiff"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d4ee5629-eb93-4fad-acd4-d8b968aeae69",
   "metadata": {},
   "outputs": [],
   "source": [
    "!save_gdal.py $mint_path/temporalCoherence.h5 --of GTIFF -o $geotiff_path/TemporalCoherence.tif"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f254355f-b1d4-48cc-90ed-c4c6f5794675",
   "metadata": {},
   "source": [
    "---\n",
    "## 7. Save the Spatial Coherence GeoTiff"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7abcb4a2-1251-454a-a33e-c71120044ec9",
   "metadata": {},
   "outputs": [],
   "source": [
    "!save_gdal.py $mint_path/avgSpatialCoh.h5 --of GTIFF -o $geotiff_path/avgSpatialCoh.tif"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7a21a0cf-8264-48e6-928e-b4e175054166",
   "metadata": {},
   "source": [
    "---\n",
    "## 8. Save the Velocity GeoTiff"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6e0bc3f5-7e2f-4fc1-bbec-c542120bb499",
   "metadata": {},
   "outputs": [],
   "source": [
    "velocity_path = mint_path / 'velocity.h5'\n",
    "era_corr_velocity_path = mint_path / 'velocityERA5.h5'\n",
    "velocity_path = era_corr_velocity_path if era_corr_velocity_path.exists() else velocity_path\n",
    "vel_tiff = geotiff_path / f'{velocity_path.stem}.tif'\n",
    "!save_gdal.py $velocity_path --of GTIFF -o $vel_tiff"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c4222472-c190-47af-89d2-82a1596cc4b4",
   "metadata": {},
   "source": [
    "---\n",
    "## 9. Compress GeoTiffs for Downloading"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c00002f5-7642-4ec9-a25a-c348102b2148",
   "metadata": {},
   "outputs": [],
   "source": [
    "zip_option = osl.select_parameter([\"Zip `GeoTiffs` directory\", \"Do not zip `GeoTiffs` directory\"], description=\"Do you wish to zip your GeoTiffs?\")\n",
    "display(zip_option)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "28751533-d26f-42c0-92c8-ca3bb0873603",
   "metadata": {},
   "outputs": [],
   "source": [
    "zip = 'Do not' not in zip_option.value\n",
    "\n",
    "if zip:\n",
    "    shutil.make_archive(mint_path/geotiff_path.stem, 'zip', geotiff_path)\n",
    "    print(f'Zip archive location: {mint_path/geotiff_path.stem}.zip')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ffe5f2b3-b4f6-42a7-aa66-01a82ffb56bd",
   "metadata": {},
   "source": [
    "---\n",
    "## 10. Delete Uncompressed Wrapped and Unwrapped Interferograms to Conserve Space\n",
    "\n",
    "If you have compressed your data to download it, you can delete the uncompressed files to conserve storage space."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c54c6180-a329-4930-8065-798fb18e5398",
   "metadata": {},
   "outputs": [],
   "source": [
    "delete_option = osl.select_parameter([\"Delete the uncompressed `displacement_maps` directory\", \"Do not delete the `displacement_maps` directory\"], description=\"Do you wish to delete the uncompressed data?\")\n",
    "display(delete_option)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c43fc621-8700-41ba-9937-d22f2961a0b1",
   "metadata": {},
   "outputs": [],
   "source": [
    "delete = 'Do not' not in delete_option.value\n",
    "\n",
    "if delete:\n",
    "    shutil.rmtree(disp_path, ignore_errors=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c1a52975-7d00-403c-b26a-971e932dac24",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "opensarlab_mintpy_recipe_book [conda env:.local-opensarlab_mintpy_recipe_book]",
   "language": "python",
   "name": "conda-env-.local-opensarlab_mintpy_recipe_book-opensarlab_mintpy_recipe_book"
  },
  "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
}