Notebook

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

This tutorial demonstrates how to query and work with the OPERA DSWx-S1 Provisional Data products locally.¶

Data Used in the Example:¶

30 meter (m) global OPERA Dynamic Surface Water Extent from Sentinel-1A/B provisional product (Version 0.1)
- This dataset contains Level-3 Dynamic OPERA surface water extent from Sentinel-1A/B provisional product version 0.1. The data are provisional non-validated surface water extent observations over selective locations and times spanning from November 2023 to May 2024. The input dataset for generating each product is the the OPERA Level-2 Radiometric Terrain-Corrected product from Sentinel-1A/B data (RTC-S1). The OPERA RTC-S1 products provide the radiometric terrain corrected γ0 backscatter coefficient in geocoded burst domain, approximately every 6 to 12 days. The OPERA RTC-S1 products are produced and distributed in the Universal Transverse Mercator (UTM) coordinate System for each single burst. The DSWx-S1 products will be produced over Military Grid Reference System (MGRS) tiles that each cover 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 3 GeoTIFF (Geographic Tagged Image File Format) files including bindary water, water classification, and associated confidence.
- Science Dataset (SDS) layers: 2024 Rio Grande do Sul flood
  - Starting in late April 2024, southern Brazil/northern Uruguay faced severe weather conditions with extreme rainfall. This led to significant and deadly flooding, which continued well into May [EO]. Here, we demonstrate how DSWx-S1 can be used to map inundation extent as a result of the flooding event.
  - B02_BWTR (Binary Water Layer)
  - B03_CONF (Confidence Layer)

Please refer to the OPERA Product Specification Document for details about the DSWx-S1 product.

1. Getting Started ¶

1.1 Import Packages ¶

In [ ]:

%load_ext autoreload
%autoreload 2

In [ ]:

import os
import sys
from datetime import datetime
import warnings

import folium
import geoviews as gv
import geopandas as gpd
import holoviews as hv
import numpy as np
from folium import plugins
import rasterio
from rioxarray.merge import merge_arrays
from shapely.geometry import Polygon, box, shape

# Holoviews and Geoviews extensions
hv.extension('bokeh')
gv.extension('bokeh', 'matplotlib')

# Local application/library-specific imports
import sys
project_root = '../../'
sys.path.append(project_root)
from src.dswx_utils import (
    colorize,
    getbasemaps,
    transform_data_for_folium,
    query_s3_bucket
)

# Suppress warnings
warnings.filterwarnings('ignore')

1.2 Set up Working Environment ¶

In [ ]:

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

2. Access data locally ¶

2.1 Initialize user-defined parameters ¶

In [ ]:

# USER-DEFINED PARAMETERS
aoi = box(-53.26679, -30.84813, -50.15623, -28.91014)  # West, South, East, North
start_date = datetime(2023, 11, 24)                    # in 2024-01-01 00:00:00 format
stop_date = datetime(2024, 5, 9)                       # in 2024-01-01 00:00:00 format
overlap_threshold = 10                                 # in percent

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

In [ ]:

# get list of BWTR layers
bwtr_ext = '_B02_BWTR.tif'

# Define the bucket name
bucket_name = 'opera-provisional-products'

# Get the list of files matching the wildcard pattern
bwtr_pattern = f'DSWx/DSWx_S1/beta/v0.1/Brazil_2024/*/*{bwtr_ext}'
bwtr_lyrs = query_s3_bucket(bucket_name, bwtr_pattern)

In [ ]:

# Create a dataframe, and pass only tiles that intersect with the user-defined bbox
search_dswx = gpd.GeoDataFrame(columns=['bwtr_name', 'TileID',
                                        'Sensor', 'Date', 'Coords', 'bbox',
                                        'SpatialOverlap', 'geometry'], crs='EPSG:4326')
