2 Reverend Bayes updates our Belief in Flood Detection

Image from wikipedia

Note

This notebooks contains interactive element. These interactive element can only be viewed on Binder by clicking on the Binder badge or 🚀 button.

%matplotlib widget

import numpy as np
import xarray as xr
import numpy as np
import datetime
import matplotlib.pyplot as plt

from scipy.stats import norm
from eomaps import Maps

sig0_dc = xr.open_dataset('../data/s1_parameters/S1_CSAR_IWGRDH/SIG0/V1M1R1/EQUI7_EU020M/E054N006T3/SIG0_20180228T043908__VV_D080_E054N006T3_EU020M_V1M1R1_S1AIWGRDH_TUWIEN.nc')

2.1 From Backscattering to Flood Mapping

This notebook explains how microwave (\(\sigma^0\)) backscattering (Figure 2.1) can be used to map the extent of a flood. We replicate in this exercise the work of Bauer-Marschallinger et al. (2022) on the TU Wien Bayesian-based flood mapping algorithm.

In the following lines we create a map with EOmaps (Quast 2024) of the \(\sigma^0\) backscattering values.

m = Maps(ax=121, crs=3857)
m.set_data(data=sig0_dc, x="x", y="y", parameter="SIG0", crs=Maps.CRS.Equi7_EU)
m.plot_map()
m.add_colorbar(label="$\sigma^0$ (dB)", orientation="vertical", hist_bins=30)
m.add_scalebar(n=5)
m2 = m.new_map(ax=122, crs=3857)
m2.set_extent(m.get_extent())
m2.add_wms.OpenStreetMap.add_layer.default()
m.apply_layout(
    {
        'figsize': [7.32, 4.59],
        '0_map': [0.05, 0.18, 0.35, 0.64],
        '1_cb': [0.8125, 0.1, 0.1, 0.8],
        '1_cb_histogram_size': 0.8,
        '2_map': [0.4375, 0.18, 0.35, 0.64]
        }
    )
m.show()

Figure 2.1: Area targeted for \(\sigma^0\) backscattering is the Greek region of Thessaly, which experienced a major flood in February of 2018.

2.2 Microwave Backscattering over Land and Water

Reverend Bayes was concerned with two events, one (the hypothesis) occurring before the other (the evidence). If we know its cause, it is easy to logically deduce the probability of an effect. However, in this case we want to deduce the probability of a cause from an observed effect, also known as “reversed probability”. In the case of flood mapping, we have \(\sigma^0\) backscatter observations over land (the effect) and we want to deduce the probability of flooding (\(F\)) and non-flooding (\(NF\)).

In other words, we want to know the probability of flooding \(P(F)\) given a pixel’s \(\sigma^0\):

\[P(F|\sigma^0)\]

and the probability of a pixel being not flooded \(P(NF)\) given a certain \(\sigma^0\):

\[P(NF|\sigma^0).\]

Bayes showed that these can be deduced from the observation that forward and reversed probability are equal, so that:

\[P(F|\sigma^0)P(\sigma^0) = P(\sigma^0|F)P(F)\]

and

\[P(NF|\sigma^0)P(\sigma^0) = P(\sigma^0|NF)P(NF).\]

The forward probability of \(\sigma^0\) given the occurrence of flooding (\(P(\sigma^0|F)\)) and \(\sigma^0\) given no flooding (\(P(\sigma^0|NF)\)) can be extracted from past information on backscattering over land and water surfaces. As seen in the sketch below (Figure 2.2), the characteristics of backscattering over land and water differ considerably.

Figure 2.2: Schematic backscattering over land and water. Image from Geological Survey Ireland

2.3 Likelihoods

The so-called likelihoods of \(P(\sigma^0|F)\) and \(P(\sigma^0|NF)\) can thus be calculated from past backscattering information. Without going into the details of how these likelihoods are calculated, you can click on a pixel of the map to plot the likelihoods of \(\sigma^0\) being governed by land or water.

RANGE = np.arange(-30, 0, 0.1)
hparam_dc = xr.open_dataset('../data/tuw_s1_harpar/S1_CSAR_IWGRDH/SIG0-HPAR/V0M2R3/EQUI7_EU020M/E054N006T3/D080.nc')
plia_dc = xr.open_dataset('../data/s1_parameters/S1_CSAR_IWGRDH/PLIA-TAG/V01R03/EQUI7_EU020M/E054N006T3/PLIA-TAG-MEAN_20200101T000000_20201231T235959__D080_E054N006T3_EU020M_V01R03_S1IWGRDH.nc')
sig0_dc['id'] = (('y', 'x'), np.arange(sig0_dc.SIG0.size).reshape(sig0_dc.SIG0.shape))
hparam_dc['id'] = (('y', 'x'), np.arange(sig0_dc.SIG0.size).reshape(sig0_dc.SIG0.shape))
plia_dc['id'] = (('y', 'x'), np.arange(sig0_dc.SIG0.size).reshape(sig0_dc.SIG0.shape))

