diff --git a/notebooks/de_dommel/01_fix_model_network.py b/notebooks/de_dommel/01_fix_model_network.py index 1c87870..f51140c 100644 --- a/notebooks/de_dommel/01_fix_model_network.py +++ b/notebooks/de_dommel/01_fix_model_network.py @@ -4,7 +4,7 @@ from ribasim import Node from ribasim.nodes import basin, level_boundary, manning_resistance, outlet from ribasim_nl import CloudStorage, Model, NetworkValidator -from shapely.geometry import Point +from shapely.geometry import Point, Polygon cloud = CloudStorage() @@ -117,13 +117,15 @@ model.level_boundary.add(boundary_node, [level_data]) model.edge.add(model.tabulated_rating_curve[614], model.level_boundary[28]) - +# %% # see: https://github.com/Deltares/Ribasim-NL/issues/102#issuecomment-2292014475 model.remove_node(node_id=1898, remove_edges=True) # see: https://github.com/Deltares/Ribasim-NL/issues/102#issuecomment-2292017813 -for node_id in [1891, 989, 1058]: - model.remove_node(node_id, remove_edges=True) +# for node_id in [1891, 989, 1058]: +# model.remove_node(node_id, remove_edges=True) +model.update_node(989, "Outlet", [outlet.Static(flow_rate=[0])]) +model.update_node(1891, "LevelBoundary", [level_data]) # see: https://github.com/Deltares/Ribasim-NL/issues/102#issuecomment-2291988317 # for from_node_id, to_node_id in [799, 1580, 625, 1123, 597, 978]: @@ -156,11 +158,50 @@ if not network_validator.node_internal_basin().empty: raise Exception("nog steeds interne basins") + +# see: https://github.com/Deltares/Ribasim-NL/issues/102#issuecomment-2367724440 +gdf = gpd.read_file( + cloud.joinpath("DeDommel", "verwerkt", "4_ribasim", "hydamo.gpkg"), + layer="stuw", + engine="pyogrio", + fid_as_index=True, +) +kst = gdf.loc[35] +geometry = Point(kst.geometry.x, kst.geometry.y) +name = kst.CODE +meta_object_type = "stuw" + +outlet_node_id = model.next_node_id + +kst_node = model.outlet.add( + Node(node_id=outlet_node_id, geometry=geometry, name=name, meta_object_type=meta_object_type), + [outlet_data], +) + + +gdf = gpd.read_file( + cloud.joinpath("DeDommel", "verwerkt", "4_ribasim", "hydamo.gpkg"), + layer="hydroobject", + engine="pyogrio", + fid_as_index=True, +) +geometry = gdf.at[2822, "geometry"].interpolate(0.5, normalized=True) +basin_node_id = model.next_node_id +basin_node = model.basin.add( + Node(node_id=basin_node_id, geometry=geometry, meta_krw_name="Witte Loop/Peelrijt", meta_krw_id="NL27_KD_3_2"), + basin_data, +) + +model.remove_edge(from_node_id=664, to_node_id=8, remove_disconnected_nodes=False) +model.edge.add(model.manning_resistance[664], basin_node) +model.edge.add(basin_node, kst_node) +model.edge.add(kst_node, model.level_boundary[8]) + +# see: https://github.com/Deltares/Ribasim-NL/issues/102#issuecomment-2293486609 df = network_validator.edge_incorrect_type_connectivity( from_node_type="ManningResistance", to_node_type="LevelBoundary" ) -# see: https://github.com/Deltares/Ribasim-NL/issues/102#issuecomment-2293486609 for node_id in df.from_node_id: model.update_node(node_id, "Outlet", [outlet.Static(flow_rate=[100])]) @@ -174,9 +215,50 @@ model.basin.area.df = model.basin.area.df[~model.basin.area.df.node_id.isin(model.unassigned_basin_area.node_id)] -# # %% write model + +# fix basin area + +# see: https://github.com/Deltares/Ribasim-NL/issues/102#issuecomment-2370880128 +basin_polygon = model.basin.area.df.union_all() +holes = [Polygon(interior) for polygon in basin_polygon.buffer(10).buffer(-10).geoms for interior in polygon.interiors] +geoseries = gpd.GeoSeries(holes, crs=28992) + +drainage_areas_df = gpd.read_file( + cloud.joinpath("DeDommel", "verwerkt", "4_ribasim", "areas.gpkg"), layer="drainage_areas" +) + +drainage_areas_df = drainage_areas_df[drainage_areas_df.buffer(-10).intersects(basin_polygon)] + +for idx, geometry in enumerate(geoseries): + # select drainage-area + drainage_area_select = drainage_areas_df[drainage_areas_df.contains(geometry.buffer(-10))] + if not drainage_area_select.empty: + if not len(drainage_area_select) == 1: + raise ValueError("hole contained by multiple drainage areas, can't fix that yet") + + drainage_area = drainage_area_select.iloc[0].geometry + + # find basin_id to merge to + selected_basins_df = model.basin.area.df[model.basin.area.df.buffer(-10).within(drainage_area)].set_index( + "node_id" + ) + intersecting_basins_df = selected_basins_df.intersection(geometry.buffer(10)) + assigned_basin_id = selected_basins_df.intersection(geometry.buffer(10)).area.idxmax() + + # clip and merge geometry + geometry = geometry.buffer(10).difference(basin_polygon) + geometry = ( + model.basin.area.df.set_index("node_id") + .at[assigned_basin_id, "geometry"] + .union(geometry) + .buffer(0.1) + .buffer(-0.1) + ) + model.basin.area.df.loc[model.basin.area.df.node_id == assigned_basin_id, "geometry"] = geometry + +# %% write model model.edge.df.reset_index(drop=True, inplace=True) model.edge.df.index.name = "edge_id" -ribasim_toml = ribasim_toml = cloud.joinpath("DeDommel", "modellen", "DeDommel_fix_model_network", "model.toml") +ribasim_toml = cloud.joinpath("DeDommel", "modellen", "DeDommel_fix_model_network", "model.toml") model.write(ribasim_toml) diff --git a/notebooks/de_dommel/03_fix_basin_area.py b/notebooks/de_dommel/03_fix_basin_area.py index 7e8dd0c..50a1636 100644 --- a/notebooks/de_dommel/03_fix_basin_area.py +++ b/notebooks/de_dommel/03_fix_basin_area.py @@ -23,13 +23,14 @@ dissolved_area_gdf.to_file(cloud.joinpath("DeDommel", "verwerkt", "water_area.gpkg")) # %% -basin_df = model.basin.node.df -basin_area_df = gpd.read_file( - cloud.joinpath("DeDommel", "verwerkt", "basin_area.gpkg"), engine="pyogrio", fid_as_index=True -) +basin_area_gpkg = cloud.joinpath("DeDommel", "verwerkt", "basin_area.gpkg") +basin_area_df = model.basin.area.df +basin_area_df.to_file(basin_area_gpkg) basin_area_df.set_index("node_id", inplace=True) +basin_df = model.basin.node.df +# %% data = [] ignore_basins = [1278, 1228, 1877, 1030] row = next(i for i in basin_df.itertuples() if i.Index == 1230) @@ -76,6 +77,8 @@ area_df = gpd.GeoDataFrame(data, crs=model.basin.node.df.crs) area_df = area_df[~area_df.is_empty] area_df.index.name = "fid" +mask = area_df.geometry.type == "Polygon" +area_df.loc[mask, "geometry"] = area_df.geometry[mask].apply(lambda x: MultiPolygon([x])) model.basin.area.df = area_df # %% ribasim_toml = cloud.joinpath("DeDommel", "modellen", "DeDommel_fix_areas", "model.toml") diff --git a/notebooks/de_dommel/05_add_berging.py b/notebooks/de_dommel/05_add_berging.py index aee0d5d..0f32eba 100644 --- a/notebooks/de_dommel/05_add_berging.py +++ b/notebooks/de_dommel/05_add_berging.py @@ -1,8 +1,11 @@ # %% + import geopandas as gpd import numpy as np +import numpy.typing as npt import pandas as pd import rasterio +from rasterio.windows import from_bounds from rasterstats import zonal_stats from ribasim import Node from ribasim.nodes import basin, tabulated_rating_curve @@ -15,64 +18,191 @@ ribasim_toml = cloud.joinpath("DeDommel", "modellen", "DeDommel_parameterized", "model.toml") model = Model.read(ribasim_toml) +banden = { + "maaiveld": 1, + "bodemhoogte_primair_winter": 2, + "bodemhoogte_primair_zomer": 3, + "bodemhoogte_secundair_winter": 4, + "bodemhoogte_secundair_zomer": 5, + "bodemhoogte_tertiair_winter": 6, + "bodemhoogte_tertiair_zomer": 7, + "ghg_2010-2019": 8, + "glg_2010-2019": 9, + "opp_primair": 10, + "opp_secundair": 11, + "opp_tertiair": 12, +} basin_area_df = gpd.read_file( cloud.joinpath("DeDommel", "verwerkt", "basin_area.gpkg"), engine="pyogrio", fid_as_index=True ) basin_area_df.set_index("node_id", inplace=True) -maaiveld_raster = cloud.joinpath("Basisgegevens", "LHM", "4.3", "GxG", "Clipped_AHN_F250_IN_MNAP.tif") -ghg_raster = cloud.joinpath("Basisgegevens", "LHM", "4.3", "GxG", "GHG_LHM431_1991-2020_CORRECTED_L_tov_MV.tif") -glg_raster = cloud.joinpath("Basisgegevens", "LHM", "4.3", "GxG", "GLG_LHM431_1991-2020_CORRECTED_L_tov_MV.tif") -ma_raster = cloud.joinpath("Rijkswaterstaat", "verwerkt", "bergingsknopen", "spafvoer1.tif") +lhm_rasters = cloud.joinpath("Basisgegevens", "LHM", "4.3", "input", "LHM_data.tif") +ma_raster = cloud.joinpath("Basisgegevens", "VanDerGaast_QH", "spafvoer1.tif") -def sample_raster(raster_file, df, all_touched=False): +def sample_raster( + raster_file, + df, + band=1, + fill_value: float | None = None, + all_touched=False, + stats="mean", + maaiveld_data: npt.ArrayLike | None = None, +): with rasterio.open(raster_file) as raster_src: - data = raster_src.read(1) + # read band + data = raster_src.read(band) + + if maaiveld_data is not None: + data = maaiveld_data - data + + # fill nodata + if fill_value is not None: + data = np.where(data == raster_src.nodata, fill_value, data) + affine = raster_src.transform - return zonal_stats(df, data, affine=affine, stats="mean", nodata=raster_src.nodata, all_touched=all_touched) + return zonal_stats(df, data, affine=affine, stats=stats, nodata=raster_src.nodata, all_touched=all_touched) -def get_rating_curve(row, area, maaiveld: None | float = None): +def get_rating_curve(row, min_level, maaiveld: None | float = None): flow_rate = np.round([0, row.ma * 0.2, row.ma * 0.33, row.ma / 2, row.ma * 2], decimals=2) - level = np.round([row.glg + 1, row.glg, row.ghg / 2, row.ghg, 0], decimals=2) + depth = np.round([row.glg + 1, row.glg, row.ghg, row.ghg / 2, 0], decimals=2) + + # set GxG < 0 to 0 + depth[depth < 0] = 0 # level relative to maaiveld if maaiveld is not None: - level = maaiveld - level + level = maaiveld - depth else: - level = row.maaiveld - level + level = row.maaiveld - depth + + # make sure level >= min_level + level[level < min_level] = min_level # flow_rate in m3/s - flow_rate = flow_rate / 1000 * area / 86400 + flow_rate = flow_rate / 1000 * row.geometry.area / 86400 + + df = pd.DataFrame({"level": np.round(level, decimals=2), "flow_rate": np.round(flow_rate, decimals=5)}) + df.drop_duplicates("level", keep="first", inplace=True) + df.drop_duplicates("flow_rate", keep="last", inplace=True) + + return tabulated_rating_curve.Static(level=df.level, flow_rate=df.flow_rate) + + +def get_basin_profile(basin_polygon, polygon, max_level, min_level): + with rasterio.open(lhm_rasters) as src: + level = np.array([], dtype=float) + area = np.array([], dtype=float) + + # Get the window and its transform + window = from_bounds(*basin_polygon.bounds, transform=src.transform) + + if (window.width) < 1 or (window.height < 1): + window = from_bounds(*basin_polygon.centroid.buffer(125).bounds, transform=src.transform) + window_transform = src.window_transform(window) + + # Primary water bottom-level + window_data = src.read(3, window=window) + + # We don't want hoofdwater / doorgaand water to be in profile + if (polygon is None) | (window_data.size == 0): + mask = ~np.isnan(window_data) + else: + mask = rasterio.features.geometry_mask( + [polygon], window_data.shape, window_transform, all_touched=True, invert=True + ) + # Include nodata as False in mask + mask[np.isnan(window_data)] = False + + # add levels + level = np.concat([level, window_data[mask].ravel()]) + + # add areas on same mask + window_data = src.read(10, window=window) + area = np.concat([area, window_data[mask].ravel()]) + + # Secondary water + window_data = src.read(5, window=window) + mask = ~np.isnan(window_data) + level = np.concat([level, window_data[mask].ravel()]) + + window_data = src.read(11, window=window) + area = np.concat([area, window_data[mask].ravel()]) + + # Tertiary water water + window_data = src.read(7, window=window) + mask = ~np.isnan(window_data) + level = np.concat([level, window_data[mask].ravel()]) - return tabulated_rating_curve.Static(level=level, flow_rate=flow_rate) + window_data = src.read(12, window=window) + area = np.concat([area, window_data[mask].ravel()]) + + # Make sure area is never larger than polygon-area + area[area > basin_polygon.area] = basin_polygon.area + + # If area is empty, we add min_level at 5% of polygon-area + if area.size == 0: + level = np.append(level, min_level) + area = np.append(area, basin_polygon.area * 0.05) + + # Add extra row with max_level at basin_polygon.area + level = np.append(level, max_level) + area = np.append(area, basin_polygon.area) + + # In pandas for magic + df = pd.DataFrame({"level": np.round(level, decimals=2), "area": np.round(area)}) + df.sort_values(by="level", inplace=True) + df = df.set_index("level").cumsum().reset_index() + df.dropna(inplace=True) + df.drop_duplicates("level", keep="last", inplace=True) + + # Return profile + return basin.Profile(area=df.area, level=df.level) # %% add columns +with rasterio.open(lhm_rasters) as raster_src: + # read band + maaiveld_data = raster_src.read(banden["maaiveld"]) + # sample rasters -ghg = sample_raster(ghg_raster, basin_area_df, all_touched=True) +ghg = sample_raster( + raster_file=lhm_rasters, + df=basin_area_df, + band=banden["ghg_2010-2019"], + all_touched=True, + maaiveld_data=maaiveld_data, +) basin_area_df.loc[:, ["ghg"]] = pd.Series(dtype=float) basin_area_df.loc[:, ["ghg"]] = [i["mean"] for i in ghg] -glg = sample_raster(glg_raster, basin_area_df, all_touched=True) +glg = sample_raster( + raster_file=lhm_rasters, + df=basin_area_df, + band=banden["glg_2010-2019"], + all_touched=True, + maaiveld_data=maaiveld_data, +) basin_area_df.loc[:, ["glg"]] = pd.Series(dtype=float) basin_area_df.loc[:, ["glg"]] = [i["mean"] for i in glg] -ma = sample_raster(ma_raster, basin_area_df, all_touched=True) +ma = sample_raster(raster_file=ma_raster, df=basin_area_df, all_touched=True, fill_value=37) # 37mm/dag is basin_area_df.loc[:, ["ma"]] = pd.Series(dtype=float) basin_area_df.loc[:, ["ma"]] = [i["mean"] for i in ma] -maaiveld = sample_raster(glg_raster, basin_area_df, all_touched=True) +maaiveld = sample_raster( + raster_file=lhm_rasters, df=basin_area_df, band=banden["maaiveld"], stats="mean min max", all_touched=True +) basin_area_df.loc[:, ["maaiveld"]] = pd.Series(dtype=float) basin_area_df.loc[:, ["maaiveld"]] = [i["mean"] if pd.isna(i["mean"]) else i["mean"] for i in maaiveld] - -# %% - -# update model +basin_area_df.loc[:, ["maaiveld_max"]] = [i["max"] if pd.isna(i["max"]) else i["max"] for i in maaiveld] +basin_area_df.loc[:, ["maaiveld_min"]] = [i["min"] if pd.isna(i["min"]) else i["min"] for i in maaiveld] +# %%update model edge_id = model.edge.df.index.max() + 1 for row in model.basin.node.df.itertuples(): # row = next(row for row in model.basin.node.df.itertuples() if row.Index == 1013) @@ -93,9 +223,20 @@ def get_rating_curve(row, area, maaiveld: None | float = None): meta_categorie="bergend", geometry=geometry, ) + + if node_id in model.basin.area.df.node_id.to_list(): + polygon = model.basin.area.df.set_index("node_id").at[node_id, "geometry"] + else: + polygon = None + + max_level = max_level = basin_area_df.at[node_id, "maaiveld_max"] + min_level = max_level = basin_area_df.at[node_id, "maaiveld_min"] + if min_level == max_level: + min_level -= 0.1 + basin_profile = get_basin_profile(basin_polygon, polygon, max_level=max_level, min_level=min_level) data = [ - basin.Profile(level=[0, 1], area=[5, round(basin_polygon.area)]), - basin.State(level=[1]), + basin_profile, + basin.State(level=[basin_profile.df.level.min() + 0.1]), basin.Area(geometry=[basin_polygon]), ] basin_node = model.basin.add(node=node, tables=data) @@ -112,9 +253,9 @@ def get_rating_curve(row, area, maaiveld: None | float = None): geometry=geometry, ) if any(pd.isna(getattr(basin_row, i)) for i in ["ghg", "glg", "ma"]): - data = [tabulated_rating_curve.Static(level=[0, 1], flow_rate=[0, 1])] + raise ValueError(f"No valid ghg, glg and/or ma for basin_id {node_id}") else: - data = [get_rating_curve(row=basin_row, area=row.geometry.area)] + data = [get_rating_curve(row=basin_row, min_level=basin_profile.df.level.min())] tbr_node = model.tabulated_rating_curve.add(node=node, tables=data) # add edges @@ -128,4 +269,10 @@ def get_rating_curve(row, area, maaiveld: None | float = None): # %% ribasim_toml = cloud.joinpath("DeDommel", "modellen", "DeDommel_bergend", "model.toml") + +# if ribasim_toml.parent.exists(): +# shutil.rmtree(ribasim_toml.parent) + model.write(ribasim_toml) + +# %% diff --git a/src/ribasim_nl/ribasim_nl/model.py b/src/ribasim_nl/ribasim_nl/model.py index e076c73..4311a35 100644 --- a/src/ribasim_nl/ribasim_nl/model.py +++ b/src/ribasim_nl/ribasim_nl/model.py @@ -154,6 +154,10 @@ def remove_node(self, node_id: int, remove_edges: bool = False): from_node_id=row.from_node_id, to_node_id=row.to_node_id, remove_disconnected_nodes=False ) + # remove from used node-ids so we can add it again in the same table + if node_id in table._parent._used_node_ids: + table._parent._used_node_ids.node_ids.remove(node_id) + def update_node(self, node_id, node_type, data, node_properties: dict = {}): existing_node_type = self.node_table().df.at[node_id, "node_type"]