Notebook

Getting Started with OPERA DSWx Product¶

This notebook serves as an introduction to the OPERA Dynamic Water eXtent (DSWx) product and visualizing relevant raster layers for reservoir monitoring applications. Note: This notebook uses provisional products, which may differ slightly from operational products. Please refer to DSWx product specification for more information.

In [1]:

# Notebook dependencies

import os
import xarray as xr
import numpy as np
import hvplot.xarray
import geoviews as gv
import pyproj
from pyproj import Proj
import holoviews as hv
from matplotlib.colors import ListedColormap

from bokeh.models import FixedTicker
hv.extension('bokeh')

import warnings
warnings.filterwarnings('ignore')

os.environ['AWS_NO_SIGN_REQUEST'] = 'YES'

Product Background¶

The DSWx product products map pixel-wise surface water detections using optical or SAR imagery. This notebook focuses on relevant layers in optical DSWx products derived from the Harmonized Landsat-8 and Sentinel-2 A/B (HLS) multispectral data.

Metadata¶

HLS products provide surface reflectance (SR) data from the Operational Land Imager (OLI) aboard the Landsat-8 remote sensing satellite and the Multi-Spectral Instrument (MSI) aboard the Sentinel-2 A/B remote sensing satellite. HLS products are distributed over projected map coordinates aligned with the Military Grid Reference System (MGRS). Each tile covers 109.8 kilometers squared divided into 3660 rows and 3660 columns at 30 meter pixel spacing. Each tile overlaps neighbors by 4900 meters in each direction.

Raster Layers¶

The DSWx product is distributed as a set of Cloud-Optimized GeoTIFF (COG) files to enable download of only particular layers of interest to a given user. All L3 DSWx layers are stored in files following GeoTIFF format specifications.

In [2]:

# User-specified input
data_dir = 's3://opera-pst-rs-pop1/products/OPERA_L3_DSWx-HLS_T11SQA_20201002T182211Z_20230118T203655Z_S2A_30_v0.0/OPERA_L3_DSWx-HLS_T11SQA_20201002T182211Z_20230118T203655Z_S2A_30_v0.0_'
bandlist = ['B01_WTR', 'B02_BWTR', 'B05_WTR-1', 'B06_WTR-2']
bandpath = f"{data_dir}%s.tif"

In [3]:

# Functions
def stack_bands(bandpath:str, bandlist:list): 
    '''
    Returns geocube with four bands stacked into one multi-dimensional array.
            Parameters:
                    bandpath (str): Path to bands that should be stacked
                    bandlist (list): Four bands that should be stacked
            Returns:
                    bandStack (xarray Dataset): Geocube with stacked bands
                    crs (int): Coordinate Reference System corresponding to bands
    '''
    bandStack = []; bandS = []; bandStack_ = [];
    for i,band in enumerate(bandlist):
        if i==0:
            bandStack_ = xr.open_rasterio(bandpath%band)
            crs = pyproj.CRS.to_epsg(pyproj.CRS.from_proj4(bandStack_.crs))
            bandStack_ = bandStack_ * bandStack_.scales[0]
            bandStack = bandStack_.squeeze(drop=True)
            bandStack = bandStack.to_dataset(name='z')
            bandStack.coords['band'] = i+1
            bandStack = bandStack.rename({'x':'longitude', 'y':'latitude', 'band':'band'})
            bandStack = bandStack.expand_dims(dim='band')  
        else:
            bandS = xr.open_rasterio(bandpath%band)
            bandS = bandS * bandS.scales[0]
            bandS = bandS.squeeze(drop=True)
            bandS = bandS.to_dataset(name='z')
            bandS.coords['band'] = i+1
            bandS = bandS.rename({'x':'longitude', 'y':'latitude', 'band':'band'})
            bandS = bandS.expand_dims(dim='band')
            bandStack = xr.concat([bandStack, bandS], dim='band')
    return bandStack, crs

In [4]:

# Creates geocube of stacked bands
dswx, crs = stack_bands(bandpath, bandlist)

