Calculate Harmonic Parameters

Calculate Harmonic Parameters#

Besides using static harmonic parameters to derive the non-flood likelihood \(P(\sigma^0|nonflood)\), we can also dynamically calculate the parameters form past backscattering information. In this notebook we show how we can extract these coefficients that describe seasonal patterns in Sentinel 1 radar backscatter variability.

import hvplot.xarray  # noqa
import numpy as np
import xarray as xr
from dask.distributed import Client
from dask_flood_mapper import flood
from dask_flood_mapper.catalog import (format_datetime_for_xarray_selection,
                                       initialize_catalog, initialize_search)
from dask_flood_mapper.harmonic_params import (
    create_harmonic_parameters, process_harmonic_parameters_datacube)
from dask_flood_mapper.processing import prepare_dc, process_sig0_dc

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

Client

Client-dc695bb4-6e14-11f0-8985-000d3a969aea

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

Cluster Info

As an example we will select only a small region of interest contained in the Zingst case study.

time_range = "2022-10-11T05:25:26"
minlon, maxlon = 12.999, 13
minlat, maxlat = 53.999, 54
bounding_box = [minlon, minlat, maxlon, maxlat]

We will now load the sigma nought datacube from EODC. In order to fit the harmonic functions to the timeseries, the temporal range is extended to include the three years before the provided time range. This is done by using the dynamic parameter of initialize_search.

eodc_catalog = initialize_catalog()
search = initialize_search(eodc_catalog, bounding_box, time_range, dynamic=True)
items_sig0 = search.item_collection()
sig0_dc = prepare_dc(items_sig0, bounding_box, bands="VV")
sig0_dc, orbit_sig0 = process_sig0_dc(sig0_dc, items_sig0, bands="VV")
sig0_dc

---------------------------------------------------------------------------
CPLE_HttpResponseError                    Traceback (most recent call last)
File rasterio/_base.pyx:310, in rasterio._base.DatasetBase.__init__()

File rasterio/_base.pyx:221, in rasterio._base.open_dataset()

File rasterio/_err.pyx:359, in rasterio._err.exc_wrap_pointer()

CPLE_HttpResponseError: CURL error: Failed to connect to data.eodc.eu port 443 after 159 ms: Couldn't connect to server

During handling of the above exception, another exception occurred:

RasterioIOError                           Traceback (most recent call last)
Cell In[4], line 5
      3 items_sig0 = search.item_collection()
      4 sig0_dc = prepare_dc(items_sig0, bounding_box, bands="VV")
----> 5 sig0_dc, orbit_sig0 = process_sig0_dc(sig0_dc, items_sig0, bands="VV")
      6 sig0_dc

File /opt/hostedtoolcache/Python/3.10.18/x64/lib/python3.10/site-packages/dask_flood_mapper/processing.py:55, in process_sig0_dc(sig0_dc, items_sig0, bands)
     48 def process_sig0_dc(
     49     sig0_dc: xr.Dataset,
     50     items_sig0: ItemCollection,
     51     bands: str | list[str],
     52 ) -> tuple[xr.Dataset, np.ndarray]:
     53     """Process the sig0 datacube."""
     54     sig0_dc: xr.Dataset = (
---> 55         post_process_eodc_cube(sig0_dc, items_sig0, bands)
     56         .rename_vars({"VV": "sig0"})
     57         .assign_coords(orbit=("time", extract_orbit_names(items_sig0)))
     58         .dropna(dim="time", how="all")
     59         .sortby("time")
     60     )
     61     orbit_sig0: np.ndarray = order_orbits(sig0_dc)
     62     sig0_dc: xr.Dataset = sig0_dc.groupby("time").mean(skipna=True)

File /opt/hostedtoolcache/Python/3.10.18/x64/lib/python3.10/site-packages/xarray/util/deprecation_helpers.py:115, in _deprecate_positional_args.<locals>._decorator.<locals>.inner(*args, **kwargs)
    111     kwargs.update({name: arg for name, arg in zip_args})
    113     return func(*args[:-n_extra_args], **kwargs)
--> 115 return func(*args, **kwargs)

File /opt/hostedtoolcache/Python/3.10.18/x64/lib/python3.10/site-packages/xarray/core/dataset.py:6356, in Dataset.dropna(self, dim, how, thresh, subset)
   6354     if dim in array.dims:
   6355         dims = [d for d in array.dims if d != dim]
