water/access_data4human.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Access Data4Human Sentinel-1 Floodmasks"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "import dask.diagnostics\n",
    "import folium\n",
    "import folium.plugins as plugins\n",
    "import geopandas as gpd\n",
    "\n",
    "from odc.stac import load\n",
    "from pystac_client import Client\n",
    "from shapely.geometry import shape\n",
    "\n",
    "import pystac"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "## Query STAC API\n",
    "The easiest way to work with the data is to use the [Geoservice STAC API](https://geoservice.dlr.de/eoc/ogc/stac). The Data4Human Sentinel-1 Floodmasks are located in the [D4H collection](https://geoservice.dlr.de/eoc/ogc/stac/collections/D4H). In the following we use [pystac-client](https://github.com/stac-utils/pystac-client) to query the items in this collections by bounding box (*bbox*) and acquisition dates (*dates_from_to*). Optionally, we can also filter the items by properties such as when they have been updated."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "# search stac api\n",
    "stacapi_endpoint = \"https://geoservice.dlr.de/eoc/ogc/stac/v1\"\n",
    "collections = [\"D4H\"]\n",
    "bbox = [33.851, -20.63, 35.859, -18.88]\n",
    "dates_from_to = [\"2019-03-01\", \"2019-03-31\"]\n",
    "filt = None #\"updated > '2022-02-20T00:00:00.000' AND updated < '2022-02-28T10:00:00.000'\"\n",
    "\n",
    "catalog = Client.open(\n",
    "    url=stacapi_endpoint,\n",
    "    ignore_conformance=True\n",
    ")\n",
    "\n",
    "search = catalog.search(\n",
    "    collections=collections, \n",
    "    bbox=bbox, \n",
    "    datetime=dates_from_to, \n",
    "    filter=filt, \n",
    "    method=\"GET\", \n",
    "    filter_lang=\"cql2-text\",\n",
    "    max_items=3000\n",
    ")\n",
    "\n",
    "items = [item for item in search.items()]\n",
    "print(f\"Search returned {len(items)} items\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "## Convert items to Geopandas Dataframe\n",
    "We can convert the pystac items to a Geopandas Dataframe for further vector analysis or simply to save them to file in a range of GIS file formats."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "items_dict = [item.to_dict() for item in items]\n",
    "items_gdf = gpd.GeoDataFrame.from_features(items_dict)\n",
    "items_gdf.head(3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "## Visualize item footprints on map\n",
    "We can use the Geopandas Dataframe to easily visualize the item footprints on an interactive Leaflet (Folium) map."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "m = folium.Map(tiles=\"Stamen Terrain\")\n",
    "layer_control = folium.LayerControl(position=\"topright\", collapsed=True)\n",
    "fullscreen = folium.plugins.Fullscreen()\n",
    "marker_cluster = folium.plugins.MarkerCluster(name=\"Item markers\").add_to(m)\n",
    "\n",
    "for item in items:\n",
    "    geom_pt = shape(item.geometry).representative_point()\n",
    "    folium.Marker(\n",
    "        [geom_pt.y, geom_pt.x], \n",
    "        popup=f\"<a href={item.links[4].target} target='_blank' rel='noopener noreferrer'>{item.id}</a><br>\"\n",
    "        f\"platform: {item.properties['platform']}<br>\"\n",
    "        f\"datetime: {item.datetime}\"\n",
    "    ).add_to(marker_cluster)\n",
    "    \n",
    "style = {\"fillColor\": '#00000000', \"color\": \"#0000ff\", \"weight\": 1}\n",
    "\n",
    "footprints = folium.GeoJson(\n",
    "    items_gdf.to_json(), \n",
    "    name=\"Item footprints\", \n",
    "    style_function=lambda x: style,\n",
    "    control=True\n",
    ")\n",
    "\n",
    "footprints.add_to(m)\n",
    "layer_control.add_to(m)\n",
    "fullscreen.add_to(m)\n",
    "m.fit_bounds(m.get_bounds())\n",
    "m "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "## Convert items to Xarray Dataset\n",
    "We can convert the pystac items to a Xarray Dataset for further raster analysis, such as timeseries analysis or mosaicing. Note that here we set a lower resolution for faster visualization (the native resolution of the products is 0.0001)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "items_ds = load(\n",
    "    items,\n",
    "    bands=[\"WATER\", \"VALID\"],\n",
    "    crs=\"EPSG:4326\",\n",
    "    resolution=0.001,\n",
    "    lon=(bbox[0], bbox[2]),\n",
    "    lat=(bbox[1], bbox[3]),\n",
    "    skip_broken_datasets=True,\n",
    "    chunks={\"latitude\": 2048, \"longitude\": 2048},\n",
    ")  \n",
    "print(items_ds)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compute coverage (number of valid observations) from valid mask assets\n",
    "Each item has several assets assigned to it. For the Data4Human Sentinel-1 Floodmasks the assets are:\n",
    "\n",
    "- \"VALID\": Valid pixel mask\n",
    "- \"WATER\": Water mask\n",
    "- \"REFER\": Reference water mask\n",
    "\n",
    "In the following we use the valid pixel mask assets to compute the number of valid observations per pixel, which provides us with an efficient way to get the satellite coverage for the given area and time of interest."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# sum valid masks along time dimension to compute number of valid observations\n",
    "coverage = items_ds[\"VALID\"].sum(dim=\"time\")\n",
    "\n",
    "# compute and plot coveragemask\n",
    "with dask.diagnostics.ProgressBar():\n",
    "    data = coverage.compute()\n",
    "    data.plot(size=10, cbar_kwargs={'label': \"Coverage (valid observations)\"})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create mosaic from water mask assets\n",
    "Here we get the maximum value along time dimension to mosaic all water masks assets."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get maximum value along time dimension to mosaic water masks\n",
    "mosaic = items_ds[\"WATER\"].max(dim=\"time\")\n",
    "\n",
    "# compute and plot mosaic\n",
    "with dask.diagnostics.ProgressBar():\n",
    "    data = mosaic.compute()\n",
    "    data.plot(size=10, cmap=\"Blues\", add_colorbar=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create water frequency and flood duration from water mask assets\n",
    "Here we do some more time-series analysis to compute a simple water frequency layer and a more advanced flood duration layer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# sum water masks along time axis and divide by the number of valid observations\n",
    "frequency = items_ds[\"WATER\"].sum(dim=\"time\") / items_ds[\"VALID\"].sum(dim=\"time\")\n",
    "\n",
    "# compute and plot frequency\n",
    "with dask.diagnostics.ProgressBar():\n",
    "    data = frequency.compute()\n",
    "    data.plot(size=10, cmap=\"Blues\", vmin=0, vmax=1, cbar_kwargs={'label': \"Relative water frequency\"})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# combine water and valid masks\n",
    "# forward fill water masks where no valid observations are available\n",
    "# set remaining nan values (at first timestamp) to 0\n",
    "water_ds = items_ds[\"WATER\"].where(items_ds[\"VALID\"] == 1).ffill(dim=\"time\").fillna(0)\n",
    "\n",
    "# resample water masks to daily temporal resolution using forward fill and sum up\n",
    "# this effectively provides us with the total flood duration in the requested reference time period\n",
    "duration = water_ds.resample(time=\"1D\").ffill().sum(dim=\"time\")\n",
    "\n",
    "# compute and plot duration\n",
    "with dask.diagnostics.ProgressBar():\n",
    "    data = duration.compute()\n",
    "    data.plot(size=10, cmap=\"Blues\", cbar_kwargs={'label': \"Flood duration (days)\"})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visualize water mask from WMS\n",
    "Besides the STAC API endpoint, we can also use common OGC webservices such as Web Map Service (WMS) to interact with the products on Geoservice. In the following we simply load a WMS layer of the water mask for a desired timestamp. Please note that this examples requires the exact timestamp of a product. There are, however, more convenient tools and Leaflet Plugins available that allow to specify dates or date ranges instead."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# NOTE: this requires the exact timestamp and accepts a list of timestamps seems\n",
    "m = folium.Map(tiles=\"Stamen Terrain\")\n",
    "watermask = folium.raster_layers.WmsTileLayer(\n",
    "    url = \"https://geoservice.dlr.de/eoc/demo/wms?\", \n",
    "    layers=\"D4H_WATER_MASK_S1\", \n",
    "    time='2022-02-21T03:02:11Z', #f\"{dates_from_to[0]}T00:00:00Z\",\n",
    "    #bounds=[[bbox[1], bbox[0]], [bbox[3], bbox[2]]],\n",
    "    fmt='image/png', \n",
    "    transparent=True, \n",
    "    version='1.1.1',\n",
    "    name=\"Water mask\",\n",
    "    control=True\n",
    ")\n",
    "watermask.add_to(m)\n",
    "fullscreen.add_to(m)\n",
    "m.fit_bounds(footprints.get_bounds())\n",
    "m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}