Local Flood Mapping

In the context of this notebook we referred to “local flood mapping” as performing the computation on you own machine, as opposed to remote processing which will be explained in the next notebook.

We will begin by loading the flood module from the dask_flood_mapper. You could use the default Dask settings, but it is often good to fine tune te settings to you machine’s specifications. We can do this by setting the Dask scheduler through the Client of dask.distributed.

from importlib.resources import files

import hvplot.xarray  # noqa
import xarray as xr
from dask.distributed import Client
from dask_flood_mapper import flood

Here we set the Dask Scheduler with the Client, where we avoid inter-worker communication which is common for working with numpy and xarray in this case.

client = Client(processes=False, threads_per_worker=2, n_workers=1, memory_limit="12GB")
client

Client

Client-dce5c823-5bfe-11f0-88ac-6045bd4c942d

Connection method: Cluster object	Cluster type: distributed.LocalCluster
Dashboard: http://10.1.0.247:8787/status

Cluster Info

LocalCluster

958352d8

Dashboard: http://10.1.0.247:8787/status	Workers: 1
Total threads: 2	Total memory: 11.18 GiB
Status: running	Using processes: False

Scheduler Info

Scheduler

Scheduler-d289547f-c196-43ae-a3d8-0293dad01821

Comm: inproc://10.1.0.247/2220/1	Workers: 1
Dashboard: http://10.1.0.247:8787/status	Total threads: 2
Started: Just now	Total memory: 11.18 GiB

Workers

Worker: 0

Comm: inproc://10.1.0.247/2220/4	Total threads: 2
Dashboard: http://10.1.0.247:46335/status	Memory: 11.18 GiB
Nanny: None
Local directory: /tmp/dask-scratch-space/worker-m8i0ys7n

No we are ready to map a flood. As an example we use here storm Babet which hit the Danish and Northern coast of Germany at the 20^th of October 2023 Wikipedia. We target an area around Zingst at the Baltic coast of Northern Germany.

time_range = "2023-10-11/2023-10-25"
bbox = [12.3, 54.3, 13.1, 54.6]
fd = flood.decision(bbox=bbox, datetime=time_range).compute()
fd

sigma naught datacube processed

harmonic parameter datacube processed

projected local incidence angle processed

<xarray.DataArray 'reproject-7e5439c8f4ad45e9bac1d4d7f67c22b1' (time: 9,
                                                                latitude: 1232,
                                                                longitude: 3283)> Size: 291MB
array([[[nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [ 0.,  0.,  0., ..., nan, nan, nan],
        ...,
        [nan, nan,  0., ..., nan, nan, nan],
        [nan, nan,  0., ..., nan, nan, nan],
        [nan, nan,  0., ..., nan, nan, nan]],

       [[nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ...,  0., nan, nan],
        ...,
        [nan, nan, nan, ...,  0.,  0.,  0.],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan]],

       [[nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ...,  0., nan, nan],
        ...,
...
        ...,
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan]],

       [[nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [ 0.,  0.,  0., ..., nan, nan, nan],
        ...,
        [nan, nan,  0., ..., nan, nan, nan],
        [nan, nan,  0., ..., nan, nan, nan],
        [nan, nan,  0., ..., nan, nan, nan]],

       [[nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ...,  0., nan, nan],
        ...,
        [nan, nan, nan, ...,  0.,  0.,  0.],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan]]])
Coordinates:
  * time         (time) datetime64[ns] 72B 2023-10-11T05:33:18 ... 2023-10-25...
  * latitude     (latitude) float64 10kB 54.6 54.6 54.6 54.6 ... 54.3 54.3 54.3
  * longitude    (longitude) float64 26kB 12.3 12.3 12.3 12.3 ... 13.1 13.1 13.1
    spatial_ref  int64 8B 0

xarray.DataArray

'reproject-7e5439c8f4ad45e9bac1d4d7f67c22b1'

• time: 9
• latitude: 1232
• longitude: 3283

• nan nan nan nan nan nan nan nan ... nan nan nan nan nan nan nan nan

  array([[[nan, nan, nan, ..., nan, nan, nan],
          [nan, nan, nan, ..., nan, nan, nan],
          [ 0.,  0.,  0., ..., nan, nan, nan],
          ...,
          [nan, nan,  0., ..., nan, nan, nan],
          [nan, nan,  0., ..., nan, nan, nan],
          [nan, nan,  0., ..., nan, nan, nan]],
  
         [[nan, nan, nan, ..., nan, nan, nan],
          [nan, nan, nan, ..., nan, nan, nan],
          [nan, nan, nan, ...,  0., nan, nan],
          ...,
          [nan, nan, nan, ...,  0.,  0.,  0.],
          [nan, nan, nan, ..., nan, nan, nan],
          [nan, nan, nan, ..., nan, nan, nan]],
  
         [[nan, nan, nan, ..., nan, nan, nan],
          [nan, nan, nan, ..., nan, nan, nan],
          [nan, nan, nan, ...,  0., nan, nan],
          ...,
  ...
          ...,
          [nan, nan, nan, ..., nan, nan, nan],
          [nan, nan, nan, ..., nan, nan, nan],
          [nan, nan, nan, ..., nan, nan, nan]],
  
         [[nan, nan, nan, ..., nan, nan, nan],
          [nan, nan, nan, ..., nan, nan, nan],
          [ 0.,  0.,  0., ..., nan, nan, nan],
          ...,
          [nan, nan,  0., ..., nan, nan, nan],
          [nan, nan,  0., ..., nan, nan, nan],
          [nan, nan,  0., ..., nan, nan, nan]],
  
         [[nan, nan, nan, ..., nan, nan, nan],
          [nan, nan, nan, ..., nan, nan, nan],
          [nan, nan, nan, ...,  0., nan, nan],
          ...,
          [nan, nan, nan, ...,  0.,  0.,  0.],
          [nan, nan, nan, ..., nan, nan, nan],
          [nan, nan, nan, ..., nan, nan, nan]]])