-> 6356         count += np.asarray(array.count(dims))
   6357         size += math.prod([self.sizes[d] for d in dims])
   6359 if thresh is not None:

File /opt/hostedtoolcache/Python/3.10.18/x64/lib/python3.10/site-packages/xarray/core/common.py:165, in AbstractArray.__array__(self, dtype)
    164 def __array__(self: Any, dtype: DTypeLike | None = None) -> np.ndarray:
--> 165     return np.asarray(self.values, dtype=dtype)

File /opt/hostedtoolcache/Python/3.10.18/x64/lib/python3.10/site-packages/xarray/core/variable.py:525, in Variable.values(self)
    522 @property
    523 def values(self):
    524     """The variable's data as a numpy.ndarray"""
--> 525     return _as_array_or_item(self._data)

File /opt/hostedtoolcache/Python/3.10.18/x64/lib/python3.10/site-packages/xarray/core/variable.py:323, in _as_array_or_item(data)
    309 def _as_array_or_item(data):
    310     """Return the given values as a numpy array, or as an individual item if
    311     it's a 0d datetime64 or timedelta64 array.
    312 
   (...)
    321     TODO: remove this (replace with np.asarray) once these issues are fixed
    322     """
--> 323     data = np.asarray(data)
    324     if data.ndim == 0:
    325         if data.dtype.kind == "M":

File /opt/hostedtoolcache/Python/3.10.18/x64/lib/python3.10/site-packages/dask/array/core.py:1693, in Array.__array__(self, dtype, **kwargs)
   1692 def __array__(self, dtype=None, **kwargs):
-> 1693     x = self.compute()
   1694     if dtype and x.dtype != dtype:
   1695         x = x.astype(dtype)

File /opt/hostedtoolcache/Python/3.10.18/x64/lib/python3.10/site-packages/dask/base.py:375, in DaskMethodsMixin.compute(self, **kwargs)
    351 def compute(self, **kwargs):
    352     """Compute this dask collection
    353 
    354     This turns a lazy Dask collection into its in-memory equivalent.
   (...)
    373     dask.compute
    374     """
--> 375     (result,) = compute(self, traverse=False, **kwargs)
    376     return result

File /opt/hostedtoolcache/Python/3.10.18/x64/lib/python3.10/site-packages/dask/base.py:661, in compute(traverse, optimize_graph, scheduler, get, *args, **kwargs)
    658     postcomputes.append(x.__dask_postcompute__())
    660 with shorten_traceback():
--> 661     results = schedule(dsk, keys, **kwargs)
    663 return repack([f(r, *a) for r, (f, a) in zip(results, postcomputes)])

File /opt/hostedtoolcache/Python/3.10.18/x64/lib/python3.10/site-packages/odc/loader/_builder.py:436, in _dask_loader_tyx()
    434 with rdr.restore_env(env, load_state):
    435     for ti, ti_srcs in enumerate(srcs):
--> 436         _fill_nd_slice(
    437             ti_srcs, gbox, cfg, chunk[ti], ydim=ydim, selection=selection
    438         )
    439     return chunk

File /opt/hostedtoolcache/Python/3.10.18/x64/lib/python3.10/site-packages/odc/loader/_builder.py:513, in _fill_nd_slice()
    510     return dst
    512 src, *rest = srcs
--> 513 yx_roi, pix = src.read(cfg, dst_gbox, dst=dst, selection=selection)
    514 assert len(yx_roi) == 2
    515 assert pix.ndim == dst.ndim

File /opt/hostedtoolcache/Python/3.10.18/x64/lib/python3.10/site-packages/odc/loader/_rio.py:115, in read()
    107 def read(
    108     self,
    109     cfg: RasterLoadParams,
   (...)
    113     selection: Optional[ReaderSubsetSelection] = None,
    114 ) -> tuple[tuple[slice, slice], np.ndarray]:
--> 115     return rio_read(self._src, cfg, dst_geobox, dst=dst, selection=selection)

File /opt/hostedtoolcache/Python/3.10.18/x64/lib/python3.10/site-packages/odc/loader/_rio.py:526, in rio_read()
    520     if cfg.fail_on_error:
    521         log.error(
    522             "Aborting load due to failure while reading: %s:%d",
    523             src.uri,
    524             src.band,
    525         )