# Creates basemap
base = gv.tile_sources.EsriImagery.opts(width=1000, height=1000, padding=0.1)

Band 1: Water classification (WTR)¶

Data Type: UInt8
Description: Masked interpreted water classification layer. This represents pixel-wise classification into one of three water classes (not water, open water, and partial surface water), cloud/cloud shadow class, or no data classes.

In [5]:

# Defines colormap for visualization
levels = [0, 0.9, 1.9, 2.9, 7.9, 8.9, 10]
color_key = {
    "Not Water": "#ffffff",
    "Open Water": "#0000ff",
    "Partial Surface Water": "#00ff00",
    "Reserved": "#000000",
    "Snow/Ice": "#00ffff",
    "Clouds/Cloud Shadow": "#7f7f7f"
}

ticks = [0.5, 1.5, 2.5, 5.5, 8.5, 9.5]
ticker = FixedTicker(ticks=ticks)
labels = dict(zip(ticks, color_key))

wtr = dswx.z.where(dswx['z']!=255).sel({'band':1})
wtr.where(wtr>0).hvplot.image(x='longitude', 
                          y='latitude', 
                          crs=crs, 
                          rasterize=True, 
                          dynamic=True, 
                          aspect='equal', 
                          frame_width=500, 
                          frame_height=500,  
                          clim=(0,10), alpha=0.8).opts(title=f"B01_WTR",xlabel='Longitude', 
                                    ylabel='Latitude',color_levels=levels,cmap=tuple(color_key.values()),
                                    colorbar_opts={'ticker':ticker,'major_label_overrides':labels}) * base

OMP: Info #276: omp_set_nested routine deprecated, please use omp_set_max_active_levels instead.

Out[5]:

Layer Values:

0: Not Water – an area with valid reflectance data that is not open water (class 1), partial surface water (class 2), or

cloud/cloud shadow (class 9). Masking can result in “not water” (class 0) where land cover masking is applied

1: Open Water – an area that is entirely water and unobstructed to the sensor, including obstructions by vegetation, terrain,

and buildings

2: Partial Surface Water – an area that is at least 50% and less than 100% open water. This may be referred to as “subpixel

inundation” when referring to a pixel’s area. Examples include inundated sinkholes, floating vegetation, and pixels bisected by coastlines

8: Snow/Ice - an area identified as snow/ice according to input HLS Fmask quality assurance (QA) data
9: Cloud/Cloud Shadow – an area identified as cloud, cloud shadow, or snow/ice according to input quality assurance (QA)

data

255: Fill value (no data)

Band 2: Binary water (BWTR)¶

Data Type: UInt8
Description: The binary water map is derived from the WTR layer as a union of water classes (open water and partial surface water) into a binary map indicating areas with and without water. This layer is meant to provide users with a quick view for water/no-water. Invalid data classes (cloud/cloud shadow and fill value) are also provided to indicate areas in which the binary classification does not provide water/no-water classification.

In [6]:

# Defines colormap for visualization
levels = [0, 0.9, 1.9, 7.9, 8.9, 10]
color_key = {
    "Not Water": "#ffffff",
    "Water": "#0000ff",
    "Reserved": "#000000",
    "Snow/Ice": "#00ffff",
    "Clouds/Cloud Shadow": "#7f7f7f"
}

ticks = [0.5, 1.5, 5.5, 8.5, 9.5]
ticker = FixedTicker(ticks=ticks)
labels = dict(zip(ticks, color_key))

bwtr = dswx.z.where(dswx['z']!=255).sel({'band':2})
bwtr.where(bwtr>0).hvplot.image(x='longitude', 
                          y='latitude', 
                          crs=crs, 
                          rasterize=True, 
                          dynamic=True, 
                          aspect='equal', 
                          frame_width=500, 
                          frame_height=500,  
                          clim=(0,10), alpha=0.8).opts(title=f"B02_BWTR",xlabel='Longitude', 
                                    ylabel='Latitude',color_levels=levels,cmap=tuple(color_key.values()),
                                    colorbar_opts={'ticker':ticker,'major_label_overrides':labels}) * base

