From af18a0465727460b886de633901ce5fec6b0ce06 Mon Sep 17 00:00:00 2001 From: hboisgon Date: Wed, 15 Mar 2023 12:18:22 +0800 Subject: [PATCH 1/9] support local dem + fix #156 --- hydromt_wflow/wflow.py | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/hydromt_wflow/wflow.py b/hydromt_wflow/wflow.py index 28877ee1..eeb0280a 100644 --- a/hydromt_wflow/wflow.py +++ b/hydromt_wflow/wflow.py @@ -168,14 +168,15 @@ def setup_basemaps( self.logger.info(f"Preparing base hydrography basemaps.") # retrieve global data (lazy!) ds_org = self.data_catalog.get_rasterdataset(hydrography_fn) - # TODO support and test (user) data from other sources with other crs! - if ds_org.raster.crs is None or ds_org.raster.crs.to_epsg() != 4326: - raise ValueError("Only EPSG:4326 base data supported.") + # TODO: add a check on resolution (degree vs meter) depending on ds_org res/crs? # get basin geometry and clip data kind, region = hydromt.workflows.parse_region(region, logger=self.logger) xy = None if kind in ["basin", "subbasin", "outlet"]: - bas_index = self.data_catalog[basin_index_fn] + if basin_index_fn is not None: + bas_index = self.data_catalog[basin_index_fn] + else: + bas_index = None geom, xy = hydromt.workflows.get_basin_geometry( ds=ds_org, kind=kind, @@ -192,6 +193,7 @@ def setup_basemaps( raise ValueError(f"wflow region argument not understood: {region}") if geom is not None and geom.crs is None: raise ValueError("wflow region geometry has no CRS") + ds_org = ds_org.raster.clip_geom(geom, align=res, buffer=10) self.logger.debug(f"Adding basins vector to staticgeoms.") self.set_staticgeoms(geom, name="basins") From b94278cc1989242c25b0d064792ef8180f9ad255 Mon Sep 17 00:00:00 2001 From: hboisgon Date: Thu, 8 Jun 2023 13:35:49 +0800 Subject: [PATCH 2/9] update toml sizeinmetres + crs/res check --- hydromt_wflow/wflow.py | 20 +++++++++++++++++++- 1 file changed, 19 insertions(+), 1 deletion(-) diff --git a/hydromt_wflow/wflow.py b/hydromt_wflow/wflow.py index eeb0280a..c79b095f 100644 --- a/hydromt_wflow/wflow.py +++ b/hydromt_wflow/wflow.py @@ -168,7 +168,7 @@ def setup_basemaps( self.logger.info(f"Preparing base hydrography basemaps.") # retrieve global data (lazy!) ds_org = self.data_catalog.get_rasterdataset(hydrography_fn) - # TODO: add a check on resolution (degree vs meter) depending on ds_org res/crs? + # get basin geometry and clip data kind, region = hydromt.workflows.parse_region(region, logger=self.logger) xy = None @@ -198,6 +198,18 @@ def setup_basemaps( self.logger.debug(f"Adding basins vector to staticgeoms.") self.set_staticgeoms(geom, name="basins") + # Check on resolution (degree vs meter) depending on ds_org res/crs + scale_ratio = int(np.round(res / ds_org.raster.res[0])) + if scale_ratio < 1: + self.logger.error( + f"The model resolution {res} should be larger than the {hydrography_fn} resolution {ds_org.raster.res[0]}" + ) + if ds_org.raster.crs.is_geographic: + if res > 1: # 111 km + self.logger.error( + "The model resolution {res} should be smaller than 1 degree (111km) for geographic coordinate systems." + "Make sure you provided res in degree rather than in meters." + ) # setup hydrography maps and set staticmap attribute with renamed maps ds_base, _ = workflows.hydrography( ds=ds_org, @@ -237,6 +249,12 @@ def setup_basemaps( self.logger.debug(f"Adding region vector to staticgeoms.") self.set_staticgeoms(self.region, name="region") + # update toml for degree/meters if needed + if ds_base.raster.crs.is_geographic: + self.set_config("model.sizeinmetres", False) + else: + self.set_config("model.sizeinmetres", True) + def setup_rivers( self, hydrography_fn, From 29ffcd4ab68f2242141cceecfafa9ea58fb161cf Mon Sep 17 00:00:00 2001 From: hboisgon Date: Fri, 9 Jun 2023 12:29:29 +0800 Subject: [PATCH 3/9] add prepare_ldd notebook --- docs/getting_started/example_index.rst | 3 +- docs/user_guide/process_analyze.rst | 8 +- examples/prepare_ldd.ipynb | 542 +++++++++++++++++++++++++ 3 files changed, 549 insertions(+), 4 deletions(-) create mode 100644 examples/prepare_ldd.ipynb diff --git a/docs/getting_started/example_index.rst b/docs/getting_started/example_index.rst index a5305baa..74cee263 100644 --- a/docs/getting_started/example_index.rst +++ b/docs/getting_started/example_index.rst @@ -26,9 +26,10 @@ For a static (non-interactive) view of the examples follow one of the links belo * `Clip Wflow model <../_examples/clip_model.ipynb>`_ -**Postprocessing and visualization** +**Pre- Post-processing and visualization** * `Convert wflow staticmaps netcdf to raster files <../_examples/convert_staticmaps_to_mapstack.ipynb>`_ +* `Prepare flow directions and related maps from a DEM <../examples/prepare_ldd.ipynb>`_ * `Plot Wflow static maps <../_examples/plot_wflow_staticmaps.ipynb>`_ * `Plot Wflow forcing data <../_examples/plot_wflow_forcing.ipynb>`_ * `Plot Wflow results data <../_examples/plot_wflow_results.ipynb>`_ diff --git a/docs/user_guide/process_analyze.rst b/docs/user_guide/process_analyze.rst index 908d9cac..4b9b3ae0 100644 --- a/docs/user_guide/process_analyze.rst +++ b/docs/user_guide/process_analyze.rst @@ -1,11 +1,12 @@ .. _process_visualize: -================================ -Postprocessing and visualization -================================ +======================================== +Pre and postprocessing and visualization +======================================== The Hydromt-Wflow plugin provides several possibilities to postprocess and visualize the model data and model results: +* `Prepare flow directions and related maps from a DEM <../examples/prepare_ldd.ipynb>`_ using the `flw methods of HydroMT `_ * `Convert Wflow staticmaps netcdf to raster files <../_examples/convert_staticmaps_to_mapstack.ipynb>`_ for further processing and analyzation * Plot `staticmaps <../_examples/plot_wflow_staticmaps.ipynb>`_, `forcing data <../_examples/plot_wflow_forcing.ipynb>`_ and `model results <../_examples/plot_wflow_results.ipynb>`_ by means of additional python packages @@ -15,6 +16,7 @@ The Hydromt-Wflow plugin provides several possibilities to postprocess and visua .. toctree:: :hidden: + Example: Prepare flow directions and related maps from a DEM <../examples/prepare_ldd.ipynb> Example: Convert wflow staticmaps netcdf to raster files <../_examples/convert_staticmaps_to_mapstack.ipynb> Example: Plot Wflow staticmaps <../_examples/plot_wflow_staticmaps.ipynb> Example: Plot Wflow forcing data <../_examples/plot_wflow_forcing.ipynb> diff --git a/examples/prepare_ldd.ipynb b/examples/prepare_ldd.ipynb new file mode 100644 index 00000000..85513af6 --- /dev/null +++ b/examples/prepare_ldd.ipynb @@ -0,0 +1,542 @@ +{ + "cells": [ + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Prepare flow directions and related data from a DEM" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "With HydroMT-Wflow, a user can choose to build a model in a geographic or projected coordinate system from an input Digital Elevation Model (DEM) and Flow Direction (flwdir) dataset.\n", + "\n", + "While DEM data are often available, this is not the always the case for the flow directions. In the plugin, we made the choice to build a Wflow model directly from user provided DEM and flwdir datasets rather than reprojecting a DEM and/or deriving flwdir on the fly. This is because there are a lot of available techniques to derive flow directions and we want the user to be sure the flow directions matches the terrain and user expectations.\n", + "\n", + "Because of this, we prefer to provide this notebook as a possible pre-processing step before calling a hydromt build wflow command. We will do this using the different [flow directions methods from HydroMT](https://deltares.github.io/hydromt/latest/api.html#flow-direction-methods) and [PyFlwDir](https://deltares.github.io/pyflwdir/latest/index.html)" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Load dependencies" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import os\n", + "import xarray as xr\n", + "import numpy as np\n", + "import geopandas as gpd\n", + "# pyflwdir\n", + "import pyflwdir\n", + "# hydromt\n", + "from hydromt import DataCatalog, flw\n", + "#plot\n", + "import matplotlib.pyplot as plt\n", + "from matplotlib import cm, colors\n", + "import cartopy.crs as ccrs" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Common plot settings\n", + "plt.style.use(\"seaborn-whitegrid\") # set nice style\n", + "# we assume the model maps are in the geographic CRS EPSG:4326\n", + "proj = ccrs.PlateCarree()\n", + "\n", + "# create nice elevation colormap\n", + "c_dem = plt.cm.terrain(np.linspace(0.25, 1, 256))\n", + "cmap = colors.LinearSegmentedColormap.from_list(\"dem\", c_dem)\n", + "norm = colors.Normalize(vmin=0, vmax=2000)\n", + "kwargs = dict(cmap=cmap, norm=norm)\n", + "\n", + "# legend settings\n", + "legend_kwargs = dict(\n", + " title=\"Legend\",\n", + " loc=\"lower right\",\n", + " frameon=True,\n", + " framealpha=0.7,\n", + " edgecolor=\"k\",\n", + " facecolor=\"white\",\n", + ")\n", + "\n", + "def add_legend_titles(ax, title, add_legend=True, projected=True):\n", + " ax.xaxis.set_visible(True)\n", + " ax.yaxis.set_visible(True)\n", + " if projected:\n", + " ax.set_xlabel(\"Easting [m]\")\n", + " ax.set_ylabel(\"Northing [m]\")\n", + " else:\n", + " ax.set_ylabel(\"latitude [degree north]\")\n", + " ax.set_xlabel(\"longitude [degree east]\")\n", + " _ = ax.set_title(title)\n", + " if add_legend:\n", + " legend = ax.legend(**legend_kwargs)\n", + "\n", + " return legend" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Elevation data and reprojection" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "In this example we will use the `merit_hydro_1k` data in the pre-defined `artifact_data` catalog of HydroMT. \n", + "\n", + "If you prefer to have a Wflow model in projected coordinate system (CRS in metres), you will need first to reproject the DEM and we will see here an example of how to do this.\n", + "\n", + "First let's read the data:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "data_catalog = DataCatalog(\"artifact_data\")\n", + "merit = data_catalog.get_rasterdataset(\"merit_hydro_1k\", bbox = [11.8,46.0,12.3,46.2])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# plot\n", + "# initialize image with geoaxes\n", + "fig = plt.figure(figsize=(10,8))\n", + "ax = fig.add_subplot(projection=proj)\n", + "\n", + "## plot elevation\\\n", + "merit['elevtn'].plot(\n", + " transform=proj, ax=ax, zorder=1, cbar_kwargs=dict(aspect=30, shrink=0.8), **kwargs\n", + ")\n", + "\n", + "# plot flwdir\n", + "flwdir = flw.flwdir_from_da(merit[\"flwdir\"], ftype=\"infer\", check_ftype=True)\n", + "feat = flwdir.streams()\n", + "gdf = gpd.GeoDataFrame.from_features(feat, crs=merit.raster.crs)\n", + "gdf.plot(ax=ax, color=\"blue\", linewidth=0.5, zorder=2, label=\"Flow directions\")\n", + "\n", + "legend = add_legend_titles(ax, \"MERIT Hydro IHU 1km\", projected=False)" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "And now let's reproject the DEM to the projected CRS EPSG 3857 using [hydromt.DataArray.raster.reproject](https://deltares.github.io/hydromt/latest/_generated/hydromt.DataArray.raster.reproject.html#hydromt.DataArray.raster.reproject) method." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "merit_elv_utm = merit[\"elevtn\"].raster.reproject(\n", + " dst_crs=3857,\n", + " dst_res = None, # keep the same resolution\n", + " method = 'bilinear',\n", + ")" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Reprojecting flow directions" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "If your original data did contain flow directions, you can decide to reproject rather than re-create them from the reprojected DEM.\n", + "\n", + "To do this, you can use the [hydromt.flw.reproject_hydrography_like](https://deltares.github.io/hydromt/latest/_generated/hydromt.flw.reproject_hydrography_like.html#hydromt.flw.reproject_hydrography_like) method from HydroMT.\n", + "\n", + "This method aims to conserve flow direction during reprojection by deriving flwdir from a reprojected grid of synthetic elevation, based on the log10 of upstream area [m2]." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "merit_utm_reproj = flw.reproject_hydrography_like(\n", + " ds_hydro = merit[[\"flwdir\", \"uparea\"]],\n", + " da_elv = merit_elv_utm,\n", + ")" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Let's compare our new reprojected flow directions to the original ones" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# plot\n", + "# initialize image with geoaxes\n", + "fig = plt.figure(figsize=(10,8))\n", + "ax = fig.add_subplot()\n", + "\n", + "## plot elevation\\\n", + "merit_elv_utm.plot(\n", + " ax=ax, zorder=1, cbar_kwargs=dict(aspect=30, shrink=0.8), **kwargs\n", + ")\n", + "\n", + "# plot flwdir\n", + "flwdir = flw.flwdir_from_da(merit[\"flwdir\"], ftype=\"infer\", check_ftype=True)\n", + "feat = flwdir.streams(min_sto=2)\n", + "gdf = gpd.GeoDataFrame.from_features(feat, crs=merit.raster.crs)\n", + "gdf.to_crs(merit_elv_utm.raster.crs).plot(ax=ax, column=\"strord\", cmap=colors.ListedColormap(cm.Blues(np.linspace(0.4, 1, 7))), label=\"Flow directions\")\n", + "\n", + "# plot flwdir reproj\n", + "flwdirr = flw.flwdir_from_da(merit_utm_reproj[\"flwdir\"], ftype=\"infer\", check_ftype=True)\n", + "featr = flwdirr.streams(min_sto=2)\n", + "gdfr = gpd.GeoDataFrame.from_features(featr, crs=merit_utm_reproj.raster.crs)\n", + "gdfr.plot(ax=ax, column=\"strord\", cmap=colors.ListedColormap(cm.Reds(np.linspace(0.4, 1, 7))), label=\"Flow directions reprojected\")\n", + "\n", + "legend = add_legend_titles(ax, \"MERIT Hydro IHU 1km Reprojected\")" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Deriving flow directions from DEM" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "To derive flow directions from a DEM, you can use the [hydromt.flw.d8_from_dem](https://deltares.github.io/hydromt/latest/_generated/hydromt.flw.d8_from_dem.html#hydromt.flw.d8_from_dem) method of HydroMT.\n", + "\n", + "This method derives D8 flow directions grid from an elevation grid and allows several options to the users:\n", + " - **outlets**: outlets can be defined at ``edge``s of the grid or force all flow to go to the minimum elevation point ``min``. Additionnally, the user can also specify the pits locations via ``idxs_pit``.\n", + " - **river burning**: it is possible to provide a river vector layer ``gdf_stream`` with ``uparea`` (km2) column which is used to burn the river in the elevation data.\n", + " - **depression filling**: local depressions are filled based on their lowest pour point level if the pour point depth is smaller than the maximum pour point depth ``max_depth``, otherwise the lowest elevation in the depression becomes a pit.\n", + "\n", + "Let's see an example:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Start a merit_utm dataset\n", + "merit_utm_flw = merit_elv_utm.to_dataset(name=\"elevtn\")\n", + "# Derive flow directions with outlets at the edges\n", + "merit_utm_flw[\"flwdir_edge\"] = flw.d8_from_dem(\n", + " da_elv = merit_elv_utm,\n", + " gdf_stream = None,\n", + " max_depth=-1,\n", + " outlets = \"edge\",\n", + " idxs_pit = None,\n", + ")\n", + "# Derive flow directions with outlet at the min elevation edge cell\n", + "merit_utm_flw[\"flwdir_min\"] = flw.d8_from_dem(\n", + " da_elv = merit_elv_utm,\n", + " gdf_stream = None,\n", + " max_depth=-1,\n", + " outlets = \"min\",\n", + " idxs_pit = None,\n", + ")" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Let's plot all the different methods" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# plot\n", + "# initialize image with geoaxes\n", + "fig = plt.figure(figsize=(10,8))\n", + "ax = fig.add_subplot()\n", + "\n", + "## plot elevation\\\n", + "merit_elv_utm.plot(\n", + " ax=ax, zorder=1, cbar_kwargs=dict(aspect=30, shrink=0.8), alpha=0.5, **kwargs\n", + ")\n", + "\n", + "# plot flwdir\n", + "flwdir = flw.flwdir_from_da(merit[\"flwdir\"], ftype=\"infer\", check_ftype=True)\n", + "feat = flwdir.streams(min_sto=2)\n", + "gdf = gpd.GeoDataFrame.from_features(feat, crs=merit.raster.crs)\n", + "gdf.to_crs(merit_elv_utm.raster.crs).plot(ax=ax, column=\"strord\", cmap=colors.ListedColormap(cm.Blues(np.linspace(0.4, 1, 7))), label=\"Flow directions\")\n", + "\n", + "# plot flwdir reproj\n", + "flwdirr = flw.flwdir_from_da(merit_utm_reproj[\"flwdir\"], ftype=\"infer\", check_ftype=True)\n", + "featr = flwdirr.streams(min_sto=2)\n", + "gdfr = gpd.GeoDataFrame.from_features(featr, crs=merit_utm_reproj.raster.crs)\n", + "gdfr.plot(ax=ax, column=\"strord\", cmap=colors.ListedColormap(cm.Reds(np.linspace(0.4, 1, 7))), label=\"Flow directions reprojected\")\n", + "\n", + "# plot flwdir edge\n", + "flwdir_edge = flw.flwdir_from_da(merit_utm_flw[\"flwdir_edge\"], ftype=\"infer\", check_ftype=True)\n", + "feate = flwdir_edge.streams(min_sto=2)\n", + "gdfe = gpd.GeoDataFrame.from_features(feate, crs=merit_utm_flw.raster.crs)\n", + "gdfe.plot(ax=ax, column=\"strord\", cmap=colors.ListedColormap(cm.Greens(np.linspace(0.4, 1, 7))), label=\"Flow directions edge\")\n", + "\n", + "# plot flwdir min\n", + "flwdir_min = flw.flwdir_from_da(merit_utm_flw[\"flwdir_min\"], ftype=\"infer\", check_ftype=True)\n", + "featm = flwdir_min.streams(min_sto=2)\n", + "gdfm = gpd.GeoDataFrame.from_features(featm, crs=merit_utm_flw.raster.crs)\n", + "gdfm.plot(ax=ax, column=\"strord\", cmap=colors.ListedColormap(cm.Purples(np.linspace(0.4, 1, 7))), label=\"Flow directions min\")\n", + "\n", + "legend = add_legend_titles(ax, \"MERIT Hydro IHU 1km Reprojected\")" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now let's see an example where we want to burn rivers.\n", + "In the data artifact, we have one river vector database : ``rivers_lin2019_v1``.\n", + "\n", + "Let's use it for burning." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Read the data from data catalog\n", + "rivers = data_catalog.get_geodataframe(\"rivers_lin2019_v1\", bbox=[11.8,46.0,12.3,46.2])\n", + "\n", + "# In this dataset, there is an Area column representing upstream area in km2\n", + "# We need to rename it to uparea\n", + "rivers = rivers.rename(columns={\"Area\": \"uparea\"})\n", + "# And finally make sure rivers is reprojected to the same CRS as the elevation\n", + "rivers = rivers.to_crs(merit_elv_utm.raster.crs)\n", + "\n", + "# Now let's use it to derive flow directions\n", + "merit_utm_flw[\"flwdir_riverburn\"] = flw.d8_from_dem(\n", + " da_elv = merit_elv_utm,\n", + " gdf_stream = rivers,\n", + " max_depth=-1,\n", + " outlets = \"edge\",\n", + " idxs_pit = None,\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# plot\n", + "# initialize image with geoaxes\n", + "fig = plt.figure(figsize=(10,8))\n", + "ax = fig.add_subplot()\n", + "\n", + "## plot elevation\\\n", + "merit_elv_utm.plot(\n", + " ax=ax, zorder=1, cbar_kwargs=dict(aspect=30, shrink=0.8), alpha=0.5, **kwargs\n", + ")\n", + "\n", + "# plot flwdir\n", + "flwdir = flw.flwdir_from_da(merit[\"flwdir\"], ftype=\"infer\", check_ftype=True)\n", + "feat = flwdir.streams(min_sto=2)\n", + "gdf = gpd.GeoDataFrame.from_features(feat, crs=merit.raster.crs)\n", + "gdf.to_crs(merit_elv_utm.raster.crs).plot(ax=ax, column=\"strord\", cmap=colors.ListedColormap(cm.Blues(np.linspace(0.4, 1, 7))), label=\"Flow directions\")\n", + "\n", + "# plot flwdir reproj\n", + "flwdirr = flw.flwdir_from_da(merit_utm_reproj[\"flwdir\"], ftype=\"infer\", check_ftype=True)\n", + "featr = flwdirr.streams(min_sto=2)\n", + "gdfr = gpd.GeoDataFrame.from_features(featr, crs=merit_utm_reproj.raster.crs)\n", + "gdfr.plot(ax=ax, column=\"strord\", cmap=colors.ListedColormap(cm.Reds(np.linspace(0.4, 1, 7))), label=\"Flow directions reprojected\")\n", + "\n", + "# plot flwdir edge\n", + "flwdir_edge = flw.flwdir_from_da(merit_utm_flw[\"flwdir_edge\"], ftype=\"infer\", check_ftype=True)\n", + "feate = flwdir_edge.streams(min_sto=2)\n", + "gdfe = gpd.GeoDataFrame.from_features(feate, crs=merit_utm_flw.raster.crs)\n", + "gdfe.plot(ax=ax, column=\"strord\", cmap=colors.ListedColormap(cm.Greens(np.linspace(0.4, 1, 7))), label=\"Flow directions edge\")\n", + "\n", + "# plot flwdir riverburn\n", + "flwdir_min = flw.flwdir_from_da(merit_utm_flw[\"flwdir_riverburn\"], ftype=\"infer\", check_ftype=True)\n", + "featm = flwdir_min.streams(min_sto=2)\n", + "gdfm = gpd.GeoDataFrame.from_features(featm, crs=merit_utm_flw.raster.crs)\n", + "gdfm.plot(ax=ax, column=\"strord\", cmap=colors.ListedColormap(cm.Purples(np.linspace(0.4, 1, 7))), label=\"Flow directions river burning\")\n", + "\n", + "rivers.plot(ax=ax, color=\"black\", label=\"Rivers Lin\")\n", + "\n", + "legend = add_legend_titles(ax, \"MERIT Hydro IHU 1km Reprojected\")" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Deriving other DEM and flow directions related data" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Once you are satisfied with your flow direction map, you can create additionnal derived variables like upstream area or streamorder that can prove useful for example to build a model based on ``subbasin`` region.\n", + "\n", + "Here are some examples how to do that using PyFlwdir methods." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Create a dataset with the riverburn flow directions\n", + "merit_utm = merit_elv_utm.to_dataset(name=\"elevtn\")\n", + "merit_utm[\"flwdir\"] = merit_utm_flw[\"flwdir_riverburn\"]\n", + "dims = merit_utm.raster.dims\n", + "\n", + "# Create a PyFlwDir object from the dataset\n", + "flwdir = flw.flwdir_from_da(merit_utm[\"flwdir\"])\n", + "\n", + "# uparea\n", + "uparea = flwdir.upstream_area(unit=\"km2\")\n", + "merit_utm[\"uparea\"] = xr.Variable(dims, uparea, attrs=dict(_FillValue=-9999))\n", + "\n", + "# stream order\n", + "strord = flwdir.stream_order()\n", + "merit_utm[\"strord\"] = xr.Variable(dims, strord)\n", + "merit_utm[\"strord\"].raster.set_nodata(255)\n", + "\n", + "# slope\n", + "slope = pyflwdir.dem.slope(\n", + " elevtn=merit_utm[\"elevtn\"].values,\n", + " nodata = merit_utm[\"elevtn\"].raster.nodata,\n", + " latlon = False, # True if geographic crs, False if projected crs\n", + " transform = merit_utm[\"elevtn\"].raster.transform,\n", + ")\n", + "merit_utm[\"slope\"] = xr.Variable(dims, slope)\n", + "merit_utm[\"slope\"].raster.set_nodata(merit_utm[\"elevtn\"].raster.nodata)\n", + "\n", + "# basin at the pits locations\n", + "basins = flwdir.basins(idxs=flwdir.idxs_pit).astype(np.int32)\n", + "merit_utm[\"basins\"] = xr.Variable(dims, basins, attrs=dict(_FillValue=0))\n", + "\n", + "merit_utm" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# plot\n", + "fig, axes = plt.subplots(3, 2, figsize=(15, 15))\n", + "\n", + "## plot elevation\n", + "merit_utm[\"elevtn\"].plot(\n", + " ax=axes[0,0], zorder=1, cbar_kwargs=dict(aspect=30, shrink=0.8), alpha=0.5, **kwargs\n", + ")\n", + "_ = add_legend_titles(axes[0,0], \"Elevation [m asl]\", add_legend=False)\n", + "\n", + "# plot flwdir riverburn\n", + "flwdir = flw.flwdir_from_da(merit_utm[\"flwdir\"], ftype=\"infer\", check_ftype=True)\n", + "feat = flwdir.streams(min_sto=2)\n", + "gdf = gpd.GeoDataFrame.from_features(feat, crs=merit_utm.raster.crs)\n", + "gdf.to_crs(merit_elv_utm.raster.crs).plot(ax=axes[0,1], column=\"strord\", cmap=colors.ListedColormap(cm.Blues(np.linspace(0.4, 1, 7))), label=\"Flow directions\")\n", + "_ = add_legend_titles(axes[0,1], \"Flow directions\", add_legend=False)\n", + "\n", + "# plot uparea\n", + "merit_utm[\"uparea\"].plot(ax=axes[1,0])\n", + "_ = add_legend_titles(axes[1,0], \"Upstream area [km2]\", add_legend=False)\n", + "\n", + "# plot strord\n", + "merit_utm[\"strord\"].plot(ax=axes[1,1], cmap=colors.ListedColormap(cm.Blues(np.linspace(0.4, 1, 7))))\n", + "_ = add_legend_titles(axes[1,1], \"Strahler Stream order\", add_legend=False)\n", + "\n", + "# plot slope\n", + "merit_utm[\"slope\"].plot(ax=axes[2,0])\n", + "_ = add_legend_titles(axes[2,0], \"Slope [m/m]\", add_legend=False)\n", + "\n", + "# plot basins\n", + "merit_utm[\"basins\"].plot(ax=axes[2,1])\n", + "_ = add_legend_titles(axes[2,1], \"Basins ID\", add_legend=False)\n" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "hydromt-wflow", + "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.11.3" + }, + "orig_nbformat": 4 + }, + "nbformat": 4, + "nbformat_minor": 2 +} From a1738ee376bc0c057fd4756fe31b34e6a4466c75 Mon Sep 17 00:00:00 2001 From: hboisgon Date: Fri, 9 Jun 2023 14:16:16 +0800 Subject: [PATCH 4/9] update docs/changelog --- docs/changelog.rst | 6 ++++++ docs/user_guide/wflow_build.rst | 5 +++++ hydromt_wflow/wflow.py | 9 ++++++++- 3 files changed, 19 insertions(+), 1 deletion(-) diff --git a/docs/changelog.rst b/docs/changelog.rst index 229cd430..0ce7dc7e 100644 --- a/docs/changelog.rst +++ b/docs/changelog.rst @@ -11,6 +11,8 @@ Unreleased Added ----- +- Support for models in CRS other than 4326. `PR #161 `_ +- Support for elevation data other than MERIT Hydro in **setup_basemaps**. - Support in toml for dir_input and dir_output options. `PR #140 `_ - Add options to calculate daily Penman-Monteith potential evaporation using the pyet package. Depending on the available variables, two options are defined ``penman-monteith_tdew`` (inputs: ['temp', 'temp_min', 'temp_max', 'wind_u', 'wind_v', 'temp_dew', 'kin', 'press_msl']) and ``penman-monteith_rh_simple`` (inputs: ['temp', 'temp_min', 'temp_max', 'wind', 'rh', 'kin']). - In **setup_reservoirs**: Global Water Watch compatibility for determining reservoir parameters. @@ -42,6 +44,10 @@ Deprecated ---------- - The **setup_hydrodem** function has been removed, and the functionality are moved to **setup_rivers** and **setup_floodplains** +Documentation +------------- +- New **prepare_ldd** example notebook to demonstrate how to prepare flow directions and other elevation related data. + v0.2.1 (22 November 2022) ========================= diff --git a/docs/user_guide/wflow_build.rst b/docs/user_guide/wflow_build.rst index aff2de4d..6764009e 100644 --- a/docs/user_guide/wflow_build.rst +++ b/docs/user_guide/wflow_build.rst @@ -21,6 +21,11 @@ for a proper implementation of this model are: - basin - subbasin +The coordinate reference system (CRS) of the model will be the same as the one of the input hydrography data. If the region +is specified using point coordinates or a bounding box, the coordinates used should match the CRS of the hydrography data. +If the user wants to use a different CRS, we advise to reproject the hydrography data to the desired CRS before building the model. +You can find some examples on how to do this in the `example notebook <../_examples/prepare_ldd.ipynb>`_. + .. _model_config: Configuration file diff --git a/hydromt_wflow/wflow.py b/hydromt_wflow/wflow.py index b83143fa..04de277f 100644 --- a/hydromt_wflow/wflow.py +++ b/hydromt_wflow/wflow.py @@ -127,7 +127,14 @@ def setup_basemaps( Iterative Hydrography Upscaling (IHU). The default ``hydrography_fn`` is "merit_hydro" (`MERIT hydro `_ at 3 arcsec resolution) Alternative sources include "merit_hydro_1k" at 30 arcsec resolution. - Users can also supply their own elevation and flow direction data. Note that only EPSG:4326 base data supported. + Users can also supply their own elevation and flow direction data in any CRS and not only EPSG:4326. + + Note that in order to define the region, using points or bounding box, the coordinates of the points / bounding box + should be in the same CRS than the hydrography data. The wflow model will then also be in the same CRS than the + hydrography data in order to avoid assumptions and reprojection errors. If the user wishes to use a different CRS, + we recommend first to reproject the hydrography data seperately because calling hydromt build. You can find examples + on how to reproject or prepare hydrography data in the + `prepare flow directions example notebok `. Adds model layers: From 81c0e3ffab7f31f2e76d8a50a8e146aaf5b7dd91 Mon Sep 17 00:00:00 2001 From: hboisgon Date: Fri, 9 Jun 2023 17:20:32 +0800 Subject: [PATCH 5/9] add tests --- hydromt_wflow/wflow.py | 10 +++--- tests/data/merit_utm/elevtn.tif | Bin 0 -> 6542 bytes tests/data/merit_utm/flwdir.tif | Bin 0 -> 1915 bytes tests/data/merit_utm/merit_utm.yml | 5 +++ tests/data/merit_utm/uparea.tif | Bin 0 -> 6542 bytes tests/test_model_methods.py | 50 +++++++++++++++++++++++++++++ 6 files changed, 59 insertions(+), 6 deletions(-) create mode 100644 tests/data/merit_utm/elevtn.tif create mode 100644 tests/data/merit_utm/flwdir.tif create mode 100644 tests/data/merit_utm/merit_utm.yml create mode 100644 tests/data/merit_utm/uparea.tif diff --git a/hydromt_wflow/wflow.py b/hydromt_wflow/wflow.py index 04de277f..efb03df6 100644 --- a/hydromt_wflow/wflow.py +++ b/hydromt_wflow/wflow.py @@ -209,13 +209,13 @@ def setup_basemaps( # Check on resolution (degree vs meter) depending on ds_org res/crs scale_ratio = int(np.round(res / ds_org.raster.res[0])) if scale_ratio < 1: - self.logger.error( + raise ValueError( f"The model resolution {res} should be larger than the {hydrography_fn} resolution {ds_org.raster.res[0]}" ) if ds_org.raster.crs.is_geographic: if res > 1: # 111 km - self.logger.error( - "The model resolution {res} should be smaller than 1 degree (111km) for geographic coordinate systems." + raise ValueError( + f"The model resolution {res} should be smaller than 1 degree (111km) for geographic coordinate systems. " "Make sure you provided res in degree rather than in meters." ) # setup hydrography maps and set staticmap attribute with renamed maps @@ -258,9 +258,7 @@ def setup_basemaps( self.set_staticgeoms(self.region, name="region") # update toml for degree/meters if needed - if ds_base.raster.crs.is_geographic: - self.set_config("model.sizeinmetres", False) - else: + if ds_base.raster.crs.is_projected: self.set_config("model.sizeinmetres", True) def setup_rivers( diff --git a/tests/data/merit_utm/elevtn.tif b/tests/data/merit_utm/elevtn.tif new file mode 100644 index 0000000000000000000000000000000000000000..7689d7dd404546fc762ae35544551632d8f31e05 GIT binary patch literal 6542 zcmciGd00+)8wc=HJuOmAiWV`|w4i9wt~&R9pOm7)G_s{gnN*4*8N8~ejrIjuk|j#n zVl0Irkq~NnEuo1PDMXAA<$dz%y8JPQ8RdDd-*tbVbMEu|oZmU;I)6M44mykk!!WW8 zBP7BI3DOY35c-hK1X=h)Rug2A54pQf?4vwDi1+^3{|I&r`e?hthul4n;)fntC4uo# zkIRQFBIuFfOEKN`ieN-UN-<)BdS5X1$A~a$f~p{xzaU1KQ5R%aK@QIre%})+$nQT3 z*VgIPr6SW)KWzR~$#u}G>*AYGg}tuw1zA_>=@WXt8gYi<39d!=_d^xt6~rhDnuH_@ zV#FC`Moe(&gqTDPMyk6G_>T9!j?kDo&PC1KSWRDTQn2?DufS1G-a&KcEDj7>@&6qw z7-p1(;KMLhLrQyN%<^5s^~OA3GKa^kaIPgU{F`};HZA6Sl^1^#Z+zx#ZbipDE`Rv^ zZ{js%P2s*Pn9i;9b^j*b?qLqxFcD|Y#>nNHc=jV$t|-BVOZT(ybMLXNF?TP{gp=~L z;2a}aZgw(pX@_k3)YVJS;}$+O;IiXQIPJ6M+-yV3KJ`??v^ZJ+;hc1vF88!vkK@-F z_L(PqZZM}gU4?V&7|M}_dY?Hx@0_>Cq`7`qCMCvGVm^;2#GSXno{1Yp6G@wDoWZ&=fg6|GbxkZqSwzu#^F$g} zm-5x#>V1pIt1Fn!I)qU05D${Cm_p$p?o=t`Pya0rCa2j8=s%VV$i*X?^g0s0%9}jY zo8~_Xpv>wZlDY3j_Tpxw73x6kIFAsobh;4YLV1g)QP=SBe{Ff0Pb8xvSF#E7 zpo3fHQChQDe0dTRvP{b)%+JCV99%}KNXsl||^hgr6C!zhrJ6)gJ4rlf=g-FTvknPg~0Vd^YZKjJ95-&hJP;3%Y=qgQ8#lJz}?R^^IO#f3A9sJ@1H5RrkV9OmJ^T^Bd3c<#2}LW8UGN(N8dI`va6c)s7>dRimzTGG19w zg^QheB%XT;`PYjvWo83fZF`8C5-+f+wH0e+o`2RF(9L|i{g*vvsqUnxuwEnXSE!W&e*{L<3 z_13>QiN5Oj*gQKQ&xIbw($^<&@#GT0x=!OplT+9xUyea86_|K_8)`4g!lVgBDC=2< zCq|#cyJ}_FGWaU$x1Yx&J4;a_`uyK|+w%9}D6e#Ut(lBw-f1|}Bp#)@5^>al4Okw# z5q-onuyDgREVz|{iknlh7?z=Iyf+3K?#5x4kKwSe-N=vFh502pxOKojG+B4-Z#|pP zRCKfVM4zo`qdVD_$ z+i@aR?AAp)w35`%( zCxS`W#L+%V0xzWrPzC5yG@Xwc>x0Cj>6;Q6QA|Wt-LVS zQw5p2E|7oJ1pUJ*Aoa%YkQI;xIT2N`-mwzcMP(o#cpS#$+yuVi6L@v20#3ia4E^hC zptWruTw0a}oqNOKnn~PeJw37ccw>Pv#>xp{*`-n_lgo#8t#GKCw*V|n)*%4O2Q67;<;_`EO9S5)$WD210sRUr@{J1M&P6249ic&K;fPppLFYgaYRZs zN1d4GFf8H(c$sa1O|?tFI&>DCxIG&pTjC(kY7Ox0=E7KB2t*fogQ$@o<(6IJdij z*PS_FuQwmsGTq@MYX^;kIXJdw6!;Gog9nZ?AvuZ#rgS(gy*UQ11%D5^6P)4A>PVPq zupU~ycmJi=^BlQ_zJOS@eeh(MEy(8?!C1F3knb@Cl3%()i1~P^kJACG)_x$Qe~+~} z@qiWmT>;wcOyJtP@i2lP0=h0gfw5lz$m^!V)6Oh-W>omMUbkE8R{`0RGa%+nD2R*s zLg|jBFrhUBjKijYZCyJ%=S?LmQgwk{Yth2Cw2g)P^)tXi+5qZx{XxHN4Fv4p0s-51 zgL~=$$n3x6vtG|yJwppY-!mChGv~r-A#4+8JZaLBpn_r+dDeIYp5 zXTI;(_A!K)HwQz|F(nvyXei7qR)NuvE(y4|)oJ6NJ?1S+?c!F4SICcgdQj+7b1#JfP#Pud_OA_FxZT`XrN z)GOY{Yuynh3v%l^*~Jt4fl>BA*x@P!x1SZTKR>Z#n+;F0Hj0hx&NFA(kB|A<{Q7PM z8?)jWYx{i@>$N4DX4JZ%FF@JW<_!l!;vqq`M~=J@vd{rh-*-8n9L!CxBx10YafAOHXW literal 0 HcmV?d00001 diff --git a/tests/data/merit_utm/flwdir.tif b/tests/data/merit_utm/flwdir.tif new file mode 100644 index 0000000000000000000000000000000000000000..1194e8af524ceb89f608692bcb5348401c0cc923 GIT binary patch literal 1915 zcmcJQziL!L6vofY4x1?!ofIZ4SWE~ODk?}s%bQBK2m}I2Et)zK6`K@SBW!0O_y!hs zvX7BQe1y2_NBcjY-JM>!!YcwLLyD-W)V4c|(tuu@aMa;S?<7KsMMTZjFNX|G$`G_K! z){3F9mB&{SwmEgNZIMaN*bAbXcDP)u7a2Zh#*UoOSjRn=i2D?7op&%>Cc?zx_=#~P zPKy2&-q^Oem6kOf<$4OiNtYPbkjvGu8_Y?${smxrcYfOG{Ibq>X;^WOqOp06WgEw>~4E_+F%;$nu*@B`TT8Vf&u X7^ljcY5dUuz*%#qkHY#<2$g_8+YpkD literal 0 HcmV?d00001 diff --git a/tests/data/merit_utm/merit_utm.yml b/tests/data/merit_utm/merit_utm.yml new file mode 100644 index 00000000..38105da6 --- /dev/null +++ b/tests/data/merit_utm/merit_utm.yml @@ -0,0 +1,5 @@ +merit_hydro_1k_utm: + path: ./{variable}.tif + data_type: RasterDataset + driver: raster + crs: 3857 \ No newline at end of file diff --git a/tests/data/merit_utm/uparea.tif b/tests/data/merit_utm/uparea.tif new file mode 100644 index 0000000000000000000000000000000000000000..b2d19b19cf124e6c88056c4e645e94031f27edc4 GIT binary patch literal 6542 zcmcgvU1(fI6rSBp(+EpVMH^!`uo6XDi3O98=A!Q1-9#JG){ zwh_~Um{JTrDN1Tm=|fymTDBPDORZ3<6r&cjNFEl%{=8IBLGgSubLQTSo4t4G?j{T~ znK^U5^PRst>FC%fs)P{LLIgr0K&uFifL^!II;huCT8H!+a+R+2+CYW8dx#G1bM3Fw zYv`%bGS&AJhD)YfuR|oWW;`q~U#4;Njj*VsIYVP^A|#?TuOt1biJ*9p);+X7aw;fg zMrkd7UY%NcZ0V#ec-Gr1O{Y5^Uix>uIsV<2M{B>S?uoNpK(bjaLkZf=ThitgAuwD0BNhHXP5gZ-}!j~rA7 z2Zd;8r5_=lh+HVf=)awg6+@379sRd0yS#Ym7~$D|Y`aq#a2regRyKZw7cZwjL&px6 z6>o8Bx@>ss+vnmPEpJimRDL2yOP!|+eSc_vz$|@TbhLyM;K-wt{I0LLCxN=-ai_KA zb$@;PrUYuQ9Q8LJX(m60DPO`SKTG9H%8SI}RDB%7$xo{!{jB1WO)--_g?Oc9o`4qy zu)ShX==qm@`w}<0y5c%!x(2f9uc(;dF>_=_=801Ss`g5eH`;fm4j9v&H--8vpFHq7 z4_uP`Xf11RdEP)KpN>M#ZgpLeb>P$;a!2J&De?ujfA!R2;`H=|cy``Mu0`$Fz`~gA ztxS%kcFS7l94Th}7)Q;PeIiqnR53$uQ8w%+{9=3`SO;S@&DSMS@4#cM$;`W^&lHji zZdvF;{Fmc3f7?KK4kz_}&+apBJ$}|2L-#6o@38EvGmQpsGE?KK3sv)Miys`8(Y``$ z$^4*MuEF>9L#^QQMV*)$c6O$YpgD++E8Ml3GyVs1^JN7GE z_v`m8a(=3Ii}5GD1(*4g{h{ph!ooS(_&;uhulFqb=ZoD2Z~pIhj6av#h&x_*@W+j= zD_Q&yGxerbTn~S+)8LKm`512{Qw!RkHu{bqEXaqB34Bjoxgytw^FyP(gs0ALQZL0F zZR&U-iP{dymi8NgV^~yJzw7)b8q>1z7R4RN<>Q7kwprJ2OvTH=o1K0-Ctl|FPYoV$ z)SDFX_x@j3W7166&$rBz@Zz)7E0C)Z>R~qr2a*faS6~OmXQx)R zLH6X Date: Tue, 25 Jul 2023 13:28:24 +0800 Subject: [PATCH 6/9] final small adjustment --- .gitignore | 1 + examples/prepare_ldd.ipynb | 108 ++++++++++++++++++++++++++++++++++++- 2 files changed, 107 insertions(+), 2 deletions(-) diff --git a/.gitignore b/.gitignore index b4b69bab..ad6854ad 100644 --- a/.gitignore +++ b/.gitignore @@ -87,6 +87,7 @@ examples/wflow_piave_forcing examples/wflow_piave_gauges examples/wflow_test_sediment examples/wflow_test_extend_sediment +examples/elevation_data # pyenv .python-version diff --git a/examples/prepare_ldd.ipynb b/examples/prepare_ldd.ipynb index 85513af6..6f1526f6 100644 --- a/examples/prepare_ldd.ipynb +++ b/examples/prepare_ldd.ipynb @@ -273,7 +273,7 @@ "merit_utm_flw[\"flwdir_edge\"] = flw.d8_from_dem(\n", " da_elv = merit_elv_utm,\n", " gdf_stream = None,\n", - " max_depth=-1,\n", + " max_depth=-1, # no local pits\n", " outlets = \"edge\",\n", " idxs_pit = None,\n", ")\n", @@ -281,7 +281,7 @@ "merit_utm_flw[\"flwdir_min\"] = flw.d8_from_dem(\n", " da_elv = merit_elv_utm,\n", " gdf_stream = None,\n", - " max_depth=-1,\n", + " max_depth=-1, # no local pits\n", " outlets = \"min\",\n", " idxs_pit = None,\n", ")" @@ -474,6 +474,9 @@ "basins = flwdir.basins(idxs=flwdir.idxs_pit).astype(np.int32)\n", "merit_utm[\"basins\"] = xr.Variable(dims, basins, attrs=dict(_FillValue=0))\n", "\n", + "# basin index file\n", + "gdf_basins = merit_utm[\"basins\"].raster.vectorize()\n", + "\n", "merit_utm" ] }, @@ -515,6 +518,107 @@ "merit_utm[\"basins\"].plot(ax=axes[2,1])\n", "_ = add_legend_titles(axes[2,1], \"Basins ID\", add_legend=False)\n" ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Exporting the newly created data and corresponding data catalog" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Finally, once we are happy with the new dataset, we can write out the data and create the corresponding data catalog so that it can be re-used to build a new wflow model." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Export the gridded data as tif files in a new folder\n", + "output_path = \"./elevation_data\"\n", + "\n", + "# export the hydrography data as tif files (one per variable)\n", + "merit_utm.raster.to_mapstack(\n", + " root = os.path.join(output_path, \"merit_utm\"),\n", + " driver = \"GTiff\",\n", + ")\n", + "\n", + "# export the basin index as geosjon\n", + "gdf_basins.to_file(os.path.join(output_path, \"merit_utm_basins.geojson\"), driver=\"GeoJSON\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now let's prepare the corresponding data catalog: (the writefile command will directly write a file using the lines in the jupyter cell)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "%%writefile ./elevation_data/data_catalog.yml\n", + "merit_utm:\n", + " data_type: RasterDataset\n", + " driver: raster\n", + " crs: 3857\n", + " path: ./merit_utm/{variable}.tif\n", + " rename:\n", + " slope: lndslp\n", + " meta:\n", + " category: topography\n", + " paper_doi: 10.5194/hess-2020-582\n", + " paper_ref: Eilander et al. (2021)\n", + " source_license: ODC-By 1.0\n", + " source_url: https://zenodo.org/record/5166932#.YVbxJ5pByUk\n", + " source_doi: 10.5281/zenodo.5166932\n", + " source_version: 1.0\n", + " processing_notes: prepared from MERIT Hydro IHU by reprojecting the DEM to EPSG:3857 and deriving flow directions with river burning from lin2019_v1 rivers.\n", + " processing_script: prepare_ldd.ipynb from hydromt_wflow repository\n", + "\n", + "merit_utm_index:\n", + " data_type: GeoDataFrame\n", + " driver: vector\n", + " crs: 3857\n", + " path: ./merit_utm_basins.geojson\n", + " meta:\n", + " category: topography\n", + " paper_doi: 10.5194/hess-2020-582\n", + " paper_ref: Eilander et al. (2021)\n", + " source_license: ODC-By 1.0\n", + " source_url: https://zenodo.org/record/5166932#.YVbxJ5pByUk\n", + " source_doi: 10.5281/zenodo.5166932\n", + " source_version: 1.0\n", + " processing_notes: prepared from MERIT Hydro IHU by reprojecting the DEM to EPSG:3857 and deriving flow directions with river burning from lin2019_v1 rivers.\n", + " processing_script: prepare_ldd.ipynb from hydromt_wflow repository\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "And now let's try to load our data again with hydromt:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "data_catalog = DataCatalog(data_libs=\"./elevation_data/data_catalog.yml\")\n", + "\n", + "merit_utm = data_catalog.get_rasterdataset(\"merit_utm\", variables=[\"elevtn\"])\n", + "merit_utm" + ] } ], "metadata": { From c6f9250ba2c7251db3386a1a3df1c4c0bbd1097d Mon Sep 17 00:00:00 2001 From: hboisgon Date: Wed, 26 Jul 2023 15:53:51 +0800 Subject: [PATCH 7/9] simplify examples/prepare_ldd.ipynb --- examples/prepare_ldd.ipynb | 233 ++++++++++--------------------------- 1 file changed, 59 insertions(+), 174 deletions(-) diff --git a/examples/prepare_ldd.ipynb b/examples/prepare_ldd.ipynb index 6f1526f6..172c7f29 100644 --- a/examples/prepare_ldd.ipynb +++ b/examples/prepare_ldd.ipynb @@ -15,7 +15,7 @@ "source": [ "With HydroMT-Wflow, a user can choose to build a model in a geographic or projected coordinate system from an input Digital Elevation Model (DEM) and Flow Direction (flwdir) dataset.\n", "\n", - "While DEM data are often available, this is not the always the case for the flow directions. In the plugin, we made the choice to build a Wflow model directly from user provided DEM and flwdir datasets rather than reprojecting a DEM and/or deriving flwdir on the fly. This is because there are a lot of available techniques to derive flow directions and we want the user to be sure the flow directions matches the terrain and user expectations.\n", + "While DEM data are often available, this is not the always the case for the flow directions (flwdir). In the plugin, we made the choice to build a Wflow model directly from user provided DEM and flwdir datasets rather than reprojecting a DEM and/or deriving flwdir on the fly. This is because there are a lot of available techniques to derive flow directions and we want the user to be sure the flow directions matches the terrain and user expectations.\n", "\n", "Because of this, we prefer to provide this notebook as a possible pre-processing step before calling a hydromt build wflow command. We will do this using the different [flow directions methods from HydroMT](https://deltares.github.io/hydromt/latest/api.html#flow-direction-methods) and [PyFlwDir](https://deltares.github.io/pyflwdir/latest/index.html)" ] @@ -96,7 +96,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "### Elevation data and reprojection" + "### Deriving flow directions from Elevation data" ] }, { @@ -104,9 +104,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "In this example we will use the `merit_hydro_1k` data in the pre-defined `artifact_data` catalog of HydroMT. \n", - "\n", - "If you prefer to have a Wflow model in projected coordinate system (CRS in metres), you will need first to reproject the DEM and we will see here an example of how to do this.\n", + "In this example we will use the `merit_hydro_1k` data in the pre-defined `artifact_data` catalog of HydroMT.\n", "\n", "First let's read the data:" ] @@ -118,7 +116,7 @@ "outputs": [], "source": [ "data_catalog = DataCatalog(\"artifact_data\")\n", - "merit = data_catalog.get_rasterdataset(\"merit_hydro_1k\", bbox = [11.8,46.0,12.3,46.2])" + "merit = data_catalog.get_rasterdataset(\"merit_hydro_1k\", variables=[\"elevtn\", \"flwdir\"], bbox = [11.8,46.0,12.3,46.2])" ] }, { @@ -146,106 +144,6 @@ "legend = add_legend_titles(ax, \"MERIT Hydro IHU 1km\", projected=False)" ] }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "And now let's reproject the DEM to the projected CRS EPSG 3857 using [hydromt.DataArray.raster.reproject](https://deltares.github.io/hydromt/latest/_generated/hydromt.DataArray.raster.reproject.html#hydromt.DataArray.raster.reproject) method." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "merit_elv_utm = merit[\"elevtn\"].raster.reproject(\n", - " dst_crs=3857,\n", - " dst_res = None, # keep the same resolution\n", - " method = 'bilinear',\n", - ")" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Reprojecting flow directions" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "If your original data did contain flow directions, you can decide to reproject rather than re-create them from the reprojected DEM.\n", - "\n", - "To do this, you can use the [hydromt.flw.reproject_hydrography_like](https://deltares.github.io/hydromt/latest/_generated/hydromt.flw.reproject_hydrography_like.html#hydromt.flw.reproject_hydrography_like) method from HydroMT.\n", - "\n", - "This method aims to conserve flow direction during reprojection by deriving flwdir from a reprojected grid of synthetic elevation, based on the log10 of upstream area [m2]." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "merit_utm_reproj = flw.reproject_hydrography_like(\n", - " ds_hydro = merit[[\"flwdir\", \"uparea\"]],\n", - " da_elv = merit_elv_utm,\n", - ")" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Let's compare our new reprojected flow directions to the original ones" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# plot\n", - "# initialize image with geoaxes\n", - "fig = plt.figure(figsize=(10,8))\n", - "ax = fig.add_subplot()\n", - "\n", - "## plot elevation\\\n", - "merit_elv_utm.plot(\n", - " ax=ax, zorder=1, cbar_kwargs=dict(aspect=30, shrink=0.8), **kwargs\n", - ")\n", - "\n", - "# plot flwdir\n", - "flwdir = flw.flwdir_from_da(merit[\"flwdir\"], ftype=\"infer\", check_ftype=True)\n", - "feat = flwdir.streams(min_sto=2)\n", - "gdf = gpd.GeoDataFrame.from_features(feat, crs=merit.raster.crs)\n", - "gdf.to_crs(merit_elv_utm.raster.crs).plot(ax=ax, column=\"strord\", cmap=colors.ListedColormap(cm.Blues(np.linspace(0.4, 1, 7))), label=\"Flow directions\")\n", - "\n", - "# plot flwdir reproj\n", - "flwdirr = flw.flwdir_from_da(merit_utm_reproj[\"flwdir\"], ftype=\"infer\", check_ftype=True)\n", - "featr = flwdirr.streams(min_sto=2)\n", - "gdfr = gpd.GeoDataFrame.from_features(featr, crs=merit_utm_reproj.raster.crs)\n", - "gdfr.plot(ax=ax, column=\"strord\", cmap=colors.ListedColormap(cm.Reds(np.linspace(0.4, 1, 7))), label=\"Flow directions reprojected\")\n", - "\n", - "legend = add_legend_titles(ax, \"MERIT Hydro IHU 1km Reprojected\")" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Deriving flow directions from DEM" - ] - }, { "attachments": {}, "cell_type": "markdown", @@ -267,19 +165,17 @@ "metadata": {}, "outputs": [], "source": [ - "# Start a merit_utm dataset\n", - "merit_utm_flw = merit_elv_utm.to_dataset(name=\"elevtn\")\n", "# Derive flow directions with outlets at the edges\n", - "merit_utm_flw[\"flwdir_edge\"] = flw.d8_from_dem(\n", - " da_elv = merit_elv_utm,\n", + "merit[\"flwdir_edge\"] = flw.d8_from_dem(\n", + " da_elv = merit[\"elevtn\"],\n", " gdf_stream = None,\n", " max_depth=-1, # no local pits\n", " outlets = \"edge\",\n", " idxs_pit = None,\n", ")\n", "# Derive flow directions with outlet at the min elevation edge cell\n", - "merit_utm_flw[\"flwdir_min\"] = flw.d8_from_dem(\n", - " da_elv = merit_elv_utm,\n", + "merit[\"flwdir_min\"] = flw.d8_from_dem(\n", + " da_elv = merit[\"elevtn\"],\n", " gdf_stream = None,\n", " max_depth=-1, # no local pits\n", " outlets = \"min\",\n", @@ -307,7 +203,7 @@ "ax = fig.add_subplot()\n", "\n", "## plot elevation\\\n", - "merit_elv_utm.plot(\n", + "merit[\"elevtn\"].plot(\n", " ax=ax, zorder=1, cbar_kwargs=dict(aspect=30, shrink=0.8), alpha=0.5, **kwargs\n", ")\n", "\n", @@ -315,24 +211,18 @@ "flwdir = flw.flwdir_from_da(merit[\"flwdir\"], ftype=\"infer\", check_ftype=True)\n", "feat = flwdir.streams(min_sto=2)\n", "gdf = gpd.GeoDataFrame.from_features(feat, crs=merit.raster.crs)\n", - "gdf.to_crs(merit_elv_utm.raster.crs).plot(ax=ax, column=\"strord\", cmap=colors.ListedColormap(cm.Blues(np.linspace(0.4, 1, 7))), label=\"Flow directions\")\n", - "\n", - "# plot flwdir reproj\n", - "flwdirr = flw.flwdir_from_da(merit_utm_reproj[\"flwdir\"], ftype=\"infer\", check_ftype=True)\n", - "featr = flwdirr.streams(min_sto=2)\n", - "gdfr = gpd.GeoDataFrame.from_features(featr, crs=merit_utm_reproj.raster.crs)\n", - "gdfr.plot(ax=ax, column=\"strord\", cmap=colors.ListedColormap(cm.Reds(np.linspace(0.4, 1, 7))), label=\"Flow directions reprojected\")\n", + "gdf.to_crs(merit.raster.crs).plot(ax=ax, column=\"strord\", cmap=colors.ListedColormap(cm.Blues(np.linspace(0.4, 1, 7))), label=\"Flow directions\")\n", "\n", "# plot flwdir edge\n", - "flwdir_edge = flw.flwdir_from_da(merit_utm_flw[\"flwdir_edge\"], ftype=\"infer\", check_ftype=True)\n", + "flwdir_edge = flw.flwdir_from_da(merit[\"flwdir_edge\"], ftype=\"infer\", check_ftype=True)\n", "feate = flwdir_edge.streams(min_sto=2)\n", - "gdfe = gpd.GeoDataFrame.from_features(feate, crs=merit_utm_flw.raster.crs)\n", + "gdfe = gpd.GeoDataFrame.from_features(feate, crs=merit.raster.crs)\n", "gdfe.plot(ax=ax, column=\"strord\", cmap=colors.ListedColormap(cm.Greens(np.linspace(0.4, 1, 7))), label=\"Flow directions edge\")\n", "\n", "# plot flwdir min\n", - "flwdir_min = flw.flwdir_from_da(merit_utm_flw[\"flwdir_min\"], ftype=\"infer\", check_ftype=True)\n", + "flwdir_min = flw.flwdir_from_da(merit[\"flwdir_min\"], ftype=\"infer\", check_ftype=True)\n", "featm = flwdir_min.streams(min_sto=2)\n", - "gdfm = gpd.GeoDataFrame.from_features(featm, crs=merit_utm_flw.raster.crs)\n", + "gdfm = gpd.GeoDataFrame.from_features(featm, crs=merit.raster.crs)\n", "gdfm.plot(ax=ax, column=\"strord\", cmap=colors.ListedColormap(cm.Purples(np.linspace(0.4, 1, 7))), label=\"Flow directions min\")\n", "\n", "legend = add_legend_titles(ax, \"MERIT Hydro IHU 1km Reprojected\")" @@ -361,12 +251,13 @@ "# In this dataset, there is an Area column representing upstream area in km2\n", "# We need to rename it to uparea\n", "rivers = rivers.rename(columns={\"Area\": \"uparea\"})\n", - "# And finally make sure rivers is reprojected to the same CRS as the elevation\n", - "rivers = rivers.to_crs(merit_elv_utm.raster.crs)\n", + "# And finally make sure rivers is reprojected to the same CRS as the elevation (if not already the case)\n", + "if rivers.crs != merit.raster.crs:\n", + " rivers = rivers.to_crs(merit.raster.crs)\n", "\n", "# Now let's use it to derive flow directions\n", - "merit_utm_flw[\"flwdir_riverburn\"] = flw.d8_from_dem(\n", - " da_elv = merit_elv_utm,\n", + "merit[\"flwdir_riverburn\"] = flw.d8_from_dem(\n", + " da_elv = merit[\"elevtn\"],\n", " gdf_stream = rivers,\n", " max_depth=-1,\n", " outlets = \"edge\",\n", @@ -386,7 +277,7 @@ "ax = fig.add_subplot()\n", "\n", "## plot elevation\\\n", - "merit_elv_utm.plot(\n", + "merit[\"elevtn\"].plot(\n", " ax=ax, zorder=1, cbar_kwargs=dict(aspect=30, shrink=0.8), alpha=0.5, **kwargs\n", ")\n", "\n", @@ -394,24 +285,18 @@ "flwdir = flw.flwdir_from_da(merit[\"flwdir\"], ftype=\"infer\", check_ftype=True)\n", "feat = flwdir.streams(min_sto=2)\n", "gdf = gpd.GeoDataFrame.from_features(feat, crs=merit.raster.crs)\n", - "gdf.to_crs(merit_elv_utm.raster.crs).plot(ax=ax, column=\"strord\", cmap=colors.ListedColormap(cm.Blues(np.linspace(0.4, 1, 7))), label=\"Flow directions\")\n", - "\n", - "# plot flwdir reproj\n", - "flwdirr = flw.flwdir_from_da(merit_utm_reproj[\"flwdir\"], ftype=\"infer\", check_ftype=True)\n", - "featr = flwdirr.streams(min_sto=2)\n", - "gdfr = gpd.GeoDataFrame.from_features(featr, crs=merit_utm_reproj.raster.crs)\n", - "gdfr.plot(ax=ax, column=\"strord\", cmap=colors.ListedColormap(cm.Reds(np.linspace(0.4, 1, 7))), label=\"Flow directions reprojected\")\n", + "gdf.to_crs(merit.raster.crs).plot(ax=ax, column=\"strord\", cmap=colors.ListedColormap(cm.Blues(np.linspace(0.4, 1, 7))), label=\"Flow directions\")\n", "\n", "# plot flwdir edge\n", - "flwdir_edge = flw.flwdir_from_da(merit_utm_flw[\"flwdir_edge\"], ftype=\"infer\", check_ftype=True)\n", + "flwdir_edge = flw.flwdir_from_da(merit[\"flwdir_edge\"], ftype=\"infer\", check_ftype=True)\n", "feate = flwdir_edge.streams(min_sto=2)\n", - "gdfe = gpd.GeoDataFrame.from_features(feate, crs=merit_utm_flw.raster.crs)\n", + "gdfe = gpd.GeoDataFrame.from_features(feate, crs=merit.raster.crs)\n", "gdfe.plot(ax=ax, column=\"strord\", cmap=colors.ListedColormap(cm.Greens(np.linspace(0.4, 1, 7))), label=\"Flow directions edge\")\n", "\n", "# plot flwdir riverburn\n", - "flwdir_min = flw.flwdir_from_da(merit_utm_flw[\"flwdir_riverburn\"], ftype=\"infer\", check_ftype=True)\n", + "flwdir_min = flw.flwdir_from_da(merit[\"flwdir_riverburn\"], ftype=\"infer\", check_ftype=True)\n", "featm = flwdir_min.streams(min_sto=2)\n", - "gdfm = gpd.GeoDataFrame.from_features(featm, crs=merit_utm_flw.raster.crs)\n", + "gdfm = gpd.GeoDataFrame.from_features(featm, crs=merit.raster.crs)\n", "gdfm.plot(ax=ax, column=\"strord\", cmap=colors.ListedColormap(cm.Purples(np.linspace(0.4, 1, 7))), label=\"Flow directions river burning\")\n", "\n", "rivers.plot(ax=ax, color=\"black\", label=\"Rivers Lin\")\n", @@ -443,41 +328,41 @@ "metadata": {}, "outputs": [], "source": [ - "# Create a dataset with the riverburn flow directions\n", - "merit_utm = merit_elv_utm.to_dataset(name=\"elevtn\")\n", - "merit_utm[\"flwdir\"] = merit_utm_flw[\"flwdir_riverburn\"]\n", - "dims = merit_utm.raster.dims\n", + "# Create a new merit_adapted dataset with the riverburn flow directions\n", + "merit_adapted = merit[\"elevtn\"].to_dataset(name=\"elevtn\")\n", + "merit_adapted[\"flwdir\"] = merit[\"flwdir_riverburn\"]\n", + "dims = merit_adapted.raster.dims\n", "\n", "# Create a PyFlwDir object from the dataset\n", - "flwdir = flw.flwdir_from_da(merit_utm[\"flwdir\"])\n", + "flwdir = flw.flwdir_from_da(merit_adapted[\"flwdir\"])\n", "\n", "# uparea\n", "uparea = flwdir.upstream_area(unit=\"km2\")\n", - "merit_utm[\"uparea\"] = xr.Variable(dims, uparea, attrs=dict(_FillValue=-9999))\n", + "merit_adapted[\"uparea\"] = xr.Variable(dims, uparea, attrs=dict(_FillValue=-9999))\n", "\n", "# stream order\n", "strord = flwdir.stream_order()\n", - "merit_utm[\"strord\"] = xr.Variable(dims, strord)\n", - "merit_utm[\"strord\"].raster.set_nodata(255)\n", + "merit_adapted[\"strord\"] = xr.Variable(dims, strord)\n", + "merit_adapted[\"strord\"].raster.set_nodata(255)\n", "\n", "# slope\n", "slope = pyflwdir.dem.slope(\n", - " elevtn=merit_utm[\"elevtn\"].values,\n", - " nodata = merit_utm[\"elevtn\"].raster.nodata,\n", + " elevtn = merit_adapted[\"elevtn\"].values,\n", + " nodata = merit_adapted[\"elevtn\"].raster.nodata,\n", " latlon = False, # True if geographic crs, False if projected crs\n", - " transform = merit_utm[\"elevtn\"].raster.transform,\n", + " transform = merit_adapted[\"elevtn\"].raster.transform,\n", ")\n", - "merit_utm[\"slope\"] = xr.Variable(dims, slope)\n", - "merit_utm[\"slope\"].raster.set_nodata(merit_utm[\"elevtn\"].raster.nodata)\n", + "merit_adapted[\"slope\"] = xr.Variable(dims, slope)\n", + "merit_adapted[\"slope\"].raster.set_nodata(merit_adapted[\"elevtn\"].raster.nodata)\n", "\n", "# basin at the pits locations\n", "basins = flwdir.basins(idxs=flwdir.idxs_pit).astype(np.int32)\n", - "merit_utm[\"basins\"] = xr.Variable(dims, basins, attrs=dict(_FillValue=0))\n", + "merit_adapted[\"basins\"] = xr.Variable(dims, basins, attrs=dict(_FillValue=0))\n", "\n", "# basin index file\n", - "gdf_basins = merit_utm[\"basins\"].raster.vectorize()\n", + "gdf_basins = merit_adapted[\"basins\"].raster.vectorize()\n", "\n", - "merit_utm" + "merit_adapted" ] }, { @@ -490,32 +375,32 @@ "fig, axes = plt.subplots(3, 2, figsize=(15, 15))\n", "\n", "## plot elevation\n", - "merit_utm[\"elevtn\"].plot(\n", + "merit_adapted[\"elevtn\"].plot(\n", " ax=axes[0,0], zorder=1, cbar_kwargs=dict(aspect=30, shrink=0.8), alpha=0.5, **kwargs\n", ")\n", "_ = add_legend_titles(axes[0,0], \"Elevation [m asl]\", add_legend=False)\n", "\n", "# plot flwdir riverburn\n", - "flwdir = flw.flwdir_from_da(merit_utm[\"flwdir\"], ftype=\"infer\", check_ftype=True)\n", + "flwdir = flw.flwdir_from_da(merit_adapted[\"flwdir\"], ftype=\"infer\", check_ftype=True)\n", "feat = flwdir.streams(min_sto=2)\n", - "gdf = gpd.GeoDataFrame.from_features(feat, crs=merit_utm.raster.crs)\n", - "gdf.to_crs(merit_elv_utm.raster.crs).plot(ax=axes[0,1], column=\"strord\", cmap=colors.ListedColormap(cm.Blues(np.linspace(0.4, 1, 7))), label=\"Flow directions\")\n", + "gdf = gpd.GeoDataFrame.from_features(feat, crs=merit_adapted.raster.crs)\n", + "gdf.to_crs(merit_adapted.raster.crs).plot(ax=axes[0,1], column=\"strord\", cmap=colors.ListedColormap(cm.Blues(np.linspace(0.4, 1, 7))), label=\"Flow directions\")\n", "_ = add_legend_titles(axes[0,1], \"Flow directions\", add_legend=False)\n", "\n", "# plot uparea\n", - "merit_utm[\"uparea\"].plot(ax=axes[1,0])\n", + "merit_adapted[\"uparea\"].plot(ax=axes[1,0])\n", "_ = add_legend_titles(axes[1,0], \"Upstream area [km2]\", add_legend=False)\n", "\n", "# plot strord\n", - "merit_utm[\"strord\"].plot(ax=axes[1,1], cmap=colors.ListedColormap(cm.Blues(np.linspace(0.4, 1, 7))))\n", + "merit_adapted[\"strord\"].plot(ax=axes[1,1], cmap=colors.ListedColormap(cm.Blues(np.linspace(0.4, 1, 7))))\n", "_ = add_legend_titles(axes[1,1], \"Strahler Stream order\", add_legend=False)\n", "\n", "# plot slope\n", - "merit_utm[\"slope\"].plot(ax=axes[2,0])\n", + "merit_adapted[\"slope\"].plot(ax=axes[2,0])\n", "_ = add_legend_titles(axes[2,0], \"Slope [m/m]\", add_legend=False)\n", "\n", "# plot basins\n", - "merit_utm[\"basins\"].plot(ax=axes[2,1])\n", + "merit_adapted[\"basins\"].plot(ax=axes[2,1])\n", "_ = add_legend_titles(axes[2,1], \"Basins ID\", add_legend=False)\n" ] }, @@ -543,13 +428,13 @@ "output_path = \"./elevation_data\"\n", "\n", "# export the hydrography data as tif files (one per variable)\n", - "merit_utm.raster.to_mapstack(\n", - " root = os.path.join(output_path, \"merit_utm\"),\n", + "merit_adapted.raster.to_mapstack(\n", + " root = os.path.join(output_path, \"merit_adapted\"),\n", " driver = \"GTiff\",\n", ")\n", "\n", "# export the basin index as geosjon\n", - "gdf_basins.to_file(os.path.join(output_path, \"merit_utm_basins.geojson\"), driver=\"GeoJSON\")" + "gdf_basins.to_file(os.path.join(output_path, \"merit_adapted_basins.geojson\"), driver=\"GeoJSON\")" ] }, { @@ -566,11 +451,11 @@ "outputs": [], "source": [ "%%writefile ./elevation_data/data_catalog.yml\n", - "merit_utm:\n", + "merit_adapted:\n", " data_type: RasterDataset\n", " driver: raster\n", " crs: 3857\n", - " path: ./merit_utm/{variable}.tif\n", + " path: ./merit_adapted/{variable}.tif\n", " rename:\n", " slope: lndslp\n", " meta:\n", @@ -581,14 +466,14 @@ " source_url: https://zenodo.org/record/5166932#.YVbxJ5pByUk\n", " source_doi: 10.5281/zenodo.5166932\n", " source_version: 1.0\n", - " processing_notes: prepared from MERIT Hydro IHU by reprojecting the DEM to EPSG:3857 and deriving flow directions with river burning from lin2019_v1 rivers.\n", + " processing_notes: prepared from MERIT Hydro IHU by deriving flow directions with river burning from lin2019_v1 rivers using pyflwdir.\n", " processing_script: prepare_ldd.ipynb from hydromt_wflow repository\n", "\n", - "merit_utm_index:\n", + "merit_adapted_index:\n", " data_type: GeoDataFrame\n", " driver: vector\n", " crs: 3857\n", - " path: ./merit_utm_basins.geojson\n", + " path: ./merit_adapted_basins.geojson\n", " meta:\n", " category: topography\n", " paper_doi: 10.5194/hess-2020-582\n", @@ -597,7 +482,7 @@ " source_url: https://zenodo.org/record/5166932#.YVbxJ5pByUk\n", " source_doi: 10.5281/zenodo.5166932\n", " source_version: 1.0\n", - " processing_notes: prepared from MERIT Hydro IHU by reprojecting the DEM to EPSG:3857 and deriving flow directions with river burning from lin2019_v1 rivers.\n", + " processing_notes: prepared from MERIT Hydro IHU by deriving flow directions with river burning from lin2019_v1 rivers using pyflwdir.\n", " processing_script: prepare_ldd.ipynb from hydromt_wflow repository\n" ] }, @@ -616,7 +501,7 @@ "source": [ "data_catalog = DataCatalog(data_libs=\"./elevation_data/data_catalog.yml\")\n", "\n", - "merit_utm = data_catalog.get_rasterdataset(\"merit_utm\", variables=[\"elevtn\"])\n", + "merit_utm = data_catalog.get_rasterdataset(\"merit_adapted\")\n", "merit_utm" ] } From 8bf1d72e3f789c649b49e5559cbd7576c482d2f7 Mon Sep 17 00:00:00 2001 From: hboisgon Date: Thu, 27 Jul 2023 09:08:18 +0800 Subject: [PATCH 8/9] run black after merge --- tests/test_model_methods.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/tests/test_model_methods.py b/tests/test_model_methods.py index e84dd1af..454de934 100644 --- a/tests/test_model_methods.py +++ b/tests/test_model_methods.py @@ -52,6 +52,7 @@ def test_setup_staticmaps(example_wflow_model): wflow_variables=["input.vertical.altitude"], ) + def test_projected_crs(tmpdir): logger = logging.getLogger(__name__) @@ -100,6 +101,7 @@ def test_projected_crs(tmpdir): assert np.quantile(mod.staticmaps["wflow_landuse"], 0.95) == 190.0 # urban assert mod.get_config("model.sizeinmetres") == True + def test_setup_lake(tmpdir, example_wflow_model): # Create dummy lake rating curves lakes = example_wflow_model.staticgeoms["lakes"] From 5563acb8c97e0e3415da1f1d87ec7ac3827dc93d Mon Sep 17 00:00:00 2001 From: Dirk Eilander Date: Thu, 27 Jul 2023 11:48:04 +0200 Subject: [PATCH 9/9] use merit hydro and remove min (not applicable) and rivers (buggy) options --- .gitignore | 3 + examples/prepare_ldd.ipynb | 230 ++++++++++--------------------------- 2 files changed, 65 insertions(+), 168 deletions(-) diff --git a/.gitignore b/.gitignore index ad6854ad..6218cf3d 100644 --- a/.gitignore +++ b/.gitignore @@ -1,5 +1,8 @@ # file based on github/gitignore +# ignore all files in this directory +sandbox + # Byte-compiled / optimized / DLL files __pycache__/ *.py[cod] diff --git a/examples/prepare_ldd.ipynb b/examples/prepare_ldd.ipynb index 172c7f29..9682f4fb 100644 --- a/examples/prepare_ldd.ipynb +++ b/examples/prepare_ldd.ipynb @@ -15,9 +15,9 @@ "source": [ "With HydroMT-Wflow, a user can choose to build a model in a geographic or projected coordinate system from an input Digital Elevation Model (DEM) and Flow Direction (flwdir) dataset.\n", "\n", - "While DEM data are often available, this is not the always the case for the flow directions (flwdir). In the plugin, we made the choice to build a Wflow model directly from user provided DEM and flwdir datasets rather than reprojecting a DEM and/or deriving flwdir on the fly. This is because there are a lot of available techniques to derive flow directions and we want the user to be sure the flow directions matches the terrain and user expectations.\n", + "While DEM data are often available, this is not the always the case for the flow directions (flwdir). We made the choice to build a Wflow model directly from user provided DEM and flwdir datasets rather than reprojecting a DEM and/or deriving flwdir on the fly. This is because deriving flow directions is often an iterative process to be sure the flow directions matche the terrain and river network. Note that for the best results the flwdir data should be derived from a high-res DEM (<100 m spatial resolution). The HydroMT-Wflow model builder will automatically resample the flwdir data to the model resolution.\n", "\n", - "Because of this, we prefer to provide this notebook as a possible pre-processing step before calling a hydromt build wflow command. We will do this using the different [flow directions methods from HydroMT](https://deltares.github.io/hydromt/latest/api.html#flow-direction-methods) and [PyFlwDir](https://deltares.github.io/pyflwdir/latest/index.html)" + "Because of this, we prefer to provide this notebook as a possible **pre-processing step** before calling a build a Wflow model with HydroMt-Wflow. Here we use the different [flow directions methods from HydroMT](https://deltares.github.io/hydromt/latest/api.html#flow-direction-methods) and [PyFlwDir](https://deltares.github.io/pyflwdir/latest/index.html)" ] }, { @@ -38,11 +38,14 @@ "import xarray as xr\n", "import numpy as np\n", "import geopandas as gpd\n", + "\n", "# pyflwdir\n", "import pyflwdir\n", + "\n", "# hydromt\n", "from hydromt import DataCatalog, flw\n", - "#plot\n", + "\n", + "# plot\n", "import matplotlib.pyplot as plt\n", "from matplotlib import cm, colors\n", "import cartopy.crs as ccrs" @@ -75,6 +78,7 @@ " facecolor=\"white\",\n", ")\n", "\n", + "\n", "def add_legend_titles(ax, title, add_legend=True, projected=True):\n", " ax.xaxis.set_visible(True)\n", " ax.yaxis.set_visible(True)\n", @@ -104,7 +108,8 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "In this example we will use the `merit_hydro_1k` data in the pre-defined `artifact_data` catalog of HydroMT.\n", + "In this example we will use the `merit_hydro` data in the pre-defined `artifact_data` catalog of HydroMT. \n", + "Note that this dataset already provides flow direction, but for sake of example we only use this as reference dataset.\n", "\n", "First let's read the data:" ] @@ -116,32 +121,9 @@ "outputs": [], "source": [ "data_catalog = DataCatalog(\"artifact_data\")\n", - "merit = data_catalog.get_rasterdataset(\"merit_hydro_1k\", variables=[\"elevtn\", \"flwdir\"], bbox = [11.8,46.0,12.3,46.2])" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# plot\n", - "# initialize image with geoaxes\n", - "fig = plt.figure(figsize=(10,8))\n", - "ax = fig.add_subplot(projection=proj)\n", - "\n", - "## plot elevation\\\n", - "merit['elevtn'].plot(\n", - " transform=proj, ax=ax, zorder=1, cbar_kwargs=dict(aspect=30, shrink=0.8), **kwargs\n", - ")\n", - "\n", - "# plot flwdir\n", - "flwdir = flw.flwdir_from_da(merit[\"flwdir\"], ftype=\"infer\", check_ftype=True)\n", - "feat = flwdir.streams()\n", - "gdf = gpd.GeoDataFrame.from_features(feat, crs=merit.raster.crs)\n", - "gdf.plot(ax=ax, color=\"blue\", linewidth=0.5, zorder=2, label=\"Flow directions\")\n", - "\n", - "legend = add_legend_titles(ax, \"MERIT Hydro IHU 1km\", projected=False)" + "merit = data_catalog.get_rasterdataset(\n", + " \"merit_hydro\", variables=[\"elevtn\", \"flwdir\"], bbox=[12, 46.0, 12.3, 46.2]\n", + ")" ] }, { @@ -152,9 +134,9 @@ "To derive flow directions from a DEM, you can use the [hydromt.flw.d8_from_dem](https://deltares.github.io/hydromt/latest/_generated/hydromt.flw.d8_from_dem.html#hydromt.flw.d8_from_dem) method of HydroMT.\n", "\n", "This method derives D8 flow directions grid from an elevation grid and allows several options to the users:\n", - " - **outlets**: outlets can be defined at ``edge``s of the grid or force all flow to go to the minimum elevation point ``min``. Additionnally, the user can also specify the pits locations via ``idxs_pit``.\n", - " - **river burning**: it is possible to provide a river vector layer ``gdf_stream`` with ``uparea`` (km2) column which is used to burn the river in the elevation data.\n", - " - **depression filling**: local depressions are filled based on their lowest pour point level if the pour point depth is smaller than the maximum pour point depth ``max_depth``, otherwise the lowest elevation in the depression becomes a pit.\n", + " - **outlets**: outlets can be defined at ``edge``s of the grid (defualt) or force all flow to go to the minimum elevation point ``min``. The latter only makes sense if your DEM only is masked to the catchment. Additionnally, the user can also force specific pits locations via ``idxs_pit``.\n", + " - **depression filling**: local depressions are filled based on their lowest pour point level if the pour point depth is smaller than the maximum pour point depth ``max_depth``, otherwise the lowest elevation in the depression becomes a pit. By default ``max_depth`` is set to -1 m filling all local depressions.\n", + "- **river burning**: while it is already possible to provide a river vector layer ``gdf_stream`` with ``uparea`` (km2) column to further guide the derivation of flow directions, this option is currently being improved and not yet fully functional. See also HydroMT core PR [305](https://github.com/Deltares/hydromt/pull/305)\n", "\n", "Let's see an example:" ] @@ -165,21 +147,13 @@ "metadata": {}, "outputs": [], "source": [ - "# Derive flow directions with outlets at the edges\n", - "merit[\"flwdir_edge\"] = flw.d8_from_dem(\n", - " da_elv = merit[\"elevtn\"],\n", - " gdf_stream = None,\n", - " max_depth=-1, # no local pits\n", - " outlets = \"edge\",\n", - " idxs_pit = None,\n", - ")\n", - "# Derive flow directions with outlet at the min elevation edge cell\n", - "merit[\"flwdir_min\"] = flw.d8_from_dem(\n", - " da_elv = merit[\"elevtn\"],\n", - " gdf_stream = None,\n", - " max_depth=-1, # no local pits\n", - " outlets = \"min\",\n", - " idxs_pit = None,\n", + "# Derive flow directions with outlets at the edges -> this is the default\n", + "merit[\"flwdir_derived\"] = flw.d8_from_dem(\n", + " da_elv=merit[\"elevtn\"],\n", + " gdf_stream=None,\n", + " max_depth=-1, # no local pits\n", + " outlets=\"edge\",\n", + " idxs_pit=None,\n", ")" ] }, @@ -188,7 +162,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Let's plot all the different methods" + "Let's compare the orignal with the derived flow direction data" ] }, { @@ -199,109 +173,38 @@ "source": [ "# plot\n", "# initialize image with geoaxes\n", - "fig = plt.figure(figsize=(10,8))\n", + "fig = plt.figure(figsize=(8, 8))\n", "ax = fig.add_subplot()\n", "\n", "## plot elevation\\\n", "merit[\"elevtn\"].plot(\n", - " ax=ax, zorder=1, cbar_kwargs=dict(aspect=30, shrink=0.8), alpha=0.5, **kwargs\n", + " ax=ax, zorder=1, cbar_kwargs=dict(aspect=30, shrink=0.5), alpha=0.5, **kwargs\n", ")\n", "\n", "# plot flwdir\n", "flwdir = flw.flwdir_from_da(merit[\"flwdir\"], ftype=\"infer\", check_ftype=True)\n", - "feat = flwdir.streams(min_sto=2)\n", + "feat = flwdir.streams(min_sto=3) # minimum stream order for plotting\n", "gdf = gpd.GeoDataFrame.from_features(feat, crs=merit.raster.crs)\n", - "gdf.to_crs(merit.raster.crs).plot(ax=ax, column=\"strord\", cmap=colors.ListedColormap(cm.Blues(np.linspace(0.4, 1, 7))), label=\"Flow directions\")\n", - "\n", - "# plot flwdir edge\n", - "flwdir_edge = flw.flwdir_from_da(merit[\"flwdir_edge\"], ftype=\"infer\", check_ftype=True)\n", - "feate = flwdir_edge.streams(min_sto=2)\n", - "gdfe = gpd.GeoDataFrame.from_features(feate, crs=merit.raster.crs)\n", - "gdfe.plot(ax=ax, column=\"strord\", cmap=colors.ListedColormap(cm.Greens(np.linspace(0.4, 1, 7))), label=\"Flow directions edge\")\n", - "\n", - "# plot flwdir min\n", - "flwdir_min = flw.flwdir_from_da(merit[\"flwdir_min\"], ftype=\"infer\", check_ftype=True)\n", - "featm = flwdir_min.streams(min_sto=2)\n", - "gdfm = gpd.GeoDataFrame.from_features(featm, crs=merit.raster.crs)\n", - "gdfm.plot(ax=ax, column=\"strord\", cmap=colors.ListedColormap(cm.Purples(np.linspace(0.4, 1, 7))), label=\"Flow directions min\")\n", - "\n", - "legend = add_legend_titles(ax, \"MERIT Hydro IHU 1km Reprojected\")" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Now let's see an example where we want to burn rivers.\n", - "In the data artifact, we have one river vector database : ``rivers_lin2019_v1``.\n", - "\n", - "Let's use it for burning." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Read the data from data catalog\n", - "rivers = data_catalog.get_geodataframe(\"rivers_lin2019_v1\", bbox=[11.8,46.0,12.3,46.2])\n", - "\n", - "# In this dataset, there is an Area column representing upstream area in km2\n", - "# We need to rename it to uparea\n", - "rivers = rivers.rename(columns={\"Area\": \"uparea\"})\n", - "# And finally make sure rivers is reprojected to the same CRS as the elevation (if not already the case)\n", - "if rivers.crs != merit.raster.crs:\n", - " rivers = rivers.to_crs(merit.raster.crs)\n", - "\n", - "# Now let's use it to derive flow directions\n", - "merit[\"flwdir_riverburn\"] = flw.d8_from_dem(\n", - " da_elv = merit[\"elevtn\"],\n", - " gdf_stream = rivers,\n", - " max_depth=-1,\n", - " outlets = \"edge\",\n", - " idxs_pit = None,\n", - ")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# plot\n", - "# initialize image with geoaxes\n", - "fig = plt.figure(figsize=(10,8))\n", - "ax = fig.add_subplot()\n", - "\n", - "## plot elevation\\\n", - "merit[\"elevtn\"].plot(\n", - " ax=ax, zorder=1, cbar_kwargs=dict(aspect=30, shrink=0.8), alpha=0.5, **kwargs\n", + "gdf.to_crs(merit.raster.crs).plot(\n", + " ax=ax, color=\"blue\", linewidth=gdf[\"strord\"] / 3, label=\"Original flow directions\"\n", ")\n", "\n", - "# plot flwdir\n", - "flwdir = flw.flwdir_from_da(merit[\"flwdir\"], ftype=\"infer\", check_ftype=True)\n", - "feat = flwdir.streams(min_sto=2)\n", - "gdf = gpd.GeoDataFrame.from_features(feat, crs=merit.raster.crs)\n", - "gdf.to_crs(merit.raster.crs).plot(ax=ax, column=\"strord\", cmap=colors.ListedColormap(cm.Blues(np.linspace(0.4, 1, 7))), label=\"Flow directions\")\n", "\n", "# plot flwdir edge\n", - "flwdir_edge = flw.flwdir_from_da(merit[\"flwdir_edge\"], ftype=\"infer\", check_ftype=True)\n", - "feate = flwdir_edge.streams(min_sto=2)\n", + "flwdir_edge = flw.flwdir_from_da(\n", + " merit[\"flwdir_derived\"], ftype=\"infer\", check_ftype=True\n", + ")\n", + "feate = flwdir_edge.streams(min_sto=3) # minimum stream order for plotting\n", "gdfe = gpd.GeoDataFrame.from_features(feate, crs=merit.raster.crs)\n", - "gdfe.plot(ax=ax, column=\"strord\", cmap=colors.ListedColormap(cm.Greens(np.linspace(0.4, 1, 7))), label=\"Flow directions edge\")\n", - "\n", - "# plot flwdir riverburn\n", - "flwdir_min = flw.flwdir_from_da(merit[\"flwdir_riverburn\"], ftype=\"infer\", check_ftype=True)\n", - "featm = flwdir_min.streams(min_sto=2)\n", - "gdfm = gpd.GeoDataFrame.from_features(featm, crs=merit.raster.crs)\n", - "gdfm.plot(ax=ax, column=\"strord\", cmap=colors.ListedColormap(cm.Purples(np.linspace(0.4, 1, 7))), label=\"Flow directions river burning\")\n", - "\n", - "rivers.plot(ax=ax, color=\"black\", label=\"Rivers Lin\")\n", + "gdfe.plot(\n", + " ax=ax,\n", + " column=\"strord\",\n", + " color=\"green\",\n", + " linewidth=gdfe[\"strord\"] / 3,\n", + " label=\"Derived flow directions\",\n", + ")\n", "\n", - "legend = add_legend_titles(ax, \"MERIT Hydro IHU 1km Reprojected\")" + "legend = add_legend_titles(ax, \"MERIT Hydro derived vs. original flow directions\")" ] }, { @@ -330,7 +233,7 @@ "source": [ "# Create a new merit_adapted dataset with the riverburn flow directions\n", "merit_adapted = merit[\"elevtn\"].to_dataset(name=\"elevtn\")\n", - "merit_adapted[\"flwdir\"] = merit[\"flwdir_riverburn\"]\n", + "merit_adapted[\"flwdir\"] = merit[\"flwdir_derived\"]\n", "dims = merit_adapted.raster.dims\n", "\n", "# Create a PyFlwDir object from the dataset\n", @@ -347,10 +250,10 @@ "\n", "# slope\n", "slope = pyflwdir.dem.slope(\n", - " elevtn = merit_adapted[\"elevtn\"].values,\n", - " nodata = merit_adapted[\"elevtn\"].raster.nodata,\n", - " latlon = False, # True if geographic crs, False if projected crs\n", - " transform = merit_adapted[\"elevtn\"].raster.transform,\n", + " elevtn=merit_adapted[\"elevtn\"].values,\n", + " nodata=merit_adapted[\"elevtn\"].raster.nodata,\n", + " latlon=False, # True if geographic crs, False if projected crs\n", + " transform=merit_adapted[\"elevtn\"].raster.transform,\n", ")\n", "merit_adapted[\"slope\"] = xr.Variable(dims, slope)\n", "merit_adapted[\"slope\"].raster.set_nodata(merit_adapted[\"elevtn\"].raster.nodata)\n", @@ -372,36 +275,25 @@ "outputs": [], "source": [ "# plot\n", - "fig, axes = plt.subplots(3, 2, figsize=(15, 15))\n", - "\n", - "## plot elevation\n", - "merit_adapted[\"elevtn\"].plot(\n", - " ax=axes[0,0], zorder=1, cbar_kwargs=dict(aspect=30, shrink=0.8), alpha=0.5, **kwargs\n", - ")\n", - "_ = add_legend_titles(axes[0,0], \"Elevation [m asl]\", add_legend=False)\n", - "\n", - "# plot flwdir riverburn\n", - "flwdir = flw.flwdir_from_da(merit_adapted[\"flwdir\"], ftype=\"infer\", check_ftype=True)\n", - "feat = flwdir.streams(min_sto=2)\n", - "gdf = gpd.GeoDataFrame.from_features(feat, crs=merit_adapted.raster.crs)\n", - "gdf.to_crs(merit_adapted.raster.crs).plot(ax=axes[0,1], column=\"strord\", cmap=colors.ListedColormap(cm.Blues(np.linspace(0.4, 1, 7))), label=\"Flow directions\")\n", - "_ = add_legend_titles(axes[0,1], \"Flow directions\", add_legend=False)\n", + "fig, axes = plt.subplots(2, 2, figsize=(12, 10))\n", "\n", "# plot uparea\n", - "merit_adapted[\"uparea\"].plot(ax=axes[1,0])\n", - "_ = add_legend_titles(axes[1,0], \"Upstream area [km2]\", add_legend=False)\n", + "merit_adapted[\"uparea\"].plot(ax=axes[0, 0])\n", + "_ = add_legend_titles(axes[0, 0], \"Upstream area [km2]\", add_legend=False)\n", "\n", "# plot strord\n", - "merit_adapted[\"strord\"].plot(ax=axes[1,1], cmap=colors.ListedColormap(cm.Blues(np.linspace(0.4, 1, 7))))\n", - "_ = add_legend_titles(axes[1,1], \"Strahler Stream order\", add_legend=False)\n", + "merit_adapted[\"strord\"].plot(\n", + " ax=axes[0, 1], cmap=colors.ListedColormap(cm.Blues(np.linspace(0.4, 1, 7)))\n", + ")\n", + "_ = add_legend_titles(axes[0, 1], \"Strahler Stream order\", add_legend=False)\n", "\n", "# plot slope\n", - "merit_adapted[\"slope\"].plot(ax=axes[2,0])\n", - "_ = add_legend_titles(axes[2,0], \"Slope [m/m]\", add_legend=False)\n", + "merit_adapted[\"slope\"].plot(ax=axes[1, 0])\n", + "_ = add_legend_titles(axes[1, 0], \"Slope [m/m]\", add_legend=False)\n", "\n", "# plot basins\n", - "merit_adapted[\"basins\"].plot(ax=axes[2,1])\n", - "_ = add_legend_titles(axes[2,1], \"Basins ID\", add_legend=False)\n" + "merit_adapted[\"basins\"].plot(ax=axes[1, 1])\n", + "_ = add_legend_titles(axes[1, 1], \"Basins ID\", add_legend=False)" ] }, { @@ -429,12 +321,14 @@ "\n", "# export the hydrography data as tif files (one per variable)\n", "merit_adapted.raster.to_mapstack(\n", - " root = os.path.join(output_path, \"merit_adapted\"),\n", - " driver = \"GTiff\",\n", + " root=os.path.join(output_path, \"merit_adapted\"),\n", + " driver=\"GTiff\",\n", ")\n", "\n", "# export the basin index as geosjon\n", - "gdf_basins.to_file(os.path.join(output_path, \"merit_adapted_basins.geojson\"), driver=\"GeoJSON\")" + "gdf_basins.to_file(\n", + " os.path.join(output_path, \"merit_adapted_basins.geojson\"), driver=\"GeoJSON\"\n", + ")" ] }, { @@ -522,7 +416,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.11.3" + "version": "3.11.4" }, "orig_nbformat": 4 },