{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Example: Working with flow direction data"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In this example we illustrate some common hydrology GIS problems based on so-called flow direction data. In HydroMT we make use of functionality from [pyflwdir](https://deltares.github.io/pyflwdir/latest/) to work with this type of data. HydroMT wraps some functionality of **pyflwdir**, to make it easier to work with [raster datasets](../guides/advanced_user/data_types.rst). However, pyflwdir has much more functionality. An overview of all the flow direction methods in HydroMT can be found in the [Reference API](../api/gis.rst#flow-direction-methods).\n",
"\n",
"Here, we will showcase the following flow direction GIS cases:\n",
"\n",
"1. Derive basin and stream geometries\n",
"2. Derive flow directions from elevation data\n",
"3. Reproject flow direction data\n",
"4. Upscale flow directions\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"from pprint import pprint\n",
"\n",
"import geopandas as gpd\n",
"\n",
"import hydromt\n",
"from hydromt.gis import utm_crs"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"First, we load some data to play with from the pre-defined artifact_data data catalog. For more information about working with data in HydroMT, see [the user guide](../guides/user_guide/data_overview.rst). As an example we will use the [MERIT Hydro](http://hydro.iis.u-tokyo.ac.jp/~yamadai/MERIT_Hydro) dataset which is set of GeoTiff files with identical grids, one for each variable of the datasets. We use the flow direction (flwdir); elevation (elevtn) and upstream area (uparea) layers."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# initialize a data catalog based on the pre-defined artifact_data catalog\n",
"data_catalog = hydromt.DataCatalog(data_libs=[\"artifact_data=v1.0.0\"])\n",
"\n",
"# we load the flow direction (flwdir); elevation (elevtn) and upstream area (uparea) layers\n",
"ds = data_catalog.get_rasterdataset(\n",
" \"merit_hydro\",\n",
" bbox=[11.7, 45.8, 12.8, 46.7],\n",
" variables=[\"flwdir\", \"elevtn\", \"uparea\"],\n",
")\n",
"ds"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Derive basin and stream geometries"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"If you have existing flow direction data from sources such as MERIT Hydro, or HydroSHEDS or similar, you can use these to delineate basins and extract streams based on a user-defined threshold. To do this we need to transform the gridded flow direction data into a `FlwdirRaster` object using the [flwdir_from_da()](../_generated/hydromt.gis.flw.flwdir_from_da.rst) method. This object is at the core of the **pyflwdir** package and creates an actionable common format from a flow direction raster which describes relations between cells."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"NOTE: that for most methods a first call might be a bit slow as the numba code is compiled just in time, a second call of the same methods (also with different arguments) will be much faster!"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# instantiate a FlwdirRaster object\n",
"flwdir = hydromt.gis.flw.flwdir_from_da(ds[\"flwdir\"], ftype=\"d8\")\n",
"print(type(flwdir))\n",
"print(flwdir)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Next, we derive streams based on a 10 km2 upstream area threshold using the pyflwdir [streams](https://deltares.github.io/pyflwdir/latest/reference.html#pyflwdir.FlwdirRaster.streams) method. Pyflwdir returns a geojson like representation of the streams per stream segment, which we parse to a GeoPandas GeoDataFrame to easily plot it."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"feats = flwdir.streams(\n",
" mask=ds[\"uparea\"].values > 10,\n",
" strord=flwdir.stream_order(), # set stream order property\n",
" uparea=ds[\"uparea\"].values, # set upstream area property\n",
")\n",
"gdf_riv = gpd.GeoDataFrame.from_features(feats, crs=ds.raster.crs)\n",
"pprint(gdf_riv.head())"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Using the [basin_map()](../_generated/hydromt.gis.flw.basin_map.rst) method we can delineate all basins in our domain."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# get the best utm zone CRS for a projected CRS\n",
"utm = utm_crs(ds.raster.bounds)\n",
"ds[\"basins\"] = hydromt.gis.flw.basin_map(\n",
" ds,\n",
" flwdir,\n",
")[0]\n",
"# use the HydroMT \"raster\" data accessor to vectorize the basin raster.\n",
"gdf_bas = ds[\"basins\"].raster.vectorize()\n",
"# calculate the area of each basin in the domain and sort the dataframe\n",
"gdf_bas[\"area\"] = gdf_bas.to_crs(utm).area / 1e6 # km2\n",
"gdf_bas = gdf_bas.sort_values(\"area\", ascending=False)\n",
"pprint(gdf_bas.head())"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# plot the results\n",
"ax = gdf_bas[:5].boundary.plot(color=\"r\", lw=1, zorder=2)\n",
"gdf_riv.plot(\n",
" zorder=2,\n",
" ax=ax,\n",
" color=\"darkblue\",\n",
" lw=gdf_riv[\"strord\"] / 8,\n",
")\n",
"ds[\"elevtn\"].plot(cmap=\"terrain\", ax=ax, vmin=-500, vmax=2000, alpha=0.7)\n",
"ax.set_title(\"Streams (darkblue) and basins (red)\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Derive flow directions from elevation data "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"If you don't have flow direction data available these can be derived from an elevation raster. HydroMT implements the algorithm proposed by [Wang & Liu (2006)](https://www.tandfonline.com/doi/abs/10.1080/13658810500433453) to do this. We use the [d8_from_dem()](../_generated/hydromt.gis.flw.d8_from_dem.rst) method which wraps the pyflwdir [fill_depressions()](https://deltares.github.io/pyflwdir/latest/reference.html#pyflwdir.dem.fill_depressions) method. \n",
"\n",
"The derivation of flow direction can be aided by a river shape file with an upstream area (\"uparea\") property. Try uncommenting the `gdf_stream` argument and compare the results."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# derive flow directions raster from elevation\n",
"da_flw = hydromt.gis.flw.d8_from_dem(\n",
" ds[\"elevtn\"],\n",
")\n",
"# parse it into a FlwdirRaster object\n",
"flwdir1 = hydromt.gis.flw.flwdir_from_da(da_flw, ftype=\"d8\")\n",
"# derive streams based on a 10 km2 threshold\n",
"feats1 = flwdir1.streams(mask=flwdir1.upstream_area(\"km2\") > 10)\n",
"gdf_riv1 = gpd.GeoDataFrame.from_features(feats1, crs=ds.raster.crs)\n",
"\n",
"# plot the new streams (red) and compare with the original (darkblue)\n",
"ax = gdf_riv.plot(zorder=2, color=\"darkblue\")\n",
"gdf_riv1.plot(zorder=2, ax=ax, color=\"r\")\n",
"ds[\"elevtn\"].plot(cmap=\"terrain\", ax=ax, vmin=-500, vmax=2000, alpha=0.7)\n",
"ax.set_title(\"Original (darkblue) and new (red) streams\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Reproject flow direction data"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Unlike continuous data such as elevation or data with discrete classes such as land use, flow direction data cannot simply be reclassified using common resampling methods. Instead, with the [reproject_hydrography_like()](../_generated/hydromt.gis.flw.reproject_hydrography_like.rst) a synthetic elevation grid is created based on an upstream area raster, this is reprojected and used to derive a new flow direction grid with the method described above. Note that this works well if we keep approximately the same resolution. For upscaling to larger grid cells different algorithms should be used, see next example."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# reproject the elevation grid first\n",
"da_elv_reproj = ds[\"elevtn\"].raster.reproject(dst_crs=utm)\n",
"# reproject the flow direction data\n",
"ds_reproj = hydromt.gis.flw.reproject_hydrography_like(\n",
" ds, # flow direction and upstream area grids\n",
" da_elv=da_elv_reproj, # destination grid\n",
")\n",
"# parse it into a FlwdirRaster object\n",
"flwdir_reproj = hydromt.gis.flw.flwdir_from_da(ds_reproj[\"flwdir\"], ftype=\"d8\")\n",
"# derive streams based on a 10 km2 threshold\n",
"feats_reproj = flwdir_reproj.streams(mask=flwdir_reproj.upstream_area(\"km2\") > 10)\n",
"gdf_riv_reproj = gpd.GeoDataFrame.from_features(feats_reproj, crs=ds_reproj.raster.crs)\n",
"\n",
"# plot the streams from the reproject data (red) and compare with the original (darkblue)\n",
"# NOTE the different coordinates on the figure axis\n",
"ax = gdf_riv_reproj.plot(zorder=3, color=\"r\")\n",
"gdf_riv.to_crs(utm).plot(ax=ax, zorder=2, color=\"darkblue\")\n",
"da_elv_reproj.raster.mask_nodata().plot(\n",
" cmap=\"terrain\", ax=ax, vmin=-500, vmax=2000, alpha=0.7\n",
")\n",
"ax.set_title(\"Original (darkblue) and new reprojected (red) streams\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Upscale flow directions"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Methods to upscale flow directions are required as models often have a coarser resolution than the elevation data used to build them. Instead of deriving flow directions from upscaled elevation data, it is better to directly upscale the flow direction data itself. The [upscale_flwdir()](../_generated/hydromt.gis.flw.upscale_flwdir.rst) method wraps a pyflwdir method that implements the recently developed Iterative Hydrography Upscaling (IHU) algorithm [(Eilander et al 2020)](https://hess.copernicus.org/articles/25/5287/2021/). Try different upscale factors and see the difference!"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# upscale flow direction with a factor \"scale_ratio\"\n",
"# this returns both a flow direction grid and a new FlwdirRaster object\n",
"da_flw_lowres, flwdir_lowres = hydromt.gis.flw.upscale_flwdir(\n",
" ds, # flow direction and upstream area grids\n",
" flwdir=flwdir, # pyflwdir FlwdirRaster object\n",
" scale_ratio=20, # upscaling factor\n",
")\n",
"\n",
"# derive streams based on a 10 km2 threshold\n",
"feats_lowres = flwdir_lowres.streams(mask=flwdir_lowres.upstream_area(\"km2\") > 10)\n",
"gdf_riv_lowres = gpd.GeoDataFrame.from_features(feats_lowres, crs=ds.raster.crs)\n",
"\n",
"# plot the streams from the upscaled flow direction (red) and compare with the original (darkblue)\n",
"ax = gdf_riv_lowres.plot(zorder=3, color=\"r\")\n",
"gdf_riv.plot(ax=ax, zorder=2, color=\"darkblue\")\n",
"ds[\"elevtn\"].raster.mask_nodata().plot(\n",
" cmap=\"terrain\", ax=ax, vmin=-500, vmax=2000, alpha=0.7\n",
")\n",
"ax.set_title(\"Original (darkblue) and new upscaled (red) streams\")"
]
}
],
"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.11.9"
},
"vscode": {
"interpreter": {
"hash": "8138b72b0e304ba55d14e4dcf8d650296065d1ee95f3f67a239b6fbf5f7328dd"
}
}
},
"nbformat": 4,
"nbformat_minor": 4
}