Out[6]:

Layer Values:

0: Not water – an area with valid reflectance data that is not water (class 1) and not cloud/cloud shadow (class 9)
1: Water – an area classified as “open water” or “partial surface water” (see WTR layer)
8: Snow/Ice - an area identified as snow/ice according to input HLS Fmask quality assurance (QA) data
9: Cloud/Cloud Shadow – an area identified as cloud, cloud shadow, or snow/ice according to input QA data
255: Fill value (no data)

Band 5: Interpretation of diagnostic layer into water classes (WTR-1)¶

Data Type: UInt8
Description: Classification of DIAG layer results into open water, partial surface water, and no-water.

In [7]:

# Defines colormap for visualization
levels = [0, 0.6, 1.3, 2]
color_key = {
    "Not Water": "#ffffff",
    "Open Water": "#0000ff",
    "Partial Surface Water": "#00ff00",
}

ticks = [0.25, 0.9, 1.6]
ticker = FixedTicker(ticks=ticks)
labels = dict(zip(ticks, color_key))

wtr1 = dswx.z.where(dswx['z']!=255).sel({'band':3})
wtr1.where(wtr1>0).hvplot.image(x='longitude', 
                          y='latitude', 
                          crs=crs, 
                          rasterize=True, 
                          dynamic=True, 
                          aspect='equal', 
                          frame_width=500, 
                          frame_height=500,  
                          clim=(0,2), alpha=0.8).opts(title=f"B05_WTR-1",xlabel='Longitude', 
                                    ylabel='Latitude',color_levels=levels,cmap=tuple(color_key.values()),
                                    colorbar_opts={'ticker':ticker,'major_label_overrides':labels}) * base

Out[7]:

Layer Values:

0: Not Water - an area with valid reflectance data that is not open water (class 1) or partial surface water (class 2)
1: Open Water – an area that is entirely water and unobstructed to the sensor, including obstructions by vegetation, terrain,

and buildings

2: Partial Surface Water – an area that is at least 50% and less than 100% open water. This may be referred to as “subpixel

inundation” when referring to a pixel’s area. Examples include wetlands, water bodies with floating vegetation, and pixels bisected by coastlines

255: Fill value (no data)

Band 6: Interpreted layer refined using land cover and terrain shadow testing (WTR-2)¶

Data Type: UInt8
Description: The WTR-2 layer is derived from the WTR-1 (Layer 5) outcome by applying additional tests based on land cover and terrain shadow information as described in [RD1][RD2] to mask (eliminate) false-positive water detections.

In [8]:

# Defines colormap for visualization
levels = [0, 0.6, 1.3, 2]
color_key = {
    "Not Water": "#ffffff",
    "Open Water": "#0000ff",
    "Partial Surface Water": "#00ff00",
}

ticks = [0.25, 0.9, 1.6]
ticker = FixedTicker(ticks=ticks)
labels = dict(zip(ticks, color_key))

wtr2 = dswx.z.where(dswx['z']!=255).sel({'band':4})
wtr2.where(wtr2>0).hvplot.image(x='longitude', 
                          y='latitude', 
                          crs=crs, 
                          rasterize=True, 
                          dynamic=True, 
                          aspect='equal', 
                          frame_width=500, 
                          frame_height=500,  
                          clim=(0,2), alpha=0.8).opts(title=f"B06_WTR-2",xlabel='Longitude', 
                                    ylabel='Latitude',color_levels=levels,cmap=tuple(color_key.values()),
                                    colorbar_opts={'ticker':ticker,'major_label_overrides':labels}) * base

Out[8]:

Layer Values:

0: Not Water - an area with valid reflectance data that is not open water (class 1) or partial surface water (class 2)
1: Open Water – an area that is entirely water and unobstructed to the sensor, including obstructions by vegetation, terrain,

and buildings

2: Partial Surface Water – an area that is at least 50% and less than 100% open water. This may be referred to as “subpixel

inundation” when referring to a pixel’s area. Examples include wetlands, water bodies with floating vegetation, and pixels bisected by coastlines

255: Fill value (no data)