Generate Flood Maps without downloading OPERA DSWx-HLS products locally¶

%load_ext autoreload
%autoreload 2

import os
import sys
import json
from datetime import datetime
from subprocess import Popen
from platform import system
from getpass import getpass

import numpy as np
import pandas as pd
import geopandas as gpd
from skimage import io

import matplotlib.pyplot as plt
from matplotlib import cm

from shapely.geometry import box, shape
from shapely.ops import transform

from osgeo import gdal
from rioxarray.merge import merge_arrays

import pyproj
from pyproj import Proj

from netrc import netrc

from pystac_client import Client, ItemSearch

import folium
from folium import plugins
import geoviews as gv
import hvplot.xarray
import holoviews as hv
hv.extension('bokeh')
gv.extension('bokeh', 'matplotlib')

import tqdm

sys.path.append('../../')
from src.dswx_utils import intersection_percent, colorize, getbasemaps, transform_data_for_folium

import warnings
warnings.filterwarnings('ignore')

inDir = os.getcwd()
os.chdir(inDir)

# Generates authentication token
# Asks for your Earthdata username and password for first time, if netrc does not exists in your home directory.

urs = 'urs.earthdata.nasa.gov'    # Earthdata URL endpoint for authentication
prompts = ['Enter NASA Earthdata Login Username: ',
           'Enter NASA Earthdata Login Password: ']

# Determine the OS (Windows machines usually use an '_netrc' file)
netrc_name = "_netrc" if system()=="Windows" else ".netrc"

# Determine if netrc file exists, and if so, if it includes NASA Earthdata Login Credentials
try:
    netrcDir = os.path.expanduser(f"~/{netrc_name}")
    netrc(netrcDir).authenticators(urs)[0]

# Below, create a netrc file and prompt user for NASA Earthdata Login Username and Password
except FileNotFoundError:
    homeDir = os.path.expanduser("~")
    Popen('touch {0}{2} | echo machine {1} >> {0}{2}'.format(homeDir + os.sep, urs, netrc_name), shell=True)
    Popen('echo login {} >> {}{}'.format(getpass(prompt=prompts[0]), homeDir + os.sep, netrc_name), shell=True)
    Popen('echo \'password {} \'>> {}{}'.format(getpass(prompt=prompts[1]), homeDir + os.sep, netrc_name), shell=True)
    # Set restrictive permissions
    Popen('chmod 0600 {0}{1}'.format(homeDir + os.sep, netrc_name), shell=True)

    # Determine OS and edit netrc file if it exists but is not set up for NASA Earthdata Login
except TypeError:
    homeDir = os.path.expanduser("~")
    Popen('echo machine {1} >> {0}{2}'.format(homeDir + os.sep, urs, netrc_name), shell=True)
    Popen('echo login {} >> {}{}'.format(getpass(prompt=prompts[0]), homeDir + os.sep, netrc_name), shell=True)
    Popen('echo \'password {} \'>> {}{}'.format(getpass(prompt=prompts[1]), homeDir + os.sep, netrc_name), shell=True)

# GDAL configurations used to successfully access PODAAC Cloud Assets via vsicurl 
gdal.SetConfigOption('GDAL_HTTP_COOKIEFILE','~/cookies.txt')
gdal.SetConfigOption('GDAL_HTTP_COOKIEJAR', '~/cookies.txt')
gdal.SetConfigOption('GDAL_DISABLE_READDIR_ON_OPEN','EMPTY_DIR')
gdal.SetConfigOption('CPL_VSIL_CURL_ALLOWED_EXTENSIONS','TIF, TIFF')

# USER-DEFINED PARAMETERS
aoi = box(67.4, 26.2, 68.0, 27.5)
start_date = datetime(2022, 1, 1)                                       # in 2022-01-01 00:00:00 format
stop_date = f"{datetime.today().strftime('%Y-%m-%d')} 23:59:59"         # in 2022-01-01 00:00:00 format
overlap_threshold = 10                                                  # in percent

print(f"Search between {start_date} and {stop_date}")
print(f"With AOI: {aoi.__geo_interface__}")

# Search data through CMR-STAC API
stac = 'https://cmr.earthdata.nasa.gov/cloudstac/'    # CMR-STAC API Endpoint
api = Client.open(f'{stac}/POCLOUD/')
collections = ['OPERA_L3_DSWX-HLS_V1']

search_params = {"collections": collections,
                 "intersects": aoi.__geo_interface__,
                 "datetime": [start_date, stop_date],
                 "max_items": 1000}
search_dswx = api.search(**search_params)

# Filter datasets based on spatial overlap 
intersects_geometry = aoi.__geo_interface__

# Check percent overlap values
print("Percent overlap before filtering: ")
print([f"{intersection_percent(i, intersects_geometry):.2f}" for i in search_dswx.items()])

# Apply spatial overlap and cloud cover filter
dswx_filtered = (
    i for i in search_dswx.items() if (intersection_percent(i, intersects_geometry) > overlap_threshold)
)

