Skip to content

Commit b43aee5

Browse files
priscavdsluispriscavdsluis
committed
Make magdalena workflow work
1 parent 93b64a1 commit b43aee5

File tree

1 file changed

+61
-35
lines changed

1 file changed

+61
-35
lines changed

‎docs/tutorials/dsd_2023/magdalena_workflow/magdalena_answers.ipynb

Lines changed: 61 additions & 35 deletions
Original file line numberDiff line numberDiff line change
@@ -5,7 +5,7 @@
55
"id": "d3a50766",
66
"metadata": {},
77
"source": [
8-
"# Tutorial Delft3D FM 1D2D fluvial flood model of Magdalena river \n",
8+
"# Creating a Delft3D FM 1D2D fluvial flood model of Magdalena river \n",
99
"\n",
1010
"This notebook is an introduction on how to use [HYDROLIB-core](https://github.com/Deltares/HYDROLIB-core) functionalities to create and adjust a Delft3D FM Suite 1D2D fluvial flood model of the Magdalena river in Columbia. First a 1D river model of the Magdalena river and Canal del Dique will be created, shown in the left figure. This model will be extended with a 2D Flexible Mesh, with a uniform roughness field. Adding a DTM on that 2D mesh is currently not part of this tutorial, but normally this is of course indispensible to mimic the fluvial flood patterns. Running the model or importing the model in the Delft3D FM Suite 1D2D is not shown in this tutorial. \n",
1111
"\n",
@@ -92,13 +92,17 @@
9292
"import numpy as np\n",
9393
"import pandas as pd\n",
9494
"import matplotlib.pyplot as plt\n",
95+
"import meshkernel as mk\n",
9596
"\n",
9697
"# I/O\n",
9798
"from pathlib import Path\n",
9899
"\n",
99100
"# vector/polygon packages\n",
100101
"import geopandas as gpd\n",
101102
"\n",
103+
"# mesh kerne\n",
104+
"from meshkernel import MeshKernel, GeometryList, MakeGridParameters\n",
105+
"\n",
102106
"# HYDROLIB-core libraries\n",
103107
"from hydrolib.core.dflowfm.mdu import FMModel, AutoStartOption\n",
104108
"from hydrolib.core.dflowfm.net.models import *\n",
@@ -303,23 +307,23 @@
303307
"outputs": [],
304308
"source": [
305309
"# Create a figure with the flowlinks of the model\n",
306-
"figure, (ax1, ax2) = plt.subplots(ncols=2)\n",
310+
"figure, ax1 = plt.subplots(ncols=1)\n",
307311
"img = plt.imread(inputdir/\"osm_background.png\")\n",
308312
"mesh1d_meshkernel = network._mesh1d._get_mesh1d()\n",
309313
"\n",
310-
"ax1.imshow(img, extent=[440000,531544.97964175534434617 ,1057000 ,1228644.01383191486820579 ])\n",
311-
"mesh1d_meshkernel.plot_edges(ax=ax1, color='red')\n",
312-
"\n",
313314
"# Add some annotation of branches:\n",
314315
"for index,branch in branches.iterrows():\n",
315316
" middleIndex = int((len(branch['geometry'].coords) - 1)/2)\n",
316317
" ax1.text(branch['geometry'].xy[0][middleIndex], branch['geometry'].xy[1][middleIndex], branch['Name'], ha='right')\n",
318+
" \n",
319+
"ax1.imshow(img, extent=[440000,\n",
320+
" 531544.97964175534434617,\n",
321+
" 1057000,\n",
322+
" 1228644.01383191486820579])\n",
323+
"mesh1d_meshkernel.plot_edges(ax=ax1, color='darkblue')\n",
324+
"\n",
325+
"\n",
317326
"\n",
318-
"# Make a second plot with a zoom-in view of the computational grid points near the junction\n",
319-
"ax2.imshow(img, extent=[440000,531544.97964175534434617 ,1057000 ,1228644.01383191486820579 ])\n",
320-
"ax2.set_xlim([500000, 525000])\n",
321-
"ax2.set_ylim([1125000 ,1150000])\n",
322-
"mesh1d_meshkernel.plot_edges(ax=ax2, color='red',marker='.',markeredgecolor='k',markersize=1)\n",
323327
"\n"
324328
]
325329
},
@@ -890,15 +894,36 @@
890894
"# Determine the boundary box\n",
891895
"extent = gdf_grid.geometry[0].bounds\n",
892896
"\n",
893-
"# Create uniform rectilinear mesh\n",
894897
"mesh2d = Mesh2d()\n",
898+
"mesh_kernel = mesh2d.meshkernel\n",
899+
"\n",
900+
"xmin, ymin, xmax, ymax = extent\n",
901+
"\n",
902+
"dx=2500.0\n",
903+
"dy=2500.0\n",
904+
"\n",
905+
"num_rows = int((ymax - ymin) / dy)\n",
906+
"num_columns = int((xmax - xmin) / dx)\n",
907+
"origin_x = xmin\n",
908+
"origin_y = ymin\n",
909+
"block_size_x = dx\n",
910+
"block_size_y = dy\n",
911+
"\n",
912+
"make_grid_parameters = MakeGridParameters()\n",
913+
"make_grid_parameters.num_columns = num_columns\n",
914+
"make_grid_parameters.num_rows = num_rows\n",
915+
"make_grid_parameters.origin_x = origin_x\n",
916+
"make_grid_parameters.origin_y = origin_y\n",
917+
"make_grid_parameters.block_size_x = dx\n",
918+
"make_grid_parameters.block_size_y = dy\n",
919+
"make_grid_parameters.upper_right_x = xmax\n",
920+
"make_grid_parameters.upper_right_y = ymax\n",
895921
"\n",
896-
"# Answer EX7\n",
897-
"mesh2d.create_rectilinear(\n",
898-
" extent=extent, \n",
899-
" dx=2500.0,\n",
900-
" dy=2500.0\n",
901-
")\n",
922+
"# Create uniform rectilinear mesh\n",
923+
"mesh_kernel.curvilinear_make_uniform(make_grid_parameters)\n",
924+
"mesh_kernel.curvilinear_convert_to_mesh2d()\n",
925+
"\n",
926+
"mesh2d._process(mesh2d.get_mesh2d())\n",
902927
"\n",
903928
"# Clip resulting mesh within polygon area of interest\n",
904929
"# Prepare a clipping polygon for MeshKernel input:\n",
@@ -957,29 +982,23 @@
957982
"metadata": {},
958983
"outputs": [],
959984
"source": [
960-
"# Create empty mesh2d object\n",
961-
"mesh2d = Mesh2d()\n",
962-
"\n",
963-
"# Create rectangular grid\n",
964-
"mesh2d.create_rectilinear(\n",
965-
" extent=extent, \n",
966-
" dx=2500.0,\n",
967-
" dy=2500.0\n",
968-
")\n",
969-
"\n",
970-
"# Clip resulting mesh within polygon area of interest\n",
971-
"# Prepare a clipping polygon for MeshKernel input:\n",
972-
"clip_polygon = np.array(gdf_grid.geometry[0].exterior.xy)\n",
973-
"clip_geom = GeometryList(clip_polygon[0],clip_polygon[1])\n",
974-
"mesh2d.clip(clip_geom)\n",
975-
"\n",
976-
"\n",
977985
"# Refine the coarse uniform grid in a buffer region around the rivers\n",
978986
"gdf_refine = gpd.read_file(inputdir / 'buffer_around_river_proj.shp')\n",
979987
"refine_polygon = np.array(gdf_refine.geometry[0].exterior.xy)\n",
980988
"refine_geom = GeometryList(refine_polygon[0], refine_polygon[1])\n",
981-
"mesh2d.refine(refine_geom, level=1)\n",
982989
"\n",
990+
"parameters = mk.MeshRefinementParameters(\n",
991+
" refine_intersected=True,\n",
992+
" use_mass_center_when_refining=False,\n",
993+
" min_edge_size=10.0,\n",
994+
" refinement_type=1,\n",
995+
" connect_hanging_nodes=True,\n",
996+
" account_for_samples_outside_face=False,\n",
997+
" max_refinement_iterations=1,\n",
998+
")\n",
999+
"\n",
1000+
"mesh_kernel.mesh2d_refine_based_on_polygon(refine_geom, parameters)\n",
1001+
"mesh2d._process(mesh2d.get_mesh2d())\n",
9831002
"fm.geometry.netfile.network._mesh2d = mesh2d\n",
9841003
"fm.geometry.netfile.filepath = \"FlowFM_1D2D_refined_net.nc\""
9851004
]
@@ -1030,6 +1049,13 @@
10301049
"fm.general.autostart = AutoStartOption.no\n",
10311050
"fm.save(recurse=True)"
10321051
]
1052+
},
1053+
{
1054+
"cell_type": "code",
1055+
"execution_count": null,
1056+
"metadata": {},
1057+
"outputs": [],
1058+
"source": []
10331059
}
10341060
],
10351061
"metadata": {

0 commit comments

Comments
 (0)