Notebook

Using Leafmap and Earthaccess to Explore OPERA DSWx-HLS Products.¶

Example below showcases a flooding event in the Indus Valley.¶

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 vizualization of OPERA data products for a user-specified area of interest (AOI).

OPERA DSWx-HLS info¶

see website https://www.jpl.nasa.gov/go/opera/products/dswx-product-suite

Import Libraries¶

In [ ]:

import earthaccess
import leafmap
import pandas as pd
import geopandas as gpd
from shapely import box
from datetime import datetime

Authentication¶

A NASA Earthdata Login 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.

In [ ]:

leafmap.nasa_data_login()

View NASA Earthdata datasets¶

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

In [ ]:

### View Earthdata datasets
earthdata_url = 'https://github.com/opengeos/NASA-Earth-Data/raw/main/nasa_earth_data.tsv'
earthdata_df = pd.read_csv(earthdata_url, sep='\t')
# earthdata_df.head()

View the available OPERA products¶

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.

In [ ]:

opera_df = earthdata_df[earthdata_df['ShortName'].str.contains('OPERA', case=False)]
# opera_df

Define an area of interest (AOI) and time period of interest (TOI)¶

Define an area of interest (AOI) for the flood event

In [ ]:

### This cell initializes the AOI and TOI.

AOI = (67.982054, 26.198739,68.543065, 28.568858) #W, S, E, N; Indus Valley, Pakistan


#Here we have selected two dates. This could expand to include date ranges but then image mosaic rules should be considered (not included here)
StartDate_PreFlood="2023-05-03T00:00:00"  #Pre-flood image start date
EndDate_PreFlood="2023-05-03T23:59:59"    #Pre-flood image end date

StartDate_SynFlood="2023-08-07T00:00:00"  #Syn-flood image start date
EndDate_SynFlood="2023-08-07T23:59:59"    #Syn-flood image end date

Query Earthdata and return metadata for OPERA products within the AOI¶

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. 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.

View OPERA Product Shortnames¶

In [ ]:

### Print the available OPERA datasets 
print('Available OPERA datasets:', opera_df['ShortName'].values)

Query the OPERA DSWx-HLS dataset for the AOI¶

In [ ]:

dswx_results_PreFlood, dswx_gdf_PreFlood = leafmap.nasa_data_search(
    short_name='OPERA_L3_DSWX-HLS_V1',
    cloud_hosted=True,
    bounding_box= AOI,
    temporal=(StartDate_PreFlood, EndDate_PreFlood),
    count=-1,  # use -1 to return all datasets
    return_gdf=True,
)

dswx_results_SynFlood, dswx_gdf_SynFlood = leafmap.nasa_data_search(
    short_name='OPERA_L3_DSWX-HLS_V1',
    cloud_hosted=True,
    bounding_box= AOI,
    temporal=(StartDate_SynFlood, EndDate_SynFlood),
    count=-1,  # use -1 to return all datasets
    return_gdf=True,
)

See the available DSWx-HLS layers¶

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.

In [ ]:

dswx_results_PreFlood[0] #Note this just shows a single MGRS/HLS tile

In [ ]:

dswx_results_SynFlood[0] #Note this just shows a single MGRS/HLS tile

View the DSWx-HLS metadata and footprints¶

In [ ]:

dswx_gdf_PreFlood.head()

In [ ]:

### Plot the location of the tiles 
dswx_gdf_PreFlood.explore(fill=False)

In [ ]:

### Plot the location of the tiles 
dswx_gdf_SynFlood.explore(fill=False)

Download data with leafmap¶

Let's download the data from one of our above queries. In the cell below we specify data from the DSWx-HLS.

Create a subdirectory¶

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.

In [ ]:

import os
from datetime import datetime

def create_data_directory():
    # Get the current date and time
    # current_datetime = datetime.now().strftime("%m_%d_%Y_%H_%M_%S")
    current_datetime = datetime.now().strftime("%m_%d_%Y")

    # Define the base directory
    base_directory = "data"

    # Create the full path for the new directory
    new_directory_path_PreFlood = os.path.join(base_directory, f"data_{current_datetime}/PreFlood")
    # Create the new directory
    os.makedirs(new_directory_path_PreFlood, exist_ok=True)

    print(f"Directory '{new_directory_path_PreFlood}' created successfully.")

    return new_directory_path_PreFlood 

directory_path_PreFlood = create_data_directory()

In [ ]:

def create_data_directory():
    # Get the current date and time
    # current_datetime = datetime.now().strftime("%m_%d_%Y_%H_%M_%S")
    current_datetime = datetime.now().strftime("%m_%d_%Y")

    # Define the base directory
    base_directory = "data"

    # Create the full path for the new directory
    new_directory_path_SynFlood = os.path.join(base_directory, f"data_{current_datetime}/SynFlood")
    # Create the new directory
    os.makedirs(new_directory_path_SynFlood, exist_ok=True)

    print(f"Directory '{new_directory_path_SynFlood}' created successfully.")

    return new_directory_path_SynFlood 

directory_path_SynFlood = create_data_directory()

Download the data¶

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.

In [ ]:

