Notebook

Using Leafmap and Earthaccess to Explore OPERA RTC-S1 Products.¶

Example below showcases a May 2024 flooding event in the Porto Algre, Brazil.¶

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

OPERA RTC-S1 info¶

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

Import Libraries¶

In [ ]:

import earthaccess
import leafmap
import pandas as pd
import numpy as np
import rasterio 
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)]

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 = (-54.215, -30.766,-50.814, -28.938) #W, S, E, N; 2024 Brazil FLoods

# RTC-S1 image dates 2023-11-25 (pre flood), 2024-05-08, 2024-05-11, 2024-05-30, 2024-06-01, 2024-06-04

#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_1="2023-11-25T00:00:00"  #Pre-flood image start date
EndDate_1="2023-11-25T23:59:59"    #Pre-flood image end date

StartDate_2="2024-05-08T00:00:00"  #Syn-flood image start date
EndDate_2="2024-05-08T23: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 RTC-S1 dataset for the AOI¶

In [ ]:

rtc_s1_results_1, rtc_s1_gdf_1 = leafmap.nasa_data_search(
    short_name='OPERA_L2_RTC-S1_V1',
    cloud_hosted=True,
    bounding_box= AOI,
    temporal=(StartDate_1, EndDate_1),
    count=-1,  # use -1 to return all datasets
    return_gdf=True,
)

rtc_s1_results_2, rtc_s1_gdf_2 = leafmap.nasa_data_search(
    short_name='OPERA_L2_RTC-S1_V1',
    cloud_hosted=True,
    bounding_box= AOI,
    temporal=(StartDate_2, EndDate_2),
    count=-1,  # use -1 to return all datasets
    return_gdf=True,
)

See the available RTC-S1 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 [ ]:

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

In [ ]:

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

View the RTC-S1 metadata and footprints¶

In [ ]:

rtc_s1_gdf_1.head()

In [ ]:

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

In [ ]:

### Plot the location of the tiles 
rtc_s1_gdf_2.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 RTC-S1

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_1 = os.path.join(base_directory, f"data_{current_datetime}/1")
    # Create the new directory
    os.makedirs(new_directory_path_1, exist_ok=True)

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

    return new_directory_path_1

directory_path_1 = 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_2 = os.path.join(base_directory, f"data_{current_datetime}/2")
    # Create the new directory
    os.makedirs(new_directory_path_2, exist_ok=True)

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

    return new_directory_path_2

directory_path_2 = 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 [ ]:

rtc_s1_data_1 = leafmap.nasa_data_download(rtc_s1_results_1, out_dir=directory_path_1)     

In [ ]:

rtc_s1_data_2 = leafmap.nasa_data_download(rtc_s1_results_2, out_dir=directory_path_2)     

View the files using Leafmap¶

Load in images from data folder¶

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. Included layers:

OPERA_L2_RTC-S1_*_VH.tif

OPERA_L2_RTC-S1_*_VV.tif

OPERA_L2_RTC-S1_*_mask.tif

OPERA_L2_RTC-S1_*.h5

In [ ]:

import os

ImageLayer='VH' #select poloarizition, either VH or VV for most places

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

# Construct the path to the data directory
data_directory_1 = os.path.join(current_directory, directory_path_1)
data_directory_2 = os.path.join(current_directory, directory_path_2)

# Create a list of file paths and a list of corresponding dates
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]
image_dates_1 = [image[25:33] for image in os.listdir(data_directory_1) if ImageLayer in image]

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]
image_dates_2 = [image[25:33] for image in os.listdir(data_directory_2) if ImageLayer in image]

Merge individual tiles into a single image¶

In [ ]:

filename_merged_1='OPERA_RTC_S1_mosaic_1.tif'
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)

In [ ]:

filename_merged_2='OPERA_RTC_S1_mosaic_2.tif'
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)

Display the merged images¶

In [ ]:

m = leafmap.Map(basemap="Esri.WorldImagery")
# m.add_raster(os.path.join(data_directory_1, filename_merged_1), opacity=1)
m.add_raster(os.path.join(data_directory_2, filename_merged_2), opacity=1,vmin=0,vmax=0.1)


m

Convert backscatter from linear scale to decibels¶

In [ ]:

raster_1_path=os.path.join(data_directory_1, filename_merged_1)
raster_2_path=os.path.join(data_directory_2, filename_merged_2)

In [ ]:

with rasterio.open(raster_1_path) as ds:
    RTC_mosaic_1 = ds.read(1)
    out_profile_1 = ds.profile

In [ ]:

with rasterio.open(raster_2_path) as ds:
    RTC_mosaic_2 = ds.read(1)
    out_profile_2 = ds.profile

In [ ]:

# Convert rasters to dB scale to help compress dynamic range
RTC_mosaic_dB_1 = 10 * np.log10(RTC_mosaic_1)
RTC_mosaic_dB_2 = 10 * np.log10(RTC_mosaic_2)

In [ ]:

output_file = os.path.join(data_directory_1, "OPERA_RTC_S1_mosaic_dB_1.tif")
with rasterio.open(output_file, 'w', **out_profile_1) as dst:
    dst.write(RTC_mosaic_dB_1,1)

In [ ]:

output_file = os.path.join(data_directory_2, "OPERA_RTC_S1_mosaic_dB_2.tif")
with rasterio.open(output_file, 'w', **out_profile_2) as dst:
    dst.write(RTC_mosaic_dB_2,1)

In [ ]:

m = leafmap.Map(basemap="Esri.WorldImagery")
# m.add_raster(os.path.join(data_directory_1, filename_merged_1), opacity=1)
m.add_raster(os.path.join(data_directory_2, "OPERA_RTC_S1_mosaic_dB_2.tif"), opacity=1,vmin=-27,vmax=-10)

m

Create a split map to show changes¶

In [ ]:

raster_path_2merged=[RTC_mosaic_dB_1,RTC_mosaic_dB_2]
#for only 2 dates - will need to update if more dates are used in merge
image_dates_merged=[image_dates_1[0],image_dates_2[0]]

In [ ]:

m = leafmap.Map(basemap="Esri.WorldImagery", zoom = 8,)
m.split_map(
     left_layer=os.path.join(data_directory_1, "OPERA_RTC_S1_mosaic_dB_1.tif"),
     right_layer=os.path.join(data_directory_2, "OPERA_RTC_S1_mosaic_dB_2.tif"),
 )

# m = leafmap.Map(basemap="Esri.WorldImagery")
m.add_raster(os.path.join(data_directory_2, "OPERA_RTC_S1_mosaic_dB_2.tif"), opacity=1,vmin=-27,vmax=-10)
m.add_raster(os.path.join(data_directory_1, "OPERA_RTC_S1_mosaic_dB_1.tif"), opacity=1,vmin=-27,vmax=-10)


m

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.