From 71697aa0f0fc8d12196f4b38578045330c1efab8 Mon Sep 17 00:00:00 2001 From: Daniel Tollenaar Date: Tue, 1 Oct 2024 13:04:43 +0200 Subject: [PATCH] Koppelen LHM 2024.9.0 (#160) Alle code voor de bouw van het LHM 2024.9.0: - Berging-code naar afzonderlijke module] - Afleiding vrij-afwaterende gebieden op het RWS-HWS - Vd Gaast ook toegepast op vrij-afwaterende gebieden op het RWS-HWS (Wilhelminakanaal) --------- Co-authored-by: Maarten Pronk --- notebooks/de_dommel/05_add_berging.py | 202 ++----------------- notebooks/koppelen/01_bergend_gebied.py | 148 ++++++++++++++ notebooks/koppelen/02_koppelen.py | 211 ++++++++++++++++++++ notebooks/koppelen/99_upload_model.py | 8 + open-vscode.bat | 1 + src/ribasim_nl/reset_index.py | 42 ++++ src/ribasim_nl/ribasim_nl/berging.py | 241 +++++++++++++++++++++++ src/ribasim_nl/ribasim_nl/concat.py | 28 ++- src/ribasim_nl/ribasim_nl/model.py | 17 +- src/ribasim_nl/ribasim_nl/reset_index.py | 34 ++-- 10 files changed, 718 insertions(+), 214 deletions(-) create mode 100644 notebooks/koppelen/01_bergend_gebied.py create mode 100644 notebooks/koppelen/02_koppelen.py create mode 100644 notebooks/koppelen/99_upload_model.py create mode 100644 open-vscode.bat create mode 100644 src/ribasim_nl/reset_index.py create mode 100644 src/ribasim_nl/ribasim_nl/berging.py diff --git a/notebooks/de_dommel/05_add_berging.py b/notebooks/de_dommel/05_add_berging.py index ea32582..bcd8f43 100644 --- a/notebooks/de_dommel/05_add_berging.py +++ b/notebooks/de_dommel/05_add_berging.py @@ -1,15 +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 +from ribasim.nodes import basin from ribasim_nl import CloudStorage, Model +from ribasim_nl.berging import add_basin_statistics, get_basin_profile, get_rating_curve from ribasim_nl.geometry import basin_to_point from shapely.geometry import LineString @@ -18,190 +14,19 @@ 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) -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, - 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: - # 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=stats, nodata=raster_src.nodata, all_touched=all_touched) - - -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) - 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 - depth - else: - 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 * 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()]) - - 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) +lhm_raster_file = cloud.joinpath("Basisgegevens", "LHM", "4.3", "input", "LHM_data.tif") +ma_raster_file = cloud.joinpath("Basisgegevens", "VanDerGaast_QH", "spafvoer1.tif") - # Return profile - return basin.Profile(area=df.area, level=df.level) +basin_area_df = add_basin_statistics(df=basin_area_df, lhm_raster_file=lhm_raster_file, ma_raster_file=ma_raster_file) -# %% add columns -with rasterio.open(lhm_rasters) as raster_src: - # read band - maaiveld_data = raster_src.read(banden["maaiveld"]) - -# sample rasters -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( - 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(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( - 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] -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(): @@ -233,7 +58,13 @@ def get_basin_profile(basin_polygon, polygon, max_level, min_level): 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) + basin_profile = get_basin_profile( + basin_polygon=basin_polygon, + polygon=polygon, + max_level=max_level, + min_level=min_level, + lhm_raster_file=lhm_raster_file, + ) data = [ basin_profile, basin.State(level=[basin_profile.df.level.min() + 0.1]), @@ -267,6 +98,15 @@ def get_basin_profile(basin_polygon, polygon, max_level, min_level): else: print(f"Geen basin-vlak voor {node_id}") +# %% +df = pd.DataFrame({"node_id": model.basin.node.df.index.to_list()}) +df.index.name = "fid" +df.loc[:, "precipitation"] = 5.787037e-08 +df.loc[:, "potential_evaporation"] = 1.157407e-08 +df.loc[:, "drainage"] = 0 +df.loc[:, "infiltration"] = 0 +model.basin.static.df = df + # %% ribasim_toml = cloud.joinpath("DeDommel", "modellen", "DeDommel_bergend", "model.toml") diff --git a/notebooks/koppelen/01_bergend_gebied.py b/notebooks/koppelen/01_bergend_gebied.py new file mode 100644 index 0000000..8038159 --- /dev/null +++ b/notebooks/koppelen/01_bergend_gebied.py @@ -0,0 +1,148 @@ +# %% +import geopandas as gpd +import pandas as pd +from ribasim import Node +from ribasim.nodes import basin +from ribasim_nl import CloudStorage, Model +from ribasim_nl.berging import add_basin_statistics, get_basin_profile, get_rating_curve +from ribasim_nl.geodataframe import split_basins +from ribasim_nl.geometry import basin_to_point +from shapely.geometry import LineString, MultiPolygon + +cloud = CloudStorage() + +# %% RWS-HWS +model_path = cloud.joinpath("Rijkswaterstaat", "modellen", "hws") +toml_file = model_path / "hws.toml" +rws_model = Model.read(toml_file) + +# %% DeDommel +model_path = cloud.joinpath("DeDommel", "modellen", "DeDommel") +toml_file = model_path / "model.toml" +model = Model.read(toml_file) +basin_polygon = model.basin.area.df[model.basin.area.df.node_id != 1228].union_all() + +drainage_area = gpd.read_file( + cloud.joinpath("DeDommel", "verwerkt", "4_ribasim", "areas.gpkg"), layer="drainage_areas" +).union_all() + + +rws_selected_areas_df = rws_model.basin.area.df[rws_model.basin.area.df.intersects(drainage_area.buffer(-10))] +rws_selected_areas = rws_selected_areas_df.union_all() + +poly = ( + rws_model.basin.area.df[rws_model.basin.area.df.intersects(drainage_area.buffer(-10))] + .buffer(0.1) + .union_all() + .buffer(3000) +) +poly = poly.difference(basin_polygon).intersection(drainage_area) + +berging_basins_df = gpd.GeoDataFrame(geometry=gpd.GeoSeries(poly.geoms, crs=28992)) + +berging_basins_df = berging_basins_df[berging_basins_df.geom_type == "Polygon"] +berging_basins_df = berging_basins_df[berging_basins_df.intersects(rws_selected_areas)] +berging_basins_df = berging_basins_df[berging_basins_df.area > 50] + +cut_lines_df = gpd.read_file(cloud.joinpath("Rijkswaterstaat", "verwerkt", "couple_user_data.gpkg"), layer="cut_lines") + +berging_basins_df = split_basins(berging_basins_df, cut_lines_df) +berging_basins_df = berging_basins_df[berging_basins_df.intersects(rws_selected_areas)] + + +rws_selected_basins_df = rws_model.basin.node.df[rws_model.basin.node.df.index.isin(rws_selected_areas_df.node_id)] + +berging_basins_df.loc[:, "node_id"] = berging_basins_df.geometry.apply( + lambda x: rws_selected_basins_df.distance(x).idxmin() +) + +berging_basins_df.to_file(cloud.joinpath("Rijkswaterstaat", "verwerkt", "bergende_basins_rws.gpkg")) + + +# %% +rws_model.update_meta_properties({"meta_categorie": "hoofdwater"}) + + +# %% toevoegen bergende gebieden +basin_area_df = gpd.read_file(cloud.joinpath("Rijkswaterstaat", "verwerkt", "bergende_basins_rws.gpkg")) +basin_area_df.set_index("node_id", inplace=True) +lhm_raster_file = cloud.joinpath("Basisgegevens", "LHM", "4.3", "input", "LHM_data.tif") +ma_raster_file = cloud.joinpath("Basisgegevens", "VanDerGaast_QH", "spafvoer1.tif") +basin_area_df = add_basin_statistics(df=basin_area_df, lhm_raster_file=lhm_raster_file, ma_raster_file=ma_raster_file) + +# %% +edge_id = rws_model.edge.df.index.max() + 1 +for row in rws_model.basin.node.df[rws_model.basin.node.df.index.isin(basin_area_df.index)].itertuples(): + # row = next(row for row in model.basin.node.df.itertuples() if row.Index == 1013) + node_id = row.Index + + if node_id in basin_area_df.index: + # basin-polygon + basin_row = basin_area_df.loc[node_id] + basin_polygon = basin_area_df.at[node_id, "geometry"] + + # add basin-node + basin_node_id = ( + rws_model.next_node_id + ) # FIXME: can be removed if issue is closed https://github.com/Deltares/Ribasim/issues/1805 + geometry = basin_to_point(basin_polygon=basin_polygon, tolerance=10) + node = Node( + node_id=basin_node_id, + meta_categorie="bergend", + geometry=geometry, + ) + + if node_id in rws_model.basin.area.df.node_id.to_list(): + polygon = rws_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=basin_polygon, + polygon=polygon, + max_level=max_level, + min_level=min_level, + lhm_raster_file=lhm_raster_file, + ) + data = [ + basin_profile, + basin.State(level=[basin_profile.df.level.min() + 0.1]), + basin.Area(geometry=[MultiPolygon([basin_polygon])]), + ] + basin_node = rws_model.basin.add(node=node, tables=data) + + # get line + line = LineString([geometry, row.geometry]) + + # add tabulated rating curve + tbr_node_id = rws_model.next_node_id + geometry = line.interpolate(0.5, normalized=True) + node = Node( + node_id=tbr_node_id, + meta_categorie="bergend", + geometry=geometry, + ) + if any(pd.isna(getattr(basin_row, i)) for i in ["ghg", "glg", "ma"]): + raise ValueError(f"No valid ghg, glg and/or ma for basin_id {node_id}") + else: + data = [get_rating_curve(row=basin_row, min_level=basin_profile.df.level.min())] + tbr_node = rws_model.tabulated_rating_curve.add(node=node, tables=data) + + # add edges + edge_id += 1 # FIXME: can be removed if issue is closed https://github.com/Deltares/Ribasim/issues/1804 + rws_model.edge.add(basin_node, tbr_node, edge_id=edge_id, meta_categorie="bergend") + edge_id += 1 + rws_model.edge.add(tbr_node, rws_model.basin[node_id], edge_id=edge_id, meta_categorie="bergend") + + else: + print(f"Geen basin-vlak voor {node_id}") + +# %% +model_path = cloud.joinpath("Rijkswaterstaat", "modellen", "hws_bergend") +toml_file = model_path / "hws.toml" + +rws_model.write(toml_file) diff --git a/notebooks/koppelen/02_koppelen.py b/notebooks/koppelen/02_koppelen.py new file mode 100644 index 0000000..08568f3 --- /dev/null +++ b/notebooks/koppelen/02_koppelen.py @@ -0,0 +1,211 @@ +# %% +from networkx import NetworkXNoPath +from ribasim_nl import CloudStorage, Model, Network, reset_index +from ribasim_nl.case_conversions import pascal_to_snake_case +from ribasim_nl.concat import concat +from shapely.geometry import LineString + +cloud = CloudStorage() + +# %% update RWS-HWS + +# RWS-HWS +model_path = cloud.joinpath("Rijkswaterstaat", "modellen", "hws_bergend") +toml_file = model_path / "hws.toml" +rws_model = Model.read(toml_file) + +# some fixes +node_id = 8413 +level = rws_model.upstream_profile(node_id).level.min() + 0.1 + +mask = (rws_model.tabulated_rating_curve.static.df.node_id == node_id) & ( + rws_model.tabulated_rating_curve.static.df.level < level +) +rws_model.tabulated_rating_curve.static.df.loc[mask, ["level"]] = level + +# reset index +rws_model = reset_index(rws_model) + +# # write model +rws_model.update_meta_properties(node_properties={"authority": "Rijkswaterstaat"}) +rws_model.write(model_path.with_name("hws_temp") / "hws.toml") + + +# %% update AGV + +# AGV +model_path = cloud.joinpath("AmstelGooienVecht", "modellen", "AmstelGooienVecht_parametrized_2024_8_47") +if not model_path.exists(): + model_url = cloud.joinurl("AmstelGooienVecht", "modellen", "AmstelGooienVecht_parametrized_2024_8_47") + cloud.download_content(model_url) +toml_file = model_path / "ribasim.toml" +# update_database(toml_file) +agv_model = Model.read(toml_file) + +# fix manning issue +agv_model.manning_resistance.static.df = agv_model.manning_resistance.static.df[ + agv_model.manning_resistance.static.df.node_id.isin(agv_model.node_table().df.index) +] + +# fix boundary-node issue +agv_model.remove_node(957, remove_edges=False) + +# reset index +agv_model = reset_index(agv_model, node_start=rws_model.next_node_id) + +# # write model +agv_model.update_meta_properties(node_properties={"authority": "AmstelGooienVecht"}) +agv_model.write(model_path.with_name("AmstelGooienVecht_temp") / "agv.toml") + +# %% update De Dommel + +# update De Dommel +model_path = cloud.joinpath("DeDommel", "modellen", "DeDommel") +toml_file = model_path / "model.toml" +# update_database(toml_file) +dommel_model = Model.read(toml_file) + +# set LevelBoundary from-to RWS +mask = dommel_model.level_boundary.static.df.node_id.isin([19, 20, 21, 22, 23, 24, 25, 27]) +dommel_model.level_boundary.static.df.loc[mask, "meta_from_authority"] = "DeDommel" +dommel_model.level_boundary.static.df.loc[mask, "meta_to_authority"] = "Rijkswaterstaat" + +mask = dommel_model.level_boundary.static.df.node_id.isin([17, 18]) +dommel_model.level_boundary.static.df.loc[mask, "meta_from_authority"] = "Rijkswaterstaat" +dommel_model.level_boundary.static.df.loc[mask, "meta_to_authority"] = "DeDommel" + +# reset index +dommel_model = reset_index(dommel_model, node_start=agv_model.next_node_id) + +# # write model +dommel_model.update_meta_properties(node_properties={"authority": "DeDommel"}) + + +dommel_model.write(model_path.with_name("DeDommel_temp") / "de_dommel.toml") + +# %% prepare coupling + +# prepare coupling +netwerk_mask_poly = agv_model.basin.area.df.union_all() + +models = [rws_model, dommel_model] +coupled_model = concat(models) + +network = Network.from_network_gpkg(cloud.joinpath("Rijkswaterstaat", "verwerkt", "netwerk.gpkg")) + +# %% coupling + +# couple boundaries +boundary_node_ids = coupled_model.level_boundary.static.df[ + (coupled_model.level_boundary.static.df.meta_to_authority == "Rijkswaterstaat") + | (coupled_model.level_boundary.static.df.meta_from_authority == "Rijkswaterstaat") +].node_id.to_list() + +mask = (coupled_model.node_table().df.meta_authority == "Rijkswaterstaat") & ( + coupled_model.node_table().df.meta_categorie == "hoofdwater" +) + +basin_ids = coupled_model.node_table().df[mask].index.to_list() +basin_areas_df = coupled_model.basin.area.df[coupled_model.basin.area.df.node_id.isin(basin_ids)].set_index("node_id") + +for boundary_node_id in boundary_node_ids: + # boundary_node_id = boundary_node_ids[0] + boundary_node = coupled_model.level_boundary[boundary_node_id] + # get upstream node to couple from + try: + # get basin-id to couple to + to_node_id = basin_areas_df.distance(boundary_node.geometry).idxmin() + to_node = coupled_model.basin[to_node_id] + listen_node_id = to_node_id + + # get to network node + to_network_node = network.nodes.distance(to_node.geometry).idxmin() + + # get node to couple from + from_node_id = coupled_model.upstream_node_id(boundary_node_id) + from_node_type = coupled_model.node_table().df.at[from_node_id, "node_type"] + from_node = getattr(coupled_model, pascal_to_snake_case(from_node_type))[from_node_id] + + # get from network node + link_idx = iter(network.links.distance(from_node.geometry).sort_values().index) + edge_geometry = None + while edge_geometry is None: + idx = next(link_idx) + try: + link_geom = network.links.at[idx, "geometry"] + if link_geom.intersects(netwerk_mask_poly): + continue + projected_point = link_geom.interpolate(link_geom.project(from_node.geometry)) + if network.nodes.distance(projected_point).min() > 10: + from_network_node = network.add_node(projected_point, max_distance=9) + else: + from_network_node = network.nodes.distance(projected_point).idxmin() + edge_geometry = network.get_line(from_network_node, to_network_node) + except NetworkXNoPath: + continue + + except KeyError: + # get basin-id to couple from + from_node_id = basin_areas_df.distance(boundary_node.geometry).idxmin() + from_node = coupled_model.basin[from_node_id] + listen_node_id = from_node_id + + # get from network node + from_network_node = network.nodes.distance(from_node.geometry).idxmin() + + # get node to couple to + to_node_id = coupled_model.downstream_node_id(boundary_node_id) + to_node_type = coupled_model.node_table().df.at[to_node_id, "node_type"] + to_node = getattr(coupled_model, pascal_to_snake_case(to_node_type))[to_node_id] + + # get edge geometry + link_idx = iter(network.links.distance(to_node.geometry).sort_values().index) + edge_geometry = None + while edge_geometry is None: + idx = next(link_idx) + try: + link_geom = network.links.at[idx, "geometry"] + if link_geom.intersects(netwerk_mask_poly): + continue + projected_point = link_geom.interpolate(link_geom.project(to_node.geometry)) + if network.nodes.distance(projected_point).min() > 10: + to_network_node = network.add_node(projected_point, max_distance=9) + else: + to_network_node = network.nodes.distance(projected_point).idxmin() + edge_geometry = network.get_line(from_network_node, to_network_node) + except NetworkXNoPath: + continue + + # remove boundary node + coupled_model.remove_node(boundary_node_id, remove_edges=True) + + # update discrete control + mask = coupled_model.discrete_control.variable.df.listen_node_id == boundary_node_id + coupled_model.discrete_control.variable.df.loc[mask, ["listen_node_id"]] = listen_node_id + + # construct edge-geometry + if edge_geometry.boundary.geoms[0].distance(from_node.geometry) > 0.001: + edge_geometry = LineString(tuple(from_node.geometry.coords) + tuple(edge_geometry.coords)) + if edge_geometry.boundary.geoms[1].distance(to_node.geometry) > 0.001: + edge_geometry = LineString(tuple(edge_geometry.coords) + tuple(to_node.geometry.coords)) + + # add edge + edge_id = coupled_model.edge.df.index.max() + 1 + coupled_model.edge.add( + edge_id=edge_id, + from_node=from_node, + to_node=to_node, + geometry=edge_geometry, + meta_from_authority="AmstelGooiEnVecht", + meta_to_authority="Rijkswaterstaat", + ) + + +# %% + +model_path = cloud.joinpath("Rijkswaterstaat", "modellen", "lhm") +toml_file = model_path / "lhm.toml" + +coupled_model.write(toml_file) + +# %% diff --git a/notebooks/koppelen/99_upload_model.py b/notebooks/koppelen/99_upload_model.py new file mode 100644 index 0000000..2d62e90 --- /dev/null +++ b/notebooks/koppelen/99_upload_model.py @@ -0,0 +1,8 @@ +# %% +from ribasim_nl import CloudStorage + +cloud = CloudStorage() + +cloud.upload_model("Rijkswaterstaat", "lhm", include_results=True, include_plots=False) + +# %% diff --git a/open-vscode.bat b/open-vscode.bat new file mode 100644 index 0000000..ae86837 --- /dev/null +++ b/open-vscode.bat @@ -0,0 +1 @@ +pixi run --environment=dev code . | exit diff --git a/src/ribasim_nl/reset_index.py b/src/ribasim_nl/reset_index.py new file mode 100644 index 0000000..cfb78b3 --- /dev/null +++ b/src/ribasim_nl/reset_index.py @@ -0,0 +1,42 @@ +# %% +import pandas as pd +from ribasim import Model +from ribasim_nl.case_conversions import pascal_to_snake_case + + +def reset_index(model: Model, node_start=1): + # only reset if we have to + node_id_min = model.node_table().df.index.min() + node_id_max = model.node_table().df.index.max() + expected_length = node_id_max - node_id_min + 1 + if not ((node_start == node_id_min) and (expected_length == len(model.node_table().df))): + # make sure column node_id == index + node_ids = model.node_table().df.index + + # create a new index for re-indexing all tables + index = pd.Series(data=[i + node_start for i in range(len(node_ids))], index=node_ids).astype("int32") + + # renumber edges + model.edge.df.loc[:, ["from_node_id"]] = model.edge.df["from_node_id"].apply(lambda x: index[x]) + + model.edge.df.loc[:, ["to_node_id"]] = model.edge.df["to_node_id"].apply(lambda x: index[x]) + + # renumber tables + for node_type in model.node_table().df.node_type.unique(): + ribasim_node = getattr(model, pascal_to_snake_case(node_type)) + for attr in ribasim_node.model_fields.keys(): + table = getattr(ribasim_node, attr) + try: + if table.df is not None: + if "node_id" in table.df.columns: + table.df.loc[:, "node_id"] = table.df["node_id"].apply(lambda x: index[x]) + table.df.index += 1 + elif table.df.index.name == "node_id": + table.df.index = table.df.reset_index("node_id")["node_id"].apply(lambda x: index.loc[x]) + # table.df.index.name = "node_id" + except KeyError as e: + raise KeyError(f"node_id {e} not in {node_type} / {attr} not in node-table") + return model + + +# %% diff --git a/src/ribasim_nl/ribasim_nl/berging.py b/src/ribasim_nl/ribasim_nl/berging.py new file mode 100644 index 0000000..fe8a3a8 --- /dev/null +++ b/src/ribasim_nl/ribasim_nl/berging.py @@ -0,0 +1,241 @@ +from pathlib import Path + +import numpy as np +import numpy.typing as npt +import pandas as pd +import rasterio +from geopandas import GeoDataFrame +from rasterio.windows import from_bounds +from rasterstats import zonal_stats +from ribasim.nodes import basin, tabulated_rating_curve +from shapely.geometry import Polygon + +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, +} + + +def sample_raster( + raster_file: Path, + df: GeoDataFrame, + band: int = 1, + fill_value: float | None = None, + all_touched: bool = False, + stats: str = "mean", + maaiveld_data: npt.ArrayLike | None = None, +): + """Sample rasters over Polygons + + Args: + raster_file (Path): Raster-file to sample + df (GeoDataFrame): GeoDataFrame with polygons + band (int, optional): Band in raster-file to sample from. Defaults to 1. + fill_value (float | None, optional): Fill-value for nodata-cells. Defaults to None. + all_touched (bool, optional): rasterize all_touched setting. Defaults to False. + stats (str, optional): rasterstats stats setting. Defaults to "mean". + maaiveld_data (npt.ArrayLike | None, optional): If numpy-array in same shape as raster, raster-data will be subtracted from it. Defaults to None. + + Returns + ------- + list[dict]: Rasterstats output + """ + with rasterio.open(raster_file) as raster_src: + # 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=stats, nodata=raster_src.nodata, all_touched=all_touched) + + +def add_basin_statistics(df: GeoDataFrame, lhm_raster_file: Path, ma_raster_file: Path) -> GeoDataFrame: + """Add Vd Gaast basin-statistics to a Polygon basin GeoDataFrame + + Args: + df (GeoDataFrame): GeoDataFrame with basins + lhm_raster_file (Path): LHM raster-file with layers + ma_raster_file (Path): Specific discharge (maatgevende afvoer) raster + + Returns + ------- + GeoDataFrame: GeoDataFrame with basins ánd statistics + """ + with rasterio.open(lhm_raster_file) as raster_src: + # read band + maaiveld_data = raster_src.read(BANDEN["maaiveld"]) + + # sample rasters + ghg = sample_raster( + raster_file=lhm_raster_file, + df=df, + band=BANDEN["ghg_2010-2019"], + all_touched=True, + maaiveld_data=maaiveld_data, + ) + df.loc[:, ["ghg"]] = pd.Series(dtype=float) + df.loc[:, ["ghg"]] = [i["mean"] for i in ghg] + + glg = sample_raster( + raster_file=lhm_raster_file, + df=df, + band=BANDEN["glg_2010-2019"], + all_touched=True, + maaiveld_data=maaiveld_data, + ) + df.loc[:, ["glg"]] = pd.Series(dtype=float) + df.loc[:, ["glg"]] = [i["mean"] for i in glg] + + ma = sample_raster(raster_file=ma_raster_file, df=df, all_touched=True, fill_value=37) # 37mm/dag is + df.loc[:, ["ma"]] = pd.Series(dtype=float) + df.loc[:, ["ma"]] = [i["mean"] for i in ma] + + maaiveld = sample_raster( + raster_file=lhm_raster_file, df=df, band=BANDEN["maaiveld"], stats="mean min max", all_touched=True + ) + df.loc[:, ["maaiveld"]] = pd.Series(dtype=float) + df.loc[:, ["maaiveld"]] = [i["mean"] if pd.isna(i["mean"]) else i["mean"] for i in maaiveld] + df.loc[:, ["maaiveld_max"]] = [i["max"] if pd.isna(i["max"]) else i["max"] for i in maaiveld] + df.loc[:, ["maaiveld_min"]] = [i["min"] if pd.isna(i["min"]) else i["min"] for i in maaiveld] + + return df + + +def get_rating_curve(row, min_level: float, maaiveld: None | float = None) -> tabulated_rating_curve.Static: + """Generate a tabulated_rating_curve.Static object from basin_statistics + + Args: + row (pd.Series): Row in a GeoDataFrame containing basin_statistics + min_level (float): minimal level in rating curve (basin.profile.level.min() of upstream basin) + maaiveld (None | float, optional): surface-level. If none, it should be in row.maaiveld. Defaults to None. + + Returns + ------- + tabulated_rating_curve.Static: Static-table for tabulated_rating_curve node + """ + flow_rate = np.round([0, row.ma * 0.2, row.ma * 0.33, row.ma / 2, row.ma * 2], 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 - depth + else: + 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 * 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, polygon: Polygon, max_level: float, min_level: float, lhm_raster_file: Path +) -> basin.Profile: + """Generate a basin.Static table for a Polygon using LHM rasters + + Args: + basin_polygon (Polygon): Polygon defining the basin + polygon (Polygon): Polygon defining all waters that are to be subtracted in primary waters + max_level (float): minimal level in basin-profile + min_level (float): maximal level in basin-profile + lhm_raster_file (Path): path to lhm-rasters + + Returns + ------- + basin.Profile: basin-profile for basin-node + """ + with rasterio.open(lhm_raster_file) 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()]) + + 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) diff --git a/src/ribasim_nl/ribasim_nl/concat.py b/src/ribasim_nl/ribasim_nl/concat.py index f08a6f9..5a9368e 100644 --- a/src/ribasim_nl/ribasim_nl/concat.py +++ b/src/ribasim_nl/ribasim_nl/concat.py @@ -1,5 +1,4 @@ import pandas as pd -import ribasim from ribasim import Model from ribasim_nl import reset_index @@ -21,24 +20,17 @@ def concat(models: list[Model]) -> Model: """ # models will be concatenated to first model. model = reset_index(models[0]) - # determine node_start of next model - node_start = model.node_table().df.node_id.max() + 1 # concat all other models into model for merge_model in models[1:]: - # reset index + # reset index of mergemodel, node_start is max node_id + node_start = model.node_table().df.index.max() + 1 merge_model = reset_index(merge_model, node_start) - # determine node_start of next model - node_start = model.node_table().df.node_id.max() + 1 - - # merge network - # model.network.node = ribasim.Node( - # df=pd.concat([model.network.node.df, merge_model.network.node.df]) - # ) - model.edge = ribasim.EdgeTable( - df=pd.concat([model.edge.df, merge_model.edge.df], ignore_index=True).reset_index(drop=True) - ) + # concat edges + edge_df = pd.concat([model.edge.df, merge_model.edge.df], ignore_index=True) + edge_df.index.name = "edge_id" + model.edge.df = edge_df # merge tables for node_type in model.node_table().df.node_type.unique(): @@ -51,7 +43,13 @@ def concat(models: list[Model]) -> Model: if merge_model_df is not None: if model_df is not None: # make sure we concat both df's into the correct ribasim-object - df = pd.concat([model_df, merge_model_df], ignore_index=True) + if "node_id" in model_df.columns: + df = pd.concat([model_df, merge_model_df], ignore_index=True) + df.index.name = "fid" + elif model_df.index.name == "node_id": + df = pd.concat([model_df, merge_model_df], ignore_index=False) + else: + raise Exception(f"{node_type} / {attr} cannot be merged") else: df = merge_model_df model_node_table.df = df diff --git a/src/ribasim_nl/ribasim_nl/model.py b/src/ribasim_nl/ribasim_nl/model.py index 4311a35..f8a4b51 100644 --- a/src/ribasim_nl/ribasim_nl/model.py +++ b/src/ribasim_nl/ribasim_nl/model.py @@ -139,7 +139,7 @@ def remove_node(self, node_id: int, remove_edges: bool = False): for attr in table.model_fields.keys(): df = getattr(table, attr).df if df is not None: - if node_id in df.columns: + if "node_id" in df.columns: getattr(table, attr).df = df[df.node_id != node_id] else: getattr(table, attr).df = df[df.index != node_id] @@ -153,11 +153,24 @@ def remove_node(self, node_id: int, remove_edges: bool = False): self.remove_edge( 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_meta_properties(self, node_properties: dict, node_types: list | None = None): + """Set properties for all, or a selection of, node-types.""" + if node_types is None: + node_types = self.node_table().df.node_type.unique() + + for node_type in node_types: + table = getattr(self, pascal_to_snake_case(node_type)) + node_df = getattr(table, "node").df + if node_df is not None: + for key, value in node_properties.items(): + if not key.startswith("meta_"): + key = f"meta_{key}" + node_df.loc[:, [key]] = value + def update_node(self, node_id, node_type, data, node_properties: dict = {}): existing_node_type = self.node_table().df.at[node_id, "node_type"] diff --git a/src/ribasim_nl/ribasim_nl/reset_index.py b/src/ribasim_nl/ribasim_nl/reset_index.py index 937dae3..8c21426 100644 --- a/src/ribasim_nl/ribasim_nl/reset_index.py +++ b/src/ribasim_nl/ribasim_nl/reset_index.py @@ -7,37 +7,39 @@ def reset_index(model: Model, node_start=1): # only reset if we have to - node_id_min = model.node_table().df.node_id.min() - node_id_max = model.node_table().df.node_id.max() + node_id_min = model.node_table().df.index.min() + node_id_max = model.node_table().df.index.max() expected_length = node_id_max - node_id_min + 1 if not ((node_start == node_id_min) and (expected_length == len(model.node_table().df))): # make sure column node_id == index - node_ids = model.node_table().df.node_id + node_ids = model.node_table().df.index # create a new index for re-indexing all tables index = pd.Series(data=[i + node_start for i in range(len(node_ids))], index=node_ids).astype("int32") - # # re-index node_id and fid - # model.network.node.df.index = model.network.node.df["fid"].apply( - # lambda x: index.loc[x] - # ) - # model.network.node.df.index.name = "fid" - # model.network.node.df.drop(columns=["fid"], inplace=True) - # model.network.node.df.loc[:, "node_id"] = model.network.node.df.index - # renumber edges - model.edge.df.loc[:, ["from_node_id"]] = model.edge.df["from_node_id"].apply(lambda x: index.loc[x]) + model.edge.df.loc[:, ["from_node_id"]] = model.edge.df["from_node_id"].apply(lambda x: index[x]) - model.edge.df.loc[:, ["to_node_id"]] = model.edge.df["to_node_id"].apply(lambda x: index.loc[x]) + model.edge.df.loc[:, ["to_node_id"]] = model.edge.df["to_node_id"].apply(lambda x: index[x]) # renumber tables for node_type in model.node_table().df.node_type.unique(): ribasim_node = getattr(model, pascal_to_snake_case(node_type)) for attr in ribasim_node.model_fields.keys(): table = getattr(ribasim_node, attr) - if table.df is not None: - table.df.loc[:, "node_id"] = table.df["node_id"].apply(lambda x: index.loc[x]) - + try: + if table.df is not None: + if "node_id" in table.df.columns: + table.df.loc[:, "node_id"] = table.df["node_id"].apply(lambda x: index[x]) + table.df.index += 1 + if "listen_node_id" in table.df.columns: + table.df.loc[:, "listen_node_id"] = table.df["listen_node_id"].apply(lambda x: index[x]) + if table.df.index.name == "node_id": + table.df.index = table.df.reset_index("node_id")["node_id"].apply(lambda x: index[x]) + + # table.df.index.name = "node_id" + except KeyError as e: + raise KeyError(f"node_id {e} not in {node_type} / {attr} not in node-table") return model