Skip to content

Commit

Permalink
Add function to compare model results with measurements (#258)
Browse files Browse the repository at this point in the history 
* fix error in modpath

* add caching function for dataframe, geodataframes and obscollections

* add functions to compare bro measurements with model data, #51

* add example for bro methods

* codacy

* fixes for comments Davíd and update to hydropandas 0.9.1

---------

Co-authored-by: Davíd Brakenhoff <[email protected]>
  • Loading branch information
@OnnoEbbens @dbrakenhoff
OnnoEbbens and dbrakenhoff authored Oct 10, 2023
1 parent 66aa48a commit 386c916
Show file tree
Hide file tree
Showing 5 changed files with 397 additions and 74 deletions.
105 changes: 32 additions & 73 deletions docs/examples/09_schoonhoven.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -29,7 +29,7 @@
"import geopandas as gpd\n",
"from nlmod.plot import DatasetCrossSection\n",
"from shapely.geometry import LineString, Point\n",
"import warnings\n"
"import warnings"
]
},
{
Expand All @@ -41,7 +41,7 @@
"source": [
"print(f\"nlmod version: {nlmod.__version__}\")\n",
"\n",
"nlmod.util.get_color_logger(\"INFO\")"
"nlmod.util.get_color_logger(\"INFO\")\n"
]
},
{
Expand All @@ -64,7 +64,7 @@
"model_ws = \"schoonhoven\"\n",
"figdir, cachedir = nlmod.util.get_model_dirs(model_ws)\n",
"extent = [116_500, 120_000, 439_000, 442_000]\n",
"time = pd.date_range(\"2020\", \"2023\", freq=\"MS\") # monthly timestep\n"
"time = pd.date_range(\"2020\", \"2023\", freq=\"MS\") # monthly timestep"
]
},
{
Expand Down Expand Up @@ -98,7 +98,7 @@
" f\"{fname_bgt} not found. Please run notebook 02_surface_water.ipynb first\"\n",
" )\n",
" )\n",
"bgt = gpd.read_file(fname_bgt)\n"
"bgt = gpd.read_file(fname_bgt)"
]
},
{
Expand All @@ -121,7 +121,7 @@
"norm = matplotlib.colors.Normalize(vmin=-3, vmax=1)\n",
"cmap = \"viridis\"\n",
"bgt.plot(\"summer_stage\", ax=ax, norm=norm, cmap=cmap)\n",
"nlmod.plot.colorbar_inside(norm=norm, cmap=cmap);"
"nlmod.plot.colorbar_inside(norm=norm, cmap=cmap)"
]
},
{
Expand All @@ -141,7 +141,7 @@
"source": [
"f, ax = nlmod.plot.get_map(extent)\n",
"bgt.plot(\"ahn_min\", ax=ax, norm=norm, cmap=cmap)\n",
"nlmod.plot.colorbar_inside(norm=norm, cmap=cmap);"
"nlmod.plot.colorbar_inside(norm=norm, cmap=cmap)"
]
},
{
Expand All @@ -167,7 +167,7 @@
" cachedir=cachedir,\n",
" cachename=\"layer_model.nc\",\n",
")\n",
"layer_model\n"
"layer_model"
]
},
{
Expand All @@ -186,7 +186,7 @@
"outputs": [],
"source": [
"ds = nlmod.to_model_ds(layer_model, model_name, model_ws, delr=100.0, delc=100.0)\n",
"ds\n"
"ds"
]
},
{
Expand All @@ -208,7 +208,7 @@
"outputs": [],
"source": [
"refinement_features = [(bgt[bgt[\"bronhouder\"] == \"L0002\"], 2)]\n",
"ds = nlmod.grid.refine(ds, refinement_features=refinement_features)\n"
"ds = nlmod.grid.refine(ds, refinement_features=refinement_features)"
]
},
{
Expand Down Expand Up @@ -245,7 +245,7 @@
"outputs": [],
"source": [
"knmi_ds = nlmod.read.knmi.get_recharge(ds, cachedir=cachedir, cachename=\"recharge.nc\")\n",
"ds.update(knmi_ds)\n"
"ds.update(knmi_ds)"
]
},
{
Expand Down Expand Up @@ -289,7 +289,7 @@
"oc = nlmod.gwf.oc(ds, gwf)\n",
"\n",
"# create storagee package\n",
"sto = nlmod.gwf.sto(ds, gwf)\n"
"sto = nlmod.gwf.sto(ds, gwf)"
]
},
{
Expand Down Expand Up @@ -341,7 +341,7 @@
"lakes[\"name\"] = \"\"\n",
"lakes.loc[lakes[\"identificatie\"].isin(ids_grote_gracht), \"name\"] = \"grote gracht\"\n",
"lakes.loc[lakes[\"identificatie\"].isin(ids_oude_haven), \"name\"] = \"oude haven\"\n",
"bgt_grid = bgt_grid[~mask]\n"
"bgt_grid = bgt_grid[~mask]"
]
},
{
Expand All @@ -363,7 +363,7 @@
"lek[\"stage\"] = 0.0\n",
"lek[\"rbot\"] = -3.0\n",
"spd = nlmod.gwf.surface_water.build_spd(lek, \"RIV\", ds)\n",
"riv = flopy.mf6.ModflowGwfriv(gwf, stress_period_data={0: spd})\n"
"riv = flopy.mf6.ModflowGwfriv(gwf, stress_period_data={0: spd})"
]
},
{
Expand All @@ -382,7 +382,7 @@
"metadata": {},
"outputs": [],
"source": [
"drn = nlmod.gwf.surface_water.gdf_to_seasonal_pkg(bgt_grid, gwf, ds)"
"drn = nlmod.gwf.surface_water.gdf_to_seasonal_pkg(bgt_grid, gwf, ds)\n"
]
},
{
Expand Down Expand Up @@ -426,7 +426,7 @@
"] = 1.0 # overstort hoogte\n",
"\n",
"# add lake to groundwaterflow model\n",
"nlmod.gwf.lake_from_gdf(gwf, lakes, ds, boundname_column=\"name\")"
"nlmod.gwf.lake_from_gdf(gwf, lakes, ds, boundname_column=\"name\")\n"
]
},
{
Expand All @@ -437,7 +437,7 @@
"outputs": [],
"source": [
"# create recharge package\n",
"rch = nlmod.gwf.rch(ds, gwf)\n"
"rch = nlmod.gwf.rch(ds, gwf)"
]
},
{
Expand All @@ -455,7 +455,7 @@
"metadata": {},
"outputs": [],
"source": [
"nlmod.sim.write_and_run(sim, ds)\n"
"nlmod.sim.write_and_run(sim, ds)"
]
},
{
Expand All @@ -474,7 +474,7 @@
"metadata": {},
"outputs": [],
"source": [
"head = nlmod.gwf.get_heads_da(ds)\n"
"head = nlmod.gwf.get_heads_da(ds)"
]
},
{
Expand Down Expand Up @@ -532,7 +532,7 @@
"cbar = nlmod.plot.colorbar_inside(pc)\n",
"dcs.plot_grid()\n",
"dcs.plot_layers(colors=\"none\", min_label_area=1000)\n",
"f.tight_layout(pad=0.0)"
"f.tight_layout(pad=0.0)\n"
]
},
{
Expand All @@ -555,7 +555,7 @@
"head_point = nlmod.gwf.get_head_at_point(head, x=x, y=y, ds=ds)\n",
"fig, ax = plt.subplots(1, 1, figsize=(10, 8))\n",
"handles = head_point.plot.line(ax=ax, hue=\"layer\")\n",
"ax.set_ylabel(\"head [m NAP]\");"
"ax.set_ylabel(\"head [m NAP]\")"
]
},
{
Expand All @@ -575,17 +575,17 @@
"source": [
"df = pd.read_csv(os.path.join(model_ws, \"lak_STAGE.csv\"), index_col=0)\n",
"df.index = ds.time.values\n",
"ax = df.plot(figsize=(10, 3))"
"ax = df.plot(figsize=(10, 3))\n"
]
},
{
"cell_type": "raw",
"cell_type": "markdown",
"id": "9355de12",
"metadata": {
"raw_mimetype": "text/x-python"
},
"source": [
"### Compare with BRO measurements"
"## Compare with BRO measurements"
]
},
{
Expand All @@ -595,35 +595,8 @@
"raw_mimetype": "text/x-python"
},
"source": [
"fname_pklz = os.path.join(ds.cachedir, 'oc_bro.pklz')\n",
"if os.path.exists(fname_pklz):\n",
" oc = pd.read_pickle(fname_pklz)\n",
"else:\n",
" oc = hpd.read_bro(extent=ds.extent, name='BRO', tmin=ds.time.values.min(), tmax=ds.time.values.max(), )\n",
" oc.to_pickle(fname_pklz)"
]
},
{
"cell_type": "raw",
"id": "06fdd4cc-a15e-485e-b12e-fda9a464d30c",
"metadata": {
"raw_mimetype": "text/x-python"
},
"source": [
"# get modellayers\n",
"oc['modellayer'] = oc.gwobs.get_modellayers(ds=ds)"
]
},
{
"cell_type": "raw",
"id": "7831eddc-3b69-4cc5-b1c0-7ee09732a2f9",
"metadata": {
"raw_mimetype": "text/x-python"
},
"source": [
"# get modelled head at measurement points\n",
"ds['heads'] = nlmod.gwf.get_heads_da(ds)\n",
"oc_modflow = hpd.read_modflow(oc, gwf, ds['heads'].values, ds.time.values)"
"oc = nlmod.read.bro.get_bro(extent, cachedir=cachedir, cachename=\"bro\")\n",
"oc_mod = nlmod.read.bro.add_modelled_head(oc, gwf, ds=ds)"
]
},
{
Expand All @@ -633,22 +606,8 @@
"raw_mimetype": "text/x-python"
},
"source": [
"# add modelled head to measured heads\n",
"obs_list_map = []\n",
"for gld in oc.index:\n",
" o = oc.loc[gld,'obs'].resample('D').last().sort_index()\n",
" modelled = oc_modflow.loc[gld, 'obs']\n",
" modelled = hpd.GroundwaterObs(modelled.rename(columns={0: 'values'}), name=f'{o.name}_mod_lay{oc.loc[gld,\"modellayer\"]}', x=o.x, y=o.y, \n",
" tube_nr=o.tube_nr+1,screen_top=o.screen_top, screen_bottom=o.screen_bottom, \n",
" tube_top=o.tube_top, monitoring_well=o.monitoring_well, source='MODFLOW', unit= 'm NAP',\n",
" ground_level=o.ground_level, metadata_available=o.metadata_available)\n",
" obs_list_map.append(o)\n",
" obs_list_map.append(modelled)\n",
"\n",
"oc_map = hpd.ObsCollection.from_list(obs_list_map, name='meting+model')\n",
"\n",
"# create interactive map\n",
"oc_map.plots.interactive_map(os.path.join(ds.figdir, 'iplots'))"
"oc_mod.plots.interactive_map(\"figures\", add_screen_to_legend=True)"
]
},
{
Expand Down Expand Up @@ -689,7 +648,7 @@
" ha=\"center\",\n",
" va=\"top\",\n",
" transform=ax.transAxes,\n",
" )\n"
" )"
]
},
{
Expand Down Expand Up @@ -722,7 +681,7 @@
"pg = nlmod.modpath.pg_from_pd(nodes, localx=0.5, localy=0.5, localz=0.5)\n",
"\n",
"# create the modpath simulation file\n",
"mpsim = nlmod.modpath.sim(mpf, pg, \"forward\", gwf=gwf)\n"
"mpsim = nlmod.modpath.sim(mpf, pg, \"forward\", gwf=gwf)"
]
},
{
Expand All @@ -733,7 +692,7 @@
"outputs": [],
"source": [
"# run modpath model\n",
"nlmod.modpath.write_and_run(mpf, script_path=\"10_modpath.ipynb\")"
"nlmod.modpath.write_and_run(mpf, script_path=\"10_modpath.ipynb\")\n"
]
},
{
Expand All @@ -743,7 +702,7 @@
"metadata": {},
"outputs": [],
"source": [
"pdata = nlmod.modpath.load_pathline_data(mpf)\n"
"pdata = nlmod.modpath.load_pathline_data(mpf)"
]
},
{
Expand Down Expand Up @@ -786,7 +745,7 @@
"line = ax.add_collection(lc)\n",
"nlmod.plot.colorbar_inside(line, label=\"Travel time (years)\")\n",
"\n",
"bgt.plot(ax=ax, edgecolor=\"k\", facecolor=\"none\")"
"bgt.plot(ax=ax, edgecolor=\"k\", facecolor=\"none\")\n"
]
},
{
Expand Down Expand Up @@ -822,7 +781,7 @@
"dcs.plot_grid()\n",
"# add labels with layer names\n",
"dcs.plot_layers(alpha=0.0, min_label_area=1000)\n",
"f.tight_layout(pad=0.0)\n"
"f.tight_layout(pad=0.0)"
]
}
],
Expand Down
Loading

0 comments on commit 386c916

Please sign in to comment.