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

import earthaccess
import leafmap
import rasterio
from rasterio.crs import CRS
from rasterio.plot import show
from rasterio.mask import mask
import pandas as pd
import geopandas as gpd
import matplotlib.pyplot as plt
import matplotlib.image as image
from PIL import Image
from shapely import box
from datetime import datetime
import numpy as np

leafmap.nasa_data_login()

### 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()

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

### This cell initializes the AOI and TOI.

AOI = f_bounds = (-118.179447, 34.182380,-118.164887,  34.208998) #W, S, E, N; This notebook looks at the Arroyo Secco Next to JPL
AOI_box = box(*AOI)

df_aoi = gpd.GeoDataFrame(geometry=[AOI_box], crs=CRS.from_epsg(4326))

# m = df_aoi.exterior.explore()

StartDate="2022-01-01T00:00:00"  #image start date
EndDate="2024-08-30T23:59:59"    #image end date

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

dswx_results, dswx_gdf = leafmap.nasa_data_search(
    short_name='OPERA_L3_DSWX-HLS_V1',
    cloud_hosted=True,
    bounding_box= AOI,
    temporal=(StartDate, EndDate),
    count=-1,  # use -1 to return all datasets
    return_gdf=True,
    # cloud_cover=10,
    # bbox_geom=AOI,
    )

dswx_gdf

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

dswx_gdf.head()

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

df_aoi.explore(fill=False)

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

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

    return new_directory_path 

directory_path = create_data_directory()

dswx_data = leafmap.nasa_data_download(dswx_results, out_dir=directory_path)     

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 = os.path.join(current_directory, directory_path)

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

images

filename_merged='Merged.tif'
merged_raster = leafmap.merge_rasters(data_directory,os.path.join(data_directory, filename_merged),input_pattern='*' + ImageLayer +'*.tif',output_format='GTiff',output_nodata=None)

m = leafmap.Map(basemap="Esri.WorldImagery")
m.add_raster(os.path.join(data_directory, filename_merged), opacity=1,nodata=0)
m.zoom_to_bounds(AOI)

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

with rasterio.open(images[0]) as src:
    # Print raster metadata
    print("Metadata:")
    print(src.meta)

from rasterio.warp import transform_bounds
from rasterio.crs import CRS

utm_bounds = transform_bounds(CRS.from_epsg(4326), src.crs, *f_bounds)
print(utm_bounds)

output_dir = os.path.join(data_directory, 'cropped')
os.makedirs(output_dir, exist_ok=True)

def crop_raster_to_roi(raster_path, utm_bounds, output_path):
    """Crop raster to the region of interest (ROI) and save it."""
    with rasterio.open(raster_path) as src:
        # Define bounding box from ROI
        bbox_geom = box(utm_bounds[0], utm_bounds[1], utm_bounds[2], utm_bounds[3])
        bbox_gdf = gpd.GeoDataFrame({'geometry': [bbox_geom]}, crs=src.crs)

        # Get the raster's transform and CRS
        transform = src.transform
        crs = src.crs

        # Reproject ROI to the raster's CRS
        bbox_gdf = bbox_gdf.to_crs(crs)

        # Mask the raster using the bounding box
        out_image, out_transform = mask(src, bbox_gdf.geometry, crop=True)
        out_meta = src.meta.copy()
        out_meta.update({
            'height': out_image.shape[1],
            'width': out_image.shape[2],
            'transform': out_transform
        })

        #  # Check for NaN values in the cropped image across all bands
        # if np.isnan(out_image).any():
        #     print(f'Skipping {os.path.basename(raster_path)} due to NaN values in the cropped area.')
        #     return

        # Check for NoData values in the cropped image
        nodata = src.nodata
        if nodata is not None:
            total_pixels = out_image.size
            nodata_pixels = (out_image == nodata).sum()
            nodata_percentage = (nodata_pixels / total_pixels) * 100

            if nodata_percentage > 90:
                print(f'Skipping {os.path.basename(raster_path)} due to more than 90% NoData values in the cropped area.')
                return

        # Check for Cloud/Cloud Shadow values in the cropped image
        cloud = 253 #cloud pixel values
        if cloud is not None:
            total_pixels = out_image.size
            cloud_pixels = (out_image == cloud).sum()
            cloud_percentage = (cloud_pixels / total_pixels) * 100

            if cloud_percentage > 50:
                print(f'Skipping {os.path.basename(raster_path)} due to more than 50% clouds in the cropped area.')
                return

        # Save the colormap
        colormap = src.colormap(1) if src.count == 1 else None

        # Save the cropped raster if no NaN values are found
        with rasterio.open(output_path, 'w', **out_meta) as dest:
            dest.write(out_image)
            if colormap:
                dest.write_colormap(1, colormap)

for raster_file in images:
    file_name = os.path.basename(raster_file)
    output_path = os.path.join(output_dir, file_name)
    crop_raster_to_roi(raster_file, utm_bounds, output_path)
    if os.path.exists(output_path):
        print(f'Cropped and saved {file_name} to {output_path}')
    else:
        print(f'Skipped {file_name}')

print('Processing complete.')

m = leafmap.Map(basemap="Esri.WorldImagery")
m.zoom_to_bounds(AOI)

# Define the directory containing your raster files and get a list of file paths
raster_files = sorted([os.path.join(output_dir, f) for f in os.listdir(output_dir) if f.endswith('.tif')])

# Optionally, you can define a dictionary of labels for the time slider
time_labels = {i: f'Image {i+1}' for i in range(len(raster_files))}

# Add the rasters with a time slider
m.add_time_slider(
    raster_files,
    time_interval=.25,  
    time_labels=time_labels,
    label_position="topcenter",
    nodata=0
)

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.add_gdf(df_aoi,fill_opacity=0, layer_name="AOI")

# Display the map
m