dswx_data_PreFlood = leafmap.nasa_data_download(dswx_results_PreFlood, out_dir=directory_path_PreFlood)     

In [ ]:

dswx_data_SynFlood = leafmap.nasa_data_download(dswx_results_SynFlood, out_dir=directory_path_SynFlood)     

View the files using Leafmap¶

Load in images from data folder¶

We load in data from only the DSWx-WTR layer below. If you'd like load data from a different layer change the B01 to suit your needs. Included layers:

OPERA_L3_DSWx-HLS_*B01_WTR.tif

OPERA_L3_DSWx-HLS_*B02_BWTR.tif

OPERA_L3_DSWx-HLS_*B03_CONF.tif

OPERA_L3_DSWx-HLS_*B04_DIAG.tif

OPERA_L3_DSWx-HLS_*B05_WTR-1.tif

OPERA_L3_DSWx-HLS_*B06_WTR-2.tif

OPERA_L3_DSWx-HLS_*B07_LAND.tif

OPERA_L3_DSWx-HLS_*B08_SHAD.tif

OPERA_L3_DSWx-HLS_*B09_CLOUD.tif

OPERA_L3_DSWx-HLS_*B10_DEM.tif

In [ ]:

import os

ImageLayer='B01' #B01 corresponds to WTR (see above)

# Get the current directory
current_directory = os.getcwd()

# Construct the path to the data directory
data_directory_PreFlood = os.path.join(current_directory, directory_path_PreFlood)
data_directory_SynFlood = os.path.join(current_directory, directory_path_SynFlood)

# Create a list of file paths and a list of corresponding dates
images_PreFlood = [os.path.join(data_directory_PreFlood, filename) for filename in os.listdir(data_directory_PreFlood) if os.path.isfile(os.path.join(data_directory_PreFlood, filename)) and ImageLayer in filename]
image_dates_PreFlood = [image[25:33] for image in os.listdir(data_directory_PreFlood) if ImageLayer in image]

images_SynFlood = [os.path.join(data_directory_SynFlood, filename) for filename in os.listdir(data_directory_SynFlood) if os.path.isfile(os.path.join(data_directory_SynFlood, filename)) and ImageLayer in filename]
image_dates_SynFlood = [image[25:33] for image in os.listdir(data_directory_SynFlood) if ImageLayer in image]

Merge individual tiles into a single image¶

In [ ]:

filename_merged_PreFlood='PreFlood_Merged.tif'
merged_raster_PreFlood = leafmap.merge_rasters(data_directory_PreFlood,os.path.join(data_directory_PreFlood, filename_merged_PreFlood),input_pattern='*' + ImageLayer +'*.tif',output_format='GTiff',output_nodata=None)

In [ ]:

filename_merged_SynFlood='SynFlood_Merged.tif'
merged_raster_SynFlood = leafmap.merge_rasters(data_directory_SynFlood,os.path.join(data_directory_SynFlood, filename_merged_SynFlood),input_pattern= '*' + ImageLayer +'*.tif',output_format='GTiff',output_nodata=None)

Display the merged images¶

In [ ]:

m = leafmap.Map(basemap="Esri.WorldImagery")
# m.add_raster(os.path.join(data_directory_PreFlood, filename_merged_PreFlood), opacity=1)
m.add_raster(os.path.join(data_directory_SynFlood, filename_merged_SynFlood), opacity=1)
legend_dict = {
    'Not Water': '##ffffff',
    'Open Surface Water': '#0000ff',
    'Partial Surface Water': '#b4d5f4',
    'HLS snow/ice mask': '#00ffff',
    'HLS cloud/cloud shadow mask': '#afafaf'
}
# Add the legend to the map
m.add_legend(legend_title="Legend Title", legend_dict=legend_dict)

m

Create a split map to show changes¶

In [ ]:

raster_PreFlood_path=os.path.join(data_directory_PreFlood, filename_merged_PreFlood)
raster_SynFlood_path=os.path.join(data_directory_SynFlood, filename_merged_SynFlood)

In [ ]:

raster_path_2merged=[raster_PreFlood_path,raster_SynFlood_path]
#for only 2 dates - will need to update if more dates are used in merge
image_dates_merged=[image_dates_PreFlood[0],image_dates_SynFlood[0]]

In [ ]:

leafmap.split_map(
    left_layer=raster_PreFlood_path,
    right_layer=raster_SynFlood_path,
    # opacity=0.5,
    left_label="First",
    right_label="Last",
    label_position="bottom",
    basemap="Esri.WorldImagery",
    zoom=10,
)

#have not yet figured out how to get legend on the split map
# legend_dict = {
#     'Not Water': '##ffffff',
#     'Open Surface Water': '#0000ff',
#     'Partial Surface Water': '#b4d5f4',
#     'HLS snow/ice mask': '#00ffff',
#     'HLS cloud/cloud shadow mask': '#afafaf'
# }

Reclassify values¶

We have found that sometimes sediment rich water is classified as snow/ice in the HLS FMask. If the user is certain there is no snow/ice in the imagery, this layer can be reclassified.

In [ ]:

#Not yet set up

Conclusions¶

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.