{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Using Leafmap and Earthaccess to Explore OPERA RTC-S1 Products. \n",
    "\n",
    "## Example below showcases a May 2024 flooding event in the Porto Algre, Brazil."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Leafmap library provides a suite of tools for interactive mapping and visualization in Jupyter Notebooks Leafmap version 0.30.0 and and later offer tools specifically for accessing NASA Earthdata by building on the newly developed NASA Earthaccess library. Earthaccess provides streamlined access to NASA Earthdata and simplifies the authentication and querying process over previously developed approaches.This notebook is designed to leverage tools within Earthaccess and Leafmap to facility easier access and visualization of OPERA data products for a user-specified area of interest (AOI). \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## OPERA RTC-S1 info\n",
    "see website https://www.jpl.nasa.gov/go/opera/products/rtc-product"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Import Libraries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import earthaccess\n",
    "import leafmap\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import rasterio \n",
    "import geopandas as gpd\n",
    "from shapely import box\n",
    "from datetime import datetime"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Authentication \n",
    "A [NASA Earthdata Login](https://urs.earthdata.nasa.gov/) account is required to download the data used in this tutorial. You can create an account at the link provided. After establishing an account, the code in the next cell will verify authentication. If this is your first time running the notebook, you will be prompted to enter your Earthdata login credentials, which will be saved in ~/.netrc."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "leafmap.nasa_data_login()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## View NASA Earthdata datasets\n",
    "A tab separated values (TSV) file, made available through the opengeos Github repository, catalogues metadata for more than 9,000 datasets available through NASA Earthdata. In the next cell we load the TSV into a pandas dataframe and view the metadata for the first five (5) Earthdata products"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### View Earthdata datasets\n",
    "earthdata_url = 'https://github.com/opengeos/NASA-Earth-Data/raw/main/nasa_earth_data.tsv'\n",
    "earthdata_df = pd.read_csv(earthdata_url, sep='\\t')\n",
    "# earthdata_df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## View the available OPERA products\n",
    "Note above that the `earthdata_df` contains a number of columns with metadata about each available product. the `ShortName` column will be used to produce a new dataframe containing only OPERA products. Let's view the available products and their metadata."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "opera_df = earthdata_df[earthdata_df['ShortName'].str.contains('OPERA', case=False)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Define an area of interest (AOI) and time period of interest (TOI)\n",
    "Define an area of interest (AOI) for the flood event"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### This cell initializes the AOI and TOI.\n",
    "\n",
    "AOI = (-54.215, -30.766,-50.814, -28.938) #W, S, E, N; 2024 Brazil FLoods\n",
    "\n",
    "# RTC-S1 image dates 2023-11-25 (pre flood), 2024-05-08, 2024-05-11, 2024-05-30, 2024-06-01, 2024-06-04\n",
    "\n",
    "#Here we have selected two dates. This could expand to include date ranges but then image mosaic rules should be considered (not included here)\n",
    "StartDate_1=\"2023-11-25T00:00:00\"  #Pre-flood image start date\n",
    "EndDate_1=\"2023-11-25T23:59:59\"    #Pre-flood image end date\n",
    "\n",
    "StartDate_2=\"2024-05-08T00:00:00\"  #Syn-flood image start date\n",
    "EndDate_2=\"2024-05-08T23:59:59\"    #Syn-flood image end date\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Query Earthdata and return metadata for OPERA products within the AOI\n",
    "The `earthaccess` library makes it simple to quickly query NASA's Common Metadata Repository (CMR) and return the associated metadata as a Geodataframe. `Leafmap` has recently added functionality that builds on `earthaccess` to enable interactive viewing of this data. \n",
    "In the next cell, the user should specify which OPERA product and the date range of interest. The AOI defined previously is used as the boundary in the query."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### View OPERA Product Shortnames"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Print the available OPERA datasets \n",
    "print('Available OPERA datasets:', opera_df['ShortName'].values)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Query the OPERA RTC-S1 dataset for the AOI\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rtc_s1_results_1, rtc_s1_gdf_1 = leafmap.nasa_data_search(\n",
    "    short_name='OPERA_L2_RTC-S1_V1',\n",
    "    cloud_hosted=True,\n",
    "    bounding_box= AOI,\n",
    "    temporal=(StartDate_1, EndDate_1),\n",
    "    count=-1,  # use -1 to return all datasets\n",
    "    return_gdf=True,\n",
    ")\n",
    "\n",
    "rtc_s1_results_2, rtc_s1_gdf_2 = leafmap.nasa_data_search(\n",
    "    short_name='OPERA_L2_RTC-S1_V1',\n",
    "    cloud_hosted=True,\n",
    "    bounding_box= AOI,\n",
    "    temporal=(StartDate_2, EndDate_2),\n",
    "    count=-1,  # use -1 to return all datasets\n",
    "    return_gdf=True,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### See the available RTC-S1 layers\n",
    "Functionality within earthaccess enables more more asthetic views of the available layers, as well as displaying the thumbnail. These links are clickable and will download in the browser when clicked. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rtc_s1_results_1[0] #Note this just shows a single MGRS/HLS tile"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rtc_s1_results_2[0] #Note this just shows a single MGRS/HLS tile"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### View the RTC-S1 metadata and footprints"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rtc_s1_gdf_1.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Plot the location of the tiles \n",
    "rtc_s1_gdf_1.explore(fill=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Plot the location of the tiles \n",
    "rtc_s1_gdf_2.explore(fill=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Download data with leafmap\n",
    "Let's download the data from one of our above queries. In the cell below we specify data from the RTC-S1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Create a subdirectory\n",
    "This will be where the files are downloaded. It will be a subdirectory inside of a directory called `data`, and the directory name will be the date that it was created."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "from datetime import datetime\n",
    "\n",
    "def create_data_directory():\n",
    "    # Get the current date and time\n",
    "    # current_datetime = datetime.now().strftime(\"%m_%d_%Y_%H_%M_%S\")\n",
    "    current_datetime = datetime.now().strftime(\"%m_%d_%Y\")\n",
    "\n",
    "    # Define the base directory\n",
    "    base_directory = \"data\"\n",
    "\n",
    "    # Create the full path for the new directory\n",
    "    new_directory_path_1 = os.path.join(base_directory, f\"data_{current_datetime}/1\")\n",
    "    # Create the new directory\n",
    "    os.makedirs(new_directory_path_1, exist_ok=True)\n",
    "\n",
    "    print(f\"Directory '{new_directory_path_1}' created successfully.\")\n",
    "\n",
    "    return new_directory_path_1\n",
    "\n",
    "directory_path_1 = create_data_directory()\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def create_data_directory():\n",
    "    # Get the current date and time\n",
    "    # current_datetime = datetime.now().strftime(\"%m_%d_%Y_%H_%M_%S\")\n",
    "    current_datetime = datetime.now().strftime(\"%m_%d_%Y\")\n",
    "\n",
    "    # Define the base directory\n",
    "    base_directory = \"data\"\n",
    "\n",
    "    # Create the full path for the new directory\n",
    "    new_directory_path_2 = os.path.join(base_directory, f\"data_{current_datetime}/2\")\n",
    "    # Create the new directory\n",
    "    os.makedirs(new_directory_path_2, exist_ok=True)\n",
    "\n",
    "    print(f\"Directory '{new_directory_path_2}' created successfully.\")\n",
    "\n",
    "    return new_directory_path_2\n",
    "\n",
    "directory_path_2 = create_data_directory()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Download the data\n",
    "The below will download the data to your newly created subdirectory. Look on your file system for a directory `/data/date/` where `date` is the date the directory was created."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rtc_s1_data_1 = leafmap.nasa_data_download(rtc_s1_results_1, out_dir=directory_path_1)     "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rtc_s1_data_2 = leafmap.nasa_data_download(rtc_s1_results_2, out_dir=directory_path_2)     "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## View the files using Leafmap"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Load in images from data folder\n",
    "We load in data from only the RTC-S1 layer below. If you'd like load data from a different layer change the to suit your needs. \n",
    "Included layers:\n",
    "\n",
    "\n",
    "OPERA_L2_RTC-S1_*_VH.tif\n",
    "\n",
    "\n",
    "OPERA_L2_RTC-S1_*_VV.tif\n",
    "\n",
    "OPERA_L2_RTC-S1_*_mask.tif\n",
    "\n",
    "OPERA_L2_RTC-S1_*.h5\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "\n",
    "ImageLayer='VH' #select poloarizition, either VH or VV for most places\n",
    "\n",
    "# Get the current directory\n",
    "current_directory = os.getcwd()\n",
    "\n",
    "# Construct the path to the data directory\n",
    "data_directory_1 = os.path.join(current_directory, directory_path_1)\n",
    "data_directory_2 = os.path.join(current_directory, directory_path_2)\n",
    "\n",
    "# Create a list of file paths and a list of corresponding dates\n",
    "images_1 = [os.path.join(data_directory_1, filename) for filename in os.listdir(data_directory_1) if os.path.isfile(os.path.join(data_directory_1, filename)) and ImageLayer in filename]\n",
    "image_dates_1 = [image[25:33] for image in os.listdir(data_directory_1) if ImageLayer in image]\n",
    "\n",
    "images_2 = [os.path.join(data_directory_2, filename) for filename in os.listdir(data_directory_2) if os.path.isfile(os.path.join(data_directory_2, filename)) and ImageLayer in filename]\n",
    "image_dates_2 = [image[25:33] for image in os.listdir(data_directory_2) if ImageLayer in image]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Merge individual tiles into a single image"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "filename_merged_1='OPERA_RTC_S1_mosaic_1.tif'\n",
    "merged_raster_1 = leafmap.merge_rasters(data_directory_1,os.path.join(data_directory_1, filename_merged_1),input_pattern='*' + ImageLayer +'*.tif',output_format='GTiff',output_nodata=None)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "filename_merged_2='OPERA_RTC_S1_mosaic_2.tif'\n",
    "merged_raster_2 = leafmap.merge_rasters(data_directory_2,os.path.join(data_directory_2, filename_merged_2),input_pattern= '*' + ImageLayer +'*.tif',output_format='GTiff',output_nodata=None)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Display the merged images"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m = leafmap.Map(basemap=\"Esri.WorldImagery\")\n",
    "# m.add_raster(os.path.join(data_directory_1, filename_merged_1), opacity=1)\n",
    "m.add_raster(os.path.join(data_directory_2, filename_merged_2), opacity=1,vmin=0,vmax=0.1)\n",
    "\n",
    "\n",
    "m"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Convert backscatter from linear scale to decibels"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "raster_1_path=os.path.join(data_directory_1, filename_merged_1)\n",
    "raster_2_path=os.path.join(data_directory_2, filename_merged_2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with rasterio.open(raster_1_path) as ds:\n",
    "    RTC_mosaic_1 = ds.read(1)\n",
    "    out_profile_1 = ds.profile"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with rasterio.open(raster_2_path) as ds:\n",
    "    RTC_mosaic_2 = ds.read(1)\n",
    "    out_profile_2 = ds.profile"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Convert rasters to dB scale to help compress dynamic range\n",
    "RTC_mosaic_dB_1 = 10 * np.log10(RTC_mosaic_1)\n",
    "RTC_mosaic_dB_2 = 10 * np.log10(RTC_mosaic_2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "output_file = os.path.join(data_directory_1, \"OPERA_RTC_S1_mosaic_dB_1.tif\")\n",
    "with rasterio.open(output_file, 'w', **out_profile_1) as dst:\n",
    "    dst.write(RTC_mosaic_dB_1,1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "output_file = os.path.join(data_directory_2, \"OPERA_RTC_S1_mosaic_dB_2.tif\")\n",
    "with rasterio.open(output_file, 'w', **out_profile_2) as dst:\n",
    "    dst.write(RTC_mosaic_dB_2,1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m = leafmap.Map(basemap=\"Esri.WorldImagery\")\n",
    "# m.add_raster(os.path.join(data_directory_1, filename_merged_1), opacity=1)\n",
    "m.add_raster(os.path.join(data_directory_2, \"OPERA_RTC_S1_mosaic_dB_2.tif\"), opacity=1,vmin=-27,vmax=-10)\n",
    "\n",
    "m\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create a split map to show changes\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "raster_path_2merged=[RTC_mosaic_dB_1,RTC_mosaic_dB_2]\n",
    "#for only 2 dates - will need to update if more dates are used in merge\n",
    "image_dates_merged=[image_dates_1[0],image_dates_2[0]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m = leafmap.Map(basemap=\"Esri.WorldImagery\", zoom = 8,)\n",
    "m.split_map(\n",
    "     left_layer=os.path.join(data_directory_1, \"OPERA_RTC_S1_mosaic_dB_1.tif\"),\n",
    "     right_layer=os.path.join(data_directory_2, \"OPERA_RTC_S1_mosaic_dB_2.tif\"),\n",
    " )\n",
    "\n",
    "# m = leafmap.Map(basemap=\"Esri.WorldImagery\")\n",
    "m.add_raster(os.path.join(data_directory_2, \"OPERA_RTC_S1_mosaic_dB_2.tif\"), opacity=1,vmin=-27,vmax=-10)\n",
    "m.add_raster(os.path.join(data_directory_1, \"OPERA_RTC_S1_mosaic_dB_1.tif\"), opacity=1,vmin=-27,vmax=-10)\n",
    "\n",
    "\n",
    "m"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Conclusions\n",
    "This is a first `earthaccess` and `leafmap` notebook for flood application. More work is needed to expand features for sophisticated filtering (cloud cover, spatial overlap) and analysis. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "opera_app",
   "language": "python",
   "name": "python3"
  },
  "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.10.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}