{ "cells": [ { "cell_type": "markdown", "id": "b809a35c-62b1-4fef-b69d-a00417343a94", "metadata": {}, "source": [ "## Accessing Deltares global flood data\n", "\n", "[Deltares](https://www.deltares.nl/en/) has produced a series of global inundation maps of flood depth using a geographic information systems (GIS)-based inundation model that takes into account water level attenuation and is forced by sea level. Multiple datasets were created using various digital elevation models (DEMs) at multiple resolutions under two different sea level rise (SLR) conditions: current (2018) and 2050. \n", "\n", "This dataset is stored in the West Europe Azure region, so this notebook will run most efficiently on Azure compute located in the same region. If you are using this data for environmental science applications, consider applying for an AI for Earth grant to support your compute requirements." ] }, { "cell_type": "markdown", "id": "f215a7ef-bdd5-497e-b9bb-53621ef8169c", "metadata": {}, "source": [ "### Environment setup" ] }, { "cell_type": "code", "execution_count": 1, "id": "30185f5a-50e3-4c9d-8644-1cb8b48aadde", "metadata": {}, "outputs": [], "source": [ "import math\n", "import warnings\n", "\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", "import xarray as xr\n", "\n", "import dask.distributed\n", "import pystac_client\n", "import planetary_computer\n", "import rioxarray # noqa: F401\n", "import contextily\n", "import shapely.geometry\n", "import cartopy.feature as cfeature\n", "import cartopy.crs as ccrs" ] }, { "cell_type": "markdown", "id": "13218229-f4ec-4d8f-8915-112ddcd6df4c", "metadata": {}, "source": [ "### Create a local Dask cluster\n", "\n", "Enable parallel reads and processing of data using Dask and xarray." ] }, { "cell_type": "code", "execution_count": 2, "id": "c2328949-b3dc-4dc9-abfa-c643ec7fbcbb", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "/user/taugspurger@microsoft.com/proxy/8787/status\n" ] } ], "source": [ "client = dask.distributed.Client(processes=False)\n", "print(client.dashboard_link)" ] }, { "cell_type": "markdown", "id": "661a2ac9-2281-42d7-b289-d17c662c9047", "metadata": {}, "source": [ "### Data access\n", "\n", "The entire dataset is made up of several dozen individual netCDF files, each representing an entire global inundation map, but derived from either a different source DEM, sea level rise condition, or return period. Return periods are occurrence probabilities for floods of a particular magnitude, often referred to as, for example, \"a 100 year flood\". Use the STAC API to query on these various properties:\n", "\n", "To start, we'll load and plot the inundation data produced from the 90m NASADEM at a 100 year return period for 2050 sea level rise conditions. " ] }, { "cell_type": "code", "execution_count": 3, "id": "f5a73cc7-1635-4cc8-bb5e-5ad4fc5c47bd", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "catalog = pystac_client.Client.open(\n", " \"https://planetarycomputer.microsoft.com/api/stac/v1\",\n", " modifier=planetary_computer.sign_inplace,\n", ")\n", "search = catalog.search(\n", " collections=[\"deltares-floods\"],\n", " query={\n", " \"deltares:dem_name\": {\"eq\": \"NASADEM\"},\n", " \"deltares:sea_level_year\": {\"eq\": 2050},\n", " \"deltares:return_period\": {\"eq\": 10},\n", " },\n", ")\n", "\n", "item = next(search.get_items())\n", "item" ] }, { "cell_type": "markdown", "id": "bd164755-8627-423b-b24d-bd595b9a8da7", "metadata": {}, "source": [ "This item has two assets: one pointing to the NetCDF file and one pointing to an index file enabling cloud-optimzed access. When accessing the data from Azure, we recommend using the index file." ] }, { "cell_type": "code", "execution_count": 4, "id": "a2e59759-e52d-46bc-9fca-9b608d4075ea", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
<xarray.Dataset>\n",
       "Dimensions:     (time: 1, lat: 216000, lon: 432000)\n",
       "Coordinates:\n",
       "  * lat         (lat) float64 -90.0 -90.0 -90.0 -90.0 ... 90.0 90.0 90.0 90.0\n",
       "  * lon         (lon) float64 -180.0 -180.0 -180.0 -180.0 ... 180.0 180.0 180.0\n",
       "  * time        (time) datetime64[ns] 2010-01-01\n",
       "Data variables:\n",
       "    inun        (time, lat, lon) float32 dask.array<chunksize=(1, 600, 600), meta=np.ndarray>\n",
       "    projection  object ...\n",
       "Attributes:\n",
       "    Conventions:  CF-1.6\n",
       "    config_file:  /mnt/globalRuns/watermask_post_NASA90m_rest/run_rp0010_slr2...\n",
       "    institution:  Deltares\n",
       "    project:      Microsoft Planetary Computer - Global Flood Maps\n",
       "    references:   https://www.deltares.nl/en/\n",
       "    source:       Global Tide and Surge Model v3.0 - ERA5\n",
       "    title:        GFM - NASA DEM 90m - 2050 slr - 0010-year return level
" ], "text/plain": [ "\n", "Dimensions: (time: 1, lat: 216000, lon: 432000)\n", "Coordinates:\n", " * lat (lat) float64 -90.0 -90.0 -90.0 -90.0 ... 90.0 90.0 90.0 90.0\n", " * lon (lon) float64 -180.0 -180.0 -180.0 -180.0 ... 180.0 180.0 180.0\n", " * time (time) datetime64[ns] 2010-01-01\n", "Data variables:\n", " inun (time, lat, lon) float32 dask.array\n", " projection object ...\n", "Attributes:\n", " Conventions: CF-1.6\n", " config_file: /mnt/globalRuns/watermask_post_NASA90m_rest/run_rp0010_slr2...\n", " institution: Deltares\n", " project: Microsoft Planetary Computer - Global Flood Maps\n", " references: https://www.deltares.nl/en/\n", " source: Global Tide and Surge Model v3.0 - ERA5\n", " title: GFM - NASA DEM 90m - 2050 slr - 0010-year return level" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "url = item.assets[\"index\"].href\n", "ds = xr.open_dataset(f\"reference::{url}\", engine=\"zarr\", consolidated=False, chunks={})\n", "ds" ] }, { "cell_type": "markdown", "id": "b90e7b14-8058-4120-9be0-12a0c95aa1eb", "metadata": {}, "source": [ "### Define an area of interest\n", "\n", "The data is 90m at a global scale, but most relevant in coastal areas. Let's zoom in a on a flood-prone region of Myanmar by defining a bounding box and clipping our xarray dataset." ] }, { "cell_type": "code", "execution_count": 5, "id": "8f7447f5-fa20-4888-9475-0a978275be00", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "AoI bounds: (93.9385986328125, 15.617746547613212, 96.96533203125, 18.37016593904468)\n" ] } ], "source": [ "myanmar_geojson = {\n", " \"type\": \"Polygon\",\n", " \"coordinates\": [\n", " [\n", " [93.9385986328125, 15.617746547613212],\n", " [96.96533203125, 15.617746547613212],\n", " [96.96533203125, 18.37016593904468],\n", " [93.9385986328125, 18.37016593904468],\n", " [93.9385986328125, 15.617746547613212],\n", " ]\n", " ],\n", "}\n", "\n", "poly = shapely.geometry.shape(myanmar_geojson)\n", "minx, miny, maxx, maxy = poly.bounds\n", "print(\"AoI bounds:\", poly.bounds)" ] }, { "cell_type": "code", "execution_count": 6, "id": "c4815ea1-1784-4a6e-8efe-f6012177699f", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
<xarray.Dataset>\n",
       "Dimensions:     (time: 1, lat: 3303, lon: 3632)\n",
       "Coordinates:\n",
       "  * lat         (lat) float64 15.62 15.62 15.62 15.62 ... 18.37 18.37 18.37\n",
       "  * lon         (lon) float64 93.94 93.94 93.94 93.94 ... 96.96 96.96 96.96\n",
       "  * time        (time) datetime64[ns] 2010-01-01\n",
       "Data variables:\n",
       "    inun        (time, lat, lon) float32 dask.array<chunksize=(1, 458, 73), meta=np.ndarray>\n",
       "    projection  object ...\n",
       "Attributes:\n",
       "    Conventions:  CF-1.6\n",
       "    config_file:  /mnt/globalRuns/watermask_post_NASA90m_rest/run_rp0010_slr2...\n",
       "    institution:  Deltares\n",
       "    project:      Microsoft Planetary Computer - Global Flood Maps\n",
       "    references:   https://www.deltares.nl/en/\n",
       "    source:       Global Tide and Surge Model v3.0 - ERA5\n",
       "    title:        GFM - NASA DEM 90m - 2050 slr - 0010-year return level
" ], "text/plain": [ "\n", "Dimensions: (time: 1, lat: 3303, lon: 3632)\n", "Coordinates:\n", " * lat (lat) float64 15.62 15.62 15.62 15.62 ... 18.37 18.37 18.37\n", " * lon (lon) float64 93.94 93.94 93.94 93.94 ... 96.96 96.96 96.96\n", " * time (time) datetime64[ns] 2010-01-01\n", "Data variables:\n", " inun (time, lat, lon) float32 dask.array\n", " projection object ...\n", "Attributes:\n", " Conventions: CF-1.6\n", " config_file: /mnt/globalRuns/watermask_post_NASA90m_rest/run_rp0010_slr2...\n", " institution: Deltares\n", " project: Microsoft Planetary Computer - Global Flood Maps\n", " references: https://www.deltares.nl/en/\n", " source: Global Tide and Surge Model v3.0 - ERA5\n", " title: GFM - NASA DEM 90m - 2050 slr - 0010-year return level" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "ds_myanmar = ds.sel(lat=slice(miny, maxy), lon=slice(minx, maxx))\n", "ds_myanmar" ] }, { "cell_type": "markdown", "id": "fe31adf8-8f30-4f59-a130-7e08b1a4b822", "metadata": {}, "source": [ "### Distribution of inundation amounts\n", "\n", "For areas with greater than zero inundation, let's bin the data in 1m increments and see how it's distributed. Counting by 90m pixels is a rough estimate of actual area." ] }, { "cell_type": "code", "execution_count": 7, "id": "e648c42b-7cfe-4575-b0c2-52929e3c3609", "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Select only flooded area\n", "warnings.filterwarnings(\"ignore\", \"All-NaN\")\n", "flooded = ds_myanmar.inun.where(ds.inun > 0)\n", "num_bins = math.ceil(flooded.max().values)\n", "flooded.plot.hist(bins=range(0, num_bins), edgecolor=\"w\")\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "628aada0-ee2c-4b54-abe5-966ceb4ffecd", "metadata": {}, "source": [ "### Plot the layer\n", "\n", "We can also look at the geographic distribution of inundation. We'll add an Esri imagery basemap for some context in our plot." ] }, { "cell_type": "code", "execution_count": 8, "id": "2463f6c9-a653-4590-89fa-b9c9f2b8173d", "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "fig = plt.figure(figsize=(18, 8), dpi=100)\n", "ax = plt.axes(projection=ccrs.PlateCarree())\n", "\n", "ax.set_extent([minx, maxx, miny, maxy])\n", "flooded.plot(\n", " ax=ax,\n", " transform=ccrs.PlateCarree(),\n", " cmap=\"PuBu\",\n", " cbar_kwargs={\"aspect\": 50},\n", " robust=True,\n", " vmin=0,\n", ")\n", "\n", "ax.set_title(\"Myanmar inundation (SLR 2050 at 250 year return period)\")\n", "\n", "contextily.add_basemap(\n", " ax,\n", " source=contextily.providers.Esri.WorldImagery,\n", " zoom=10,\n", " crs=ds.projection.attrs[\"EPSG_code\"],\n", ")" ] }, { "cell_type": "markdown", "id": "b76a73ec-6f4d-4f4c-9253-ebeaa4550ae6", "metadata": {}, "source": [ "### Working with two sea level rise conditions\n", "\n", "We've been looking at the 2050 sea level rise conditions. Let's compare with the current SLR conditions under the same DEM, resolution, and return period.\n", "\n", "First we'll read in the other dataset and concatenate to a single xarray dataset." ] }, { "cell_type": "code", "execution_count": 9, "id": "e6f6007c-e847-4ce5-aa3a-c89103ea75a5", "metadata": {}, "outputs": [], "source": [ "search = catalog.search(\n", " collections=[\"deltares-floods\"],\n", " query={\n", " \"deltares:dem_name\": {\"eq\": \"NASADEM\"},\n", " \"deltares:sea_level_year\": {\"eq\": 2018},\n", " \"deltares:return_period\": {\"eq\": 10},\n", " },\n", ")\n", "\n", "item = next(search.get_items())\n", "item\n", "\n", "ds_2018 = xr.open_dataset(\n", " f\"reference::{item.assets['index'].href}\",\n", " engine=\"zarr\",\n", " chunks={},\n", " consolidated=False,\n", ")\n", "\n", "ds_2018_myanmar = ds_2018.sel(lat=slice(miny, maxy), lon=slice(minx, maxx))" ] }, { "cell_type": "code", "execution_count": 10, "id": "d6a0b1f9-cd6d-46cb-ba29-d29c27ee69a1", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
<xarray.Dataset>\n",
       "Dimensions:     (time: 2, lat: 3303, lon: 3632)\n",
       "Coordinates:\n",
       "  * lat         (lat) float64 15.62 15.62 15.62 15.62 ... 18.37 18.37 18.37\n",
       "  * lon         (lon) float64 93.94 93.94 93.94 93.94 ... 96.96 96.96 96.96\n",
       "  * time        (time) datetime64[ns] 2018-01-01 2050-01-01\n",
       "Data variables:\n",
       "    inun        (time, lat, lon) float32 dask.array<chunksize=(1, 458, 73), meta=np.ndarray>\n",
       "    projection  (time) object nan nan\n",
       "Attributes:\n",
       "    Conventions:  CF-1.6\n",
       "    config_file:  /mnt/globalRuns/watermask_post_NASA/run_rp0010_slr2018/coas...\n",
       "    institution:  Deltares\n",
       "    project:      Microsoft Planetary Computer - Global Flood Maps\n",
       "    references:   https://www.deltares.nl/en/\n",
       "    source:       Global Tide and Surge Model v3.0 - ERA5\n",
       "    title:        GFM - NASA DEM 90m - 2018 slr - 0010-year return level
" ], "text/plain": [ "\n", "Dimensions: (time: 2, lat: 3303, lon: 3632)\n", "Coordinates:\n", " * lat (lat) float64 15.62 15.62 15.62 15.62 ... 18.37 18.37 18.37\n", " * lon (lon) float64 93.94 93.94 93.94 93.94 ... 96.96 96.96 96.96\n", " * time (time) datetime64[ns] 2018-01-01 2050-01-01\n", "Data variables:\n", " inun (time, lat, lon) float32 dask.array\n", " projection (time) object nan nan\n", "Attributes:\n", " Conventions: CF-1.6\n", " config_file: /mnt/globalRuns/watermask_post_NASA/run_rp0010_slr2018/coas...\n", " institution: Deltares\n", " project: Microsoft Planetary Computer - Global Flood Maps\n", " references: https://www.deltares.nl/en/\n", " source: Global Tide and Surge Model v3.0 - ERA5\n", " title: GFM - NASA DEM 90m - 2018 slr - 0010-year return level" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Concat the two datasets along the time dimension\n", "mds = xr.concat([ds_2018_myanmar, ds_myanmar], dim=\"time\")\n", "\n", "# Time coordinates are not set in the data files. Set them correctly\n", "# to allow selecting by label.\n", "mds = mds.assign_coords(\n", " time=np.array([np.datetime64(\"2018-01-01\"), np.datetime64(\"2050-01-01\")])\n", ")\n", "mds" ] }, { "cell_type": "markdown", "id": "f8085a9e-8c31-4085-ba1b-755bbbc99997", "metadata": {}, "source": [ "#### Plot the two layers" ] }, { "cell_type": "code", "execution_count": 11, "id": "fd024fb4-8a19-4e84-9bcc-25afb7b5d7db", "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "flooded = mds.where(mds.inun > 0).load()\n", "\n", "prj = ccrs.PlateCarree()\n", "\n", "p = flooded.inun.plot(\n", " col=\"time\",\n", " col_wrap=None,\n", " transform=prj,\n", " subplot_kws={\"projection\": prj},\n", " cmap=\"PuBu\",\n", " size=10,\n", " cbar_kwargs={\"aspect\": 30, \"shrink\": 0.6},\n", " robust=True,\n", " vmin=0,\n", ")\n", "\n", "for ax in p.axes.flat:\n", " ax.coastlines(linewidth=0.5)\n", " ax.set_extent([minx, maxx, miny, maxy])\n", " ax.add_feature(cfeature.LAND, zorder=0, linewidth=0.5, facecolor=\"#f5f5ed\")\n", "\n", "plt.draw()" ] }, { "cell_type": "markdown", "id": "fd19cf50-eeec-453b-8841-ec850f8c9aef", "metadata": {}, "source": [ "#### Compare the distribution of inundation depth by area" ] }, { "cell_type": "code", "execution_count": 12, "id": "534fa37e-48b8-4ec1-abcf-dcda8d215b5c", "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "fig, axes = plt.subplots(ncols=2, figsize=(15, 5), sharey=True, sharex=True)\n", "\n", "for i in range(2):\n", " num_bins = math.ceil(flooded.inun.max().values)\n", " year = str(flooded.isel(time=i).time.values)[:4]\n", " axes[i].set_ylabel(f\"Pixel count ({year} SLR conditions)\")\n", " flooded.inun.isel(time=i).plot.hist(\n", " ax=axes[i], edgecolor=\"white\", bins=range(0, num_bins)\n", " )\n", "\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "1f63293e-3163-42ae-bba7-21663d08f0f2", "metadata": {}, "source": [ "### Plot areas of increased inundation under 2050 sea level rise conditions\n", "\n", "It was a bit hard to see predicted changes in inundation in the side-by-side map plots, but the histograms showed a definite increase in the higher bins. Let's calculate the difference and plot just the increase." ] }, { "cell_type": "code", "execution_count": 13, "id": "44667dd0-5442-4da5-9bf3-ad64f04e4e01", "metadata": {}, "outputs": [], "source": [ "diff = mds[\"inun\"].diff(\"time\")\n", "diff = diff.where(diff > 0)" ] }, { "cell_type": "code", "execution_count": 14, "id": "bff2a355-b5bf-47fa-a615-6bd123591a8b", "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "fig, axis = plt.subplots(\n", " subplot_kw=dict(projection=ccrs.PlateCarree()), figsize=(12, 8), dpi=100\n", ")\n", "\n", "diff.plot(\n", " ax=axis,\n", " transform=ccrs.PlateCarree(),\n", " cmap=\"plasma\",\n", " cbar_kwargs={\"aspect\": 50},\n", ")\n", "\n", "plt.title(\"Increased inundation under 2050 SLR conditions\")\n", "axis.coastlines(linewidth=0.5)\n", "axis.set_extent([minx, maxx, miny, maxy])\n", "axis.add_feature(cfeature.LAND, zorder=0, linewidth=0, facecolor=\"#f5f5ed\")\n", "\n", "plt.draw()" ] } ], "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.8.13" }, "widgets": { "application/vnd.jupyter.widget-state+json": { "state": {}, "version_major": 2, "version_minor": 0 } } }, "nbformat": 4, "nbformat_minor": 5 }