for i in bwtr_lyrs:
    # get iterative names
    fname_base = os.path.basename(i.split(bwtr_ext)[0])
    fn = os.path.basename(i).split('_')
    ID = fn[3]
    sensor = fn[6]
    dat = datetime.strptime(fn[4][:8], '%Y%m%d')

    ## get geometry
    ds, _ = transform_data_for_folium(i)
    ds_bounds = ds.rio.bounds()
    x_max = ds_bounds[2]
    x_min = ds_bounds[0]
    y_max = ds_bounds[3]
    y_min = ds_bounds[1]
    # Create a Polygon geometry
    ds_polygon = Polygon([(x_min, y_min), (x_min, y_max),
                         (x_max, y_max), (x_max, y_min)])
    coords = [i for i in ds_polygon.exterior.coords]
    bbox = ds_polygon.bounds

    # get intersection
    intersected_geom = ds_polygon.intersection(aoi)
    spatial_overlap = (intersected_geom.area * 100) / ds_polygon.area

    # Query DSWx-S1 tiles based on spatial overlap with respect to defined AOI
    if spatial_overlap > overlap_threshold:
        # Create a dictionary for the new row
        new_row = {'bwtr_name': i,
                   'TileID': ID,
                   'Sensor': sensor,
                   'Date': dat,
                   'Coords': coords,
                   'bbox': bbox,
                   'SpatialOverlap': spatial_overlap,
                   'geometry': ds_polygon}
        # Convert the dictionary to a GeoDataFrame
        new_gdf = gpd.GeoDataFrame([new_row], crs='EPSG:4326')
        # Append the new GeoDataFrame to the existing one
        search_dswx = search_dswx.append(new_gdf, ignore_index=True)

In [ ]:

## View snippet of dataframe
search_dswx.head()

In [ ]:

## Print search information
# Total granules
print(f"{len(search_dswx)} granules after search filter, out of {len(bwtr_lyrs)} total")

# Check percent overlap values
min_index = search_dswx['SpatialOverlap'].idxmin()
print(f"tile {search_dswx.loc[min_index]['bwtr_name']} \n" + \
      f"has the minimum spatial overlap in the set of {min(search_dswx['SpatialOverlap']):.2f} \n\n")

max_index = search_dswx['SpatialOverlap'].idxmax()
print(f"tile {search_dswx.loc[max_index]['bwtr_name']} \n" + \
      f"has the maximum spatial overlap in the set of {max(search_dswx['SpatialOverlap']):.2f}")

In [ ]:

# Visualize the DSWx tile boundary and the user-defined bbox
granules_poly = gv.Polygons(search_dswx['geometry'], 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
intersects_geometry = aoi.__geo_interface__
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

3. Load and visualize the flood extent ¶

3.1.a Merge tiles for images captured before the flood event ¶

In [ ]:

# calculate difference in time between tile acquisition and specified start date
search_dswx['TimeDiff'] = (search_dswx['Date'] - start_date).abs()
index_of_closest_date = search_dswx['TimeDiff'].idxmin()

# determine closest acquisition date
capture_date = search_dswx.loc[index_of_closest_date]['Date']
pre_flood_tiles = search_dswx[search_dswx['Date'] == capture_date]

# capture relevant images to merge
pre_flood_bwrt = pre_flood_tiles['bwtr_name'].to_list()

# drop time delta column
search_dswx.drop(columns=['TimeDiff'], inplace=True)

In [ ]:

# merge all BWTR layers
transformed_BWTR_arr = []
for i in pre_flood_bwrt:
    transformed_arr, cm_BWTR = transform_data_for_folium(i)
    transformed_BWTR_arr.append(transformed_arr)
merged_BWTR_pre = merge_arrays(transformed_BWTR_arr)

# Check one of the DataArrays
merged_BWTR_pre

3.1.b Visualize merged "pre-flood" tiles using Folium ¶

In [ ]:

# get different expected input derivatives of the bounds
WSEN = merged_BWTR_pre.rio.bounds()
ESWN_bounds = [WSEN[2], WSEN[1], WSEN[0], WSEN[3]]
SWNE_bounds = [[WSEN[1], WSEN[0]], [WSEN[3], WSEN[2]]]

In [ ]:

# save pre-flooding image to file, with static water masked out
merged_BWTR_pre.rio.to_raster("merged_BWTR_preflood.tif")

In [ ]:

# Colorize the map using predefined colors from DSWx for Folium display
colored_BWTR_pre,cmap_BWTR = colorize(merged_BWTR_pre[0], cmap=cm_BWTR)

In [ ]:

# Initialize Folium basemap
xmid = (SWNE_bounds[0][1]+SWNE_bounds[1][1])/2
ymid = (SWNE_bounds[0][0]+SWNE_bounds[1][0])/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_BWTR_pre, 
                                  opacity=0.6, 
                                  bounds=SWNE_bounds,
                                  mercator_project=True, # need to add to fix pixel offset
                                  name='Pre-flood water extent',
                                  show=True).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

