OPERA_RTC_download_reproject_mosaic_sample_bursts.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "79aef48e-10f1-4684-8b63-0478576fd71b",
   "metadata": {},
   "source": [
    "# **OPERA RTC CalVal Mosaics**\n",
    "\n",
    "Author: Alex Lewandowski\n",
    "\n",
    "**This notebook prepares RTC-ADT CalVal v0.4 release mosaics from samples bursts located in `s3://opera-pst-rs-pop1/products/RTC_S1/`**\n",
    "\n",
    "1. Input a Sentinel-1 scene ID\n",
    "1. Identify burst IDs using ASF burst extractor\n",
    "1. Download bursts\n",
    "1. Reproject bursts to predominant UTM\n",
    "1. Mosaic bursts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b696fc01-2fa4-49a5-9015-d76e4815fac6",
   "metadata": {},
   "outputs": [],
   "source": [
    "from collections import Counter\n",
    "import boto3\n",
    "from datetime import datetime\n",
    "from ipyfilechooser import FileChooser\n",
    "from pathlib import Path\n",
    "import re\n",
    "import requests\n",
    "import shutil\n",
    "from tqdm.auto import tqdm \n",
    "from typing import Union, List, Dict\n",
    "\n",
    "from ipywidgets import Layout\n",
    "\n",
    "import numpy as np\n",
    "from osgeo import gdal\n",
    "gdal.UseExceptions()\n",
    "\n",
    "import asf_search as disco"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e423a955-3de9-47d0-ba69-f14e7f6efda0",
   "metadata": {},
   "source": [
    "### **1. Input S1 scene name and identify bursts**\n",
    "\n",
    "Currently (08-16-2023) v0.4 processed sample burst data are available for scenes in this spreadsheet: https://docs.google.com/spreadsheets/d/1tdj59dZpyx_BFO1NDc46CnW1Bt9UhNs2OK-hTUdhLh0/edit?usp=sharing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a1b4e0c8-5148-490e-84e9-a9f50ba56b03",
   "metadata": {
    "tags": [
     "parameters"
    ]
   },
   "outputs": [],
   "source": [
    "session = requests.Session()\n",
    "\n",
    "# try/except for papermill\n",
    "try:\n",
    "    scene = input(\"Enter a Sentinel-1 granule ID\")\n",
    "except:\n",
    "    pass"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0eabbd73-44bf-4c52-998d-3357433bf6f3",
   "metadata": {
    "scrolled": true,
    "tags": []
   },
   "outputs": [],
   "source": [
    "polarizations = ['VV', 'VH']\n",
    "beam_modes = ['IW1', 'IW2', 'IW3']\n",
    "bursts = range(15)\n",
    "\n",
    "orbital_path = f\"T{disco.granule_search(scene)[1].properties['pathNumber']:03d}\"\n",
    "\n",
    "run = True\n",
    "while run:\n",
    "    opera_bursts = []\n",
    "    run = False\n",
    "    for beam_mode in beam_modes:\n",
    "        for polarization in polarizations:\n",
    "            for burst in bursts:\n",
    "                url = f'https://sentinel1-burst.asf.alaska.edu/{scene}/{beam_mode}/{polarization}/{burst}.xml'\n",
    "                response = session.get(url)\n",
    "                status = response.status_code\n",
    "                if status != 404:\n",
    "                    print(f'{url} {status}')\n",
    "                    if status == 200:\n",
    "                        print(response.headers['content-disposition'])\n",
    "                        burst_swath = '-'.join(response.headers['content-disposition'].split(\"filename=\")[-1].split('_')[1:3])\n",
    "                        opera_bursts.append(f\"OPERA_L2_RTC-S1_{orbital_path}-{burst_swath}\")\n",
    "                    elif status == 202:\n",
    "                        run = True"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c72694c6-01ca-4c87-8721-df7fc71c8f09",
   "metadata": {},
   "outputs": [],
   "source": [
    "len(opera_bursts)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4d90197a-c3ba-40d5-8ceb-423c39436305",
   "metadata": {},
   "source": [
    "### **2. Create data directories**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7775a579-abdd-4e75-b72b-d3b036476590",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Select the directory in which to store your OPERA data\")\n",
    "fc = FileChooser(Path.cwd(), layout=Layout(width='1000px'))\n",
    "display(fc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1675d859-8aa3-4c29-a571-696253f868c4",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# try/except for papermill\n",
    "try:\n",
    "    opera_dir = Path(fc.selected)\n",
    "except TypeError:\n",
    "    pass\n",
    "opera_dir = Path(opera_dir)\n",
    "opera_dir.mkdir(exist_ok=True)\n",
    "\n",
    "scene_dir = opera_dir/f\"OPERA_RTC_{scene}\"\n",
    "vh_dir = scene_dir/\"vh\"\n",
    "vv_dir = scene_dir/\"vv\"\n",
    "inc_angle_dir = scene_dir/\"inc_angle\"\n",
    "local_inc_angle_dir = scene_dir/\"local_inc_angle_dir\"\n",
    "ls_mask_dir = scene_dir/\"ls_mask_dir\"\n",
    "\n",
    "for d in [scene_dir, vh_dir, vv_dir, inc_angle_dir, local_inc_angle_dir, ls_mask_dir]:\n",
    "    d.mkdir(exist_ok=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f0f1c9fc-1ec2-4e28-8449-7ec52e0a2652",
   "metadata": {},
   "source": [
    "### **3. Build dictionary of burst level backscatter and static file S3 URIs**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fda18afd-731f-46d7-b5e5-97736bbcb583",
   "metadata": {},
   "outputs": [],
   "source": [
    "s1_date_regex = '(?<=_)\\d{8}(?=T\\d{6}_\\d{8}T\\d{6})'\n",
    "s1_date_str = re.search(s1_date_regex, scene)\n",
    "if s1_date_str:\n",
    "    s1_date_str = s1_date_str.group(0)\n",
    "else:\n",
    "    raise Exception(\"date string not found in Sentinel-1 scene name\")\n",
    "s1_date_str"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c4d670c8-2915-455b-838b-3fa352b540d5",
   "metadata": {
    "scrolled": true,
    "tags": []
   },
   "outputs": [],
   "source": [
    "bucket_name = 'opera-pst-rs-pop1'\n",
    "s3 = boto3.resource('s3')\n",
    "my_bucket = s3.Bucket(bucket_name)\n",
    "\n",
    "opera_ds = {\n",
    "    'vv': {'dst_dir': vv_dir, 's3_uri': []},\n",
    "    'vh': {'dst_dir': vh_dir, 's3_uri': []},\n",
    "    'inc_angle': {'dst_dir': inc_angle_dir, 's3_uri': []},\n",
    "    'local_inc_angle': {'dst_dir': local_inc_angle_dir, 's3_uri': []},\n",
    "    'ls_mask': {'dst_dir': ls_mask_dir, 's3_uri': []}\n",
    "}\n",
    "\n",
    "vh_regex = '.*_VH.tif$'\n",
    "vv_regex = '.*_VV.tif$'\n",
    "ls_regex = '.*_static_layers/.*_static_layover_shadow_mask.tif$'\n",
    "inc_angle_regex = '.*_static_layers/.*_static_incidence_angle.tif$'\n",
    "local_inc_angle_regex = '.*_static_layers/.*_static_local_incidence_angle.tif$'\n",
    "\n",
    "for burst in tqdm(opera_bursts):  \n",
    "    for object_summary in my_bucket.objects.filter(Prefix=f\"products/RTC_S1/{burst}\"):\n",
    "        if s1_date_str in object_summary.key and re.search(vh_regex, object_summary.key):\n",
    "            opera_ds['vh']['s3_uri'].append(f\"s3://{bucket_name}/{object_summary.key}\")\n",
    "        elif s1_date_str in object_summary.key and re.search(vv_regex, object_summary.key):\n",
    "            opera_ds['vv']['s3_uri'].append(f\"s3://{bucket_name}/{object_summary.key}\")\n",
    "        elif re.search(ls_regex, object_summary.key):\n",
    "            opera_ds['ls_mask']['s3_uri'].append(f\"s3://{bucket_name}/{object_summary.key}\")\n",
    "        elif re.search(inc_angle_regex, object_summary.key):\n",
    "            opera_ds['inc_angle']['s3_uri'].append(f\"s3://{bucket_name}/{object_summary.key}\")\n",
    "        elif re.search(local_inc_angle_regex, object_summary.key):\n",
    "            opera_ds['local_inc_angle']['s3_uri'].append(f\"s3://{bucket_name}/{object_summary.key}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e33621fa-ac6b-4db5-8541-bb4852da0ad8",
   "metadata": {},
   "source": [
    "### **4. Filter duplicated bursts for those most recently processed**\n",
    "\n",
    "The bucket contains bursts from multiple processing runs, some of which used different orbit data. The August 3, 2023 processing run used precise orbits."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7814ae5b-b1d9-4686-a20b-938873e1e054",
   "metadata": {},
   "outputs": [],
   "source": [
    "def burst_id_from_path(pth):\n",
    "    burst_id_regex = '(?<=S1_)T\\d{3}-\\d{6}-IW(1|2|3)'\n",
    "    results = re.search(burst_id_regex, pth)\n",
    "    if results:\n",
    "        return results.group(0)\n",
    "    else:\n",
    "        raise Exception(f\"No burst ID found in path: {pth}\")\n",
    "        \n",
    "def acquisition_timestamp_from_path(pth):\n",
    "    s1_date_regex = '(?<=_)\\d{8}T\\d{6}Z(?=_\\d{8}T\\d{6})'\n",
    "    timestamp_str = re.search(s1_date_regex, pth)\n",
    "    if timestamp_str:\n",
    "        return timestamp_str.group(0)\n",
    "    else:\n",
    "        raise Exception(f\"acquisition timestamp string not found in Sentinel-1 path: {pth}\")\n",
    "        \n",
    "def burst_id_acquisition_timestamp_from_path(pth):\n",
    "    burst_id_regex = '(?<=S1_)T\\d{3}-\\d{6}-IW(1|2|3)_\\d{8}T\\d{6}Z(?=_\\d{8}T\\d{6})'\n",
    "    results = re.search(burst_id_regex, pth)\n",
    "    if results:\n",
    "        return results.group(0)\n",
    "    else:\n",
    "        raise Exception(f'No \"burst-ID_timestamp\"  found in path: {pth}')\n",
    "        \n",
    "\n",
    "\n",
    "def processing_timestamp_from_path(pth):\n",
    "    s1_date_regex = '(?<=_\\d{8}T\\d{6}Z_)\\d{8}T\\d{6}Z(?=_S1)'\n",
    "    timestamp_str = re.search(s1_date_regex, pth)\n",
    "    if timestamp_str:\n",
    "        return timestamp_str.group(0)\n",
    "    else:\n",
    "        raise Exception(f\"processing timestamp string not found in Sentinel-1 path: {pth}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ae865bf7-f70b-4171-a328-d40afeaee3a6",
   "metadata": {},
   "outputs": [],
   "source": [
    "def most_recently_processed_burst_filter(path_list, keep_date_index=-1):\n",
    "    burst_dict = {}\n",
    "    for p in path_list:\n",
    "        burst_id = burst_id_from_path(p)\n",
    "        processing_ts = processing_timestamp_from_path(p)\n",
    "        \n",
    "        if burst_id not in burst_dict.keys():\n",
    "            burst_dict[burst_id] = {}\n",
    "        \n",
    "        \n",
    "        burst_dict[burst_id][processing_ts] = p\n",
    "        \n",
    "\n",
    "    keep_bursts = []\n",
    "    for burst in burst_dict:\n",
    "        date = [datetime.strptime(i, \"%Y%m%dT%H%M%SZ\") for i in burst_dict[burst]]\n",
    "        date.sort()\n",
    "        keep_bursts.append(burst_dict[burst][date[keep_date_index].strftime(\"%Y%m%dT%H%M%SZ\")])\n",
    "\n",
    "    return keep_bursts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "da1c516a-0a9e-4ffd-bcbb-307afd5c0249",
   "metadata": {},
   "outputs": [],
   "source": [
    "# try/except for Papermill\n",
    "try:\n",
    "    if type(keep_date_index) == int:\n",
    "        pass\n",
    "except NameError:\n",
    "    keep_date_index = -1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eaa99252-9ee3-4022-9958-fc67365ec852",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "for ds in tqdm(opera_ds):\n",
    "    opera_ds[ds]['s3_uri'] = most_recently_processed_burst_filter(opera_ds[ds]['s3_uri'], keep_date_index=keep_date_index)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f7241913-a456-41df-8a1d-285208b56515",
   "metadata": {},
   "outputs": [],
   "source": [
    "for ds in opera_ds:\n",
    "    ds_len = len(opera_ds[ds]['s3_uri'])\n",
    "    print(f\"{ds} length: {ds_len}\")\n",
    "    if 26 > ds_len > 30:\n",
    "        raise Exception(f\"Unexpected number of bursts identified for {ds}: {ds_len}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4fad1c95-4d90-4248-a550-0e85f7277306",
   "metadata": {},
   "source": [
    "### **5. Download the data**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "52adeddb-8bb8-4953-b914-09ec7948da0a",
   "metadata": {
    "scrolled": true,
    "tags": []
   },
   "outputs": [],
   "source": [
    "for ds in tqdm(opera_ds):\n",
    "    for burst in opera_ds[ds]['s3_uri']:\n",
    "        !aws s3 cp {burst} {opera_ds[ds]['dst_dir']}/{burst.split('/')[-1]}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d6c70b32-8b88-4fde-9870-ddab6fe4f3e9",
   "metadata": {
    "tags": []
   },
   "source": [
    "### **6. Reproject all bursts (backscatter and static files) to predominant UTM**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4a6af25e-232d-4e99-a3a9-ad26afd10eca",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_projection(img_path: Union[Path, str]) -> Union[str, None]:\n",
    "    \"\"\"\n",
    "    Takes: a string or posix path to a product in a UTM projection\n",
    "\n",
    "    Returns: the projection (as a string) or None if none found\n",
    "    \"\"\"\n",
    "    img_path = str(img_path)\n",
    "    try:\n",
    "        info = gdal.Info(img_path, format='json')['coordinateSystem']['wkt']\n",
    "    except KeyError:\n",
    "        return None\n",
    "    except TypeError:\n",
    "        raise FileNotFoundError\n",
    "\n",
    "    regex = 'ID\\[\"EPSG\",[0-9]{4,5}\\]\\]$'\n",
    "    results = re.search(regex, info)\n",
    "    if results:\n",
    "        return results.group(0).split(',')[1][:-2]\n",
    "    else:\n",
    "        return None\n",
    "    \n",
    "def get_projections(tiff_paths: List[Union[Path, str]]) -> Dict:\n",
    "    \"\"\"\n",
    "    Takes: List of string or posix paths to geotiffs\n",
    "    \n",
    "    Returns: Dictionary key: epsg, value: number of tiffs in that epsg \n",
    "    \"\"\"\n",
    "    epsgs = []\n",
    "    for p in tiff_paths:\n",
    "        epsgs.append(get_projection(p))\n",
    "\n",
    "    epsgs = dict(Counter(epsgs))\n",
    "    return epsgs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6786f87d-a1ae-4941-bcdc-e17c4e97edac",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "vv_paths = list(vv_dir.glob(\"*VV.tif\"))\n",
    "vh_paths = list(vh_dir.glob(\"*VH.tif\"))\n",
    "local_inc_angle_paths = list(local_inc_angle_dir.glob(\"*static_local_incidence_angle.tif\"))\n",
    "inc_angle_paths = list(inc_angle_dir.glob(\"*static_incidence_angle.tif\"))\n",
    "ls_mask_paths = list(ls_mask_dir.glob(\"*static_layover_shadow_mask.tif\"))\n",
    "\n",
    "epsgs = get_projections(vh_paths)\n",
    "predominant_epsg = None if len(epsgs) == 1 else max(epsgs, key=epsgs.get)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "76dada1f-9f4f-468a-8789-d7c2ae6a88e2",
   "metadata": {},
   "outputs": [],
   "source": [
    "epsgs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d3b5cab4-2636-4631-9abd-5748c5667654",
   "metadata": {},
   "outputs": [],
   "source": [
    "predominant_epsg"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8d1014d0-a5b4-4fb5-8509-3db2ecfcdf4d",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_res(tiff):\n",
    "    tiff = str(tiff)\n",
    "    f =  gdal.Open(tiff)\n",
    "    return f.GetGeoTransform()[1] \n",
    "\n",
    "def get_no_data_val(pth):\n",
    "    pth = str(pth)\n",
    "    f = gdal.Open(str(pth))\n",
    "    return np.nan if not f.GetRasterBand(1).GetNoDataValue() else f.GetRasterBand(1).GetNoDataValue()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "77ddeff4-7cd8-48e7-a967-f8a164286b09",
   "metadata": {},
   "outputs": [],
   "source": [
    "if predominant_epsg:\n",
    "    for ds in [vh_paths, vv_paths, local_inc_angle_paths, inc_angle_paths, ls_mask_paths]:\n",
    "        for pth in ds:\n",
    "            src_SRS = get_projection(str(pth))\n",
    "            if src_SRS != predominant_epsg:\n",
    "                res = get_res(pth)\n",
    "                \n",
    "                no_data_val = get_no_data_val(pth)\n",
    "                \n",
    "                temp = pth.parent/f\"temp_{pth.stem}.tif\"\n",
    "                pth.rename(temp)\n",
    "\n",
    "                warp_options = {\n",
    "                    \"dstSRS\":f\"EPSG:{predominant_epsg}\", \"srcSRS\":f\"EPSG:{src_SRS}\",\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()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "99ff9dae-cb7b-42d7-8bc1-9280bdb915d2",
   "metadata": {},
   "source": [
    "## **7. Merge bursts to create mosaics**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c6377fa6-ffa2-4ab5-8c72-f7e89532e2a8",
   "metadata": {
    "scrolled": true,
    "tags": []
   },
   "outputs": [],
   "source": [
    "ds_pth_dict = {\n",
    "    'vh': vh_paths,\n",
    "    'vv': vv_paths,\n",
    "    'inc_angle': inc_angle_paths,\n",
    "    'local_inc_angle': local_inc_angle_paths,\n",
    "    'ls_mask': ls_mask_paths\n",
    "}\n",
    "\n",
    "for ds in tqdm(ds_pth_dict):\n",
    "    merge_str = ''\n",
    "    for pth in ds_pth_dict[ds]:\n",
    "        merge_str = f\"{merge_str} {str(pth)}\"\n",
    "\n",
    "    no_data_val = get_no_data_val(ds_pth_dict[ds][0])\n",
    "\n",
    "\n",
    "    output = scene_dir/f\"OPERA_RTC_v0.4_{ds}_{scene}_mosaic.tif\"\n",
    "    merge_command = f\"gdal_merge.py -n {no_data_val} -a_nodata {no_data_val} -o {output} {merge_str}\"\n",
    "    print(merge_command)\n",
    "    !{merge_command}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a9148de6-d4e1-4287-abd6-72fe4864bc02",
   "metadata": {},
   "source": [
    "## **8. Delete Bursts**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "946ce5cb-f9ec-4df2-af97-4eaabbacfdaa",
   "metadata": {},
   "outputs": [],
   "source": [
    "for d in [vh_dir, vv_dir, inc_angle_dir, local_inc_angle_dir, ls_mask_dir]:\n",
    "    shutil.rmtree(d)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "opera_calval_rtc [conda env:.local-opera_calval_rtc]",
   "language": "python",
   "name": "conda-env-.local-opera_calval_rtc-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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}