--> 526         raise e
    527 except rasterio.errors.RasterioError as e:
    528     if cfg.fail_on_error:

File /opt/hostedtoolcache/Python/3.10.18/x64/lib/python3.10/site-packages/odc/loader/_rio.py:512, in rio_read()
    508     return roi, out.transpose([1, 2, 0])
    510 try:
    511     return fixup_out(
--> 512         _rio_read(src, cfg, dst_geobox, prep_dst(dst), selection=selection)
    513     )
    514 except (
    515     rasterio.errors.RasterioIOError,
    516     rasterio.errors.RasterBlockError,
    517     rasterio.errors.WarpOperationError,
    518     rasterio.errors.WindowEvaluationError,
    519 ) as e:
    520     if cfg.fail_on_error:

File /opt/hostedtoolcache/Python/3.10.18/x64/lib/python3.10/site-packages/odc/loader/_rio.py:560, in _rio_read()
    549 def _rio_read(
    550     src: RasterSource,
    551     cfg: RasterLoadParams,
   (...)
    556     # if resampling is `nearest` then ignore sub-pixel translation when deciding
    557     # whether we can just paste source into destination
    558     ttol = 0.9 if cfg.nearest else 0.05
--> 560     with rasterio.open(src.uri, "r", sharing=False) as rdr:
    561         assert isinstance(rdr, rasterio.DatasetReader)
    562         ovr_idx: Optional[int] = None

File /opt/hostedtoolcache/Python/3.10.18/x64/lib/python3.10/site-packages/rasterio/env.py:463, in wrapper()
    460     session = DummySession()
    462 with env_ctor(session=session):
--> 463     return f(*args, **kwds)

File /opt/hostedtoolcache/Python/3.10.18/x64/lib/python3.10/site-packages/rasterio/__init__.py:356, in open()
    353     path = _parse_path(raw_dataset_path)
    355 if mode == "r":
--> 356     dataset = DatasetReader(path, driver=driver, sharing=sharing, **kwargs)
    357 elif mode == "r+":
    358     dataset = get_writer_for_path(path, driver=driver)(
    359         path, mode, driver=driver, sharing=sharing, **kwargs
    360     )

File rasterio/_base.pyx:312, in rasterio._base.DatasetBase.__init__()

RasterioIOError: CURL error: Failed to connect to data.eodc.eu port 443 after 159 ms: Couldn't connect to server

Calculate Harmonic Parameters#

This function fits sine and cosine functions known as harmonic oscillators to each pixel of the Sentinel 1 \(\sigma^0\) datacube. These seasonally varying curves can then be extracted from time series. What is left is the noise or transient events, for example flood events, superimposed on the seasonal trend.

# Function to reformat time range for selecting from the Xarray
datetime_xr = format_datetime_for_xarray_selection(search, time_range)
hpar_list = create_harmonic_parameters(sig0_dc)
__, hpar_dc, __ = process_harmonic_parameters_datacube(sig0_dc, datetime_xr, hpar_list)
hpar_dc

t = sig0_dc.time.dt.dayofyear
n = 365
y = (
    hpar_dc.M0
    + hpar_dc.C1 * np.cos(2 * np.pi * t / n)
    + hpar_dc.S1 * np.sin(2 * np.pi * t / n)
    + hpar_dc.C2 * np.cos(2 * np.pi * t / n)
    + hpar_dc.S2 * np.sin(2 * np.pi * t / n)
    + hpar_dc.C3 * np.cos(2 * np.pi * t / n)
    + hpar_dc.S3 * np.sin(2 * np.pi * t / n)
)

Fit Harmonic Function to Original Data#

Finally, we merge the two datasets and superimpose the fitted harmonic function on the raw sigma nought data.

xr.merge([y.rename("pred"), sig0_dc]).squeeze().assign_coords(
    x=range(len(sig0_dc.x)), y=range(len(sig0_dc.y))
).hvplot(x="time")

Calculating of harmonic parameters can be selected during flood mapping by setting the parameter to True.

flood.decision(bbox=bounding_box, datetime=time_range, dynamic=True)

This setting theoretically frees the user from the need for a precomputed harmonic parameters dataset, thereby only requiring sentinel-1 sigma nought and PLIA datasets as static inputs. Do mind, however, that calculating harmonic parameters can be a computationally expensive operation.

Contents