{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot Delft3D FM model using HydroMT"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**HydroMT** provides a simple interface to model schematization from which we can make beautiful plots:\n",
    "\n",
    "- Mesh layers are saved to the model `mesh.data` component as a  `xugrid.UgridDataset`\n",
    "- Vector layers are saved to the model `geoms.data` component as a `geopandas.GeoDataFrame`. Note that in case of Delft3D FM these are not used by the model kernel, but only for analysis and visualization purposes.\n",
    "- Gridded data like bedlevels or infiltration capacity are saved to the model `inifield.data` component as a `xarray.DataArray`. Here the maps are regular grid in the same CRS as the Delft3D FM model but not necessarily the same resolution or grid, as Delft3D FM kernel can do the interpolation.\n",
    "\n",
    "We use the [cartopy](https://scitools.org.uk/cartopy/docs/latest/) package to plot maps. This packages provides a simple interface to plot geographic data and add background satellite imagery."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Load dependencies"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from hydromt_delft3dfm import DFlowFMModel"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot maps dependencies\n",
    "import matplotlib.pyplot as plt\n",
    "import cartopy.crs as ccrs\n",
    "import cartopy.io.img_tiles as cimgt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Read the model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "root = \"dflowfm_piave\"\n",
    "mod = DFlowFMModel(root, mode=\"r\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# read the mesh to get the grids as geodataframes for plotting\n",
    "mesh1d = mod.mesh.mesh_gdf[\"mesh1d\"]\n",
    "mesh2d = mod.mesh.mesh_gdf[\"mesh2d\"]\n",
    "\n",
    "# Get the different types of branches in mesh1d\n",
    "rivers = mod.rivers\n",
    "pipes = mod.pipes\n",
    "# Additional geometry and structures\n",
    "manholes = mod.geoms.data[\"manholes\"]\n",
    "crosssections = mod.geoms.data[\"crosssections\"]\n",
    "\n",
    "# Read the elevation from maps\n",
    "elv = mod.inifield.data[\"elevtn\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot model schematization base maps"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we plot the model base mesh information as well as the topography map. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# we assume the model maps are in the projected CRS EPSG:3857\n",
    "proj = ccrs.epsg(3857)\n",
    "# adjust zoomlevel and figure size to your basis size & aspect\n",
    "zoom_level = 12\n",
    "figsize = (10, 8)\n",
    "\n",
    "# initialize image with geoaxes\n",
    "fig = plt.figure(figsize=figsize)\n",
    "ax = fig.add_subplot(projection=proj)\n",
    "bbox = elv.raster.box.to_crs(3857).buffer(3e3).to_crs(elv.raster.crs).total_bounds\n",
    "extent = np.array(bbox)[[0, 2, 1, 3]]\n",
    "ax.set_extent(extent, crs=proj)\n",
    "\n",
    "# add sat background image\n",
    "ax.add_image(cimgt.QuadtreeTiles(), zoom_level, alpha=0.5)\n",
    "\n",
    "## plot elevation\\\n",
    "elv.plot(transform=proj, ax=ax, zorder=0.5, cmap=\"BrBG\", add_colorbar=True)\n",
    "\n",
    "# plot rivers\n",
    "rivers.plot(ax=ax, linewidth= 1, color=\"blue\", zorder=3, label=\"rivers\")\n",
    "# plot pipes\n",
    "pipes.plot(ax=ax, color=\"k\", linewidth=1, zorder=3, label=\"pipes\")\n",
    "\n",
    "## plot mesh\n",
    "mesh1d.plot(ax=ax, color=\"r\", zorder=2, label=\"mesh1d\")\n",
    "mesh2d.plot(ax=ax, facecolor=\"none\", edgecolor=\"gray\", linewidth=0.5, zorder=2, label=\"mesh2d\")\n",
    "\n",
    "ax.xaxis.set_visible(True)\n",
    "ax.yaxis.set_visible(True)\n",
    "ax.set_ylabel(f\"latitude [degree north]\")\n",
    "ax.set_xlabel(f\"longitude [degree east]\")\n",
    "_ = ax.set_title(f\"dflowfm base map\")\n",
    "legend = ax.legend(\n",
    "    title=\"Legend\",\n",
    "    loc=\"lower right\",\n",
    "    frameon=True,\n",
    "    framealpha=0.7,\n",
    "    edgecolor=\"k\",\n",
    "    facecolor=\"white\",\n",
    ")\n",
    "\n",
    "# save figure\n",
    "# NOTE create figs folder in model root if it does not exist\n",
    "# fn_out = join(mod.root, \"figs\", \"basemap.png\")\n",
    "# plt.savefig(fn_out, dpi=225, bbox_inches=\"tight\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And zoom in to see better some of the structures:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# we assume the model maps are in the projected CRS EPSG:3857\n",
    "proj = ccrs.epsg(3857)\n",
    "# adjust zoomlevel and figure size to your basis size & aspect\n",
    "zoom_level = 14\n",
    "figsize = (10, 8)\n",
    "#extent = np.array([1385000, 1390000, 5861000, 5863000])\n",
    "extent = np.array([1385000, 1390000, 5857000, 5860000])\n",
    "\n",
    "# initialize image with geoaxes\n",
    "fig = plt.figure(figsize=figsize)\n",
    "ax = fig.add_subplot(projection=proj)\n",
    "ax.set_extent(extent, crs=proj)\n",
    "\n",
    "# add sat background image\n",
    "ax.add_image(cimgt.QuadtreeTiles(), zoom_level, alpha=0.5)\n",
    "\n",
    "# plot rivers\n",
    "rivers.plot(ax=ax, linewidth= 1, color=\"blue\", zorder=2, label=\"rivers\")\n",
    "# plot pipes\n",
    "pipes.plot(ax=ax, color=\"k\", linewidth=1, zorder=2, label=\"pipes\")\n",
    "\n",
    "# plot structures\n",
    "manholes.plot(ax=ax, facecolor=\"y\", markersize=2, zorder=4, label=\"manholes\")\n",
    "crosssections.plot(ax=ax, facecolor=\"grey\", marker = '|',  markersize=15, zorder=3, label=\"cross-sections\")\n",
    "\n",
    "ax.xaxis.set_visible(True)\n",
    "ax.yaxis.set_visible(True)\n",
    "ax.set_ylabel(f\"latitude [degree north]\")\n",
    "ax.set_xlabel(f\"longitude [degree east]\")\n",
    "_ = ax.set_title(f\"dflowfm base map\")\n",
    "legend = ax.legend(\n",
    "    title=\"Legend\",\n",
    "    loc=\"lower right\",\n",
    "    frameon=True,\n",
    "    framealpha=0.7,\n",
    "    edgecolor=\"k\",\n",
    "    facecolor=\"white\",\n",
    ")"
   ]
  }
 ],
 "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.18"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}