# Inspect the items inside the filtered query
dswx_data = list(dswx_filtered)
# Inspect one data
dswx_data[0].to_dict()

## Print search information
# Total granules
print(f"Total granules after search filter: {len(dswx_data)}")

# Check percent overlap values
print("Percent overlap after filtering: ")
print([f"{intersection_percent(i, intersects_geometry):.2f}" for i in dswx_data])

# Visualize the DSWx tile boundary and the user-defined bbox
geom_df = []
for d,_ in enumerate(dswx_data):
    geom_df.append(shape(dswx_data[d].geometry))

geom_granules = gpd.GeoDataFrame({'geometry':geom_df})
granules_poly = gv.Polygons(geom_granules, label='DSWx tile boundary').opts(line_color='blue', color=None, show_legend=True)

# Use geoviews to combine a basemap with the shapely polygon of our Region of Interest (ROI)
base = gv.tile_sources.EsriImagery.opts(width=1000, height=1000)

# Get the user-specified aoi
geom_aoi = shape(intersects_geometry)
aoi_poly = gv.Polygons(geom_aoi, label='User-specified bbox').opts(line_color='yellow', color=None, show_legend=True)

# Plot using geoviews wrapper
granules_poly*base*aoi_poly

# Create table of search results
dswx_data_df = []
for item in dswx_data:
    item.to_dict()
    fn = item.id.split('_')
    ID = fn[3]
    sensor = fn[6]
    dat = item.datetime.strftime('%Y-%m-%d')
    spatial_overlap = intersection_percent(item, intersects_geometry)
    geom = item.geometry
    bbox = item.bbox

    # Take all the band href information 
    band_links = [item.assets[links].href for links in item.assets.keys()]
    dswx_data_df.append([ID,sensor,dat,geom,bbox,spatial_overlap,band_links])

dswx_data_df = pd.DataFrame(dswx_data_df, columns = ['TileID', 'Sensor', 'Date', 'Coords', 'bbox', 'SpatialOverlap', 'BandLinks'])
dswx_data_df

# Take one of the flooded dataset and check what files are included
dswx_data_df.iloc[43].BandLinks

# Get B02-BWTR layer for tiles acquired on 2022-09-30, project to folium's projection and merge tiles
T42RUR_B02, T42RUR_B02_cm = transform_data_for_folium(dswx_data_df.iloc[42].BandLinks[1])
T42RUQ_B02, T42RUQ_B02_cm = transform_data_for_folium(dswx_data_df.iloc[43].BandLinks[1])
merged_B02 = merge_arrays([T42RUR_B02, T42RUQ_B02])

# Get B03-CONF layer for tiles acquired on 2022-09-30, project to folium's projection and merge tiles
T42RUR_B03, T42RUR_B03_cm = transform_data_for_folium(dswx_data_df.iloc[42].BandLinks[2])
T42RUQ_B03, T42RUQ_B03_cm = transform_data_for_folium(dswx_data_df.iloc[43].BandLinks[2])
merged_B03 = merge_arrays([T42RUR_B03, T42RUQ_B03])

# Check one of the DataArrays
merged_B02

# Colorize the map using predefined colors from DSWx for Folium display
colored_B02,cmap_B02 = colorize(merged_B02[0], cmap=T42RUR_B02_cm)
colored_B03,cmap_B03 = colorize(merged_B03[0], cmap=T42RUR_B03_cm) 

# Initialize Folium basemap
xmid =(merged_B02.x.values.min()+merged_B02.x.values.max())/2 ; ymid = (merged_B02.y.values.min()+merged_B02.y.values.max())/2
m = folium.Map(location=[ymid, xmid], zoom_start=9, tiles='CartoDB positron', show=True)

# Add custom basemaps
basemaps = getbasemaps()
for basemap in basemaps:
    basemaps[basemap].add_to(m)

# Overlay B02 and B03 layers
folium.raster_layers.ImageOverlay(colored_B02, 
                                        opacity=0.6, 
                                        bounds=[[merged_B02.y.values.min(),merged_B02.x.values.min()],[merged_B02.y.values.max(),merged_B02.x.values.max()]],
                                        name='Flooded Area',
                                        show=True).add_to(m)

folium.raster_layers.ImageOverlay(colored_B03, 
                                        opacity=0.8, 
                                        bounds=[[merged_B03.y.values.min(),merged_B03.x.values.min()],[merged_B03.y.values.max(),merged_B03.x.values.max()]],
                                        name='Confidence Layer', 
                                        show=False).add_to(m)

# layer control
m.add_child(folium.LayerControl())

# Add fullscreen button
plugins.Fullscreen().add_to(m)

# Add inset minimap image
minimap = plugins.MiniMap(width=300, height=300)
m.add_child(minimap)

# Mouse Position
fmtr = "function(num) {return L.Util.formatNum(num, 3) + ' º ';};"
plugins.MousePosition(position='bottomright', separator=' | ', prefix="Lat/Lon:",
                     lat_formatter=fmtr, lng_formatter=fmtr).add_to(m)

# Display
m