Reverend Bayes updates our Belief in Flood Detection

How an 275 year old idea helps map the extent of floods

Image from wikipedia

Note

This notebook contains interactive elements. The full interactive elements can only be viewed on Binder by clicking on the Binder badge or 🚀 button.

Overview

This notebook explains how microwave (\(\sigma^0\)) backscattering can be used to map the extent of a flood. We replicate in this exercise the work of [1] on the TU Wien Bayesian-based flood mapping algorithm.

Prerequisites

Concepts

Importance

Notes

Intro to xarray

Necessary

Intro to Harmonic parameters

Necessary

Documentation hvPlot

Helpful

Interactive plotting

Documentation odc-stac

Helpful

Data access

  • Time to learn: 10 min

Imports

import datetime

import holoviews as hv
import hvplot.pandas
import hvplot.xarray
import numpy as np
import pandas as pd
import panel as pn
import pystac_client
import rioxarray  # noqa: F401
import xarray as xr
from bokeh.models import FixedTicker
from odc import stac as odc_stac
from scipy.stats import norm

pn.extension()
hv.extension("bokeh")

Greece Flooding 2018

In this exercise we will replicate the case study of the above mentioned paper, the February 2018 flooding of the Greek region of Thessaly.

time_range = "2018-02-28T04:00:00Z/2018-02-28T05:00:00Z"
minlon, maxlon = 21.93, 22.23
minlat, maxlat = 39.47, 39.64
bounding_box = [minlon, minlat, maxlon, maxlat]

EODC STAC Catalog

The data required for TU Wien flood mapping algorithm consists of terrain corrected sigma naught backscatter data \(\sigma^{0}\), the projected local incidence angle (PLIA) values of those measurements, and the harmonic parameters (HPAR) of a model fit on the pixel’s backscatter time series. The latter two datasets will needed to calculate the probability density functions over land and water for. We will be getting the required data from the EODC STAC Catalog. Specifically the collections: SENTINEL_SIG0_20M, SENTINEL1_MPLIA and SENTINEL1_HPAR. We use the pystac-client and odc_stac packages to, respectively, discover and fetch the data.

Due to the way the data is acquired and stored, some items include “no data” areas. In our case, no data has the value -9999, but this can vary from data provider to data provider. This information can usually be found in the metadata. Furthermore, to save memory, data is often stored as integer (e.g. 25) and not in float (e.g. 2.5) format. For this reason, the backscatter values are often multiplied by a scale factor. Hence we define the function post_process_eodc_cube to correct for these factors as obtained from the STAC metadata.

Sigma naught

eodc_catalog = pystac_client.Client.open("https://stac.eodc.eu/api/v1")
search = eodc_catalog.search(
    collections="SENTINEL1_SIG0_20M",
    bbox=bounding_box,
    datetime=time_range,
)
items_sig0 = search.item_collection()


def post_process_eodc_cube(dc, items, bands):
    """
    Postprocessing of EODC data cubes.

    Parameters
    ----------
    x : xarray.Dataset
    items: pystac.item_collection.ItemCollection
        STAC items that concern the Xarray Dataset
    bands: array
        Selected bands

    Returns
    -------
    xarray.Dataset
    """
    if not isinstance(bands, tuple):
        bands = tuple([bands])
    for i in bands:
        dc[i] = post_process_eodc_cube_(dc[i], items, i)
    return dc


def post_process_eodc_cube_(dc, items, band):
    fields = items[0].assets[band].extra_fields
    scale = fields.get("raster:bands")[0]["scale"]
    nodata = fields.get("raster:bands")[0]["nodata"]
    return dc.where(dc != nodata) / scale


bands = "VV"
sig0_dc = odc_stac.load(items_sig0, bands=bands, bbox=bounding_box)
sig0_dc = (
    post_process_eodc_cube(sig0_dc, items_sig0, bands)
    .rename_vars({"VV": "sig0"})
    .dropna(dim="time", how="all")
    .median("time")
)

sig0_dc

<xarray.Dataset> Size: 5MB
Dimensions:      (y: 977, x: 1324)
Coordinates:
  * y            (y) float64 8kB 6.388e+05 6.388e+05 ... 6.193e+05 6.193e+05
  * x            (x) float64 11kB 5.658e+06 5.658e+06 ... 5.684e+06 5.684e+06
    spatial_ref  int32 4B 27704
Data variables:
    sig0         (y, x) float32 5MB -9.6 -9.2 -8.3 -8.7 ... -12.3 -11.6 -9.7

Harmonic Parameters

search = eodc_catalog.search(
    collections="SENTINEL1_HPAR",
    bbox=bounding_box,
    query=["sat:relative_orbit=80"],
)