3.2.a Merge tiles for images captured during the flood event ¶

In [ ]:

# calculate difference in time between tile acquisition and specified start date
search_dswx['TimeDiff'] = (search_dswx['Date'] - stop_date).abs()
index_of_closest_date = search_dswx['TimeDiff'].idxmin()

# determine closest acquisition date
capture_date = search_dswx.loc[index_of_closest_date]['Date']
flood_tiles = search_dswx[search_dswx['Date'] == capture_date]

# capture relevant images to merge
flood_bwrt = flood_tiles['bwtr_name'].to_list()

# drop time delta column
search_dswx.drop(columns=['TimeDiff'], inplace=True)

In [ ]:

# merge all BWTR layers
transformed_BWTR_arr = []
for i in flood_bwrt:
    transformed_arr, cm_BWTR = transform_data_for_folium(i)
    transformed_BWTR_arr.append(transformed_arr)
merged_BWTR_flood = merge_arrays(transformed_BWTR_arr)

# Check one of the DataArrays
merged_BWTR_flood

3.2.b Visualize merged flood tiles using Folium ¶

In [ ]:

WSEN = merged_BWTR_flood.rio.bounds()
ESWN_bounds = [WSEN[2], WSEN[1], WSEN[0], WSEN[3]]
SWNE_bounds = [[WSEN[1], WSEN[0]], [WSEN[3], WSEN[2]]]

In [ ]:

# save late flooding image to file, with static water masked out
merged_BWTR_flood.data[0][merged_BWTR_pre[0][:,:] == 255] = 255
merged_BWTR_flood.rio.to_raster("merged_BWTR_floodimage.tif")

In [ ]:

# Colorize the map using predefined colors from DSWx for Folium display
colored_BWTR_flood,cmap_BWTR = colorize(merged_BWTR_flood[0], cmap=cm_BWTR)

In [ ]:

# Compute change in flooding between early and late-flooding tiles
colored_BWTR_change = np.logical_and(colored_BWTR_flood == 1, \
                                     colored_BWTR_pre == 0).astype(int)

In [ ]:

# delete some variables to free up space
del merged_BWTR_flood, merged_BWTR_pre, colored_BWTR_pre

In [ ]:

# Initialize Folium basemap
xmid = (SWNE_bounds[0][1]+SWNE_bounds[1][1])/2
ymid = (SWNE_bounds[0][0]+SWNE_bounds[1][0])/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_BWTR_change, 
                                  opacity=0.6, 
                                  bounds=SWNE_bounds,
                                  mercator_project=True, # need to add to fix pixel offset
                                  name='Flood water extent',
                                  show=True).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

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

This tutorial demonstrates how to query and work with the OPERA DSWx-S1 Provisional Data products locally.¶

Data Used in the Example:¶

1. Getting Started ¶

1.1 Import Packages ¶

1.2 Set up Working Environment ¶

2. Access data locally ¶

2.1 Initialize user-defined parameters ¶

3. Load and visualize the flood extent ¶

3.1.a Merge tiles for images captured before the flood event ¶

3.1.b Visualize merged "pre-flood" tiles using Folium ¶

3.2.a Merge tiles for images captured during the flood event¶

3.2.b Visualize merged flood tiles using Folium ¶

3.2.a Merge tiles for images captured during the flood event ¶