def calc_water_likelihood(id, x):
    point = plia_dc.where(plia_dc.id == id, drop=True)
    wbsc_mean = point.PLIA * -0.394181 + -4.142015
    wbsc_std = 2.754041
    return norm.pdf(x, wbsc_mean.to_numpy(), wbsc_std).flatten()

def expected_land_backscatter(data, dtime_str):
    w = np.pi * 2 / 365
    dt = datetime.datetime.strptime(dtime_str, "%Y-%m-%d")
    t = dt.timetuple().tm_yday
    wt = w * t

    M0 = data.M0
    S1 = data.S1
    S2 = data.S2
    S3 = data.S3
    C1 = data.C1
    C2 = data.C2
    C3 = data.C3
    hm_c1 = (M0 + S1 * np.sin(wt)) + (C1 * np.cos(wt))
    hm_c2 = ((hm_c1 + S2 * np.sin(2 * wt)) + C2 * np.cos(2 * wt))
    hm_c3 = ((hm_c2 + S3 * np.sin(3 * wt)) + C3 * np.cos(3 * wt))
    return hm_c3

def calc_land_likelihood(id, x):
    point = hparam_dc.where(hparam_dc.id == id, drop=True)
    lbsc_mean = expected_land_backscatter(point, '2018-02-01')
    lbsc_std = point.STD
    return norm.pdf(x, lbsc_mean.to_numpy(), lbsc_std.to_numpy()).flatten()

def calc_likelihoods(id, x):
    if isinstance(x, list):
        x = np.arange(x[0], x[1], 0.1)
    water_likelihood, land_likelihood = calc_water_likelihood(id=id, x=x), calc_land_likelihood(id=id, x=x)
    return water_likelihood, land_likelihood

def view_bayes_flood(sig0_dc, calc_posteriors=None, bayesian_flood_decision=None):

    # initialize a map on top
    m = Maps(ax=122, layer="data", crs=Maps.CRS.Equi7_EU)

    # initialize 2 matplotlib plot-axes next to the map
    ax_upper = m.f.add_subplot(221)
    ax_upper.set_ylabel("likelihood")
    ax_upper.set_xlabel("$\sigma^0 (dB)$")

    ax_lower = m.f.add_subplot(223)
    ax_lower.set_ylabel("probability")
    ax_lower.set_xlabel("$\sigma^0 (dB)$")

    # -------- assign data to the map and plot it
    if bayesian_flood_decision is not None:
        # add map
        m2 = m.new_layer(layer="map")
        m2.add_wms.OpenStreetMap.add_layer.default()
        flood_classification =  bayesian_flood_decision(sig0_dc.id, sig0_dc.SIG0)
        sig0_dc["decision"] = (('y', 'x'), flood_classification.reshape(sig0_dc.SIG0.shape))
        sig0_dc["decision"] = sig0_dc.decision.where(sig0_dc.SIG0.notnull())
        sig0_dc["decision"] = sig0_dc.decision.where(sig0_dc.decision==0)
        m.set_data(data=sig0_dc, x="x", y="y", parameter="decision", crs=Maps.CRS.Equi7_EU)
        m.plot_map()
        m.show_layer("map", ("data", 0.5))
        m.apply_layout(
            {
                'figsize': [7.32, 4.59],
                '0_map': [0.44573, 0.11961, 0.3375, 0.75237],
                '1_': [0.10625, 0.5781, 0.3125, 0.29902],
                '2_': [0.10625, 0.11961, 0.3125, 0.29902],
            }
        )

    else:
        m.set_data(data=sig0_dc, x="x", y="y", parameter="SIG0", crs=Maps.CRS.Equi7_EU)
        m.plot_map()
        m.add_colorbar(label="$\sigma^0$ (dB)", orientation="vertical", hist_bins=30)
        m.apply_layout(
            {
                'figsize': [7.32, 4.59],
                '0_map': [0.44573, 0.11961, 0.3375, 0.75237],
                '1_': [0.10625, 0.5781, 0.3125, 0.29902],
                '2_': [0.10625, 0.11961, 0.3125, 0.29902],
                '3_cb': [0.8, 0.09034, 0.1, 0.85],
                '3_cb_histogram_size': 0.8
            }
        )

    def update_plots(ID, **kwargs):
        
        # get the data
        value = sig0_dc.where(sig0_dc.id == ID, drop=True).SIG0.to_numpy()
        y1_pdf, y2_pdf = calc_water_likelihood(ID, RANGE), calc_land_likelihood(ID, RANGE)

        # plot the lines and vline
        (water,) = ax_upper.plot(RANGE, y1_pdf, 'k-', lw=2, label="water")
        (land,) = ax_upper.plot(RANGE, y2_pdf,'r-', lw=5, alpha=0.6, label="land")
        value_left = ax_upper.vlines(x=value, ymin=0, ymax=np.max((y1_pdf, y2_pdf)), lw=3, label="observed")
        ax_upper.legend(loc="upper left")

        # add all artists as "temporary pick artists" so that they
        # are removed when the next datapoint is selected
        for a in [water, land, value_left]:
            m.cb.pick.add_temporary_artist(a)

        if calc_posteriors is not None:
            f_post, nf_post = calc_posteriors(y1_pdf, y2_pdf)
            (f,) = ax_lower.plot(RANGE, f_post, 'k-', lw=2, label="flood")
            (nf,) = ax_lower.plot(RANGE, nf_post,'r-', lw=5, alpha=0.6, label="non-flood")
            value_right = ax_lower.vlines(x=value, ymin=-0.1, ymax=1.1, lw=3, label="observed")
            ax_lower.legend(loc="upper left")
            for a in [f, nf, value_right]:
                m.cb.pick.add_temporary_artist(a)

        # re-compute axis limits based on the new artists
        ax_upper.relim()
        ax_upper.autoscale()

    m.cb.pick.attach(update_plots)
    m.cb.pick.attach.mark(permanent=False, buffer=1, fc="none", ec="r")
    m.cb.pick.attach.mark(permanent=False, buffer=2, fc="none", ec="r", ls=":")
    m.show()