items_hpar = search.item_collection()
bands = ("C1", "C2", "C3", "M0", "S1", "S2", "S3", "STD")
hpar_dc = odc_stac.load(
    items_hpar,
    bands=bands,
    bbox=bounding_box,
    groupby=None,
)
hpar_dc = post_process_eodc_cube(hpar_dc, items_hpar, bands).median("time")
hpar_dc

<xarray.Dataset> Size: 41MB
Dimensions:      (y: 977, x: 1324)
Coordinates:
  * y            (y) float64 8kB 6.388e+05 6.388e+05 ... 6.193e+05 6.193e+05
  * x            (x) float64 11kB 5.658e+06 5.658e+06 ... 5.684e+06 5.684e+06
    spatial_ref  int32 4B 27704
Data variables:
    C1           (y, x) float32 5MB -0.1 -0.1 0.0 0.1 0.3 ... 1.2 1.6 1.8 1.4
    C2           (y, x) float32 5MB -0.1 -0.2 -0.2 0.0 -0.1 ... 0.2 0.2 0.6 0.6
    C3           (y, x) float32 5MB 0.2 0.1 0.0 0.0 0.1 ... -0.4 -0.6 -0.5 -0.6
    M0           (y, x) float32 5MB -9.0 -9.7 -10.0 -9.7 ... -11.8 -11.3 -11.5
    S1           (y, x) float32 5MB -0.3 -0.2 -0.2 -0.1 ... -0.3 -0.2 -0.7 -1.1
    S2           (y, x) float32 5MB -0.2 0.0 0.0 -0.2 ... -0.2 -0.3 -0.4 -0.2
    S3           (y, x) float32 5MB -0.1 0.0 0.0 -0.1 -0.1 ... 0.0 0.1 0.1 0.4
    STD          (y, x) float32 5MB 1.3 1.2 1.1 1.0 1.2 ... 1.9 1.9 1.8 1.8 1.9

Projected Local Incidence Angles

search = eodc_catalog.search(
    collections="SENTINEL1_MPLIA",
    bbox=bounding_box,
    query=["sat:relative_orbit=80"],
)

items_plia = search.item_collection()

bands = "MPLIA"
plia_dc = odc_stac.load(
    items_plia,
    bands=bands,
    bbox=bounding_box,
)

plia_dc = post_process_eodc_cube(plia_dc, items_plia, bands).median("time")
plia_dc

<xarray.Dataset> Size: 5MB
Dimensions:      (y: 977, x: 1324)
Coordinates:
  * y            (y) float64 8kB 6.388e+05 6.388e+05 ... 6.193e+05 6.193e+05
  * x            (x) float64 11kB 5.658e+06 5.658e+06 ... 5.684e+06 5.684e+06
    spatial_ref  int32 4B 27704
Data variables:
    MPLIA        (y, x) float32 5MB 27.32 29.22 32.16 ... 33.79 34.02 34.27

Finally, we merged the datasets as one big dataset and reproject the data in EPSG 4326 for easier visualizing of the data.

flood_dc = xr.merge([sig0_dc, plia_dc, hpar_dc])
flood_dc = flood_dc.rio.reproject("EPSG:4326").rio.write_crs("EPSG:4326")
flood_dc

<xarray.Dataset> Size: 49MB
Dimensions:      (x: 1443, y: 846)
Coordinates:
  * x            (x) float64 12kB 21.92 21.92 21.92 21.93 ... 22.23 22.23 22.23
  * y            (y) float64 7kB 39.65 39.65 39.65 39.65 ... 39.46 39.46 39.46
    spatial_ref  int64 8B 0
Data variables:
    sig0         (y, x) float32 5MB nan nan nan nan nan ... nan nan nan nan nan
    MPLIA        (y, x) float32 5MB nan nan nan nan nan ... nan nan nan nan nan
    C1           (y, x) float32 5MB nan nan nan nan nan ... nan nan nan nan nan
    C2           (y, x) float32 5MB nan nan nan nan nan ... nan nan nan nan nan
    C3           (y, x) float32 5MB nan nan nan nan nan ... nan nan nan nan nan
    M0           (y, x) float32 5MB nan nan nan nan nan ... nan nan nan nan nan
    S1           (y, x) float32 5MB nan nan nan nan nan ... nan nan nan nan nan
    S2           (y, x) float32 5MB nan nan nan nan nan ... nan nan nan nan nan
    S3           (y, x) float32 5MB nan nan nan nan nan ... nan nan nan nan nan
    STD          (y, x) float32 5MB nan nan nan nan nan ... nan nan nan nan nan

From Backscattering to Flood Mapping

In the following lines we create a map with microwave backscattering values.

# | label: fig-area
# | fig-cap: Area targeted for $\sigma^0$ backscattering is the Greek region of Thessaly, which experienced a major flood in February of 2018.
mrs_view = flood_dc.sig0.hvplot.image(
    x="x", y="y", cmap="viridis", geo=True, tiles=True
).opts(frame_height=400)
mrs_view