• Coordinates: (4)
  • time
    (time)
    datetime64[ns]
    2023-10-11T05:33:18 ... 2023-10-...

    array(['2023-10-11T05:33:18.000000000', '2023-10-13T05:16:57.000000000',
           '2023-10-16T16:53:04.000000000', '2023-10-16T16:53:29.000000000',
           '2023-10-18T05:25:06.000000000', '2023-10-21T17:01:14.000000000',
           '2023-10-21T17:01:39.000000000', '2023-10-23T05:33:19.000000000',
           '2023-10-25T05:16:57.000000000'], dtype='datetime64[ns]')

  • latitude
    (latitude)
    float64
    54.6 54.6 54.6 ... 54.3 54.3 54.3

    units :

    degrees_north

    resolution :

    -0.00024378434985427075

    crs :

    EPSG:4326

    axis :

    Y

    long_name :

    latitude

    standard_name :

    latitude

    array([54.600015, 54.599771, 54.599528, ..., 54.300404, 54.30016 , 54.299917])

  • longitude
    (longitude)
    float64
    12.3 12.3 12.3 ... 13.1 13.1 13.1

    units :

    degrees_east

    resolution :

    0.00024378434985427075

    crs :

    EPSG:4326

    axis :

    X

    long_name :

    longitude

    standard_name :

    longitude

    array([12.300017, 12.300261, 12.300505, ..., 13.09963 , 13.099874, 13.100118])

  • spatial_ref
    ()
    int64
    0

    crs_wkt :

    GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]

    semi_major_axis :

    6378137.0

    semi_minor_axis :

    6356752.314245179

    inverse_flattening :

    298.257223563

    reference_ellipsoid_name :

    WGS 84

    longitude_of_prime_meridian :

    0.0

    prime_meridian_name :

    Greenwich

    geographic_crs_name :

    WGS 84

    horizontal_datum_name :

    World Geodetic System 1984

    grid_mapping_name :

    latitude_longitude

    spatial_ref :

    GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]

    GeoTransform :

    12.299895587547377 0.0002437843498542707 0.0 54.60013705251117 0.0 -0.0002437843498542731

    array(0)

• Indexes: (3)
  • time
    PandasIndex

    PandasIndex(DatetimeIndex(['2023-10-11 05:33:18', '2023-10-13 05:16:57',
                   '2023-10-16 16:53:04', '2023-10-16 16:53:29',
                   '2023-10-18 05:25:06', '2023-10-21 17:01:14',
                   '2023-10-21 17:01:39', '2023-10-23 05:33:19',
                   '2023-10-25 05:16:57'],
                  dtype='datetime64[ns]', name='time', freq=None))

  • latitude
    PandasIndex

    PandasIndex(Index([ 54.60001516033625,  54.59977137598639,  54.59952759163654,
            54.59928380728668, 54.599040022936826,  54.59879623858698,
            54.59855245423712,  54.59830866988727,  54.59806488553741,
           54.597821101187556,
           ...
            54.30211068481433, 54.301866900464475,  54.30162311611462,
           54.301379331764764,  54.30113554741491,  54.30089176306505,
           54.300647978715205,  54.30040419436535,  54.30016041001549,
            54.29991662566564],
          dtype='float64', name='latitude', length=1232))

  • longitude
    PandasIndex

    PandasIndex(Index([12.300017479722305, 12.300261264072159, 12.300505048422014,
           12.300748832771868, 12.300992617121722, 12.301236401471577,
            12.30148018582143, 12.301723970171285,  12.30196775452114,
           12.302211538870994,
           ...
           13.097923656795333, 13.098167441145188, 13.098411225495042,
           13.098655009844896, 13.098898794194751, 13.099142578544605,
           13.099386362894458, 13.099630147244314, 13.099873931594168,
           13.100117715944021],
          dtype='float64', name='longitude', length=3283))

• Attributes: (0)

Since flood.decision does not make a distinction between Sentinel-1 observations over land or over sea, we need to mask pixels over water. For the example here we will load a mask distributed along with this package. But in general this step is left to the user’s own discretion.

data_text = files("dask_flood_mapper.data").joinpath("wcover.nc")
wcover = xr.open_dataarray(data_text, decode_coords="all")

Then we can mask water bodies, like so:

fd = fd.where(wcover != 80)

Now we are ready to view our results. We use here hvplot to create an animation of the flood’s extent over time.

fd.hvplot.image(
    x="longitude",
    y="latitude",
    rasterize=True,
    geo=True,
    tiles=True,
    project=True,
    cmap=["rgba(0, 0, 1, 0.1)", "darkred"],
    cticks=[(0, "non-flood"), (1, "flood")],
    frame_height=400,
)

client.close()