view_bayes_flood(sig0_dc)

Figure 2.3: Likelihoods for \(\sigma^0\) being associated with land or water for 1 pixel in the Greek area of Thessaly. Likelihoods are calculated over a range of \(\sigma^0\). The pixel’s observed \(\sigma^0\) is given with a vertical line. Click on the map to re-calculate and update this figure for another pixel in the study area. Map created with EOmaps (Quast 2024).

2.4 Posteriors

Having calculated the likelihoods, we can now move on to calculate the probability of (non-)flooding given a pixel’s \(\sigma^0\). These so-called posteriors need one more piece of information, as can be seen in the equation above. We need the probability that a pixel is flooded \(P(F)\) or not flooded \(P(NF)\). Of course, these are the figures we’ve been trying to find this whole time. We don’t actually have them yet, so what can we do? In Bayesian statistics, we can just start with our best guess. These guesses are called our “priors”, because they are the beliefs we hold prior to looking at the data. This subjective prior belief is the foundation Bayesian statistics, and we use the likelihoods we just calculated to update our belief in this particular hypothesis. This updated belief is called the “posterior”.

Let’s say that our best estimate for the chance of flooding versus non-flooding of a pixel is 50-50: a coin flip. We now can also calculate the probability of backscattering \(P(\sigma^0)\), as the weighted average of the water and land likelihoods, ensuring that our posteriors range between 0 to 1.

The following code block shows how we calculate the priors.

def calc_posteriors(water_likelihood, land_likelihood):
    evidence = (water_likelihood * 0.5) + (land_likelihood * 0.5)
    return (water_likelihood * 0.5) / evidence,  (land_likelihood * 0.5) / evidence

We can plot the posterior probabilities of flooding and non-flooding again and compare these to pixel’s measured \(\sigma^0\). Click on a pixel to calculate the posterior probability.

view_bayes_flood(sig0_dc, calc_posteriors)

Figure 2.4: Posterior probabilities for \(\sigma^0\) of 1 pixel being associated with land for water in the Greek area of Thessaly. Click on the map to re-calculate and update this figure for another pixel in the study area. Map created with EOmaps (Quast 2024).

2.5 Flood Classification

We are now ready to combine all this information and classify the pixels according to the probability of flooding given the backscatter value of each pixel. Here we just look whether the probability of flooding is higher than non-flooding:

def bayesian_flood_decision(id, sig0_dc):
    nf_post_prob, f_post_prob = calc_posteriors(*calc_likelihoods(id,  sig0_dc))
    return np.greater(f_post_prob, nf_post_prob)

Click on a point in the below map to see the likelihoods and posterior distributions (in the left-hand subplots).

view_bayes_flood(sig0_dc, calc_posteriors, bayesian_flood_decision)

Figure 2.5: Flood extent of the Greek region of Thessaly based on Bayesian probabilities are shown on the map superimposed on an open street map. Click on a pixel to generate the point’s water and land likelihoods as well as the posterior probabilities. Map created with EOmaps (Quast 2024).

Bauer-Marschallinger, Bernhard, Senmao Cao, Mark Edwin Tupas, Florian Roth, Claudio Navacchi, Thomas Melzer, Vahid Freeman, and Wolfgang Wagner. 2022. “Satellite-Based Flood Mapping Through Bayesian Inference from a Sentinel-1 SAR Datacube.” Remote Sensing 14 (15): 3673. https://doi.org/10.3390/rs14153673.

Quast, Raphael. 2024. “EOmaps: A Python Package to Visualize and Analyze Geographical Datasets.” https://doi.org/10.5281/zenodo.6459598.