From 8d13795967e14bc6ebeb4dbc66080cdc9414834f Mon Sep 17 00:00:00 2001 From: Martijn Visser Date: Tue, 12 Nov 2024 12:18:17 +0100 Subject: [PATCH 1/2] Remove stash folder --- stash/5_model_netwerk_old.py | 621 ---------------------------------- stash/5b_upgrade_to_main.py | 69 ---- stash/6_model_sturing copy.py | 322 ------------------ 3 files changed, 1012 deletions(-) delete mode 100644 stash/5_model_netwerk_old.py delete mode 100644 stash/5b_upgrade_to_main.py delete mode 100644 stash/6_model_sturing copy.py diff --git a/stash/5_model_netwerk_old.py b/stash/5_model_netwerk_old.py deleted file mode 100644 index 2b77342..0000000 --- a/stash/5_model_netwerk_old.py +++ /dev/null @@ -1,621 +0,0 @@ -# %% init -import geopandas as gpd -import pandas as pd -import ribasim - -from ribasim_nl import CloudStorage, Network, reset_index -from ribasim_nl.rating_curve import read_rating_curve -from ribasim_nl.verdeelsleutels import ( - read_verdeelsleutel, - verdeelsleutel_to_fractions, -) - -cloud = CloudStorage() -node_list = [] -edge_list = [] - -boundaries_passed = [] -PRECIPITATION = 0.005 / 86400 # m/s -EVAPORATION = 0.001 / 86400 # m/s - -network = Network.from_network_gpkg(cloud.joinpath("Rijkswaterstaat", "verwerkt", "netwerk.gpkg")) -boundary_gdf = gpd.read_file( - cloud.joinpath("Rijkswaterstaat", "verwerkt", "model_user_data.gpkg"), - engine="pyogrio", - layer="boundary", - fid_as_index=True, -) - -structures_gdf = gpd.read_file( - cloud.joinpath("Rijkswaterstaat", "verwerkt", "hydamo.gpkg"), - layer="kunstwerken", - engine="pyogrio", - fid_as_index=True, -) - -basin_poly_gdf = gpd.read_file( - cloud.joinpath("Rijkswaterstaat", "verwerkt", "basins.gpkg"), - layer="ribasim_basins", - engine="pyogrio", - fid_as_index=True, -) - -structures_df = pd.read_excel(cloud.joinpath("Rijkswaterstaat", "verwerkt", "kunstwerk_complexen.xlsx")) -structures_df = structures_df[structures_df.in_model] -structures_df["code"] = structures_df["code"].astype(str) - -level_area_df = pd.read_csv(cloud.joinpath("Rijkswaterstaat", "verwerkt", "basins_level_area.csv")) - -rating_curves_gdf = gpd.read_file( - cloud.joinpath("Rijkswaterstaat", "verwerkt", "rating_curves.gpkg"), - engine="pyogrio", - fid_as_index=True, -) - -verdeelsleutel_gdf = gpd.read_file( - cloud.joinpath("Rijkswaterstaat", "verwerkt", "verdeelsleutel.gpkg"), - engine="pyogrio", - fid_as_index=True, -) - -verdeelsleutel_df = read_verdeelsleutel(cloud.joinpath("Rijkswaterstaat", "verwerkt", "verdeelsleutel_driel.xlsx")) - - -# %% inlezen netwerk nodes en links -# define nodes and links -nodes_gdf = network.nodes -links_gdf = network.links -network_union_lines = links_gdf.unary_union -network.overlay(basin_poly_gdf[["basin_id", "geometry"]]) - -basin_admin: dict[str, list] = {i: {"nodes_from": [], "nodes_to": [], "neighbors": []} for i in basin_poly_gdf.index} - -# %% toevoegen reating-curves aan basin-admin -# toevoegen reating-curves aan basin-admin -for row in rating_curves_gdf.itertuples(): - basin_id = basin_poly_gdf[basin_poly_gdf.contains(row.geometry)].index[0] - basin_admin[basin_id]["rating_curve"] = row.curve_id - -# %% definieren boundaries -# definieren boundaries -boundary_gdf["node_id"] = boundary_gdf["geometry"].apply( - lambda x: network.move_node( - x, - max_distance=10, - align_distance=10, - node_types=["downstream_boundary", "connection", "upstream_boundary"], - ) -) - - -def get_structure_codes(structures_gdf, complex_codes, line_string): - """Return list of nearest structure_codes from a list of complex_codes and a linestring.""" - structure_codes = [] - for complex_code, gdf in structures_gdf[structures_gdf.complex_code.isin(complex_codes)].groupby("complex_code"): - gdf_select = gdf[ - gdf.soort.isin( - [ - "Stuwen", - "Spuisluizen", - "Stormvloedkeringen", - "keersluizen", - "Gemalen", - ] - ) - ] - if gdf_select.empty: - gdf_select = gdf[gdf.soort == "Schutsluizen"] - if gdf_select.empty: - raise Exception(f"kan geen kunstwerk vinden voor {complex_code}") - structure_codes += [gdf_select.at[gdf_select.distance(line_string).sort_values().index[0], "code"]] - return structure_codes - - -def get_type_and_value(code, structures_gdf, structures_df): - if code not in structures_df.code.to_numpy(): - complex_code = structures_gdf.set_index("code").loc[code].complex_code - row = structures_df.set_index("code").loc[complex_code] - else: - row = structures_df.set_index("code").loc[code] - return row.ribasim_type, row.ribasim_waarde - - -# %% kunstwerken toekennen aan basins -# kunstwerken toekennen aan basins -# maak lijst met kunstwerk-codes -structure_codes = structures_df[structures_df.code.isin(structures_gdf.code)]["code"].to_list() - -complex_codes = structures_df[~structures_df.code.isin(structures_gdf.code)]["code"].to_list() -structure_codes += get_structure_codes( - structures_gdf, - complex_codes=complex_codes, - line_string=network_union_lines, -) - - -# % itereer over kunstwerken -kwk_select_gdf = structures_gdf[structures_gdf.code.isin(structure_codes)] -condition = kwk_select_gdf.complex_code.duplicated(keep=False) -condition = condition | kwk_select_gdf.complex_naam.isna() -kwk_select_gdf.loc[condition, ["model_naam"]] = kwk_select_gdf[condition].naam -kwk_select_gdf.loc[~condition, ["model_naam"]] = kwk_select_gdf[~condition].complex_naam - -for kwk in kwk_select_gdf.itertuples(): - # get network_node_id - node_id = network.move_node( - kwk.geometry, - max_distance=150, - align_distance=100, - node_types=["connection", "upstream_boundary"], - ) - if node_id is None: - node_id = network.add_node(kwk.geometry, max_distance=150) - - # add to node_list - node_type, node_value = get_type_and_value(kwk.code, structures_gdf, structures_df) - - if network.has_upstream_nodes(node_id) and network.has_downstream_nodes( - node_id - ): # normal case, structure is between basins - upstream, downstream = network.get_upstream_downstream(node_id, "basin_id") - if upstream is not None: - basin_admin[upstream]["nodes_to"] += [node_id] - if downstream is not None: - basin_admin[upstream]["neighbors"] += [downstream] - else: - raise ValueError(f"Kan geen boven en benedenstroomse knoop vinden voor {kwk.naam}") - if downstream is not None: - basin_admin[downstream]["nodes_from"] += [node_id] - if upstream is not None: - basin_admin[downstream]["neighbors"] += [upstream] - else: - print(f"{kwk.naam} is een outlet") - boundary = boundary_gdf.loc[ - boundary_gdf[boundary_gdf["type"] == "LevelBoundary"].distance(kwk.geometry).sort_values().index[0] - ] - boundaries_passed += [boundary.name] - node_list += [ - { - "type": "LevelBoundary", - "value": boundary.level, - "is_structure": False, - "code_waterbeheerder": boundary.code, - "waterbeheerder": "Rijkswaterstaat", - "name": boundary.naam, - "node_id": boundary.node_id, - "geometry": network.nodes.at[boundary.node_id, "geometry"], - } - ] - print(f"Processing node_id: {node_id}") - edge_list += [ - { - "from_node_id": node_id, - "to_node_id": boundary.node_id, - "edge_type": "flow", - "geometry": network.get_line(node_id, boundary.node_id), - } - ] - - elif network.has_downstream_nodes(node_id): - first, second = network.get_downstream(node_id, "basin_id") - if first is not None: - basin_admin[first]["nodes_from"] += [node_id] - if second is not None: - basin_admin[first]["neighbors"] += [second] - else: - raise ValueError(f"{kwk.naam} heeft problemen") - if second is not None: - basin_admin[second]["nodes_from"] += [node_id] - if first is not None: - basin_admin[second]["neighbors"] += [first] - else: - raise ValueError(f"{kwk.naam} heeft problemen") - - node_list += [ - { - "type": node_type, - "value": node_value, - "is_structure": True, - "code_waterbeheerder": kwk.code, - "waterbeheerder": "Rijkswaterstaat", - "name": kwk.model_naam, - "node_id": node_id, - "geometry": network.nodes.at[node_id, "geometry"], - } - ] - kwk_select_gdf.loc[kwk.Index, ["node_id", "upstream", "downstream"]] = ( - node_id, - upstream, - downstream, - ) -# ignore edges with "circular structures combinations" -ignore_links = [] -for kwk in kwk_select_gdf.itertuples(): - condition = kwk_select_gdf["upstream"] == kwk.downstream - condition = condition & (kwk_select_gdf["downstream"] == kwk.upstream) - df = kwk_select_gdf[condition] - for row in df.itertuples(): - ignore_links += [(row.node_id, kwk.node_id)] - -# %% netwerk opbouwen per basin -# netwerk opbouwen per basin -for basin in basin_poly_gdf.itertuples(): - # for basin in basin_poly_gdf.loc[[55, 70]].itertuples(): - print(f"{basin.Index} {basin.naam}") - # get upstream an downstream basin boundary nodes - boundary_nodes = network.nodes[network.nodes.distance(basin.geometry.boundary) < 0.01] - us_ds = [network.get_upstream_downstream(i, "basin_id", max_iters=3, max_length=1000) for i in boundary_nodes.index] - boundary_nodes["upstream"] = [i[0] for i in us_ds] - boundary_nodes["downstream"] = [i[1] for i in us_ds] - boundary_nodes = boundary_nodes.dropna(subset=["upstream", "downstream"], how="any") - boundary_nodes.loc[:, ["count"]] = [ - len(network.upstream_nodes(i)) + len(network.downstream_nodes(i)) for i in boundary_nodes.index - ] - - boundary_nodes.sort_values(by="count", ascending=False, inplace=True) - boundary_nodes.loc[:, ["sum"]] = boundary_nodes["upstream"] + boundary_nodes["downstream"] - boundary_nodes.drop_duplicates("sum", inplace=True) - connected = basin_admin[basin.Index]["neighbors"] - boundary_nodes = boundary_nodes[~boundary_nodes["upstream"].isin(connected)] - boundary_nodes = boundary_nodes[~boundary_nodes["downstream"].isin(connected)] - - # add boundary nodes - for node in boundary_nodes.itertuples(): - neighbor = next(i for i in [node.upstream, node.downstream] if not i == basin.Index) - if basin.Index not in basin_admin[neighbor]["neighbors"]: - ds_node = node.upstream == basin.Index - node_type = "ManningResistance" - - # in case basin has rating curve we use FractionalFlow - if ds_node and ("rating_curve" in basin_admin[basin.Index].keys()): - node_type = "FractionalFlow" - # in us_case basin has rating curve we use FractionalFlow - elif (not ds_node) and ("rating_curve" in basin_admin[neighbor].keys()): - node_type = "FractionalFlow" - - if node_type == "FractionalFlow": - code = verdeelsleutel_gdf.at[ - verdeelsleutel_gdf.distance(network.nodes.at[node.Index, "geometry"]).sort_values().index[0], - "fractie", - ] - else: - code = None - - node_list += [ - { - "type": node_type, - "is_structure": False, - "node_id": node.Index, - "code_waterbeheerder": code, - "geometry": network.nodes.at[node.Index, "geometry"], - } - ] - if node.upstream != basin.Index: - basin_admin[basin.Index]["nodes_from"] += [node.Index] - basin_admin[basin.Index]["neighbors"] += [node.upstream] - elif node.downstream != basin.Index: - basin_admin[basin.Index]["nodes_to"] += [node.Index] - basin_admin[basin.Index]["neighbors"] += [node.downstream] - else: - raise Exception(f"iets gaat fout bij {basin.Index}") - - nodes_from = basin_admin[basin.Index]["nodes_from"] - nodes_to = basin_admin[basin.Index]["nodes_to"] - - # add basin-node - if (len(nodes_from) > 0) and (len(nodes_to) > 0): - gdf = network.subset_nodes( - nodes_from, - nodes_to, - duplicated_nodes=True, - directed=True, - inclusive=False, - ignore_links=ignore_links, - ) - if gdf.empty: - gdf = network.subset_nodes( - nodes_from, - nodes_to, - duplicated_nodes=True, - directed=False, - inclusive=False, - ignore_links=ignore_links, - ) - elif (len(nodes_from) > 0) or (len(nodes_to) > 0): - gdf = network.nodes[network.nodes.basin_id == basin.Index] - - gdf = gdf[gdf["basin_id"] == basin.Index] - - basin_node_id = gdf.distance(basin.geometry.centroid).sort_values().index[0] - basin_admin[basin.Index]["node_id"] = basin_node_id - node_list += [ - { - "type": "Basin", - "is_structure": False, - "waterbeheerder": "Rijkswaterstaat", - "name": basin.naam, - "node_id": basin_node_id, - "basin_id": basin.Index, - "geometry": network.nodes.at[basin_node_id, "geometry"], - } - ] - - # add rating-curve and extra edge - if "rating_curve" in basin_admin[basin.Index].keys(): - rc_node_id = next(i for i in network.downstream_nodes(basin_node_id) if i in gdf.index) - node_list += [ - { - "type": "TabulatedRatingCurve", - "is_structure": False, - "waterbeheerder": "Rijkswaterstaat", - "code_waterbeheerder": basin_admin[basin.Index]["rating_curve"], - "name": basin_admin[basin.Index]["rating_curve"], - "node_id": rc_node_id, - "basin_id": basin.Index, - "geometry": network.nodes.at[rc_node_id, "geometry"], - } - ] - edge_list += [ - { - "from_node_id": basin_node_id, - "to_node_id": rc_node_id, - "name": basin.naam, - "edge_type": "flow", - "geometry": network.get_line(basin_node_id, rc_node_id, directed=True), - } - ] - from_id = rc_node_id - else: - from_id = basin_node_id - - # add edge to basin - network.add_weight("basin_id", basin.basin_id) - for node_from in nodes_from: - edge_list += [ - { - "from_node_id": node_from, - "to_node_id": basin_node_id, - "name": basin.naam, - "edge_type": "flow", - "geometry": network.get_line(node_from, basin_node_id, directed=False, weight="weight"), - } - ] - - # add from basin to edge - for node_to in nodes_to: - edge_list += [ - { - "from_node_id": from_id, - "to_node_id": node_to, - "name": basin.naam, - "edge_type": "flow", - "geometry": network.get_line(from_id, node_to, directed=False, weight="weight"), - } - ] -# %% afronden boundaries -# afronden boundaries - -for boundary in boundary_gdf[~boundary_gdf.index.isin(boundaries_passed)].itertuples(): - if boundary.type == "FlowBoundary": - basin_id = network.find_downstream(boundary.node_id, "basin_id") - value = boundary.flow_rate - edge_list += [ - { - "from_node_id": boundary.node_id, - "to_node_id": basin_admin[basin_id]["node_id"], - "name": basin_poly_gdf.loc[basin_id].naam, - "edge_type": "flow", - "geometry": network.get_line(boundary.node_id, basin_admin[basin_id]["node_id"], directed=True), - } - ] - basin_admin[basin_id]["nodes_from"] += [boundary.node_id] - else: - basin_id = network.find_upstream(boundary.node_id, "basin_id", max_iters=50) - basin_node_id = basin_admin[basin_id]["node_id"] - - gdf = network.nodes[network.nodes["basin_id"] == basin_id] - manning_node_id = ( - gdf.distance(network.get_line(basin_node_id, boundary.node_id).interpolate(0.5, normalized=True)) - .sort_values() - .index[0] - ) # get manning node, closest to half-way on the line between basin and boundary - - value = boundary.level - node_list += [ - { - "type": "ManningResistance", - "is_structure": False, - "node_id": manning_node_id, - "geometry": network.nodes.at[manning_node_id, "geometry"], - } - ] - edge_list += [ - { - "from_node_id": basin_node_id, - "to_node_id": manning_node_id, - "name": basin_poly_gdf.loc[basin_id].naam, - "edge_type": "flow", - "geometry": network.get_line(basin_node_id, manning_node_id, directed=True), - }, - { - "from_node_id": manning_node_id, - "to_node_id": boundary.node_id, - "name": basin_poly_gdf.loc[basin_id].naam, - "edge_type": "flow", - "geometry": network.get_line(manning_node_id, boundary.node_id, directed=True), - }, - ] - basin_admin[basin_id]["nodes_to"] += [boundary.node_id] - - boundaries_passed += [boundary.Index] - node_list += [ - { - "type": boundary.type, - "value": value, - "is_structure": False, - "code_waterbeheerder": boundary.code, - "waterbeheerder": "Rijkswaterstaat", - "name": boundary.naam, - "node_id": boundary.node_id, - "geometry": network.nodes.at[boundary.node_id, "geometry"], - } - ] - -# %% opslaan interim-netwerk -# opslaan interim-netwerk -gpd.GeoDataFrame(node_list, crs=28992).to_file( - cloud.joinpath("Rijkswaterstaat", "verwerkt", "ribasim_intermediate.gpkg"), - layer="Node", -) -gpd.GeoDataFrame(edge_list, crs=28992).to_file( - cloud.joinpath("Rijkswaterstaat", "verwerkt", "ribasim_intermediate.gpkg"), - layer="Edge", -) # %% - -# %% conversie naar ribasim -# conversie naar ribasim -node_df = gpd.GeoDataFrame(node_list, crs=28992).drop_duplicates() -node_df.set_index("node_id", drop=False, inplace=True) -node_df.index.name = "fid" -node = ribasim.Node( - df=node_df[["node_id", "name", "type", "waterbeheerder", "code_waterbeheerder", "geometry"]].rename( - columns={ - "type": "node_type", - "code_waterbeheerder": "meta_code_waterbeheerder", - "waterbeheerder": "meta_waterbeheerder", - } - ) -) - -edge_df = gpd.GeoDataFrame(edge_list, crs=28992) -edge_df.loc[:, "from_node_type"] = edge_df.from_node_id.apply(lambda x: node_df.at[x, "type"]) - -edge_df.loc[:, "to_node_type"] = edge_df.to_node_id.apply(lambda x: node_df.at[x, "type"]) - -edge = ribasim.Edge(df=edge_df) -network = ribasim.Network(node=node, edge=edge) - -# define Basin -static_df = node_df[node_df["type"] == "Basin"][["node_id", "basin_id"]].set_index("basin_id") - -# Area -area_df = node_df[node_df["type"] == "Basin"][["basin_id", "node_id"]] -area_df.loc[:, "geometry"] = area_df.basin_id.apply(lambda x: basin_poly_gdf.at[x, "geometry"]) -area_df = gpd.GeoDataFrame(area_df[["node_id", "geometry"]], crs=28992) - -level_area_df.drop_duplicates(["level", "area", "id"], inplace=True) -profile_df = level_area_df[level_area_df["id"].isin(static_df.index)] -profile_df["node_id"] = profile_df["id"].apply(lambda x: static_df.at[x, "node_id"]) -profile_df = profile_df[["node_id", "area", "level"]] - -static_df["precipitation"] = PRECIPITATION -static_df["potential_evaporation"] = EVAPORATION -static_df["drainage"] = 0 -static_df["infiltration"] = 0 -static_df["urban_runoff"] = 0 - -state_df = profile_df.groupby("node_id").min()["level"].reset_index() -state_df.loc[:, ["level"]] = state_df["level"].apply(lambda x: max(x + 1, 0)) - -basin = ribasim.Basin(profile=profile_df, static=static_df, state=state_df, area=area_df) - -resistance_df = node_df[node_df["type"] == "LinearResistance"][["node_id", "value"]] -resistance_df.rename(columns={"value": "resistance"}, inplace=True) -linear_resistance = ribasim.LinearResistance(static=resistance_df) - -# % define Pump -pump_df = node_df[node_df["type"] == "Pump"][["node_id", "value"]] -pump_df.rename(columns={"value": "flow_rate"}, inplace=True) -pump = ribasim.Pump(static=pump_df) - - -# % define Outlet -outlet_df = node_df[node_df["type"] == "Outlet"][["node_id", "value"]] -outlet_df.rename(columns={"value": "flow_rate"}, inplace=True) -outlet = ribasim.Outlet(static=outlet_df) - -# % define fraction -node_index = node_df[node_df["type"] == "FractionalFlow"][["code_waterbeheerder", "node_id"]].set_index( - "code_waterbeheerder" -)["node_id"] - -fractional_flow_df = verdeelsleutel_to_fractions(verdeelsleutel_df, node_index) -fractional_flow = ribasim.FractionalFlow(static=fractional_flow_df) - -# % -node_index = node_df[node_df["type"] == "TabulatedRatingCurve"][["code_waterbeheerder", "node_id"]].set_index( - "code_waterbeheerder" -)["node_id"] -tabulated_rating_curve_df = read_rating_curve( - cloud.joinpath("Rijkswaterstaat", "verwerkt", "rating_curves.xlsx"), - node_index, -) -tabulated_rating_curve = ribasim.TabulatedRatingCurve( - static=tabulated_rating_curve_df.rename(columns={"code_waterbeheerder": "meta_code_waterbeheerder"}) -) - -# % define Manning -manning_df = node_df[node_df["type"] == "ManningResistance"][["node_id"]] -manning_df["length"] = 10000 -manning_df["manning_n"] = 0.04 -manning_df["profile_width"] = 10000 -manning_df["profile_slope"] = 1 - -manning_resistance = ribasim.ManningResistance(static=manning_df) - -# % define FlowBoundary -flow_boundary_df = node_df[node_df["type"] == "FlowBoundary"][["node_id", "value"]] -flow_boundary_df.rename(columns={"value": "flow_rate"}, inplace=True) -flow_boundary = ribasim.FlowBoundary(static=flow_boundary_df) - -# % define LevelBoundary -level_boundary_df = node_df[node_df["type"] == "LevelBoundary"][["node_id", "value"]] -level_boundary_df.rename(columns={"value": "level"}, inplace=True) -level_boundary = ribasim.LevelBoundary(static=level_boundary_df) - - -# % write model -model = ribasim.Model( - network=network, - basin=basin, - flow_boundary=flow_boundary, - level_boundary=level_boundary, - linear_resistance=linear_resistance, - manning_resistance=manning_resistance, - tabulated_rating_curve=tabulated_rating_curve, - fractional_flow=fractional_flow, - pump=pump, - outlet=outlet, - # terminal=terminal, - starttime="2020-01-01 00:00:00", - endtime="2021-01-01 00:00:00", -) - -# -model = reset_index(model) - -# Nijkerkersluis naar inactive -nijkerk_idx = model.network.node.df[model.network.node.df["meta_code_waterbeheerder"] == "32E-001-04"].index - - -model.linear_resistance.static.df.loc[model.linear_resistance.static.df["node_id"].isin(nijkerk_idx), ["active"]] = ( - False -) - -model.linear_resistance.static.df.loc[~model.linear_resistance.static.df["node_id"].isin(nijkerk_idx), ["active"]] = ( - True -) - - -# -model.solver.algorithm = "RK4" -model.solver.dt = 10.0 -model.solver.saveat = 3600 - -# %% wegschrijven model -# wegschrijven model -print("write ribasim model") -ribasim_toml = cloud.joinpath("Rijkswaterstaat", "modellen", "hws_network", "hws.toml") -model.write(ribasim_toml) - -# %% diff --git a/stash/5b_upgrade_to_main.py b/stash/5b_upgrade_to_main.py deleted file mode 100644 index b5a43c5..0000000 --- a/stash/5b_upgrade_to_main.py +++ /dev/null @@ -1,69 +0,0 @@ -# %% -import sqlite3 - -import pandas as pd -import ribasim -import tomli -import tomli_w -from ribasim import Node -from ribasim.nodes import linear_resistance - -from ribasim_nl import CloudStorage - - -def upgrade_feature(database_path): - conn = sqlite3.connect(database_path) - - table = "DiscreteControl / condition" - df = pd.read_sql_query(f"SELECT * FROM '{table}'", conn) - df_renamed = df.rename( - columns={ - "listen_feature_type": "listen_node_type", - "listen_feature_id": "listen_node_id", - } - ) - df_renamed.to_sql(table, conn, if_exists="replace", index=False) - conn.close() - - -def upgrade_crs(toml_path): - with open(toml_path, "rb") as f: - config = tomli.load(f) - - config["crs"] = "EPSG:28992" - - with open(toml_path, "wb") as f: - tomli_w.dump(config, f) - - -cloud = CloudStorage() - -toml_file = cloud.joinpath("rijkswaterstaat", "modellen", "hws_network_upgraded", "hws.toml") -upgrade_crs(toml_file) -# upgrade_feature(toml_file.with_name("database").with_suffix(".gpkg")) - -# %% -model = ribasim.Model.read(toml_file) - -mask = (model.edge.df.to_node_type == "LevelBoundary") & (model.edge.df.from_node_type == "ManningResistance") - -node_ids = model.edge.df[mask].from_node_id.to_list() - - -model.edge.df.loc[model.edge.df.from_node_id.isin(node_ids), ["from_node_type"]] = "LinearResistance" -model.edge.df.loc[model.edge.df.to_node_id.isin(node_ids), ["to_node_type"]] = "LinearResistance" - -# %% change - -nodes = model.node_table().df[model.node_table().df["node_id"].isin(node_ids)] - - -for attr in model.manning_resistance.__fields__.keys(): - table = getattr(model.manning_resistance, attr) - table.df = table.df[~table.df["node_id"].isin(node_ids)] - - -for _, node in nodes.iterrows(): - model.linear_resistance.add(Node(**node.to_dict()), [linear_resistance.Static(resistance=[0.005])]) - -model.write(toml_file) diff --git a/stash/6_model_sturing copy.py b/stash/6_model_sturing copy.py deleted file mode 100644 index 0edbfb0..0000000 --- a/stash/6_model_sturing copy.py +++ /dev/null @@ -1,322 +0,0 @@ -# %% -from datetime import datetime - -import pandas as pd -import ribasim -from ribasim import nodes - -from ribasim_nl import CloudStorage, discrete_control -from ribasim_nl.case_conversions import pascal_to_snake_case -from ribasim_nl.model import add_control_node_to_network -from ribasim_nl.verdeelsleutels import read_verdeelsleutel, verdeelsleutel_to_control - -cloud = CloudStorage() -waterbeheerder = "Rijkswaterstaat" - -ribasim_toml = cloud.joinpath("Rijkswaterstaat", "modellen", "hws_network_upgraded", "hws.toml") -model = ribasim.Model.read(ribasim_toml) - -verdeelsleutel_df = read_verdeelsleutel(cloud.joinpath("Rijkswaterstaat", "verwerkt", "verdeelsleutel_driel.xlsx")) - - -def add_pid(node_id, node_type, ctrl_node, pid): - kwargs = {k: [v] for k, v in pid.items() if k != "listen_node_type"} - model.pid_control.add( - ctrl_node, - [nodes.pid_control.Static(listen_node_type=pid["listen_node_type"], **kwargs)], - ) - - model.edge.add( - model.pid_control[ctrl_node.node_id], - getattr(model, pascal_to_snake_case(node_type))[node_id], - ) - - -# %% toevoegen Driel -# Driel toevoegen met een control op Fractional Nodes -name = "Driel" -code_waterbeheerder = "40A-004-02" - -offset_node_id = ( - model.node_table().df[model.node_table().df["meta_code_waterbeheerder"] == code_waterbeheerder] -).node_id.to_numpy()[0] - -model = verdeelsleutel_to_control( - verdeelsleutel_df, - model, - code_waterbeheerder=code_waterbeheerder, - offset_node_id=offset_node_id, - waterbeheerder="Rijkswaterstaat", -) - -# %% add Haringvliet -# Haringvliet toevoegen met control op outlet -name = "Haringvlietsluizen" -code_waterbeheerder = "37C-350-09" -flow_rate_haringvliet = [0, 0, 50, 50, 400, 2500, 3800, 5200, 6800, 8000, 9000] -flow_rate_lobith = [0, 1100, 1100, 1700, 2000, 4000, 6000, 8000, 10000, 12000, 15000] - -node_ids = ( - model.node_table().df[model.node_table().df["meta_code_waterbeheerder"] == code_waterbeheerder].node_id.to_numpy() -) - - -ctrl_node = add_control_node_to_network( - model, - node_ids, - ctrl_node_geom=(62430, 427430), - meta_waterbeheerder="Rijkswaterstaat", - name=name, - meta_code_waterbeheerder=name, -) - -listen_node_id = ( - model.node_table().df[model.node_table().df["meta_code_waterbeheerder"] == "LOBH"].node_id.to_numpy()[0] -) - -listen_node_type = model.node_table().df.set_index("node_id").at[listen_node_id, "node_type"] - - -condition_df = discrete_control.condition( - values=flow_rate_lobith, - node_id=ctrl_node.node_id, - listen_feature_id=listen_node_id, - name=name, -) - -logic_df = discrete_control.logic( - node_id=ctrl_node.node_id, - length=len(flow_rate_haringvliet), - name=name, -) - -model.discrete_control.add( - ctrl_node, - [ - nodes.discrete_control.Condition( - listen_node_id=listen_node_id, - listen_node_type=listen_node_type, - variable="flow_rate", - greater_than=condition_df["greater_than"].to_list(), - ), - nodes.discrete_control.Logic( - truth_state=logic_df.truth_state.to_list(), - control_state=logic_df.control_state.to_list(), - ), - ], -) - -model.edge.add( - model.discrete_control[ctrl_node.node_id], - model.outlet[node_ids[0]], -) - - -outlet_df = discrete_control.node_table(values=flow_rate_lobith, variable="flow_rate", name=name, node_id=node_ids[0]) -for col, dtype in model.outlet.static.df.dtypes.to_dict().items(): - if col not in outlet_df.columns: - outlet_df.loc[:, [col]] = None - outlet_df[col] = outlet_df[col].astype(dtype) - - -model.outlet.static.df = model.outlet.static.df[model.outlet.static.df.node_id != node_ids[0]] -model.outlet.static.df = pd.concat([model.outlet.static.df, outlet_df]) - -# %% toevoegen IJsselmeer - -# toevoegen peil IJsselmeer door opgave peil waddenzee -node_ids = ( - model.node_table().df[model.node_table().df["meta_code_waterbeheerder"].isin(["KOBU", "OEBU"])].node_id.to_numpy() -) - -time = pd.date_range(model.starttime, model.endtime) - -day_of_year = [ - "01-01", - "03-01", - "03-11", - "03-21", - "04-01", - "04-10", - "04-15", - "08-11", - "08-15", - "08-21", - "08-31", - "09-11", - "09-15", - "09-21", - "10-01", - "10-11", - "10-21", - "10-31", - "12-31", -] - -level = [ - -0.4, - -0.1, - -0.1, - -0.1, - -0.2, - -0.2, - -0.2, - -0.2, - -0.3, - -0.3, - -0.3, - -0.3, - -0.3, - -0.3, - -0.4, - -0.4, - -0.4, - -0.4, - -0.4, -] -level_cycle_df = pd.DataFrame( - { - "dayofyear": [datetime.strptime(i, "%m-%d").timetuple().tm_yday for i in day_of_year], - "level": level, - } -).set_index("dayofyear") - - -def get_level(timestamp, level_cycle_df): - return level_cycle_df.at[level_cycle_df.index[level_cycle_df.index <= timestamp.dayofyear].max(), "level"] - - -time = pd.date_range(model.starttime, model.endtime) -level_df = pd.DataFrame({"time": time, "level": [get_level(i, level_cycle_df) for i in time]}) - -level_df = pd.concat( - [pd.concat([level_df, pd.DataFrame({"node_id": [node_id] * len(level_df)})], axis=1) for node_id in node_ids], - ignore_index=True, -) -model.level_boundary.time.df = level_df -model.level_boundary.static.df = model.level_boundary.static.df[~model.level_boundary.static.df.node_id.isin(node_ids)] - -# %% toevoegen Irenesluizen - -# PID control op ARK volgens peilbesluit op -0.4m NAP: https://www.helpdeskwater.nl/publish/pages/138359/41_peilbesluit_boezem_noordzeekanaal_amsterdam-rijnkanaal.pdf -name = "Irenesluizen" -code_waterbeheerder = "39B-002-02" - -node_id, node_type = ( - model.node_table() - .df[model.node_table().df["meta_code_waterbeheerder"] == code_waterbeheerder] - .iloc[0][["node_id", "node_type"]] -) - -condition = model.outlet.static.df.node_id == node_id -model.outlet.static.df.loc[condition, ["min_flow_rate", "max_flow_rate"]] = 0, 60 - -gdf = model.edge.df[model.edge.df.from_node_id == node_id] -gdf = gdf[gdf.to_node_type == "Basin"] -listen_node_id = gdf.iloc[0].to_node_id - - -ctrl_node = add_control_node_to_network( - model, - [node_id], - meta_waterbeheerder="Rijkswaterstaat", - name=name, - meta_code_waterbeheerder=name, - ctrl_type="PidControl", -) - -pid = { - "listen_node_type": "Basin", - "listen_node_id": listen_node_id, - "target": -0.4, - "proportional": -50000, - "integral": -1e-7, - "derivative": 0, -} - -add_pid(node_id, node_type, ctrl_node, pid) - -# %% toevoegen Gemaal Panheel - -# PID control op Kanaal Wessem-Nederweerd volgens peilbesluit op 28.65m NAP: https://www.helpdeskwater.nl/publish/pages/188982/watersystemen-midden-limburg-en-noord-brabantse-kanalen.pdf -name = "Gemaal Panheel" -code_waterbeheerder = "58C-001-06" - -node_id, node_type = ( - model.node_table() - .df[model.node_table().df["meta_code_waterbeheerder"] == code_waterbeheerder] - .iloc[0][["node_id", "node_type"]] -) - -condition = model.pump.static.df.node_id == node_id -model.pump.static.df.loc[condition, ["min_flow_rate", "max_flow_rate"]] = 0, 3 - -gdf = model.edge.df[model.edge.df.from_node_id == node_id] -gdf = gdf[gdf.to_node_type == "Basin"] -listen_node_id = gdf.iloc[0].to_node_id - - -ctrl_node = add_control_node_to_network( - model, - [node_id], - meta_waterbeheerder="Rijkswaterstaat", - name=name, - meta_code_waterbeheerder=name, - ctrl_type="PidControl", - offset=20, -) - -pid = { - "listen_node_type": "Basin", - "listen_node_id": listen_node_id, - "target": 28.60, - "proportional": -50000, - "integral": -1e-7, - "derivative": 0, -} - -add_pid(node_id, node_type, ctrl_node, pid) - -name = "Sluis Panheel" -code_waterbeheerder = "58C-001-04" - -node_id, node_type = ( - model.node_table() - .df[model.node_table().df["meta_code_waterbeheerder"] == code_waterbeheerder] - .iloc[0][["node_id", "node_type"]] -) - -condition = model.outlet.static.df.node_id == node_id -model.outlet.static.df.loc[condition, ["min_flow_rate", "max_flow_rate"]] = 2.5, 10 - -gdf = model.edge.df[model.edge.df.to_node_id == node_id] -gdf = gdf[gdf.from_node_type == "Basin"] -listen_node_id = gdf.iloc[0].from_node_id - - -ctrl_node = add_control_node_to_network( - model, - [node_id], - meta_waterbeheerder="Rijkswaterstaat", - name=name, - meta_code_waterbeheerder=name, - ctrl_type="PidControl", - offset=20, -) -pid = { - "listen_node_type": "Basin", - "listen_node_id": listen_node_id, - "target": 28.70, - "proportional": -50000, - "integral": -1e-7, - "derivative": 0, -} - -add_pid(node_id, node_type, ctrl_node, pid) - - -# %% wegschrijven model -ribasim_toml = cloud.joinpath("Rijkswaterstaat", "modellen", "hws_sturing_upgraded", "hws.toml") -model.write(ribasim_toml) - -# %% From 4aa3ee10e56b3a3264a04ecaff632f60f2eb20f1 Mon Sep 17 00:00:00 2001 From: Martijn Visser Date: Tue, 12 Nov 2024 12:24:33 +0100 Subject: [PATCH 2/2] Convert remaining notebooks to scripts --- notebooks/ijsselmeermodel/brongegevens.ipynb | 44 -- notebooks/ijsselmeermodel/brongegevens.py | 11 + notebooks/nl-kunstwerken.ipynb | 186 ------- notebooks/nl-kunstwerken.py | 102 ++++ notebooks/uitlaten_inlaten.ipynb | 254 ---------- notebooks/uitlaten_inlaten.py | 171 +++++++ .../hydamo_0_analyse_data_waterboard.ipynb | 281 ----------- .../hydamo_0_analyse_data_waterboard.py | 133 +++++ .../hydamo_1_process_hydamo_data.ipynb | 246 ---------- .../notebooks/hydamo_1_process_hydamo_data.py | 136 ++++++ ...amo_2_run_ribasim_lumping_waterboard.ipynb | 104 ---- ...hydamo_2_run_ribasim_lumping_waterboard.py | 38 ++ ...damo_3_create_ribasim_model_networks.ipynb | 313 ------------ .../hydamo_3_create_ribasim_model_networks.py | 232 +++++++++ ...hydamo_4_create_dummy_ribasim_models.ipynb | 174 ------- .../hydamo_4_create_dummy_ribasim_models.py | 88 ++++ .../xxxx_combine_waterschap_layers.ipynb | 462 ------------------ .../xxxx_combine_waterschap_layers.py | 314 ++++++++++++ 18 files changed, 1225 insertions(+), 2064 deletions(-) delete mode 100644 notebooks/ijsselmeermodel/brongegevens.ipynb create mode 100644 notebooks/ijsselmeermodel/brongegevens.py delete mode 100644 notebooks/nl-kunstwerken.ipynb create mode 100644 notebooks/nl-kunstwerken.py delete mode 100644 notebooks/uitlaten_inlaten.ipynb create mode 100644 notebooks/uitlaten_inlaten.py delete mode 100644 scripts/notebooks/hydamo_0_analyse_data_waterboard.ipynb create mode 100644 scripts/notebooks/hydamo_0_analyse_data_waterboard.py delete mode 100644 scripts/notebooks/hydamo_1_process_hydamo_data.ipynb create mode 100644 scripts/notebooks/hydamo_1_process_hydamo_data.py delete mode 100644 scripts/notebooks/hydamo_2_run_ribasim_lumping_waterboard.ipynb create mode 100644 scripts/notebooks/hydamo_2_run_ribasim_lumping_waterboard.py delete mode 100644 scripts/notebooks/hydamo_3_create_ribasim_model_networks.ipynb create mode 100644 scripts/notebooks/hydamo_3_create_ribasim_model_networks.py delete mode 100644 scripts/notebooks/hydamo_4_create_dummy_ribasim_models.ipynb create mode 100644 scripts/notebooks/hydamo_4_create_dummy_ribasim_models.py delete mode 100644 scripts/notebooks/xxxx_combine_waterschap_layers.ipynb create mode 100644 scripts/notebooks/xxxx_combine_waterschap_layers.py diff --git a/notebooks/ijsselmeermodel/brongegevens.ipynb b/notebooks/ijsselmeermodel/brongegevens.ipynb deleted file mode 100644 index 29917f0..0000000 --- a/notebooks/ijsselmeermodel/brongegevens.ipynb +++ /dev/null @@ -1,44 +0,0 @@ -{ - "cells": [ - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "import os\n", - "from pathlib import Path\n", - "\n", - "import geopandas as gpd" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "DATA_DIR = os.getenv(\"RIBASIM_NL_DATA_DIR\")\n", - "\n", - "# file-paths\n", - "kunstwerken_gpkg = Path(DATA_DIR) / \"nl_kunstwerken.gpkg\"" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "kunstwerken_gdf = gpd.read_file(kunstwerken_gpkg)" - ] - } - ], - "metadata": { - "language_info": { - "name": "python" - } - }, - "nbformat": 4, - "nbformat_minor": 2 -} diff --git a/notebooks/ijsselmeermodel/brongegevens.py b/notebooks/ijsselmeermodel/brongegevens.py new file mode 100644 index 0000000..69935a4 --- /dev/null +++ b/notebooks/ijsselmeermodel/brongegevens.py @@ -0,0 +1,11 @@ +import os +from pathlib import Path + +import geopandas as gpd + +DATA_DIR = os.getenv("RIBASIM_NL_DATA_DIR") + +# file-paths +kunstwerken_gpkg = Path(DATA_DIR) / "nl_kunstwerken.gpkg" + +kunstwerken_gdf = gpd.read_file(kunstwerken_gpkg) diff --git a/notebooks/nl-kunstwerken.ipynb b/notebooks/nl-kunstwerken.ipynb deleted file mode 100644 index 6dd7d44..0000000 --- a/notebooks/nl-kunstwerken.ipynb +++ /dev/null @@ -1,186 +0,0 @@ -{ - "cells": [ - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "import os\n", - "from pathlib import Path\n", - "\n", - "import geopandas as gpd\n", - "import pandas as pd\n", - "import requests\n", - "from shapely.geometry import Point" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# Voorbereiding\n", - "\n", - "Globale variabelen\n", - "`DATA_DIR`: De locale directory waar project-data staat opgeslagen\n", - "`EXCEL_FILE`: Het Excel-bestand dat moet worden ingelezen\n", - "`CRS`: De projectile waarin de ruimtelijke data moet worden opgeslagen (28992 = Rijksdriehoekstelsel) " - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# environmnt variables\n", - "DATA_DIR = os.getenv(\"RIBASIM_NL_DATA_DIR\")\n", - "RIBASIM_NL_CLOUD_PASS = os.getenv(\"RIBASIM_NL_CLOUD_PASS\")\n", - "assert DATA_DIR is not None\n", - "assert RIBASIM_NL_CLOUD_PASS is not None\n", - "\n", - "EXCEL_FILE = r\"# Overzicht kunstwerken primaire keringen waterschappen_ET.xlsx\"\n", - "CRS = 28992\n", - "RIBASIM_NL_CLOUD_USER = \"nhi_api\"\n", - "WEBDAV_URL = \"https://deltares.thegood.cloud/remote.php/dav\"\n", - "BASE_URL = f\"{WEBDAV_URL}/files/{RIBASIM_NL_CLOUD_USER}/D-HYDRO modeldata\"\n", - "\n", - "# file-paths\n", - "kunstwerken_xlsx = Path(DATA_DIR) / EXCEL_FILE\n", - "kunstwerken_gpkg = kunstwerken_xlsx.parent / \"nl_kunstwerken.gpkg\"\n", - "\n", - "\n", - "def upload_file(url, path):\n", - " with open(path, \"rb\") as f:\n", - " r = requests.put(url, data=f, auth=(RIBASIM_NL_CLOUD_USER, RIBASIM_NL_CLOUD_PASS))\n", - " r.raise_for_status()" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Inlezen NL kunstwerken vanuit data dir\n", - "We lezen de geleverde excel in en:\n", - "- skippen de eerste 6 regels boven de header\n", - "- gooien, voor dit project, irrelevante kolommen weg\n", - "- hernoemen de kolom met organisatie-naam, x en y coordinaat\n", - "- transformeren de x en y coordinaat; wordt NaN wanneer data dat niet toelaat (text of missend)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# skip first rows\n", - "kunstwerken_df = pd.read_excel(kunstwerken_xlsx, skiprows=6)\n", - "\n", - "# drop irrelevant columns\n", - "columns = kunstwerken_df.columns[1:13]\n", - "kunstwerken_df = kunstwerken_df.loc[:, columns]\n", - "\n", - "# rename columns into our liking\n", - "kunstwerken_df.rename(\n", - " columns={\n", - " \"Unnamed: 1\": \"organisatie\",\n", - " \"Y coördinaat RD\": \"y\",\n", - " \"X coördinaat RD\": \"x\",\n", - " },\n", - " inplace=True,\n", - ")\n", - "\n", - "# drop no-data rows\n", - "kunstwerken_df = kunstwerken_df[~kunstwerken_df[\"organisatie\"].isna()]\n", - "\n", - "# convert x/y to numeric\n", - "kunstwerken_df[\"x\"] = pd.to_numeric(kunstwerken_df[\"x\"], errors=\"coerce\")\n", - "kunstwerken_df[\"y\"] = pd.to_numeric(kunstwerken_df[\"y\"], errors=\"coerce\")" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Aanmaken niet-ruimtelijke table\n", - "Waar x/y coordinaten mizzen maken we een niet-ruimtelijke table die we wegschrijven in geen GeoPackage" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# make a non-spatial GeoDataFrame where x/y are missing\n", - "kunstwerken_non_spatial_df = kunstwerken_df[kunstwerken_df[\"x\"].isna() | kunstwerken_df[\"y\"].isna()]\n", - "kunstwerken_non_spatial_gdf = gpd.GeoDataFrame(kunstwerken_non_spatial_df, geometry=gpd.GeoSeries(), crs=28992)\n", - "\n", - "# writ to GeoPackage\n", - "kunstwerken_non_spatial_gdf.to_file(kunstwerken_gpkg, layer=\"kunstwerken (geen coordinaten)\")" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Aanmaken ruimtelijke table\n", - "Waar x/y coordinaten beschikbaar zijn maken we een ruimtelijke table die we wegschrijven in geen GeoPackage" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# make a spatial GeoDataFrame where x/y exist\n", - "kunstwerken_spatial_df = kunstwerken_df[~kunstwerken_df[\"x\"].isna() & ~kunstwerken_df[\"y\"].isna()]\n", - "geometry_series = gpd.GeoSeries(kunstwerken_spatial_df.apply((lambda x: Point(x.x, x.y)), axis=1))\n", - "kunstwerken_spatial_gdf = gpd.GeoDataFrame(kunstwerken_spatial_df, geometry=geometry_series, crs=CRS)\n", - "\n", - "# write to GeoPackage\n", - "kunstwerken_spatial_gdf.to_file(kunstwerken_gpkg, layer=\"kunstwerken\")" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Upload geopackage" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "to_url = f\"{BASE_URL}/Rijkswaterstaat/{kunstwerken_gpkg.name}\"\n", - "upload_file(to_url, kunstwerken_gpkg)" - ] - } - ], - "metadata": { - "kernelspec": { - "display_name": "Python 3", - "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.12.1" - } - }, - "nbformat": 4, - "nbformat_minor": 2 -} diff --git a/notebooks/nl-kunstwerken.py b/notebooks/nl-kunstwerken.py new file mode 100644 index 0000000..c41f6c8 --- /dev/null +++ b/notebooks/nl-kunstwerken.py @@ -0,0 +1,102 @@ +import os +from pathlib import Path + +import geopandas as gpd +import pandas as pd +import requests +from shapely.geometry import Point + +# # Voorbereiding +# +# Globale variabelen +# `DATA_DIR`: De locale directory waar project-data staat opgeslagen +# `EXCEL_FILE`: Het Excel-bestand dat moet worden ingelezen +# `CRS`: De projectile waarin de ruimtelijke data moet worden opgeslagen (28992 = Rijksdriehoekstelsel) + + +# environmnt variables +DATA_DIR = os.getenv("RIBASIM_NL_DATA_DIR") +RIBASIM_NL_CLOUD_PASS = os.getenv("RIBASIM_NL_CLOUD_PASS") +assert DATA_DIR is not None +assert RIBASIM_NL_CLOUD_PASS is not None + +EXCEL_FILE = r"# Overzicht kunstwerken primaire keringen waterschappen_ET.xlsx" +CRS = 28992 +RIBASIM_NL_CLOUD_USER = "nhi_api" +WEBDAV_URL = "https://deltares.thegood.cloud/remote.php/dav" +BASE_URL = f"{WEBDAV_URL}/files/{RIBASIM_NL_CLOUD_USER}/D-HYDRO modeldata" + +# file-paths +kunstwerken_xlsx = Path(DATA_DIR) / EXCEL_FILE +kunstwerken_gpkg = kunstwerken_xlsx.parent / "nl_kunstwerken.gpkg" + + +def upload_file(url, path): + with open(path, "rb") as f: + r = requests.put(url, data=f, auth=(RIBASIM_NL_CLOUD_USER, RIBASIM_NL_CLOUD_PASS)) + r.raise_for_status() + + +# ## Inlezen NL kunstwerken vanuit data dir +# We lezen de geleverde excel in en: +# - skippen de eerste 6 regels boven de header +# - gooien, voor dit project, irrelevante kolommen weg +# - hernoemen de kolom met organisatie-naam, x en y coordinaat +# - transformeren de x en y coordinaat; wordt NaN wanneer data dat niet toelaat (text of missend) + + +# skip first rows +kunstwerken_df = pd.read_excel(kunstwerken_xlsx, skiprows=6) + +# drop irrelevant columns +columns = kunstwerken_df.columns[1:13] +kunstwerken_df = kunstwerken_df.loc[:, columns] + +# rename columns into our liking +kunstwerken_df.rename( + columns={ + "Unnamed: 1": "organisatie", + "Y coördinaat RD": "y", + "X coördinaat RD": "x", + }, + inplace=True, +) + +# drop no-data rows +kunstwerken_df = kunstwerken_df[~kunstwerken_df["organisatie"].isna()] + +# convert x/y to numeric +kunstwerken_df["x"] = pd.to_numeric(kunstwerken_df["x"], errors="coerce") +kunstwerken_df["y"] = pd.to_numeric(kunstwerken_df["y"], errors="coerce") + + +# ## Aanmaken niet-ruimtelijke table +# Waar x/y coordinaten mizzen maken we een niet-ruimtelijke table die we wegschrijven in geen GeoPackage + + +# make a non-spatial GeoDataFrame where x/y are missing +kunstwerken_non_spatial_df = kunstwerken_df[kunstwerken_df["x"].isna() | kunstwerken_df["y"].isna()] +kunstwerken_non_spatial_gdf = gpd.GeoDataFrame(kunstwerken_non_spatial_df, geometry=gpd.GeoSeries(), crs=28992) + +# writ to GeoPackage +kunstwerken_non_spatial_gdf.to_file(kunstwerken_gpkg, layer="kunstwerken (geen coordinaten)") + + +# ## Aanmaken ruimtelijke table +# Waar x/y coordinaten beschikbaar zijn maken we een ruimtelijke table die we wegschrijven in geen GeoPackage + + +# make a spatial GeoDataFrame where x/y exist +kunstwerken_spatial_df = kunstwerken_df[~kunstwerken_df["x"].isna() & ~kunstwerken_df["y"].isna()] +geometry_series = gpd.GeoSeries(kunstwerken_spatial_df.apply((lambda x: Point(x.x, x.y)), axis=1)) +kunstwerken_spatial_gdf = gpd.GeoDataFrame(kunstwerken_spatial_df, geometry=geometry_series, crs=CRS) + +# write to GeoPackage +kunstwerken_spatial_gdf.to_file(kunstwerken_gpkg, layer="kunstwerken") + + +# ## Upload geopackage + + +to_url = f"{BASE_URL}/Rijkswaterstaat/{kunstwerken_gpkg.name}" +upload_file(to_url, kunstwerken_gpkg) diff --git a/notebooks/uitlaten_inlaten.ipynb b/notebooks/uitlaten_inlaten.ipynb deleted file mode 100644 index c914db2..0000000 --- a/notebooks/uitlaten_inlaten.ipynb +++ /dev/null @@ -1,254 +0,0 @@ -{ - "cells": [ - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "import os\n", - "import warnings\n", - "from pathlib import Path\n", - "\n", - "import fiona\n", - "import geopandas as gpd\n", - "import pandas as pd\n", - "import requests\n", - "from shapely.geometry import Point\n", - "\n", - "from hydamo import HyDAMO, code_utils\n", - "\n", - "warnings.simplefilter(\"ignore\", UserWarning)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# Voorbereiding\n", - "\n", - "Globale variabelen\n", - "- `DATA_DIR`: De locale directory waar project-data staat opgeslagen\n", - "- `EXCEL_FILE`: Het Excel-bestand dat moet worden ingelezen\n", - "- `CRS`: De projectile waarin de ruimtelijke data moet worden opgeslagen (28992 = Rijksdriehoekstelsel) " - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# environmnt variables\n", - "DATA_DIR = os.getenv(\"RIBASIM_NL_DATA_DIR\")\n", - "RIBASIM_NL_CLOUD_PASS = os.getenv(\"RIBASIM_NL_CLOUD_PASS\")\n", - "assert DATA_DIR is not None\n", - "assert RIBASIM_NL_CLOUD_PASS is not None\n", - "\n", - "DATA_DIR = Path(DATA_DIR)\n", - "EXCEL_FILE = r\"uitlaten_inlaten.xlsx\"\n", - "CRS = 28992\n", - "RIBASIM_NL_CLOUD_USER = \"nhi_api\"\n", - "WEBDAV_URL = \"https://deltares.thegood.cloud/remote.php/dav\"\n", - "BASE_URL = f\"{WEBDAV_URL}/files/{RIBASIM_NL_CLOUD_USER}/D-HYDRO modeldata\"\n", - "\n", - "# file-paths\n", - "kunstwerken_xlsx = Path(DATA_DIR) / EXCEL_FILE\n", - "kunstwerken_gpkg = kunstwerken_xlsx.parent / f\"{kunstwerken_xlsx.stem}.gpkg\"\n", - "\n", - "\n", - "def upload_file(url, path):\n", - " with open(path, \"rb\") as f:\n", - " r = requests.put(url, data=f, auth=(RIBASIM_NL_CLOUD_USER, RIBASIM_NL_CLOUD_PASS))\n", - " r.raise_for_status()" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Inlezen NL kunstwerken vanuit data dir\n", - "- Inlezen Excel\n", - "- Nu masken we (nog) op kunstwerken die we uit files kunnen trekken" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "kunstwerken_df = pd.read_excel(kunstwerken_xlsx)\n", - "kunstwerken_df[\"user_id\"] = None\n", - "files_mask = ~kunstwerken_df[\"damo_bestand\"].isna()\n", - "data = {}" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Toevoegen kunstwerken op XY-locatie" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "for kwk_row in kunstwerken_df[~files_mask].itertuples():\n", - " layer = kwk_row.hydamo_object\n", - " if layer not in data.keys():\n", - " data[layer] = []\n", - " name = kwk_row.dm_naam\n", - " geometry = Point(kwk_row.x, kwk_row.y)\n", - " code = code_utils.code_from_geometry(geometry)\n", - " result = {\"code\": code, \"naam\": name, \"geometry\": geometry}\n", - "\n", - " kunstwerken_df.loc[kwk_row.Index, [\"user_id\"]] = code_utils.generate_model_id(\n", - " result[\"code\"], layer, bgt_code=kwk_row.bgt_code\n", - " )\n", - "\n", - " data[layer] += [result]" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Ophalen kunstwerken\n", - "- Aanmaken lege data-dict\n", - "- loopen over kunstwerken/lagen en daar de relevante kunstwerken uit halen" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# create data-dict for every layer\n", - "for layer in kunstwerken_df[\"hydamo_object\"].unique():\n", - " if layer not in data.keys():\n", - " data[layer] = []\n", - "\n", - "# group by file and check if exists\n", - "for file, file_df in kunstwerken_df[files_mask].groupby(\"damo_bestand\"):\n", - " file = DATA_DIR.joinpath(file)\n", - " assert file.exists()\n", - "\n", - " # open per layer, and check if specified layer xists\n", - " for layer, layer_df in file_df.groupby(\"damo_laag\"):\n", - " file_layers = fiona.listlayers(file)\n", - " if (\n", - " len(file_layers) == 1\n", - " ): # in case single-layer files, users don't understand a `layer-property` and make mistakes\n", - " layer = file_layers[0]\n", - " if layer not in fiona.listlayers(file):\n", - " raise ValueError(f\"layer '{layer}' not a layer in '{file}'. Specify one of {fiona.listlayers(file)}\")\n", - " print(f\"reading {file.name}, layer {layer}\")\n", - " gdf = gpd.read_file(file, layer=layer)\n", - " gdf.columns = [i.lower() for i in gdf.columns]\n", - " gdf.to_crs(CRS, inplace=True)\n", - "\n", - " # read every row this file-layer group and get the source-info\n", - " for kwk_row in layer_df.itertuples():\n", - " # get the index from the used code or name column\n", - " damo_index = kwk_row.damo_ident_kolom\n", - " src_index = getattr(kwk_row, f\"damo_{damo_index}_kolom\").strip().lower()\n", - " index_value = str(kwk_row.damo_waarde)\n", - "\n", - " # read the source\n", - " if src_index in gdf.columns:\n", - " if index_value in gdf[src_index].to_numpy():\n", - " src_row = gdf.set_index(src_index).loc[index_value]\n", - " else:\n", - " raise KeyError(f\"'{index_value}' not a value in '{file}', layer '{layer}', column '{src_index}'\")\n", - " else:\n", - " raise KeyError(\n", - " f\"'{src_index}' not a column in '{file}', layer '{layer}' (searching value '{index_value}')\"\n", - " )\n", - "\n", - " # populate the result\n", - " result = {}\n", - "\n", - " # populate code and naam fields\n", - " for damo_att in [\"code\", \"naam\"]:\n", - " if damo_att == damo_index:\n", - " result[damo_att] = index_value\n", - " else:\n", - " column = getattr(kwk_row, f\"damo_{damo_att}_kolom\")\n", - " if pd.isna(column) and (damo_att == \"code\"):\n", - " result[damo_att] = code_utils.code_from_geometry(src_row.geometry)\n", - " else:\n", - " column = getattr(kwk_row, f\"damo_{damo_att}_kolom\").strip().lower()\n", - " result[damo_att] = str(getattr(src_row, column))\n", - "\n", - " # get the geometry. We get the centroid to avoid flatten all kinds of mult-features\n", - " result[\"geometry\"] = src_row.geometry.centroid\n", - "\n", - " # add it to our data dictionary\n", - " data[kwk_row.hydamo_object] += [result]\n", - "\n", - " # update Escel with user-id\n", - " kunstwerken_df.loc[kwk_row.Index, [\"user_id\"]] = code_utils.generate_model_id(\n", - " result[\"code\"], kwk_row.hydamo_object, bgt_code=kwk_row.bgt_code\n", - " )\n", - " kunstwerken_df.loc[kwk_row.Index, [\"x\"]] = result[\"geometry\"].x\n", - " kunstwerken_df.loc[kwk_row.Index, [\"y\"]] = result[\"geometry\"].y" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Wegschrijven HyDAMO\n", - "- lokaal\n", - "- op TheGoodCloud" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "hydamo = HyDAMO(\"2.2.1\")\n", - "for layer in data.keys():\n", - " if layer != \"duikersifonhevel\":\n", - " gdf = gpd.GeoDataFrame(data[layer], crs=CRS)\n", - " getattr(hydamo, layer).set_data(gdf, check_columns=False)\n", - "\n", - "hydamo.to_geopackage(kunstwerken_gpkg)\n", - "\n", - "kunstwerken_df.to_excel(kunstwerken_xlsx, index=False, sheet_name=\"kunstwerken\")\n", - "\n", - "for file in [kunstwerken_xlsx, kunstwerken_gpkg]:\n", - " to_url = f\"{BASE_URL}/HyDAMO_geconstrueerd/{file.name}\"\n", - " upload_file(to_url, file)" - ] - } - ], - "metadata": { - "kernelspec": { - "display_name": "Python 3", - "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.6" - } - }, - "nbformat": 4, - "nbformat_minor": 2 -} diff --git a/notebooks/uitlaten_inlaten.py b/notebooks/uitlaten_inlaten.py new file mode 100644 index 0000000..c9e4e91 --- /dev/null +++ b/notebooks/uitlaten_inlaten.py @@ -0,0 +1,171 @@ +import os +import warnings +from pathlib import Path + +import fiona +import geopandas as gpd +import pandas as pd +import requests +from shapely.geometry import Point + +from hydamo import HyDAMO, code_utils + +warnings.simplefilter("ignore", UserWarning) + + +# # Voorbereiding +# +# Globale variabelen +# - `DATA_DIR`: De locale directory waar project-data staat opgeslagen +# - `EXCEL_FILE`: Het Excel-bestand dat moet worden ingelezen +# - `CRS`: De projectile waarin de ruimtelijke data moet worden opgeslagen (28992 = Rijksdriehoekstelsel) + + +# environmnt variables +DATA_DIR = os.getenv("RIBASIM_NL_DATA_DIR") +RIBASIM_NL_CLOUD_PASS = os.getenv("RIBASIM_NL_CLOUD_PASS") +assert DATA_DIR is not None +assert RIBASIM_NL_CLOUD_PASS is not None + +DATA_DIR = Path(DATA_DIR) +EXCEL_FILE = r"uitlaten_inlaten.xlsx" +CRS = 28992 +RIBASIM_NL_CLOUD_USER = "nhi_api" +WEBDAV_URL = "https://deltares.thegood.cloud/remote.php/dav" +BASE_URL = f"{WEBDAV_URL}/files/{RIBASIM_NL_CLOUD_USER}/D-HYDRO modeldata" + +# file-paths +kunstwerken_xlsx = Path(DATA_DIR) / EXCEL_FILE +kunstwerken_gpkg = kunstwerken_xlsx.parent / f"{kunstwerken_xlsx.stem}.gpkg" + + +def upload_file(url, path): + with open(path, "rb") as f: + r = requests.put(url, data=f, auth=(RIBASIM_NL_CLOUD_USER, RIBASIM_NL_CLOUD_PASS)) + r.raise_for_status() + + +# ## Inlezen NL kunstwerken vanuit data dir +# - Inlezen Excel +# - Nu masken we (nog) op kunstwerken die we uit files kunnen trekken + + +kunstwerken_df = pd.read_excel(kunstwerken_xlsx) +kunstwerken_df["user_id"] = None +files_mask = ~kunstwerken_df["damo_bestand"].isna() +data = {} + + +# ## Toevoegen kunstwerken op XY-locatie + + +for kwk_row in kunstwerken_df[~files_mask].itertuples(): + layer = kwk_row.hydamo_object + if layer not in data.keys(): + data[layer] = [] + name = kwk_row.dm_naam + geometry = Point(kwk_row.x, kwk_row.y) + code = code_utils.code_from_geometry(geometry) + result = {"code": code, "naam": name, "geometry": geometry} + + kunstwerken_df.loc[kwk_row.Index, ["user_id"]] = code_utils.generate_model_id( + result["code"], layer, bgt_code=kwk_row.bgt_code + ) + + data[layer] += [result] + + +# ## Ophalen kunstwerken +# - Aanmaken lege data-dict +# - loopen over kunstwerken/lagen en daar de relevante kunstwerken uit halen + + +# create data-dict for every layer +for layer in kunstwerken_df["hydamo_object"].unique(): + if layer not in data.keys(): + data[layer] = [] + +# group by file and check if exists +for file, file_df in kunstwerken_df[files_mask].groupby("damo_bestand"): + file = DATA_DIR.joinpath(file) + assert file.exists() + + # open per layer, and check if specified layer xists + for layer, layer_df in file_df.groupby("damo_laag"): + file_layers = fiona.listlayers(file) + if ( + len(file_layers) == 1 + ): # in case single-layer files, users don't understand a `layer-property` and make mistakes + layer = file_layers[0] + if layer not in fiona.listlayers(file): + raise ValueError(f"layer '{layer}' not a layer in '{file}'. Specify one of {fiona.listlayers(file)}") + print(f"reading {file.name}, layer {layer}") + gdf = gpd.read_file(file, layer=layer) + gdf.columns = [i.lower() for i in gdf.columns] + gdf.to_crs(CRS, inplace=True) + + # read every row this file-layer group and get the source-info + for kwk_row in layer_df.itertuples(): + # get the index from the used code or name column + damo_index = kwk_row.damo_ident_kolom + src_index = getattr(kwk_row, f"damo_{damo_index}_kolom").strip().lower() + index_value = str(kwk_row.damo_waarde) + + # read the source + if src_index in gdf.columns: + if index_value in gdf[src_index].to_numpy(): + src_row = gdf.set_index(src_index).loc[index_value] + else: + raise KeyError(f"'{index_value}' not a value in '{file}', layer '{layer}', column '{src_index}'") + else: + raise KeyError( + f"'{src_index}' not a column in '{file}', layer '{layer}' (searching value '{index_value}')" + ) + + # populate the result + result = {} + + # populate code and naam fields + for damo_att in ["code", "naam"]: + if damo_att == damo_index: + result[damo_att] = index_value + else: + column = getattr(kwk_row, f"damo_{damo_att}_kolom") + if pd.isna(column) and (damo_att == "code"): + result[damo_att] = code_utils.code_from_geometry(src_row.geometry) + else: + column = getattr(kwk_row, f"damo_{damo_att}_kolom").strip().lower() + result[damo_att] = str(getattr(src_row, column)) + + # get the geometry. We get the centroid to avoid flatten all kinds of mult-features + result["geometry"] = src_row.geometry.centroid + + # add it to our data dictionary + data[kwk_row.hydamo_object] += [result] + + # update Escel with user-id + kunstwerken_df.loc[kwk_row.Index, ["user_id"]] = code_utils.generate_model_id( + result["code"], kwk_row.hydamo_object, bgt_code=kwk_row.bgt_code + ) + kunstwerken_df.loc[kwk_row.Index, ["x"]] = result["geometry"].x + kunstwerken_df.loc[kwk_row.Index, ["y"]] = result["geometry"].y + + +# ## Wegschrijven HyDAMO +# - lokaal +# - op TheGoodCloud + + +hydamo = HyDAMO("2.2.1") +for layer in data.keys(): + if layer != "duikersifonhevel": + gdf = gpd.GeoDataFrame(data[layer], crs=CRS) + getattr(hydamo, layer).set_data(gdf, check_columns=False) + +hydamo.to_geopackage(kunstwerken_gpkg) + +kunstwerken_df.to_excel(kunstwerken_xlsx, index=False, sheet_name="kunstwerken") + +for file in [kunstwerken_xlsx, kunstwerken_gpkg]: + to_url = f"{BASE_URL}/HyDAMO_geconstrueerd/{file.name}" + upload_file(to_url, file) diff --git a/scripts/notebooks/hydamo_0_analyse_data_waterboard.ipynb b/scripts/notebooks/hydamo_0_analyse_data_waterboard.ipynb deleted file mode 100644 index 7628456..0000000 --- a/scripts/notebooks/hydamo_0_analyse_data_waterboard.ipynb +++ /dev/null @@ -1,281 +0,0 @@ -{ - "cells": [ - { - "cell_type": "code", - "execution_count": null, - "id": "0", - "metadata": { - "tags": [] - }, - "outputs": [], - "source": [ - "from pathlib import Path\n", - "\n", - "import pandas as pd\n", - "from pandas_xlsx_tables import xlsx_tables_to_dfs\n", - "from ribasim_lumping_tools.LHM_data_bewerking_analyse_utils import (\n", - " check_ids_hydamo_data,\n", - " check_if_object_on_hydroobject,\n", - " read_original_data,\n", - " translate_data_to_hydamo_format,\n", - ")\n", - "\n", - "from hydamo.datamodel import HyDAMO" - ] - }, - { - "cell_type": "markdown", - "id": "1", - "metadata": {}, - "source": [ - "Vertaal originele data naar Hydamo data zoals gedefinieerd in de tabel hydamo_data_format.xlsx" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "2", - "metadata": { - "tags": [] - }, - "outputs": [], - "source": [ - "base_dir = \"..\\\\\"\n", - "\n", - "waterboard = \"AAenMaas\"\n", - "waterboard_code = 1" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "3", - "metadata": {}, - "outputs": [], - "source": [ - "waterboard_dir = Path(base_dir, waterboard, \"verwerkt\")\n", - "path_hydamo_format = Path(waterboard_dir, \"HyDAMO_format_AAenMaas.xlsx\")\n", - "hydamo_format = xlsx_tables_to_dfs(path_hydamo_format)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "4", - "metadata": { - "tags": [] - }, - "outputs": [], - "source": [ - "# eerst inlezen hydroobject, vertalen naar hydamo\n", - "hydamo_object = \"hydroobject\"\n", - "hydamo_translate_table, data_original = read_original_data(waterboard_dir, hydamo_format, hydamo_object, waterboard)\n", - "hydroobject = translate_data_to_hydamo_format(hydamo_translate_table, data_original)\n", - "\n", - "# maak een created_date aan indien nodig\n", - "if \"created_date\" not in data_original.columns:\n", - " hydroobject[\"created_date\"] = pd.NaT\n", - "# transformeer created_date waardes indien nodig\n", - "hydroobject[\"created_date\"] = hydroobject[\"created_date\"].replace(\"\", pd.NaT)\n", - "\n", - "# hydroobject.loc[hydroobject['code'].duplicated(keep=False), 'data_issue'] = 'duplicate_id'\n", - "data_hydamo_dict = {\"hydroobject\": hydroobject.set_crs(28992)}\n", - "\n", - "# geometry hydroobject bufferen met 10 cm voor de spatial join\n", - "hydroobject[\"buffer\"] = hydroobject.copy().buffer(5) # 5 meter buffer omdat anders relevante gemalen wegvallen\n", - "hydroobject_buffered = hydroobject.set_geometry(\"buffer\").set_crs(28992)" - ] - }, - { - "cell_type": "markdown", - "id": "5", - "metadata": {}, - "source": [ - "Specificeer welke HyDAMO data je wilt omzetten" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "6", - "metadata": { - "tags": [] - }, - "outputs": [], - "source": [ - "hydamo_objects = [\n", - " \"stuw\",\n", - " \"gemaal\",\n", - " \"afvoergebiedaanvoergebied\",\n", - " \"pomp\",\n", - " ##'peilgebiedvigerend',\n", - " ##'peilgebiedpraktijk',\n", - " ##'streefpeil',\n", - " \"duikersifonhevel\",\n", - " ##'afsluiter',\n", - " ##'sluis',\n", - "]" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "7", - "metadata": { - "tags": [] - }, - "outputs": [], - "source": [ - "for hydamo_object in hydamo_objects:\n", - " # lees aangeleverde data en hydamo tabel voor gegeven kunstwerk en waterschap\n", - " table_hydamo, data_original = read_original_data(waterboard_dir, hydamo_format, hydamo_object, waterboard)\n", - " if data_original is None:\n", - " data_hydamo_dict[hydamo_object] = None\n", - " else:\n", - " # vertaal data naar hydamo-ribasim format\n", - " data_hydamo = translate_data_to_hydamo_format(table_hydamo, data_original)\n", - "\n", - " # maak een created_date aan indien nodig\n", - " if \"created_date\" not in data_original.columns and hydamo_object != \"sluis\":\n", - " hydroobject[\"created_date\"] = pd.NaT\n", - " if \"last_edited_date\" not in data_original.columns and hydamo_object == \"afsluiter\":\n", - " hydroobject[\"last_edited_date\"] = pd.NaT\n", - " if \"lvpublicatiedatum\" not in data_original.columns and hydamo_object == \"afsluiter\":\n", - " hydroobject[\"lvpublicatiedatum\"] = pd.NaT\n", - "\n", - " # transformeer created_date waardes indien nodig\n", - " if hydamo_object != \"sluis\":\n", - " data_hydamo[\"created_date\"] = data_hydamo[\"created_date\"].replace(\"\", pd.NaT)\n", - " if hydamo_object == \"afsluiter\":\n", - " data_hydamo[\"last_edited_date\"] = data_hydamo[\"last_edited_date\"].replace(\"\", pd.NaT)\n", - " data_hydamo[\"lvpublicatiedatum\"] = data_hydamo[\"lvpublicatiedatum\"].replace(\"\", pd.NaT)\n", - "\n", - " # check dubbele id's\n", - " if hydamo_object not in [\"streefpeil\"]: # streefpeil heeft geen code, alleen globalid etc\n", - " data_hydamo.loc[data_hydamo[\"code\"].duplicated(keep=False), \"data_issue\"] = \"duplicate_id\"\n", - " # TODO check op 'code' lijkt met logischer want die kolom wordt vaker gebruikt. Maar bij WDOD bijv. is die niet ingevuld. Toch op globalid?\n", - " # check of kuntstwerk op hydroobject ligt\n", - " if hydamo_object in [\"stuw\", \"gemaal\", \"duikersifonhevel\", \"sluis\"]:\n", - " data_hydamo = check_if_object_on_hydroobject(\n", - " data_hydamo=data_hydamo, hydroobject_buffered=hydroobject_buffered\n", - " )\n", - " # verwijder kunstwerken die niet op hydroobject liggen\n", - " data_hydamo = data_hydamo[data_hydamo[\"code_hydroobject\"] != \"niet op hydroobject\"]\n", - " data_hydamo = data_hydamo.reset_index()\n", - " # voeg toe aan de hydamo dataset\n", - " data_hydamo_dict[hydamo_object] = data_hydamo" - ] - }, - { - "cell_type": "markdown", - "id": "8", - "metadata": {}, - "source": [ - "Waterschap specifieke acties" - ] - }, - { - "cell_type": "markdown", - "id": "9", - "metadata": {}, - "source": [ - "Export normal" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "10", - "metadata": { - "tags": [] - }, - "outputs": [], - "source": [ - "# for hydamo_object in ['hydroobject'] + hydamo_objects:\n", - "# # export to geopackage\n", - "# export_to_geopackage(\n", - "# data_hydamo=data_hydamo_dict[hydamo_object],\n", - "# hydamo_format=hydamo_format,\n", - "# waterboard=waterboard,\n", - "# hydamo_object=hydamo_object\n", - "# )" - ] - }, - { - "cell_type": "markdown", - "id": "11", - "metadata": {}, - "source": [ - "### ribasim-nl hydamo" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "12", - "metadata": { - "tags": [] - }, - "outputs": [], - "source": [ - "hydamo = HyDAMO(version=\"2.2.1_sweco\")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "13", - "metadata": {}, - "outputs": [], - "source": [ - "for hydamo_object in [\"hydroobject\"] + hydamo_objects:\n", - " data_hydamo = data_hydamo_dict[hydamo_object]\n", - " if hydamo_object == \"stuw\":\n", - " data_hydamo = data_hydamo.drop(columns=[\"code_hydroobject\", \"data_issue\"]) # ,'index_right'\n", - " data_hydamo = check_ids_hydamo_data(data_hydamo, waterboard_code, hydamo_object)\n", - " setattr(hydamo, hydamo_object, data_hydamo)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "14", - "metadata": { - "tags": [] - }, - "outputs": [], - "source": [ - "hydamo.to_geopackage(\"..\\\\hydamo.gpkg\")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "15", - "metadata": {}, - "outputs": [], - "source": [] - } - ], - "metadata": { - "kernelspec": { - "display_name": "ribasim_lumping_venv", - "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" - } - }, - "nbformat": 4, - "nbformat_minor": 5 -} diff --git a/scripts/notebooks/hydamo_0_analyse_data_waterboard.py b/scripts/notebooks/hydamo_0_analyse_data_waterboard.py new file mode 100644 index 0000000..7c717af --- /dev/null +++ b/scripts/notebooks/hydamo_0_analyse_data_waterboard.py @@ -0,0 +1,133 @@ +from pathlib import Path + +import pandas as pd +from pandas_xlsx_tables import xlsx_tables_to_dfs +from ribasim_lumping_tools.LHM_data_bewerking_analyse_utils import ( + check_ids_hydamo_data, + check_if_object_on_hydroobject, + read_original_data, + translate_data_to_hydamo_format, +) + +from hydamo.datamodel import HyDAMO + +# Vertaal originele data naar Hydamo data zoals gedefinieerd in de tabel hydamo_data_format.xlsx + + +base_dir = "..\\" + +waterboard = "AAenMaas" +waterboard_code = 1 + + +waterboard_dir = Path(base_dir, waterboard, "verwerkt") +path_hydamo_format = Path(waterboard_dir, "HyDAMO_format_AAenMaas.xlsx") +hydamo_format = xlsx_tables_to_dfs(path_hydamo_format) + + +# eerst inlezen hydroobject, vertalen naar hydamo +hydamo_object = "hydroobject" +hydamo_translate_table, data_original = read_original_data(waterboard_dir, hydamo_format, hydamo_object, waterboard) +hydroobject = translate_data_to_hydamo_format(hydamo_translate_table, data_original) + +# maak een created_date aan indien nodig +if "created_date" not in data_original.columns: + hydroobject["created_date"] = pd.NaT +# transformeer created_date waardes indien nodig +hydroobject["created_date"] = hydroobject["created_date"].replace("", pd.NaT) + +# hydroobject.loc[hydroobject['code'].duplicated(keep=False), 'data_issue'] = 'duplicate_id' +data_hydamo_dict = {"hydroobject": hydroobject.set_crs(28992)} + +# geometry hydroobject bufferen met 10 cm voor de spatial join +hydroobject["buffer"] = hydroobject.copy().buffer(5) # 5 meter buffer omdat anders relevante gemalen wegvallen +hydroobject_buffered = hydroobject.set_geometry("buffer").set_crs(28992) + + +# Specificeer welke HyDAMO data je wilt omzetten + + +hydamo_objects = [ + "stuw", + "gemaal", + "afvoergebiedaanvoergebied", + "pomp", + ##'peilgebiedvigerend', + ##'peilgebiedpraktijk', + ##'streefpeil', + "duikersifonhevel", + ##'afsluiter', + ##'sluis', +] + + +for hydamo_object in hydamo_objects: + # lees aangeleverde data en hydamo tabel voor gegeven kunstwerk en waterschap + table_hydamo, data_original = read_original_data(waterboard_dir, hydamo_format, hydamo_object, waterboard) + if data_original is None: + data_hydamo_dict[hydamo_object] = None + else: + # vertaal data naar hydamo-ribasim format + data_hydamo = translate_data_to_hydamo_format(table_hydamo, data_original) + + # maak een created_date aan indien nodig + if "created_date" not in data_original.columns and hydamo_object != "sluis": + hydroobject["created_date"] = pd.NaT + if "last_edited_date" not in data_original.columns and hydamo_object == "afsluiter": + hydroobject["last_edited_date"] = pd.NaT + if "lvpublicatiedatum" not in data_original.columns and hydamo_object == "afsluiter": + hydroobject["lvpublicatiedatum"] = pd.NaT + + # transformeer created_date waardes indien nodig + if hydamo_object != "sluis": + data_hydamo["created_date"] = data_hydamo["created_date"].replace("", pd.NaT) + if hydamo_object == "afsluiter": + data_hydamo["last_edited_date"] = data_hydamo["last_edited_date"].replace("", pd.NaT) + data_hydamo["lvpublicatiedatum"] = data_hydamo["lvpublicatiedatum"].replace("", pd.NaT) + + # check dubbele id's + if hydamo_object not in ["streefpeil"]: # streefpeil heeft geen code, alleen globalid etc + data_hydamo.loc[data_hydamo["code"].duplicated(keep=False), "data_issue"] = "duplicate_id" + # TODO check op 'code' lijkt met logischer want die kolom wordt vaker gebruikt. Maar bij WDOD bijv. is die niet ingevuld. Toch op globalid? + # check of kuntstwerk op hydroobject ligt + if hydamo_object in ["stuw", "gemaal", "duikersifonhevel", "sluis"]: + data_hydamo = check_if_object_on_hydroobject( + data_hydamo=data_hydamo, hydroobject_buffered=hydroobject_buffered + ) + # verwijder kunstwerken die niet op hydroobject liggen + data_hydamo = data_hydamo[data_hydamo["code_hydroobject"] != "niet op hydroobject"] + data_hydamo = data_hydamo.reset_index() + # voeg toe aan de hydamo dataset + data_hydamo_dict[hydamo_object] = data_hydamo + + +# Waterschap specifieke acties + +# Export normal + + +# for hydamo_object in ['hydroobject'] + hydamo_objects: +# # export to geopackage +# export_to_geopackage( +# data_hydamo=data_hydamo_dict[hydamo_object], +# hydamo_format=hydamo_format, +# waterboard=waterboard, +# hydamo_object=hydamo_object +# ) + + +# ### ribasim-nl hydamo + + +hydamo = HyDAMO(version="2.2.1_sweco") + + +for hydamo_object in ["hydroobject"] + hydamo_objects: + data_hydamo = data_hydamo_dict[hydamo_object] + if hydamo_object == "stuw": + data_hydamo = data_hydamo.drop(columns=["code_hydroobject", "data_issue"]) # ,'index_right' + data_hydamo = check_ids_hydamo_data(data_hydamo, waterboard_code, hydamo_object) + setattr(hydamo, hydamo_object, data_hydamo) + + +hydamo.to_geopackage("..\\hydamo.gpkg") diff --git a/scripts/notebooks/hydamo_1_process_hydamo_data.ipynb b/scripts/notebooks/hydamo_1_process_hydamo_data.ipynb deleted file mode 100644 index ccf6472..0000000 --- a/scripts/notebooks/hydamo_1_process_hydamo_data.ipynb +++ /dev/null @@ -1,246 +0,0 @@ -{ - "cells": [ - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Verwerk wijzigingen aan HyDAMO geopackages in nieuwe HyDAMO input geopackage\n", - "\n", - "Contactpersoon: Harm Nomden (Sweco)\n", - "\n", - "Laatste update: 15-03-2024" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "import warnings\n", - "from pathlib import Path\n", - "\n", - "import fiona\n", - "import geopandas as gpd\n", - "import numpy as np\n", - "import pandas as pd\n", - "from hydamo_preprocessing.preprocessing import preprocess_hydamo_hydroobjects\n", - "\n", - "warnings.simplefilter(action=\"ignore\", category=UserWarning)\n", - "warnings.simplefilter(action=\"ignore\", category=FutureWarning)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "def process_hydamo_changes(\n", - " dir_waterschap, dir_hydamo_preprocess, dir_hydamo_changes, dir_hydamo_processed, sel_layers=None\n", - "):\n", - " # process hydamo changes (toevoegen en verwijderen) to new hydamo geopackage\n", - " path_hydamo_gpkg_preprocess = Path(dir_hydamo_preprocess, \"hydamo.gpkg\")\n", - " path_hydamo_gpkg_processed = Path(dir_hydamo_processed, \"hydamo.gpkg\")\n", - " path_hydamo_gpkg_remove = Path(dir_hydamo_changes, \"hydamo_verwijderen.gpkg\")\n", - " path_hydamo_gpkg_add = Path(dir_hydamo_changes, \"hydamo_toevoegen.gpkg\")\n", - "\n", - " if sel_layers is None or sel_layers == []:\n", - " sel_layers = fiona.listlayers(path_hydamo_gpkg_preprocess)\n", - " print(sel_layers)\n", - " for layer in sel_layers:\n", - " if layer == \"layer_styles\":\n", - " continue\n", - " print(f\" - {layer}\")\n", - " # read original hydamo gpkg (from specified region)\n", - " gdf = gpd.read_file(str(path_hydamo_gpkg_preprocess), layer=layer, crs=28992)\n", - "\n", - " # remove objects\n", - " if layer in fiona.listlayers(path_hydamo_gpkg_remove):\n", - " gdf_remove = gpd.read_file(path_hydamo_gpkg_remove, layer=layer, crs=28992)\n", - " try:\n", - " gdf = gdf.loc[~np.isin(gdf[\"code\"], gdf_remove[\"code\"])]\n", - " except KeyError:\n", - " gdf = gdf.loc[~np.isin(gdf[\"globalid\"], gdf_remove[\"globalid\"])]\n", - " # add new objects\n", - " if layer in fiona.listlayers(path_hydamo_gpkg_add):\n", - " gdf_add = gpd.read_file(path_hydamo_gpkg_add, layer=layer, crs=28992)\n", - " gdf_add = gdf_add.to_crs(28992)\n", - " gdf = gdf.to_crs(28992)\n", - " gdf = gpd.GeoDataFrame(pd.concat([gdf, gdf_add])).reset_index()\n", - "\n", - " # save to new hydamo gpkg\n", - " layer_options = \"ASPATIAL_VARIANT=GPKG_ATTRIBUTES\"\n", - " if gdf.geometry.isnull().all():\n", - " gdf.to_file(str(path_hydamo_gpkg_processed), layer=layer, driver=\"GPKG\", layer_options=layer_options)\n", - " else:\n", - " gdf.to_file(str(path_hydamo_gpkg_processed), layer=layer, driver=\"GPKG\")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "main_dir = \"..\\\\Ribasim modeldata\"" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "waterschappen = [\n", - " # \"AaenMaas\",\n", - " # \"BrabantseDelta\",\n", - " \"DeDommel\",\n", - " # \"DrentsOverijsselseDelta\",\n", - " # \"HunzeenAas\",\n", - " # \"Limburg\",\n", - " # \"RijnenIJssel\",\n", - " # \"StichtseRijnlanden\",\n", - " # \"ValleienVeluwe\",\n", - " # \"Vechtstromen\"\n", - "]" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "optional: preprocess the hydro-objects (check and adapt endpoints)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "preprocess_hydroobjects = False" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "if preprocess_hydroobjects:\n", - " for waterschap in waterschappen:\n", - " dir_waterschap = main_dir / waterschap / \"verwerkt\"\n", - " dir_hydamo_preprocess = dir_waterschap / \"2_voorbewerking\"\n", - " dir_hydamo_processed = dir_waterschap / \"4_ribasim\"\n", - "\n", - " hydroobjects = gpd.read_file(Path(dir_hydamo_preprocess, \"hydamo.gpkg\"), layer=\"hydroobject\")\n", - " wfd_lines = gpd.read_file(Path(dir_hydamo_processed, \"krw.gpkg\"), layer=\"krw_line\")\n", - " wfd_polygons = gpd.read_file(Path(dir_hydamo_processed, \"krw.gpkg\"), layer=\"krw_polygon\")\n", - "\n", - " hydroobject_new = preprocess_hydamo_hydroobjects(\n", - " hydroobjects,\n", - " wfd_lines=wfd_lines,\n", - " wfd_polygons=wfd_polygons,\n", - " buffer_distance_endpoints=0.5,\n", - " wfd_id_column=\"owmident\",\n", - " buffer_distance_wfd=10,\n", - " overlap_ratio_wfd=0.9,\n", - " )\n", - "\n", - " hydroobject_new.to_file(Path(dir_hydamo_preprocess, \"hydamo.gpkg\"), layer=\"hydroobject\", driver=\"GPKG\")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "sel_layers = [\n", - " \"hydroobject\",\n", - " # 'stuw',\n", - " # 'gemaal',\n", - " # 'afvoergebiedaanvoergebied',\n", - " # 'pomp',\n", - " # 'peilgebiedvigerend',\n", - " # 'peilgebiedpraktijk',\n", - " # 'streefpeil',\n", - " # 'duikersifonhevel',\n", - " # 'afsluiter',\n", - " # 'sluis',\n", - "]" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "for waterschap in waterschappen:\n", - " print(f\"Waterschap {waterschap}\")\n", - " dir_waterschap = Path(main_dir, waterschap, \"verwerkt\")\n", - " dir_hydamo_preprocess = Path(dir_waterschap, \"2_voorbewerking\")\n", - " dir_hydamo_changes = Path(dir_waterschap, \"3_input\")\n", - " dir_hydamo_processed = Path(dir_waterschap, \"4_ribasim\")\n", - "\n", - " process_hydamo_changes(\n", - " dir_waterschap=dir_waterschap,\n", - " dir_hydamo_preprocess=dir_hydamo_preprocess,\n", - " dir_hydamo_changes=dir_hydamo_changes,\n", - " dir_hydamo_processed=dir_hydamo_processed,\n", - " sel_layers=sel_layers,\n", - " )" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] - } - ], - "metadata": { - "kernelspec": { - "display_name": "ribasim_lumping", - "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": "a036bb1803af6fe22f064fcf42d66cd9fc5247b5d3b121167c30abfc8c1c6b18" - } - } - }, - "nbformat": 4, - "nbformat_minor": 2 -} diff --git a/scripts/notebooks/hydamo_1_process_hydamo_data.py b/scripts/notebooks/hydamo_1_process_hydamo_data.py new file mode 100644 index 0000000..a0503dc --- /dev/null +++ b/scripts/notebooks/hydamo_1_process_hydamo_data.py @@ -0,0 +1,136 @@ +# Verwerk wijzigingen aan HyDAMO geopackages in nieuwe HyDAMO input geopackage +# +# Contactpersoon: Harm Nomden (Sweco) +# +# Laatste update: 15-03-2024 + + +import warnings +from pathlib import Path + +import fiona +import geopandas as gpd +import numpy as np +import pandas as pd +from hydamo_preprocessing.preprocessing import preprocess_hydamo_hydroobjects + +warnings.simplefilter(action="ignore", category=UserWarning) +warnings.simplefilter(action="ignore", category=FutureWarning) + + +def process_hydamo_changes( + dir_waterschap, dir_hydamo_preprocess, dir_hydamo_changes, dir_hydamo_processed, sel_layers=None +): + # process hydamo changes (toevoegen en verwijderen) to new hydamo geopackage + path_hydamo_gpkg_preprocess = Path(dir_hydamo_preprocess, "hydamo.gpkg") + path_hydamo_gpkg_processed = Path(dir_hydamo_processed, "hydamo.gpkg") + path_hydamo_gpkg_remove = Path(dir_hydamo_changes, "hydamo_verwijderen.gpkg") + path_hydamo_gpkg_add = Path(dir_hydamo_changes, "hydamo_toevoegen.gpkg") + + if sel_layers is None or sel_layers == []: + sel_layers = fiona.listlayers(path_hydamo_gpkg_preprocess) + print(sel_layers) + for layer in sel_layers: + if layer == "layer_styles": + continue + print(f" - {layer}") + # read original hydamo gpkg (from specified region) + gdf = gpd.read_file(str(path_hydamo_gpkg_preprocess), layer=layer, crs=28992) + + # remove objects + if layer in fiona.listlayers(path_hydamo_gpkg_remove): + gdf_remove = gpd.read_file(path_hydamo_gpkg_remove, layer=layer, crs=28992) + try: + gdf = gdf.loc[~np.isin(gdf["code"], gdf_remove["code"])] + except KeyError: + gdf = gdf.loc[~np.isin(gdf["globalid"], gdf_remove["globalid"])] + # add new objects + if layer in fiona.listlayers(path_hydamo_gpkg_add): + gdf_add = gpd.read_file(path_hydamo_gpkg_add, layer=layer, crs=28992) + gdf_add = gdf_add.to_crs(28992) + gdf = gdf.to_crs(28992) + gdf = gpd.GeoDataFrame(pd.concat([gdf, gdf_add])).reset_index() + + # save to new hydamo gpkg + layer_options = "ASPATIAL_VARIANT=GPKG_ATTRIBUTES" + if gdf.geometry.isnull().all(): + gdf.to_file(str(path_hydamo_gpkg_processed), layer=layer, driver="GPKG", layer_options=layer_options) + else: + gdf.to_file(str(path_hydamo_gpkg_processed), layer=layer, driver="GPKG") + + +main_dir = "..\\Ribasim modeldata" + + +waterschappen = [ + # "AaenMaas", + # "BrabantseDelta", + "DeDommel", + # "DrentsOverijsselseDelta", + # "HunzeenAas", + # "Limburg", + # "RijnenIJssel", + # "StichtseRijnlanden", + # "ValleienVeluwe", + # "Vechtstromen" +] + + +# optional: preprocess the hydro-objects (check and adapt endpoints) + + +preprocess_hydroobjects = False + + +if preprocess_hydroobjects: + for waterschap in waterschappen: + dir_waterschap = main_dir / waterschap / "verwerkt" + dir_hydamo_preprocess = dir_waterschap / "2_voorbewerking" + dir_hydamo_processed = dir_waterschap / "4_ribasim" + + hydroobjects = gpd.read_file(Path(dir_hydamo_preprocess, "hydamo.gpkg"), layer="hydroobject") + wfd_lines = gpd.read_file(Path(dir_hydamo_processed, "krw.gpkg"), layer="krw_line") + wfd_polygons = gpd.read_file(Path(dir_hydamo_processed, "krw.gpkg"), layer="krw_polygon") + + hydroobject_new = preprocess_hydamo_hydroobjects( + hydroobjects, + wfd_lines=wfd_lines, + wfd_polygons=wfd_polygons, + buffer_distance_endpoints=0.5, + wfd_id_column="owmident", + buffer_distance_wfd=10, + overlap_ratio_wfd=0.9, + ) + + hydroobject_new.to_file(Path(dir_hydamo_preprocess, "hydamo.gpkg"), layer="hydroobject", driver="GPKG") + + +sel_layers = [ + "hydroobject", + # 'stuw', + # 'gemaal', + # 'afvoergebiedaanvoergebied', + # 'pomp', + # 'peilgebiedvigerend', + # 'peilgebiedpraktijk', + # 'streefpeil', + # 'duikersifonhevel', + # 'afsluiter', + # 'sluis', +] + + +for waterschap in waterschappen: + print(f"Waterschap {waterschap}") + dir_waterschap = Path(main_dir, waterschap, "verwerkt") + dir_hydamo_preprocess = Path(dir_waterschap, "2_voorbewerking") + dir_hydamo_changes = Path(dir_waterschap, "3_input") + dir_hydamo_processed = Path(dir_waterschap, "4_ribasim") + + process_hydamo_changes( + dir_waterschap=dir_waterschap, + dir_hydamo_preprocess=dir_hydamo_preprocess, + dir_hydamo_changes=dir_hydamo_changes, + dir_hydamo_processed=dir_hydamo_processed, + sel_layers=sel_layers, + ) diff --git a/scripts/notebooks/hydamo_2_run_ribasim_lumping_waterboard.ipynb b/scripts/notebooks/hydamo_2_run_ribasim_lumping_waterboard.ipynb deleted file mode 100644 index 64c0795..0000000 --- a/scripts/notebooks/hydamo_2_run_ribasim_lumping_waterboard.ipynb +++ /dev/null @@ -1,104 +0,0 @@ -{ - "cells": [ - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "import warnings\n", - "from pathlib import Path\n", - "\n", - "import pandas as pd\n", - "from ribasim_lumping_tools.run_ribasim_lumping_waterboard import run_ribasim_lumping_for_waterboard\n", - "\n", - "warnings.simplefilter(\"ignore\")\n", - "pd.options.mode.chained_assignment = None\n", - "warnings.simplefilter(action=\"ignore\", category=UserWarning)\n", - "warnings.simplefilter(action=\"ignore\", category=FutureWarning)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "base_dir = Path(\"..\\\\..\\\\Ribasim modeldata\")\n", - "dx = 250.0\n", - "\n", - "waterschappen = [\n", - " \"HunzeenAas\",\n", - " \"DrentsOverijsselseDelta\",\n", - " \"Vechtstromen\",\n", - " \"RijnenIJssel\",\n", - " \"ValleienVeluwe\",\n", - " \"StichtseRijnlanden\",\n", - " \"BrabantseDelta\",\n", - " \"DeDommel\",\n", - " \"AaenMaas\",\n", - " \"Limburg\",\n", - "]" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "for waterschap in waterschappen:\n", - " run_ribasim_lumping_for_waterboard(\n", - " base_dir=base_dir,\n", - " waterschap=waterschap,\n", - " dx=dx,\n", - " buffer_distance=1.0,\n", - " assign_unassigned_areas_to_basins=False if waterschap == \"ValleienVeluwe\" else True,\n", - " remove_isolated_basins=True,\n", - " )" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] - } - ], - "metadata": { - "kernelspec": { - "display_name": "ribasim_lumping_venv", - "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" - } - }, - "nbformat": 4, - "nbformat_minor": 2 -} diff --git a/scripts/notebooks/hydamo_2_run_ribasim_lumping_waterboard.py b/scripts/notebooks/hydamo_2_run_ribasim_lumping_waterboard.py new file mode 100644 index 0000000..0a2579e --- /dev/null +++ b/scripts/notebooks/hydamo_2_run_ribasim_lumping_waterboard.py @@ -0,0 +1,38 @@ +import warnings +from pathlib import Path + +import pandas as pd +from ribasim_lumping_tools.run_ribasim_lumping_waterboard import run_ribasim_lumping_for_waterboard + +warnings.simplefilter("ignore") +pd.options.mode.chained_assignment = None +warnings.simplefilter(action="ignore", category=UserWarning) +warnings.simplefilter(action="ignore", category=FutureWarning) + + +base_dir = Path("..\\..\\Ribasim modeldata") +dx = 250.0 + +waterschappen = [ + "HunzeenAas", + "DrentsOverijsselseDelta", + "Vechtstromen", + "RijnenIJssel", + "ValleienVeluwe", + "StichtseRijnlanden", + "BrabantseDelta", + "DeDommel", + "AaenMaas", + "Limburg", +] + + +for waterschap in waterschappen: + run_ribasim_lumping_for_waterboard( + base_dir=base_dir, + waterschap=waterschap, + dx=dx, + buffer_distance=1.0, + assign_unassigned_areas_to_basins=False if waterschap == "ValleienVeluwe" else True, + remove_isolated_basins=True, + ) diff --git a/scripts/notebooks/hydamo_3_create_ribasim_model_networks.ipynb b/scripts/notebooks/hydamo_3_create_ribasim_model_networks.ipynb deleted file mode 100644 index e161cdc..0000000 --- a/scripts/notebooks/hydamo_3_create_ribasim_model_networks.ipynb +++ /dev/null @@ -1,313 +0,0 @@ -{ - "cells": [ - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "import warnings\n", - "from pathlib import Path\n", - "\n", - "import geopandas as gpd\n", - "import pandas as pd\n", - "\n", - "warnings.simplefilter(action=\"ignore\", category=UserWarning)\n", - "warnings.simplefilter(action=\"ignore\", category=FutureWarning)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "%load_ext autoreload\n", - "%autoreload 2" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "def generate_ribasim_network(waterschap_path, split_nodes_type_conversion, split_nodes_id_conversion=None):\n", - " ribasim_network_dir = Path(waterschap_path, \"4_ribasim\")\n", - " ribasim_network_path = Path(ribasim_network_dir, \"ribasim_network.gpkg\")\n", - " Path(ribasim_network_dir, \"foreign_input.gpkg\")\n", - " ribasim_model_path = Path(ribasim_network_dir, \"ribasim_model.gpkg\")\n", - "\n", - " # Ribasim basins\n", - " basins = gpd.read_file(ribasim_network_path, layer=\"basins\")\n", - " basins = basins[[\"basin\", \"basin_node_id\", \"geometry\"]].rename(\n", - " columns={\"basin\": \"meta_id\", \"basin_node_id\": \"node_id\"}\n", - " )\n", - " basins[\"node_type\"] = \"Basin\"\n", - " basin_areas = gpd.read_file(ribasim_network_path, layer=\"basin_areas\")\n", - " basin_areas.basin_node_id = basin_areas.basin_node_id.astype(int)\n", - " basin_areas[\"meta_area\"] = basin_areas.geometry.area\n", - " basins = basins.merge(\n", - " basin_areas[[\"basin_node_id\", \"meta_area\"]].rename(columns={\"basin_node_id\": \"node_id\"}),\n", - " how=\"left\",\n", - " on=\"node_id\",\n", - " )\n", - " basin_areas = basin_areas[[\"basin_node_id\", \"meta_area\", \"geometry\"]]\n", - " basins[\"name\"] = \"Basin\"\n", - "\n", - " # Ribasim boundaries\n", - " boundaries = gpd.read_file(ribasim_network_path, layer=\"boundaries\")\n", - " boundaries[\"node_id\"] = boundaries[\"boundary_id\"]\n", - " boundaries = boundaries[[\"node_id\", \"boundary_id\", \"boundary_name\", \"boundary_type\", \"geometry\"]].rename(\n", - " columns={\"boundary_id\": \"meta_id\", \"boundary_name\": \"meta_name\", \"boundary_type\": \"node_type\"}\n", - " )\n", - " level_boundaries = boundaries.loc[boundaries[\"node_type\"] == \"LevelBoundary\", :].copy()\n", - " level_boundaries[\"meta_water_level\"] = -10.0\n", - "\n", - " boundaries.loc[boundaries[\"node_type\"] == \"FlowBoundary\", :].copy()\n", - "\n", - " # RIBASIM LUMPING split nodes\n", - " split_nodes = gpd.read_file(ribasim_network_path, layer=\"split_nodes\")\n", - " split_nodes = gpd.read_file(ribasim_network_path, layer=\"split_nodes\")\n", - " split_nodes = split_nodes.rename(\n", - " columns={\"object_type\": \"meta_object_type\", \"object_function\": \"meta_object_function\"}\n", - " )\n", - "\n", - " # Conversion split_node from object-type + object-function to node_type\n", - " if \"meta_object_function\" in split_nodes.columns:\n", - " split_nodes = split_nodes.merge(\n", - " split_nodes_type_conversion, how=\"left\", on=[\"meta_object_type\", \"meta_object_function\"]\n", - " )\n", - " else:\n", - " split_nodes = split_nodes.merge(split_nodes_type_conversion_dhydro, how=\"left\", on=\"meta_object_type\")\n", - "\n", - " if isinstance(split_nodes_id_conversion, dict):\n", - " for key, value in split_nodes_id_conversion.items():\n", - " if len(split_nodes[split_nodes[\"split_node_id\"] == key]) == 0:\n", - " print(f\" * split_node type conversion id={key} (type={value}) does not exist\")\n", - " split_nodes.loc[split_nodes[\"split_node_id\"] == key, \"meta_object_function\"] = value\n", - "\n", - " # define final split_nodes\n", - " split_nodes = split_nodes.loc[split_nodes.status & (split_nodes.meta_object_function != \"harde_knip\"), :].copy()\n", - " if \"meta_object_function\" not in split_nodes:\n", - " split_nodes[\"meta_object_function\"] = \"\"\n", - " split_nodes = split_nodes[\n", - " [\n", - " \"split_node\",\n", - " \"split_node_id\",\n", - " \"split_node_node_id\",\n", - " \"node_type\",\n", - " \"meta_object_type\",\n", - " \"meta_object_function\",\n", - " \"geometry\",\n", - " ]\n", - " ]\n", - " split_nodes = split_nodes.rename(\n", - " columns={\"split_node\": \"meta_split_node_id\", \"split_node_id\": \"name\", \"split_node_node_id\": \"node_id\"}\n", - " )\n", - "\n", - " # Combine all nodes\n", - " nodes = gpd.GeoDataFrame(pd.concat([boundaries, split_nodes, basins]), crs=28992).reset_index(drop=True)\n", - "\n", - " # Combine all edges\n", - " basin_connections = gpd.read_file(ribasim_network_path, layer=\"basin_connections\")\n", - " boundary_connections = gpd.read_file(ribasim_network_path, layer=\"boundary_connections\")\n", - "\n", - " edges_columns = [\"from_node_id\", \"to_node_id\", \"connection\", \"geometry\"]\n", - " edges = gpd.GeoDataFrame(\n", - " pd.concat([basin_connections[edges_columns], boundary_connections[edges_columns]]),\n", - " geometry=\"geometry\",\n", - " crs=28992,\n", - " )\n", - "\n", - " edges = edges[[\"from_node_id\", \"to_node_id\", \"geometry\"]]\n", - "\n", - " edges = edges.merge(\n", - " nodes[[\"node_id\", \"node_type\"]].rename(columns={\"node_id\": \"from_node_id\", \"node_type\": \"from_node_type\"}),\n", - " how=\"left\",\n", - " on=\"from_node_id\",\n", - " )\n", - " edges = edges.merge(\n", - " nodes[[\"node_id\", \"node_type\"]].rename(columns={\"node_id\": \"to_node_id\", \"node_type\": \"to_node_type\"}),\n", - " how=\"left\",\n", - " on=\"to_node_id\",\n", - " )\n", - " edges[\"edge_type\"] = \"flow\"\n", - "\n", - " # Export nodes and edges\n", - " nodes.drop_duplicates(keep=\"first\").to_file(ribasim_model_path, layer=\"Node\")\n", - " edges.drop_duplicates(keep=\"first\").to_file(ribasim_model_path, layer=\"Edge\")\n", - "\n", - " print(f\" - no of nodes: {len(nodes)}\")\n", - " print(f\" - no of edges: {len(edges)}\")\n", - " return nodes, edges, split_nodes" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "base_dir = Path(\"..\\\\Ribasim modeldata\")\n", - "\n", - "waterschappen = [\n", - " \"Noorderzijlvest\",\n", - " \"HunzeenAas\",\n", - " \"DrentsOverijsselseDelta\",\n", - " \"Vechtstromen\",\n", - " \"RijnenIJssel\",\n", - " \"ValleienVeluwe\",\n", - " \"StichtseRijnlanden\",\n", - " \"BrabantseDelta\",\n", - " \"DeDommel\",\n", - " \"AaenMaas\",\n", - " \"Limburg\",\n", - "]" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# HYDAMO (10 waterschappen)\n", - "split_nodes_type_conversion_hydamo = pd.DataFrame(\n", - " columns=[\"meta_object_type\", \"meta_object_function\", \"node_type\"],\n", - " data=[\n", - " [\"stuw\", \"\", \"TabulatedRatingCurve\"],\n", - " [\"stuw\", \"afwaterend\", \"TabulatedRatingCurve\"],\n", - " [\"stuw\", \"inlaat\", \"Outlet\"],\n", - " [\"afsluitmiddel\", \"\", \"Outlet\"],\n", - " [\"afsluitmiddel\", \"inlaat\", \"Outlet\"],\n", - " [\"afsluitmiddel\", \"uitlaat\", \"Outlet\"],\n", - " [\"duikersifonhevel\", \"\", \"ManningResistance\"],\n", - " [\"duikersifonhevel\", \"inlaat\", \"Outlet\"],\n", - " [\"duikersifonhevel\", \"afwaterend\", \"TabulatedRatingCurve\"],\n", - " [\"duikersifonhevel\", \"open verbinding\", \"ManningResistance\"],\n", - " [\"openwater\", \"\", \"ManningResistance\"],\n", - " [\"openwater\", \"open verbinding\", \"ManningResistance\"],\n", - " [\"openwater\", \"afwaterend\", \"TabulatedRatingCurve\"],\n", - " [\"gemaal\", \"\", \"Pump\"],\n", - " [\"gemaal\", \"afvoer\", \"Pump\"],\n", - " [\"gemaal\", \"aanvoer\", \"Pump\"],\n", - " [\"gemaal\", \"aanvoer/afvoer\", \"Pump\"],\n", - " [\"sluis\", \"\", \"Outlet\"],\n", - " [\"sluis\", \"schut- en lekverlies\", \"Outlet\"],\n", - " [\"sluis\", \"spui\", \"Outlet\"],\n", - " [\"sluis\", \"keersluis\", \"Outlet\"],\n", - " ],\n", - ")\n", - "\n", - "# DHYDRO (NOORDERZIJLVEST)\n", - "split_nodes_type_conversion_dhydro = pd.DataFrame(\n", - " columns=[\"meta_object_type\", \"node_type\"],\n", - " data=[\n", - " [\"weir\", \"TabulatedRatingCurve\"],\n", - " [\"uniweir\", \"TabulatedRatingCurve\"],\n", - " [\"universalWeir\", \"TabulatedRatingCurve\"],\n", - " [\"pump\", \"Pump\"],\n", - " [\"culvert\", \"ManningResistance\"],\n", - " [\"openwater\", \"ManningResistance\"],\n", - " [\"orifice\", \"Outlet\"],\n", - " ],\n", - ")\n", - "\n", - "te_verwijderen_aanvoergemalen = [\n", - " \"iKGM004\",\n", - " \"iKGM036\",\n", - " \"iKGM069\",\n", - " \"iKGM073\",\n", - " \"iKGM086\",\n", - " \"iKGM101\",\n", - " \"iKGM102\",\n", - " \"iKGM129\",\n", - " \"iKGM157\",\n", - " \"iKGM163\",\n", - " \"iKGM165\",\n", - " \"iKGM189\",\n", - " \"iKGM190\",\n", - " \"iKGM192\",\n", - " \"iKGM194\",\n", - " \"iKGM198\",\n", - " \"iKGM206\",\n", - " \"iKGM214\",\n", - " \"iKGM226\",\n", - " \"iKGM241\",\n", - " \"iKGM248\",\n", - " \"iKGM260\",\n", - " \"iKGM265\",\n", - " \"iKGM295\",\n", - " \"iKGM302\",\n", - " \"iKST0163\",\n", - " \"iKST0470\",\n", - " \"iKST0569\",\n", - " \"iKST0572\",\n", - " \"iKST0624\",\n", - " \"iKST0707\",\n", - " \"iKST6330\",\n", - " \"iKST6352\",\n", - " \"iKST6386\",\n", - " \"iKST6388\",\n", - " \"iKST6415\",\n", - " \"iKST6622\",\n", - " \"iKST9950\",\n", - "]\n", - "split_nodes_id_conversion_dhydro = {g: \"harde_knip\" for g in te_verwijderen_aanvoergemalen}" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "for waterschap in waterschappen:\n", - " print(f\"Waterschap: {waterschap}\")\n", - " waterschap_path = Path(base_dir, waterschap, \"verwerkt\")\n", - " if waterschap == \"Noorderzijlvest\":\n", - " nodes, edges, split_nodes = generate_ribasim_network(\n", - " waterschap_path, split_nodes_type_conversion_dhydro, split_nodes_id_conversion_dhydro\n", - " )\n", - " else:\n", - " nodes, edges, split_nodes = generate_ribasim_network(waterschap_path, split_nodes_type_conversion_hydamo)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] - } - ], - "metadata": { - "kernelspec": { - "display_name": "ribasim_lumping_venv", - "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" - } - }, - "nbformat": 4, - "nbformat_minor": 2 -} diff --git a/scripts/notebooks/hydamo_3_create_ribasim_model_networks.py b/scripts/notebooks/hydamo_3_create_ribasim_model_networks.py new file mode 100644 index 0000000..63bc50a --- /dev/null +++ b/scripts/notebooks/hydamo_3_create_ribasim_model_networks.py @@ -0,0 +1,232 @@ +import warnings +from pathlib import Path + +import geopandas as gpd +import pandas as pd + +warnings.simplefilter(action="ignore", category=UserWarning) +warnings.simplefilter(action="ignore", category=FutureWarning) + + +def generate_ribasim_network(waterschap_path, split_nodes_type_conversion, split_nodes_id_conversion=None): + ribasim_network_dir = Path(waterschap_path, "4_ribasim") + ribasim_network_path = Path(ribasim_network_dir, "ribasim_network.gpkg") + Path(ribasim_network_dir, "foreign_input.gpkg") + ribasim_model_path = Path(ribasim_network_dir, "ribasim_model.gpkg") + + # Ribasim basins + basins = gpd.read_file(ribasim_network_path, layer="basins") + basins = basins[["basin", "basin_node_id", "geometry"]].rename( + columns={"basin": "meta_id", "basin_node_id": "node_id"} + ) + basins["node_type"] = "Basin" + basin_areas = gpd.read_file(ribasim_network_path, layer="basin_areas") + basin_areas.basin_node_id = basin_areas.basin_node_id.astype(int) + basin_areas["meta_area"] = basin_areas.geometry.area + basins = basins.merge( + basin_areas[["basin_node_id", "meta_area"]].rename(columns={"basin_node_id": "node_id"}), + how="left", + on="node_id", + ) + basin_areas = basin_areas[["basin_node_id", "meta_area", "geometry"]] + basins["name"] = "Basin" + + # Ribasim boundaries + boundaries = gpd.read_file(ribasim_network_path, layer="boundaries") + boundaries["node_id"] = boundaries["boundary_id"] + boundaries = boundaries[["node_id", "boundary_id", "boundary_name", "boundary_type", "geometry"]].rename( + columns={"boundary_id": "meta_id", "boundary_name": "meta_name", "boundary_type": "node_type"} + ) + level_boundaries = boundaries.loc[boundaries["node_type"] == "LevelBoundary", :].copy() + level_boundaries["meta_water_level"] = -10.0 + + boundaries.loc[boundaries["node_type"] == "FlowBoundary", :].copy() + + # RIBASIM LUMPING split nodes + split_nodes = gpd.read_file(ribasim_network_path, layer="split_nodes") + split_nodes = gpd.read_file(ribasim_network_path, layer="split_nodes") + split_nodes = split_nodes.rename( + columns={"object_type": "meta_object_type", "object_function": "meta_object_function"} + ) + + # Conversion split_node from object-type + object-function to node_type + if "meta_object_function" in split_nodes.columns: + split_nodes = split_nodes.merge( + split_nodes_type_conversion, how="left", on=["meta_object_type", "meta_object_function"] + ) + else: + split_nodes = split_nodes.merge(split_nodes_type_conversion_dhydro, how="left", on="meta_object_type") + + if isinstance(split_nodes_id_conversion, dict): + for key, value in split_nodes_id_conversion.items(): + if len(split_nodes[split_nodes["split_node_id"] == key]) == 0: + print(f" * split_node type conversion id={key} (type={value}) does not exist") + split_nodes.loc[split_nodes["split_node_id"] == key, "meta_object_function"] = value + + # define final split_nodes + split_nodes = split_nodes.loc[split_nodes.status & (split_nodes.meta_object_function != "harde_knip"), :].copy() + if "meta_object_function" not in split_nodes: + split_nodes["meta_object_function"] = "" + split_nodes = split_nodes[ + [ + "split_node", + "split_node_id", + "split_node_node_id", + "node_type", + "meta_object_type", + "meta_object_function", + "geometry", + ] + ] + split_nodes = split_nodes.rename( + columns={"split_node": "meta_split_node_id", "split_node_id": "name", "split_node_node_id": "node_id"} + ) + + # Combine all nodes + nodes = gpd.GeoDataFrame(pd.concat([boundaries, split_nodes, basins]), crs=28992).reset_index(drop=True) + + # Combine all edges + basin_connections = gpd.read_file(ribasim_network_path, layer="basin_connections") + boundary_connections = gpd.read_file(ribasim_network_path, layer="boundary_connections") + + edges_columns = ["from_node_id", "to_node_id", "connection", "geometry"] + edges = gpd.GeoDataFrame( + pd.concat([basin_connections[edges_columns], boundary_connections[edges_columns]]), + geometry="geometry", + crs=28992, + ) + + edges = edges[["from_node_id", "to_node_id", "geometry"]] + + edges = edges.merge( + nodes[["node_id", "node_type"]].rename(columns={"node_id": "from_node_id", "node_type": "from_node_type"}), + how="left", + on="from_node_id", + ) + edges = edges.merge( + nodes[["node_id", "node_type"]].rename(columns={"node_id": "to_node_id", "node_type": "to_node_type"}), + how="left", + on="to_node_id", + ) + edges["edge_type"] = "flow" + + # Export nodes and edges + nodes.drop_duplicates(keep="first").to_file(ribasim_model_path, layer="Node") + edges.drop_duplicates(keep="first").to_file(ribasim_model_path, layer="Edge") + + print(f" - no of nodes: {len(nodes)}") + print(f" - no of edges: {len(edges)}") + return nodes, edges, split_nodes + + +base_dir = Path("..\\Ribasim modeldata") + +waterschappen = [ + "Noorderzijlvest", + "HunzeenAas", + "DrentsOverijsselseDelta", + "Vechtstromen", + "RijnenIJssel", + "ValleienVeluwe", + "StichtseRijnlanden", + "BrabantseDelta", + "DeDommel", + "AaenMaas", + "Limburg", +] + + +# HYDAMO (10 waterschappen) +split_nodes_type_conversion_hydamo = pd.DataFrame( + columns=["meta_object_type", "meta_object_function", "node_type"], + data=[ + ["stuw", "", "TabulatedRatingCurve"], + ["stuw", "afwaterend", "TabulatedRatingCurve"], + ["stuw", "inlaat", "Outlet"], + ["afsluitmiddel", "", "Outlet"], + ["afsluitmiddel", "inlaat", "Outlet"], + ["afsluitmiddel", "uitlaat", "Outlet"], + ["duikersifonhevel", "", "ManningResistance"], + ["duikersifonhevel", "inlaat", "Outlet"], + ["duikersifonhevel", "afwaterend", "TabulatedRatingCurve"], + ["duikersifonhevel", "open verbinding", "ManningResistance"], + ["openwater", "", "ManningResistance"], + ["openwater", "open verbinding", "ManningResistance"], + ["openwater", "afwaterend", "TabulatedRatingCurve"], + ["gemaal", "", "Pump"], + ["gemaal", "afvoer", "Pump"], + ["gemaal", "aanvoer", "Pump"], + ["gemaal", "aanvoer/afvoer", "Pump"], + ["sluis", "", "Outlet"], + ["sluis", "schut- en lekverlies", "Outlet"], + ["sluis", "spui", "Outlet"], + ["sluis", "keersluis", "Outlet"], + ], +) + +# DHYDRO (NOORDERZIJLVEST) +split_nodes_type_conversion_dhydro = pd.DataFrame( + columns=["meta_object_type", "node_type"], + data=[ + ["weir", "TabulatedRatingCurve"], + ["uniweir", "TabulatedRatingCurve"], + ["universalWeir", "TabulatedRatingCurve"], + ["pump", "Pump"], + ["culvert", "ManningResistance"], + ["openwater", "ManningResistance"], + ["orifice", "Outlet"], + ], +) + +te_verwijderen_aanvoergemalen = [ + "iKGM004", + "iKGM036", + "iKGM069", + "iKGM073", + "iKGM086", + "iKGM101", + "iKGM102", + "iKGM129", + "iKGM157", + "iKGM163", + "iKGM165", + "iKGM189", + "iKGM190", + "iKGM192", + "iKGM194", + "iKGM198", + "iKGM206", + "iKGM214", + "iKGM226", + "iKGM241", + "iKGM248", + "iKGM260", + "iKGM265", + "iKGM295", + "iKGM302", + "iKST0163", + "iKST0470", + "iKST0569", + "iKST0572", + "iKST0624", + "iKST0707", + "iKST6330", + "iKST6352", + "iKST6386", + "iKST6388", + "iKST6415", + "iKST6622", + "iKST9950", +] +split_nodes_id_conversion_dhydro = {g: "harde_knip" for g in te_verwijderen_aanvoergemalen} + + +for waterschap in waterschappen: + print(f"Waterschap: {waterschap}") + waterschap_path = Path(base_dir, waterschap, "verwerkt") + if waterschap == "Noorderzijlvest": + nodes, edges, split_nodes = generate_ribasim_network( + waterschap_path, split_nodes_type_conversion_dhydro, split_nodes_id_conversion_dhydro + ) + else: + nodes, edges, split_nodes = generate_ribasim_network(waterschap_path, split_nodes_type_conversion_hydamo) diff --git a/scripts/notebooks/hydamo_4_create_dummy_ribasim_models.ipynb b/scripts/notebooks/hydamo_4_create_dummy_ribasim_models.ipynb deleted file mode 100644 index 0e1693f..0000000 --- a/scripts/notebooks/hydamo_4_create_dummy_ribasim_models.ipynb +++ /dev/null @@ -1,174 +0,0 @@ -{ - "cells": [ - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "import os\n", - "import warnings\n", - "from pathlib import Path\n", - "\n", - "import fiona\n", - "import geopandas as gpd\n", - "import pandas as pd\n", - "from ribasim_lumping_tools.default_model import DEFAULTS, default_model\n", - "\n", - "warnings.simplefilter(action=\"ignore\", category=UserWarning)\n", - "warnings.simplefilter(action=\"ignore\", category=FutureWarning)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "%load_ext autoreload\n", - "%autoreload 2" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "#### Read RIBASIM model" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "base_dir = Path(\"..\\\\Ribasim modeldata\\\\\")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "waterschappen = {\n", - " \"AaenMaas\": \"NL38\",\n", - " \"BrabantseDelta\": \"NL25\",\n", - " \"DeDommel\": \"NL27\",\n", - " \"DrentsOverijsselseDelta\": \"NL59\",\n", - " \"HunzeenAas\": \"NL33\",\n", - " \"Limburg\": \"NL60\",\n", - " \"Noorderzijlvest\": \"NL34\",\n", - " \"RijnenIJssel\": \"NL7\",\n", - " \"StichtseRijnlanden\": \"NL14\",\n", - " \"ValleienVeluwe\": \"NL8\",\n", - " \"Vechtstromen\": \"NL44\",\n", - "}" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "new_model_dir_string = \"..\\\\modellen\\\\WATERBOARD\\\\modellen\\\\WATERBOARD_2024_6_3\"\n", - "\n", - "for waterschap, waterschap_code in waterschappen.items():\n", - " print(waterschap)\n", - " new_model_dir = Path(new_model_dir_string.replace(\"WATERBOARD\", waterschap))\n", - " print(new_model_dir)\n", - "\n", - " if not new_model_dir.exists():\n", - " os.makedirs(new_model_dir)\n", - "\n", - " # gpkg\n", - " old_ribasim_model_gpkg = Path(base_dir, waterschap, \"verwerkt\", \"4_ribasim\", \"ribasim_model.gpkg\")\n", - " old_krw_gpkg = Path(base_dir, waterschap, \"verwerkt\", \"4_ribasim\", \"krw.gpkg\")\n", - "\n", - " # read nodes\n", - " node_df = gpd.read_file(old_ribasim_model_gpkg, layer=\"Node\", engine=\"pyogrio\", fid_as_index=True)\n", - " node_df = node_df.rename(columns={\"type\": \"node_type\"})\n", - " node_df[\"meta_code\"] = waterschap_code\n", - "\n", - " # read edges\n", - " edge_df = gpd.read_file(old_ribasim_model_gpkg, layer=\"Edge\", engine=\"pyogrio\", fid_as_index=True)\n", - "\n", - " # read basin areas\n", - " basin_areas = gpd.read_file(\n", - " str(old_ribasim_model_gpkg).replace(\"ribasim_model.gpkg\", \"ribasim_network.gpkg\"), layer=\"basin_areas\"\n", - " )\n", - " basin_areas = basin_areas[[\"basin_node_id\", \"geometry\"]].rename(columns={\"basin_node_id\": \"node_id\"})\n", - " basin_areas.node_id = basin_areas.node_id.astype(int)\n", - "\n", - " # read krw\n", - " krw = gpd.GeoDataFrame()\n", - " krw_layers = fiona.listlayers(str(old_krw_gpkg))\n", - " if \"krw_line\" in krw_layers:\n", - " krw_line = gpd.read_file(str(old_krw_gpkg), layer=\"krw_line\").explode(index_parts=True)\n", - " krw_line.geometry = krw_line.geometry.buffer(10, join_style=\"bevel\")\n", - " krw = pd.concat([krw, krw_line])\n", - " if \"krw_vlak\" in krw_layers:\n", - " krw_vlak = gpd.read_file(str(old_krw_gpkg), layer=\"krw_vlak\").explode(index_parts=True)\n", - " krw = pd.concat([krw, krw_vlak])\n", - " krw = krw[[\"owmident\", \"owmnaam\", \"owmtype\", \"geometry\"]].reset_index(drop=True)\n", - " krw.columns = [\"meta_krw_id\", \"meta_krw_name\", \"meta_krw_type\", \"geometry\"]\n", - "\n", - " node_df = (\n", - " node_df.sjoin(krw, how=\"left\").drop(columns=[\"index_right\"]).drop_duplicates(subset=\"node_id\", keep=\"first\")\n", - " )\n", - " node_df[\"meta_categorie\"] = \"doorgaand\"\n", - " node_df.loc[~node_df.meta_krw_id.isna(), \"meta_categorie\"] = \"hoofdwater\"\n", - "\n", - " # create default model\n", - " model = default_model(node_df, edge_df, basin_areas, **DEFAULTS)\n", - "\n", - " # write model to disk\n", - " ribasim_toml = Path(new_model_dir, \"model.toml\")\n", - " model.write(str(ribasim_toml))" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] - } - ], - "metadata": { - "kernelspec": { - "display_name": "ribasim_lumping_venv", - "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" - } - }, - "nbformat": 4, - "nbformat_minor": 2 -} diff --git a/scripts/notebooks/hydamo_4_create_dummy_ribasim_models.py b/scripts/notebooks/hydamo_4_create_dummy_ribasim_models.py new file mode 100644 index 0000000..3f9be52 --- /dev/null +++ b/scripts/notebooks/hydamo_4_create_dummy_ribasim_models.py @@ -0,0 +1,88 @@ +import os +import warnings +from pathlib import Path + +import fiona +import geopandas as gpd +import pandas as pd +from ribasim_lumping_tools.default_model import DEFAULTS, default_model + +warnings.simplefilter(action="ignore", category=UserWarning) +warnings.simplefilter(action="ignore", category=FutureWarning) + + +# #### Read RIBASIM model + + +base_dir = Path("..\\Ribasim modeldata\\") + + +waterschappen = { + "AaenMaas": "NL38", + "BrabantseDelta": "NL25", + "DeDommel": "NL27", + "DrentsOverijsselseDelta": "NL59", + "HunzeenAas": "NL33", + "Limburg": "NL60", + "Noorderzijlvest": "NL34", + "RijnenIJssel": "NL7", + "StichtseRijnlanden": "NL14", + "ValleienVeluwe": "NL8", + "Vechtstromen": "NL44", +} + + +new_model_dir_string = "..\\modellen\\WATERBOARD\\modellen\\WATERBOARD_2024_6_3" + +for waterschap, waterschap_code in waterschappen.items(): + print(waterschap) + new_model_dir = Path(new_model_dir_string.replace("WATERBOARD", waterschap)) + print(new_model_dir) + + if not new_model_dir.exists(): + os.makedirs(new_model_dir) + + # gpkg + old_ribasim_model_gpkg = Path(base_dir, waterschap, "verwerkt", "4_ribasim", "ribasim_model.gpkg") + old_krw_gpkg = Path(base_dir, waterschap, "verwerkt", "4_ribasim", "krw.gpkg") + + # read nodes + node_df = gpd.read_file(old_ribasim_model_gpkg, layer="Node", engine="pyogrio", fid_as_index=True) + node_df = node_df.rename(columns={"type": "node_type"}) + node_df["meta_code"] = waterschap_code + + # read edges + edge_df = gpd.read_file(old_ribasim_model_gpkg, layer="Edge", engine="pyogrio", fid_as_index=True) + + # read basin areas + basin_areas = gpd.read_file( + str(old_ribasim_model_gpkg).replace("ribasim_model.gpkg", "ribasim_network.gpkg"), layer="basin_areas" + ) + basin_areas = basin_areas[["basin_node_id", "geometry"]].rename(columns={"basin_node_id": "node_id"}) + basin_areas.node_id = basin_areas.node_id.astype(int) + + # read krw + krw = gpd.GeoDataFrame() + krw_layers = fiona.listlayers(str(old_krw_gpkg)) + if "krw_line" in krw_layers: + krw_line = gpd.read_file(str(old_krw_gpkg), layer="krw_line").explode(index_parts=True) + krw_line.geometry = krw_line.geometry.buffer(10, join_style="bevel") + krw = pd.concat([krw, krw_line]) + if "krw_vlak" in krw_layers: + krw_vlak = gpd.read_file(str(old_krw_gpkg), layer="krw_vlak").explode(index_parts=True) + krw = pd.concat([krw, krw_vlak]) + krw = krw[["owmident", "owmnaam", "owmtype", "geometry"]].reset_index(drop=True) + krw.columns = ["meta_krw_id", "meta_krw_name", "meta_krw_type", "geometry"] + + node_df = ( + node_df.sjoin(krw, how="left").drop(columns=["index_right"]).drop_duplicates(subset="node_id", keep="first") + ) + node_df["meta_categorie"] = "doorgaand" + node_df.loc[~node_df.meta_krw_id.isna(), "meta_categorie"] = "hoofdwater" + + # create default model + model = default_model(node_df, edge_df, basin_areas, **DEFAULTS) + + # write model to disk + ribasim_toml = Path(new_model_dir, "model.toml") + model.write(str(ribasim_toml)) diff --git a/scripts/notebooks/xxxx_combine_waterschap_layers.ipynb b/scripts/notebooks/xxxx_combine_waterschap_layers.ipynb deleted file mode 100644 index 8e6ce6c..0000000 --- a/scripts/notebooks/xxxx_combine_waterschap_layers.ipynb +++ /dev/null @@ -1,462 +0,0 @@ -{ - "cells": [ - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "from pathlib import Path\n", - "\n", - "import fiona\n", - "import geopandas as gpd\n", - "import matplotlib.pyplot as plt\n", - "import numpy as np\n", - "import pandas as pd\n", - "\n", - "from ribasim_lumping.utils.general_functions import remove_holes_from_polygons" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# locatie van de waterschapsmappen\n", - "base_dir = \"..\\\\Ribasim modeldata\\\\\"\n", - "\n", - "# creeer een lijst met alle namen van de waterschappen\n", - "waterschappen = {\n", - " \"Noorderzijlvest\": \"Noorderzijlvest\",\n", - " \"HunzeenAas\": \"Hunze en Aa's\",\n", - " \"DrentsOverijsselseDelta\": \"Drents Overijsselse Delta\",\n", - " \"Vechtstromen\": \"Vechtstromen\",\n", - " \"RijnenIJssel\": \"Rijn en IJssel\",\n", - " \"ValleienVeluwe\": \"Vallei en Veluwe\",\n", - " \"StichtseRijnlanden\": \"De Stichtse Rijnlanden\",\n", - " \"BrabantseDelta\": \"Brabantse Delta\",\n", - " \"DeDommel\": \"De Dommel\",\n", - " \"AaenMaas\": \"Aa en Maas\",\n", - " \"Limburg\": \"Limburg\",\n", - "}\n", - "\n", - "# lijst met de benodigde layers\n", - "layers = {\n", - " \"basins\": \"ribasim_network.gpkg\",\n", - " \"basin_areas\": \"ribasim_network.gpkg\",\n", - " \"split_nodes\": \"ribasim_network.gpkg\",\n", - " \"boundaries\": \"ribasim_network.gpkg\",\n", - " \"boundary_connections\": \"ribasim_network.gpkg\",\n", - " \"basin_connections\": \"ribasim_network.gpkg\",\n", - " \"areas\": \"areas.gpkg\",\n", - " \"drainage_areas\": \"areas.gpkg\",\n", - " \"foreign_drainage_areas\": \"foreign_input.gpkg\",\n", - " # \"gemaal\": \"hydamo.gpkg\",\n", - " # \"stuw\": \"hydamo.gpkg\",\n", - " # \"onderdoorlaat\": \"hydamo.gpkg\",\n", - " # \"afsluitmiddel\": \"hydamo.gpkg\",\n", - " # \"duikersifonhevel\": \"hydamo.gpkg\",\n", - " # \"hydroobject\": \"hydamo.gpkg\",\n", - "}\n", - "\n", - "output_gpkg = \"data//alle_waterschappen.gpkg\"\n", - "# output_gpkg = \"data//foreign_input.gpkg\"" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# waterschappen_geoms = gpd.read_file(\"data_oud//waterschappen.gpkg\").to_crs(28992)\n", - "waterschappen_labels = list(waterschappen.keys())" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "split_nodes = gpd.read_file(\n", - " Path(base_dir, list(waterschappen.keys())[1], \"verwerkt\", \"4_ribasim\", layers[\"split_nodes\"])\n", - ")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# loop door de verschillende shapefiles die je wilt hebben per waterschap\n", - "for layer, gpkg_file in layers.items():\n", - " print(layer)\n", - " layer_totaal = None\n", - " # loop door de directories van de waterschappen\n", - " print(\" - \", end=\"\")\n", - " for i, waterschap in enumerate(waterschappen):\n", - " print(waterschap[:3], end=\" \")\n", - " # maak de directory\n", - " locatie_gpkg = Path(base_dir, waterschap, \"verwerkt\", \"4_ribasim\", gpkg_file)\n", - " if not locatie_gpkg.exists():\n", - " continue\n", - " if layer not in fiona.listlayers(locatie_gpkg):\n", - " continue\n", - "\n", - " # read the shapefile layers\n", - " layer_data = gpd.read_file(locatie_gpkg, layer=layer, engine=\"pyogrio\")\n", - " if layer == \"areas\":\n", - " layer_data = layer_data[[\"code\", \"geometry\"]]\n", - " if layer == \"foreign_drainage_areas\":\n", - " layer_data = layer_data[[\"name\", \"boundary_name\", \"geometry\"]]\n", - " if layer in [\n", - " \"drainage_areas\",\n", - " \"gemaal\",\n", - " \"stuw\",\n", - " \"afsluitmiddel\",\n", - " \"onderdoorlaat\",\n", - " \"duikersifonhevel\",\n", - " \"hydroobject\",\n", - " ]:\n", - " if \"code\" not in layer_data.columns:\n", - " layer_data[\"code\"] = None\n", - " layer_data = layer_data[[\"code\", \"geometry\"]]\n", - "\n", - " # add waterschap name\n", - " layer_data[\"waterschap\"] = waterschap\n", - "\n", - " layer_data = layer_data.set_crs(28992, allow_override=True)\n", - "\n", - " if layer_totaal is None:\n", - " layer_totaal = layer_data.copy()\n", - " else:\n", - " layer_totaal = pd.concat([layer_totaal, layer_data])\n", - "\n", - " if layer_totaal is not None:\n", - " layer_totaal.to_file(output_gpkg, layer=layer, driver=\"GPKG\")\n", - " print(\" -> saved\")" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Plots" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# load the data\n", - "areas = gpd.read_file(output_gpkg, layer=\"areas\")\n", - "basins = gpd.read_file(output_gpkg, layer=\"basins\")\n", - "basin_areas = gpd.read_file(output_gpkg, layer=\"basin_areas\")\n", - "split_nodes = gpd.read_file(output_gpkg, layer=\"split_nodes\")\n", - "boundaries = gpd.read_file(output_gpkg, layer=\"boundaries\")\n", - "\n", - "drainage_areas = gpd.read_file(output_gpkg, layer=\"drainage_areas\")\n", - "foreign_drainage_areas = gpd.read_file(output_gpkg, layer=\"foreign_drainage_areas\")\n", - "gemaal = gpd.read_file(output_gpkg, layer=\"gemaal\")\n", - "stuw = gpd.read_file(output_gpkg, layer=\"stuw\")\n", - "onderdoorlaat = gpd.read_file(output_gpkg, layer=\"onderdoorlaat\")\n", - "afsluitmiddel = gpd.read_file(output_gpkg, layer=\"afsluitmiddel\")\n", - "duikersifonhevel = gpd.read_file(output_gpkg, layer=\"duikersifonhevel\")\n", - "hydroobject = gpd.read_file(output_gpkg, layer=\"hydroobject\")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# BOUNDARIES: FILL TYPE\n", - "boundaries[\"Type\"] = (\n", - " boundaries[\"Type\"]\n", - " .fillna(boundaries[\"quantity\"])\n", - " .replace({\"dischargebnd\": \"FlowBoundary\", \"waterlevelbnd\": \"LevelBoundary\"})\n", - ")\n", - "# CHECK BOUNDARIES\n", - "boundaries[[\"Type\", \"quantity\", \"waterschap\"]].fillna(\"\").groupby(\n", - " [\"Type\", \"quantity\", \"waterschap\"]\n", - ").size().reset_index() # .rename(columns={0:'count'})\n", - "boundaries.to_file(output_gpkg, layer=\"boundaries\")\n", - "# SEPARATE FLOW AND LEVEL BOUNDARIES\n", - "flow_boundaries = boundaries[boundaries[\"Type\"] == \"FlowBoundary\"]\n", - "level_boundaries = boundaries[boundaries[\"Type\"] == \"LevelBoundary\"]" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# BASIN AREAS\n", - "basin_areas_waterschap = areas.dissolve(by=[\"waterschap\", \"basin_node_id\"])\n", - "basin_areas_waterschap.area = basin_areas_waterschap.geometry.area\n", - "rng = np.random.default_rng()\n", - "basin_areas_waterschap[\"color_no\"] = rng.choice(np.arange(50), size=len(basin_areas_waterschap))" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "basin_areas_waterschap = remove_holes_from_polygons(basin_areas_waterschap.explode(), 100_000)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "basin_areas_waterschap.to_file(output_gpkg, layer=\"basin_areas\")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "basin_areas_waterschap.reset_index().to_file(output_gpkg, layer=\"basin_areas\")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "rng = np.random.default_rng()\n", - "basin_areas_waterschap[\"color_no\"] = rng.choice(np.arange(50), size=len(basin_areas_waterschap))" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# BASIN AREAS\n", - "fig, ax = plt.subplots()\n", - "basin_areas_waterschap.reset_index(drop=True).plot(ax=ax, column=\"color_no\")\n", - "waterschappen.plot(ax=ax, facecolor=\"None\")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# CALCULATE SURFACE AREA OF WATER BOARDS\n", - "areas[\"area\"] = areas.geometry.area / 1_000_000\n", - "areas[[\"area\", \"waterschap\"]].groupby(\"waterschap\").sum()" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# PLOT FOR SURFACE AREA, BOUNDARIES, SPLIT NODES, BASINS, BASIN AREAS\n", - "\n", - "\n", - "def addlabels(ax, x, y):\n", - " for _x, _y in zip(x, y):\n", - " ax.text(_x, _y, _y, ha=\"center\", va=\"bottom\", fontsize=7)\n", - "\n", - "\n", - "# make the plots\n", - "fig, axs = plt.subplots(4, 1, figsize=(5, 7), sharex=True, gridspec_kw={\"hspace\": 0.25, \"wspace\": 0.3})\n", - "# fig.tight_layout()\n", - "\n", - "data_sets = [boundaries, split_nodes, basins, basin_areas]\n", - "columns = [\"Boundaries\", \"Split nodes\", \"Basins\", \"Basin areas\"]\n", - "data_labels = [\"Boundaries\", \"Split nodes\", \"Basins\", \"Basin areas\"]\n", - "\n", - "for data_set, data_label, ax in zip(data_sets, data_labels, axs.flatten()):\n", - " labels, counts = np.unique(data_set.waterschap, return_counts=True)\n", - " counts_def = []\n", - " for w_lab in waterschappen.keys():\n", - " counts_new = 0\n", - " for label, count in zip(labels, counts):\n", - " if label == w_lab:\n", - " counts_new = count\n", - " counts_def += [counts_new]\n", - " ax.bar(waterschappen.values, counts_def, align=\"center\")\n", - " addlabels(ax, waterschappen.values, counts_def)\n", - " ax.set_ylim([0, max(counts_def) * 1.2])\n", - " ax.set_title(data_label, fontsize=10, ha=\"left\", x=-0.1, fontweight=\"bold\")\n", - " ax.tick_params(axis=\"x\", which=\"major\", labelsize=10)\n", - " ax.tick_params(axis=\"y\", which=\"major\", labelsize=9)\n", - "\n", - "basin_areas.area = basin_areas.geometry.area\n", - "basin_areas[\"area_km2\"] = basin_areas.geometry.area / 1000000\n", - "# basin_areas[basin_areas.waterschap==\"Noorderzijlvest\", \"color_no\"] =\n", - "\n", - "ax = axs[-1] # [-1]\n", - "# basin_areas_km2 = basin_areas[[\"waterschap\", \"area_km2\"]].groupby(\"waterschap\").sum().area_km2\n", - "# ax.bar(basin_areas_km2.index, basin_areas_km2.values, align='center')\n", - "# addlabels(ax, basin_areas_km2.index, basin_areas_km2.round(0).values)#((basin_areas_km2/1000).round(0)*1000.0).values)\n", - "# ax.set_ylim([0, basin_areas_km2.max()*1.2])\n", - "# ax.set_ylabel(\"area [km2]\")\n", - "ax.tick_params(axis=\"x\", labelrotation=90)\n", - "ax.set_xticklabels(waterschappen.values);" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# PLOT FOR PUMPS, WEIRS, CULVERTS, HYDROOBJECTS\n", - "\n", - "# make the plots\n", - "fig, axs = plt.subplots(4, 1, figsize=(5, 7), sharex=True, gridspec_kw={\"hspace\": 0.25, \"wspace\": 0.3})\n", - "fig.tight_layout()\n", - "\n", - "waterschap_areas = areas[[\"area\", \"waterschap\"]].groupby(\"waterschap\").sum()\n", - "counts_def = []\n", - "for w_lab in waterschappen.keys():\n", - " counts_new = 0\n", - " for label, count in zip(waterschap_areas.index, waterschap_areas.area.round(0).values):\n", - " if label == w_lab:\n", - " counts_new = count\n", - " counts_def += [int(counts_new)]\n", - "axs[0].bar(waterschappen_labels, counts_def, align=\"center\")\n", - "addlabels(axs[0], waterschappen_labels, counts_def)\n", - "axs[0].set_ylim([0, max(counts_def) * 1.2])\n", - "axs[0].set_title(\"Surface area [km2]\", fontsize=10, ha=\"left\", x=-0.1, fontweight=\"bold\")\n", - "axs[0].tick_params(axis=\"x\", which=\"major\", labelsize=10)\n", - "axs[0].tick_params(axis=\"y\", which=\"major\", labelsize=9)\n", - "\n", - "hydroobject[\"length\"] = hydroobject.geometry.length / 1000\n", - "hydroobject_length = hydroobject[[\"length\", \"waterschap\"]].groupby(\"waterschap\").sum()\n", - "counts_def = []\n", - "for w_lab in waterschappen.keys():\n", - " counts_new = 0\n", - " for label, count in zip(hydroobject_length.index, hydroobject_length.length.round(0).values):\n", - " if label == w_lab:\n", - " counts_new = count\n", - " counts_def += [int(counts_new)]\n", - "axs[1].bar(waterschappen_labels, counts_def, align=\"center\")\n", - "addlabels(axs[1], waterschappen_labels, counts_def)\n", - "axs[1].set_ylim([0, max(counts_def) * 1.2])\n", - "axs[1].set_title(\"Hydro-objects [km]\", fontsize=10, ha=\"left\", x=-0.1, fontweight=\"bold\")\n", - "axs[1].tick_params(axis=\"x\", which=\"major\", labelsize=10)\n", - "axs[1].tick_params(axis=\"y\", which=\"major\", labelsize=9)\n", - "\n", - "afsluitmiddel = pd.concat([afsluitmiddel, onderdoorlaat])\n", - "\n", - "data_sets = [gemaal, stuw]\n", - "columns = [\"Gemaal\", \"Stuw\"]\n", - "data_labels = [\"Pumping stations\", \"Weirs\"]\n", - "\n", - "\n", - "def addlabels(ax, x, y):\n", - " for _x, _y in zip(x, y):\n", - " ax.text(_x, _y, _y, ha=\"center\", va=\"bottom\", fontsize=7)\n", - "\n", - "\n", - "for data_set, data_label, ax in zip(data_sets, data_labels, axs.flatten()[2:]):\n", - " labels, counts = np.unique(data_set.waterschap, return_counts=True)\n", - " counts_def = []\n", - " for w_lab in waterschappen.keys():\n", - " counts_new = 0\n", - " for label, count in zip(labels, counts):\n", - " if label == w_lab:\n", - " counts_new = count\n", - " counts_def += [int(counts_new)]\n", - " ax.bar(waterschappen_labels, counts_def, align=\"center\")\n", - " addlabels(ax, waterschappen_labels, counts_def)\n", - " ax.set_ylim([0, max(counts_def) * 1.2])\n", - " ax.set_title(data_label, fontsize=10, ha=\"left\", x=-0.1, fontweight=\"bold\")\n", - " ax.tick_params(axis=\"x\", which=\"major\", labelsize=10)\n", - " ax.tick_params(axis=\"y\", which=\"major\", labelsize=9)\n", - "\n", - "basin_areas.area = basin_areas.geometry.area\n", - "basin_areas[\"area_km2\"] = basin_areas.geometry.area / 1000000\n", - "\n", - "ax = axs[-1] # [-1]\n", - "ax.tick_params(axis=\"x\", labelrotation=90)\n", - "ax.set_xticklabels(waterschappen_labels);" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "from shapely.geometry import Polygon\n", - "\n", - "\n", - "def remove_small_holes_from_areas(gdf, min_area):\n", - " list_geometry = []\n", - " for polygon in gdf.geometry:\n", - " list_interiors = []\n", - " for interior in polygon.interiors:\n", - " p = Polygon(interior)\n", - " if p.area > min_area:\n", - " list_interiors.append(interior)\n", - " temp_pol = Polygon(polygon.exterior.coords, holes=list_interiors)\n", - " list_geometry.append(temp_pol)\n", - " gdf.geometry = list_geometry\n", - " return gdf\n", - "\n", - "\n", - "drainage_areas = remove_small_holes_from_areas(drainage_areas, 1000.0)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "drainage_areas.to_file(Path(base_dir, \"areas.gpkg\"), layer=\"drainage_areas\", driver=\"GPKG\")" - ] - } - ], - "metadata": { - "kernelspec": { - "display_name": "ribasim_lumping_venv", - "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": "5bd738815a0acbb3ad0a69908638385386c988630a46c6a41055953a8964d49b" - } - } - }, - "nbformat": 4, - "nbformat_minor": 2 -} diff --git a/scripts/notebooks/xxxx_combine_waterschap_layers.py b/scripts/notebooks/xxxx_combine_waterschap_layers.py new file mode 100644 index 0000000..da2d2c0 --- /dev/null +++ b/scripts/notebooks/xxxx_combine_waterschap_layers.py @@ -0,0 +1,314 @@ +from pathlib import Path + +import fiona +import geopandas as gpd +import matplotlib.pyplot as plt +import numpy as np +import pandas as pd +from shapely.geometry import Polygon + +from ribasim_lumping.utils.general_functions import remove_holes_from_polygons + +# locatie van de waterschapsmappen +base_dir = "..\\Ribasim modeldata\\" + +# creeer een lijst met alle namen van de waterschappen +waterschappen = { + "Noorderzijlvest": "Noorderzijlvest", + "HunzeenAas": "Hunze en Aa's", + "DrentsOverijsselseDelta": "Drents Overijsselse Delta", + "Vechtstromen": "Vechtstromen", + "RijnenIJssel": "Rijn en IJssel", + "ValleienVeluwe": "Vallei en Veluwe", + "StichtseRijnlanden": "De Stichtse Rijnlanden", + "BrabantseDelta": "Brabantse Delta", + "DeDommel": "De Dommel", + "AaenMaas": "Aa en Maas", + "Limburg": "Limburg", +} + +# lijst met de benodigde layers +layers = { + "basins": "ribasim_network.gpkg", + "basin_areas": "ribasim_network.gpkg", + "split_nodes": "ribasim_network.gpkg", + "boundaries": "ribasim_network.gpkg", + "boundary_connections": "ribasim_network.gpkg", + "basin_connections": "ribasim_network.gpkg", + "areas": "areas.gpkg", + "drainage_areas": "areas.gpkg", + "foreign_drainage_areas": "foreign_input.gpkg", + # "gemaal": "hydamo.gpkg", + # "stuw": "hydamo.gpkg", + # "onderdoorlaat": "hydamo.gpkg", + # "afsluitmiddel": "hydamo.gpkg", + # "duikersifonhevel": "hydamo.gpkg", + # "hydroobject": "hydamo.gpkg", +} + +output_gpkg = "data//alle_waterschappen.gpkg" +# output_gpkg = "data//foreign_input.gpkg" + + +# waterschappen_geoms = gpd.read_file("data_oud//waterschappen.gpkg").to_crs(28992) +waterschappen_labels = list(waterschappen.keys()) + + +split_nodes = gpd.read_file( + Path(base_dir, list(waterschappen.keys())[1], "verwerkt", "4_ribasim", layers["split_nodes"]) +) + + +# loop door de verschillende shapefiles die je wilt hebben per waterschap +for layer, gpkg_file in layers.items(): + print(layer) + layer_totaal = None + # loop door de directories van de waterschappen + print(" - ", end="") + for i, waterschap in enumerate(waterschappen): + print(waterschap[:3], end=" ") + # maak de directory + locatie_gpkg = Path(base_dir, waterschap, "verwerkt", "4_ribasim", gpkg_file) + if not locatie_gpkg.exists(): + continue + if layer not in fiona.listlayers(locatie_gpkg): + continue + + # read the shapefile layers + layer_data = gpd.read_file(locatie_gpkg, layer=layer, engine="pyogrio") + if layer == "areas": + layer_data = layer_data[["code", "geometry"]] + if layer == "foreign_drainage_areas": + layer_data = layer_data[["name", "boundary_name", "geometry"]] + if layer in [ + "drainage_areas", + "gemaal", + "stuw", + "afsluitmiddel", + "onderdoorlaat", + "duikersifonhevel", + "hydroobject", + ]: + if "code" not in layer_data.columns: + layer_data["code"] = None + layer_data = layer_data[["code", "geometry"]] + + # add waterschap name + layer_data["waterschap"] = waterschap + + layer_data = layer_data.set_crs(28992, allow_override=True) + + if layer_totaal is None: + layer_totaal = layer_data.copy() + else: + layer_totaal = pd.concat([layer_totaal, layer_data]) + + if layer_totaal is not None: + layer_totaal.to_file(output_gpkg, layer=layer, driver="GPKG") + print(" -> saved") + + +# ### Plots + + +# load the data +areas = gpd.read_file(output_gpkg, layer="areas") +basins = gpd.read_file(output_gpkg, layer="basins") +basin_areas = gpd.read_file(output_gpkg, layer="basin_areas") +split_nodes = gpd.read_file(output_gpkg, layer="split_nodes") +boundaries = gpd.read_file(output_gpkg, layer="boundaries") + +drainage_areas = gpd.read_file(output_gpkg, layer="drainage_areas") +foreign_drainage_areas = gpd.read_file(output_gpkg, layer="foreign_drainage_areas") +gemaal = gpd.read_file(output_gpkg, layer="gemaal") +stuw = gpd.read_file(output_gpkg, layer="stuw") +onderdoorlaat = gpd.read_file(output_gpkg, layer="onderdoorlaat") +afsluitmiddel = gpd.read_file(output_gpkg, layer="afsluitmiddel") +duikersifonhevel = gpd.read_file(output_gpkg, layer="duikersifonhevel") +hydroobject = gpd.read_file(output_gpkg, layer="hydroobject") + + +# BOUNDARIES: FILL TYPE +boundaries["Type"] = ( + boundaries["Type"] + .fillna(boundaries["quantity"]) + .replace({"dischargebnd": "FlowBoundary", "waterlevelbnd": "LevelBoundary"}) +) +# CHECK BOUNDARIES +boundaries[["Type", "quantity", "waterschap"]].fillna("").groupby( + ["Type", "quantity", "waterschap"] +).size().reset_index() # .rename(columns={0:'count'}) +boundaries.to_file(output_gpkg, layer="boundaries") +# SEPARATE FLOW AND LEVEL BOUNDARIES +flow_boundaries = boundaries[boundaries["Type"] == "FlowBoundary"] +level_boundaries = boundaries[boundaries["Type"] == "LevelBoundary"] + + +# BASIN AREAS +basin_areas_waterschap = areas.dissolve(by=["waterschap", "basin_node_id"]) +basin_areas_waterschap.area = basin_areas_waterschap.geometry.area +rng = np.random.default_rng() +basin_areas_waterschap["color_no"] = rng.choice(np.arange(50), size=len(basin_areas_waterschap)) + + +basin_areas_waterschap = remove_holes_from_polygons(basin_areas_waterschap.explode(), 100_000) + + +basin_areas_waterschap.to_file(output_gpkg, layer="basin_areas") + + +basin_areas_waterschap.reset_index().to_file(output_gpkg, layer="basin_areas") + + +rng = np.random.default_rng() +basin_areas_waterschap["color_no"] = rng.choice(np.arange(50), size=len(basin_areas_waterschap)) + + +# BASIN AREAS +fig, ax = plt.subplots() +basin_areas_waterschap.reset_index(drop=True).plot(ax=ax, column="color_no") +waterschappen.plot(ax=ax, facecolor="None") + + +# CALCULATE SURFACE AREA OF WATER BOARDS +areas["area"] = areas.geometry.area / 1_000_000 +areas[["area", "waterschap"]].groupby("waterschap").sum() + + +# PLOT FOR SURFACE AREA, BOUNDARIES, SPLIT NODES, BASINS, BASIN AREAS + + +def addlabels(ax, x, y): + for _x, _y in zip(x, y): + ax.text(_x, _y, _y, ha="center", va="bottom", fontsize=7) + + +# make the plots +fig, axs = plt.subplots(4, 1, figsize=(5, 7), sharex=True, gridspec_kw={"hspace": 0.25, "wspace": 0.3}) +# fig.tight_layout() + +data_sets = [boundaries, split_nodes, basins, basin_areas] +columns = ["Boundaries", "Split nodes", "Basins", "Basin areas"] +data_labels = ["Boundaries", "Split nodes", "Basins", "Basin areas"] + +for data_set, data_label, ax in zip(data_sets, data_labels, axs.flatten()): + labels, counts = np.unique(data_set.waterschap, return_counts=True) + counts_def = [] + for w_lab in waterschappen.keys(): + counts_new = 0 + for label, count in zip(labels, counts): + if label == w_lab: + counts_new = count + counts_def += [counts_new] + ax.bar(waterschappen.values, counts_def, align="center") + addlabels(ax, waterschappen.values, counts_def) + ax.set_ylim([0, max(counts_def) * 1.2]) + ax.set_title(data_label, fontsize=10, ha="left", x=-0.1, fontweight="bold") + ax.tick_params(axis="x", which="major", labelsize=10) + ax.tick_params(axis="y", which="major", labelsize=9) + +basin_areas.area = basin_areas.geometry.area +basin_areas["area_km2"] = basin_areas.geometry.area / 1000000 +# basin_areas[basin_areas.waterschap=="Noorderzijlvest", "color_no"] = + +ax = axs[-1] # [-1] +# basin_areas_km2 = basin_areas[["waterschap", "area_km2"]].groupby("waterschap").sum().area_km2 +# ax.bar(basin_areas_km2.index, basin_areas_km2.values, align='center') +# addlabels(ax, basin_areas_km2.index, basin_areas_km2.round(0).values)#((basin_areas_km2/1000).round(0)*1000.0).values) +# ax.set_ylim([0, basin_areas_km2.max()*1.2]) +# ax.set_ylabel("area [km2]") +ax.tick_params(axis="x", labelrotation=90) +ax.set_xticklabels(waterschappen.values) + + +# PLOT FOR PUMPS, WEIRS, CULVERTS, HYDROOBJECTS + +# make the plots +fig, axs = plt.subplots(4, 1, figsize=(5, 7), sharex=True, gridspec_kw={"hspace": 0.25, "wspace": 0.3}) +fig.tight_layout() + +waterschap_areas = areas[["area", "waterschap"]].groupby("waterschap").sum() +counts_def = [] +for w_lab in waterschappen.keys(): + counts_new = 0 + for label, count in zip(waterschap_areas.index, waterschap_areas.area.round(0).values): + if label == w_lab: + counts_new = count + counts_def += [int(counts_new)] +axs[0].bar(waterschappen_labels, counts_def, align="center") +addlabels(axs[0], waterschappen_labels, counts_def) +axs[0].set_ylim([0, max(counts_def) * 1.2]) +axs[0].set_title("Surface area [km2]", fontsize=10, ha="left", x=-0.1, fontweight="bold") +axs[0].tick_params(axis="x", which="major", labelsize=10) +axs[0].tick_params(axis="y", which="major", labelsize=9) + +hydroobject["length"] = hydroobject.geometry.length / 1000 +hydroobject_length = hydroobject[["length", "waterschap"]].groupby("waterschap").sum() +counts_def = [] +for w_lab in waterschappen.keys(): + counts_new = 0 + for label, count in zip(hydroobject_length.index, hydroobject_length.length.round(0).values): + if label == w_lab: + counts_new = count + counts_def += [int(counts_new)] +axs[1].bar(waterschappen_labels, counts_def, align="center") +addlabels(axs[1], waterschappen_labels, counts_def) +axs[1].set_ylim([0, max(counts_def) * 1.2]) +axs[1].set_title("Hydro-objects [km]", fontsize=10, ha="left", x=-0.1, fontweight="bold") +axs[1].tick_params(axis="x", which="major", labelsize=10) +axs[1].tick_params(axis="y", which="major", labelsize=9) + +afsluitmiddel = pd.concat([afsluitmiddel, onderdoorlaat]) + +data_sets = [gemaal, stuw] +columns = ["Gemaal", "Stuw"] +data_labels = ["Pumping stations", "Weirs"] + + +def addlabels(ax, x, y): + for _x, _y in zip(x, y): + ax.text(_x, _y, _y, ha="center", va="bottom", fontsize=7) + + +for data_set, data_label, ax in zip(data_sets, data_labels, axs.flatten()[2:]): + labels, counts = np.unique(data_set.waterschap, return_counts=True) + counts_def = [] + for w_lab in waterschappen.keys(): + counts_new = 0 + for label, count in zip(labels, counts): + if label == w_lab: + counts_new = count + counts_def += [int(counts_new)] + ax.bar(waterschappen_labels, counts_def, align="center") + addlabels(ax, waterschappen_labels, counts_def) + ax.set_ylim([0, max(counts_def) * 1.2]) + ax.set_title(data_label, fontsize=10, ha="left", x=-0.1, fontweight="bold") + ax.tick_params(axis="x", which="major", labelsize=10) + ax.tick_params(axis="y", which="major", labelsize=9) + +basin_areas.area = basin_areas.geometry.area +basin_areas["area_km2"] = basin_areas.geometry.area / 1000000 + +ax = axs[-1] # [-1] +ax.tick_params(axis="x", labelrotation=90) +ax.set_xticklabels(waterschappen_labels) + + +def remove_small_holes_from_areas(gdf, min_area): + list_geometry = [] + for polygon in gdf.geometry: + list_interiors = [] + for interior in polygon.interiors: + p = Polygon(interior) + if p.area > min_area: + list_interiors.append(interior) + temp_pol = Polygon(polygon.exterior.coords, holes=list_interiors) + list_geometry.append(temp_pol) + gdf.geometry = list_geometry + return gdf + + +drainage_areas = remove_small_holes_from_areas(drainage_areas, 1000.0) + + +drainage_areas.to_file(Path(base_dir, "areas.gpkg"), layer="drainage_areas", driver="GPKG")