{ "cells": [ { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "# Generate Flood Maps without downloading OPERA DSWx-HLS products locally\n", "## This tutorial demonstrates how to query and work with the OPERA DSWx-HLS products from the cloud ([OPERA_L3_DSWX-HLS_V1](https://dx.doi.org/10.5067/OPDSW-PL3V0)).\n", "\n", "--- \n", "\n", "### Data Used in the Example: \n", "\n", "- **30 meter (m) global OPERA Dynamic Surface Water Extent from Harmonized Landsat Sentinel-2A/B product (Version 1) - [OPERA_L3_DSWX-HLS_V1](https://dx.doi.org/10.5067/OPDSW-PL3V0)**\n", " - This dataset contains OPERA Level-3 Dynamic Surface Water Extent product version 1. The input dataset for generating each product is the Harmonized Landsat-8 and Sentinel-2A/B (HLS) product version 2.0. HLS products provide surface reflectance (SR) data from the Operational Land Imager (OLI) aboard the Landsat 8 satellite and the MultiSpectral Instrument (MSI) aboard the Sentinel-2A/B satellite. The surface water extent products are distributed over projected map coordinates using the Universal Transverse Mercator (UTM) projection. Each UTM tile covers an area of 109.8 km × 109.8 km. This area is divided into 3,660 rows and 3,660 columns at 30-m pixel spacing. Each product is distributed as a set of 10 GeoTIFF (Geographic Tagged Image File Format) files including water classification, associated confidence, land cover classification, terrain shadow layer, cloud/cloud-shadow classification, Digital elevation model (DEM), and Diagnostic layer.\n", " - **Science Dataset (SDS) layers: Pakistan Floods**\n", " - In 2022, Pakistan’s monsoon season produced significant rainfall, devastating floods and landslides, affecting all four of the country's provinces and ~14% of its population [[CDP](https://disasterphilanthropy.org/disasters/2022-pakistan-floods/)]. Here, we demonstrate how DSWx-HLS can be used to map inundation extent as a result of the monsoon event in Sep 2022.\n", " - B02_BWTR (Binary Water Layer) \n", " - B03_CONF (Confidence Layer) \n", "\n", "Please refer to the [OPERA Product Specification Document](https://d2pn8kiwq2w21t.cloudfront.net/documents/ProductSpec_DSWX_URS309746.pdf) for details about the DSWx-HLS product.\n", "\n", "---\n", "\n", "## Before Starting this Tutorial \n", "\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." ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "## 1. Getting Started " ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "### 1.1 Import Packages " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%load_ext autoreload\n", "%autoreload 2" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import os\n", "import sys\n", "import json\n", "from datetime import datetime\n", "from subprocess import Popen\n", "from platform import system\n", "from getpass import getpass\n", "\n", "import numpy as np\n", "import pandas as pd\n", "import geopandas as gpd\n", "from skimage import io\n", "\n", "import matplotlib.pyplot as plt\n", "from matplotlib import cm\n", "\n", "from shapely.geometry import box, shape\n", "from shapely.ops import transform\n", "\n", "from osgeo import gdal\n", "from rioxarray.merge import merge_arrays\n", "\n", "import pyproj\n", "from pyproj import Proj\n", "\n", "from netrc import netrc\n", "\n", "from pystac_client import Client, ItemSearch\n", "\n", "import folium\n", "from folium import plugins\n", "import geoviews as gv\n", "import hvplot.xarray\n", "import holoviews as hv\n", "hv.extension('bokeh')\n", "gv.extension('bokeh', 'matplotlib')\n", "\n", "import tqdm\n", "\n", "sys.path.append('../../')\n", "from src.dswx_utils import intersection_percent, colorize, getbasemaps, transform_data_for_folium\n", "\n", "import warnings\n", "warnings.filterwarnings('ignore')" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "### 1.2 Set up Working Environment " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "inDir = os.getcwd()\n", "os.chdir(inDir)" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "### 1.3 Generate Authentication Token " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Generates authentication token\n", "# Asks for your Earthdata username and password for first time, if netrc does not exists in your home directory.\n", "\n", "urs = 'urs.earthdata.nasa.gov' # Earthdata URL endpoint for authentication\n", "prompts = ['Enter NASA Earthdata Login Username: ',\n", " 'Enter NASA Earthdata Login Password: ']\n", "\n", "# Determine the OS (Windows machines usually use an '_netrc' file)\n", "netrc_name = \"_netrc\" if system()==\"Windows\" else \".netrc\"\n", "\n", "# Determine if netrc file exists, and if so, if it includes NASA Earthdata Login Credentials\n", "try:\n", " netrcDir = os.path.expanduser(f\"~/{netrc_name}\")\n", " netrc(netrcDir).authenticators(urs)[0]\n", "\n", "# Below, create a netrc file and prompt user for NASA Earthdata Login Username and Password\n", "except FileNotFoundError:\n", " homeDir = os.path.expanduser(\"~\")\n", " Popen('touch {0}{2} | echo machine {1} >> {0}{2}'.format(homeDir + os.sep, urs, netrc_name), shell=True)\n", " Popen('echo login {} >> {}{}'.format(getpass(prompt=prompts[0]), homeDir + os.sep, netrc_name), shell=True)\n", " Popen('echo \\'password {} \\'>> {}{}'.format(getpass(prompt=prompts[1]), homeDir + os.sep, netrc_name), shell=True)\n", " # Set restrictive permissions\n", " Popen('chmod 0600 {0}{1}'.format(homeDir + os.sep, netrc_name), shell=True)\n", "\n", " # Determine OS and edit netrc file if it exists but is not set up for NASA Earthdata Login\n", "except TypeError:\n", " homeDir = os.path.expanduser(\"~\")\n", " Popen('echo machine {1} >> {0}{2}'.format(homeDir + os.sep, urs, netrc_name), shell=True)\n", " Popen('echo login {} >> {}{}'.format(getpass(prompt=prompts[0]), homeDir + os.sep, netrc_name), shell=True)\n", " Popen('echo \\'password {} \\'>> {}{}'.format(getpass(prompt=prompts[1]), homeDir + os.sep, netrc_name), shell=True)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# GDAL configurations used to successfully access PODAAC Cloud Assets via vsicurl \n", "gdal.SetConfigOption('GDAL_HTTP_COOKIEFILE','~/cookies.txt')\n", "gdal.SetConfigOption('GDAL_HTTP_COOKIEJAR', '~/cookies.txt')\n", "gdal.SetConfigOption('GDAL_DISABLE_READDIR_ON_OPEN','EMPTY_DIR')\n", "gdal.SetConfigOption('CPL_VSIL_CURL_ALLOWED_EXTENSIONS','TIF, TIFF')" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "## 2. CMR-STAC API: Search for data based on spatial query " ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "### 2.1 Initialize user-defined parameters " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# USER-DEFINED PARAMETERS\n", "aoi = box(67.4, 26.2, 68.0, 27.5)\n", "start_date = datetime(2022, 1, 1) # in 2022-01-01 00:00:00 format\n", "stop_date = f\"{datetime.today().strftime('%Y-%m-%d')} 23:59:59\" # in 2022-01-01 00:00:00 format\n", "overlap_threshold = 10 # in percent\n", "\n", "print(f\"Search between {start_date} and {stop_date}\")\n", "print(f\"With AOI: {aoi.__geo_interface__}\")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Search data through CMR-STAC API\n", "stac = 'https://cmr.earthdata.nasa.gov/cloudstac/' # CMR-STAC API Endpoint\n", "api = Client.open(f'{stac}/POCLOUD/')\n", "collections = ['OPERA_L3_DSWX-HLS_V1']\n", "\n", "search_params = {\"collections\": collections,\n", " \"intersects\": aoi.__geo_interface__,\n", " \"datetime\": [start_date, stop_date],\n", " \"max_items\": 1000}\n", "search_dswx = api.search(**search_params)" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "### 2.2 Query DSWx-HLS tiles based on spatial overlap with respect to defined AOI " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Filter datasets based on spatial overlap \n", "intersects_geometry = aoi.__geo_interface__\n", "\n", "# Check percent overlap values\n", "print(\"Percent overlap before filtering: \")\n", "print([f\"{intersection_percent(i, intersects_geometry):.2f}\" for i in search_dswx.items()])\n", "\n", "# Apply spatial overlap and cloud cover filter\n", "dswx_filtered = (\n", " i for i in search_dswx.items() if (intersection_percent(i, intersects_geometry) > overlap_threshold)\n", ")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Inspect the items inside the filtered query\n", "dswx_data = list(dswx_filtered)\n", "# Inspect one data\n", "dswx_data[0].to_dict()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "## Print search information\n", "# Total granules\n", "print(f\"Total granules after search filter: {len(dswx_data)}\")\n", "\n", "# Check percent overlap values\n", "print(\"Percent overlap after filtering: \")\n", "print([f\"{intersection_percent(i, intersects_geometry):.2f}\" for i in dswx_data])" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Visualize the DSWx tile boundary and the user-defined bbox\n", "geom_df = []\n", "for d,_ in enumerate(dswx_data):\n", " geom_df.append(shape(dswx_data[d].geometry))\n", "\n", "geom_granules = gpd.GeoDataFrame({'geometry':geom_df})\n", "granules_poly = gv.Polygons(geom_granules, label='DSWx tile boundary').opts(line_color='blue', color=None, show_legend=True)\n", "\n", "# Use geoviews to combine a basemap with the shapely polygon of our Region of Interest (ROI)\n", "base = gv.tile_sources.EsriImagery.opts(width=1000, height=1000)\n", "\n", "# Get the user-specified aoi\n", "geom_aoi = shape(intersects_geometry)\n", "aoi_poly = gv.Polygons(geom_aoi, label='User-specified bbox').opts(line_color='yellow', color=None, show_legend=True)\n", "\n", "# Plot using geoviews wrapper\n", "granules_poly*base*aoi_poly" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "### 2.3 Create a table of search results " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Create table of search results\n", "dswx_data_df = []\n", "for item in dswx_data:\n", " item.to_dict()\n", " fn = item.id.split('_')\n", " ID = fn[3]\n", " sensor = fn[6]\n", " dat = item.datetime.strftime('%Y-%m-%d')\n", " spatial_overlap = intersection_percent(item, intersects_geometry)\n", " geom = item.geometry\n", " bbox = item.bbox\n", "\n", " # Take all the band href information \n", " band_links = [item.assets[links].href for links in item.assets.keys()]\n", " dswx_data_df.append([ID,sensor,dat,geom,bbox,spatial_overlap,band_links])\n", "\n", "dswx_data_df = pd.DataFrame(dswx_data_df, columns = ['TileID', 'Sensor', 'Date', 'Coords', 'bbox', 'SpatialOverlap', 'BandLinks'])\n", "dswx_data_df" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "## 3. Load and visualize the flood extent " ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "### 3.1 Load B02-BWTR (Binary Water Layer) and B03-CONF (Confidence Layer) " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Take one of the flooded dataset and check what files are included\n", "dswx_data_df.iloc[43].BandLinks" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "### 3.2 Merge tiles " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Get B02-BWTR layer for tiles acquired on 2022-09-30, project to folium's projection and merge tiles\n", "T42RUR_B02, T42RUR_B02_cm = transform_data_for_folium(dswx_data_df.iloc[42].BandLinks[1])\n", "T42RUQ_B02, T42RUQ_B02_cm = transform_data_for_folium(dswx_data_df.iloc[43].BandLinks[1])\n", "merged_B02 = merge_arrays([T42RUR_B02, T42RUQ_B02])\n", "\n", "# Get B03-CONF layer for tiles acquired on 2022-09-30, project to folium's projection and merge tiles\n", "T42RUR_B03, T42RUR_B03_cm = transform_data_for_folium(dswx_data_df.iloc[42].BandLinks[2])\n", "T42RUQ_B03, T42RUQ_B03_cm = transform_data_for_folium(dswx_data_df.iloc[43].BandLinks[2])\n", "merged_B03 = merge_arrays([T42RUR_B03, T42RUQ_B03])\n", "\n", "# Check one of the DataArrays\n", "merged_B02" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "### 3.3 Visualize merged tiles using Folium " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Colorize the map using predefined colors from DSWx for Folium display\n", "colored_B02,cmap_B02 = colorize(merged_B02[0], cmap=T42RUR_B02_cm)\n", "colored_B03,cmap_B03 = colorize(merged_B03[0], cmap=T42RUR_B03_cm) " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Initialize Folium basemap\n", "xmid =(merged_B02.x.values.min()+merged_B02.x.values.max())/2 ; ymid = (merged_B02.y.values.min()+merged_B02.y.values.max())/2\n", "m = folium.Map(location=[ymid, xmid], zoom_start=9, tiles='CartoDB positron', show=True)\n", "\n", "# Add custom basemaps\n", "basemaps = getbasemaps()\n", "for basemap in basemaps:\n", " basemaps[basemap].add_to(m)\n", "\n", "# Overlay B02 and B03 layers\n", "folium.raster_layers.ImageOverlay(colored_B02, \n", " opacity=0.6, \n", " bounds=[[merged_B02.y.values.min(),merged_B02.x.values.min()],[merged_B02.y.values.max(),merged_B02.x.values.max()]],\n", " name='Flooded Area',\n", " show=True).add_to(m)\n", "\n", "folium.raster_layers.ImageOverlay(colored_B03, \n", " opacity=0.8, \n", " bounds=[[merged_B03.y.values.min(),merged_B03.x.values.min()],[merged_B03.y.values.max(),merged_B03.x.values.max()]],\n", " name='Confidence Layer', \n", " show=False).add_to(m)\n", "\n", "# layer control\n", "m.add_child(folium.LayerControl())\n", "\n", "# Add fullscreen button\n", "plugins.Fullscreen().add_to(m)\n", "\n", "# Add inset minimap image\n", "minimap = plugins.MiniMap(width=300, height=300)\n", "m.add_child(minimap)\n", "\n", "# Mouse Position\n", "fmtr = \"function(num) {return L.Util.formatNum(num, 3) + ' º ';};\"\n", "plugins.MousePosition(position='bottomright', separator=' | ', prefix=\"Lat/Lon:\",\n", " lat_formatter=fmtr, lng_formatter=fmtr).add_to(m)\n", "\n", "# Display\n", "m" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "opera_app_dev", "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.1.-1" }, "orig_nbformat": 4, "vscode": { "interpreter": { "hash": "bff12bf762c2135bafcb19b3b536001906af3097d6f26a86554c7fc4e262c650" } } }, "nbformat": 4, "nbformat_minor": 2 }