Bouw Ribasim schematisatie¶
Met deze Notebook bouwen we een Ribasim schematisatie op basis van het netwerk van het Landelijk Hydrologisch Model (LHM).
De werkwijze vastgelegd in hoofdstukken:
- Inlezen van LSW en DM knopen, en profielen, als Ribasim Basins
- Het bouwen van het Ribasim netwerk (nodes en edges) tussen Ribasim Basins
- Het leggen van boundaries aan de rand van het Ribasim netwerk
- Het wegschrijven van het Ribasim netwerk
In [1]:
Copied!
from config import LHM_DIR, DATA_DIR, MODEL_DIR, load_src
import geopandas as gpd
import pandas as pd
import xarray as xr
from shapely.geometry import Point, LineString
import ribasim
load_src()
from lhm.network import Network
from lhm.utils import next_index
from lhm.read import read_dm_nds
from lhm.utils import report_progress
from config import LHM_DIR, DATA_DIR, MODEL_DIR, load_src
import geopandas as gpd
import pandas as pd
import xarray as xr
from shapely.geometry import Point, LineString
import ribasim
load_src()
from lhm.network import Network
from lhm.utils import next_index
from lhm.read import read_dm_nds
from lhm.utils import report_progress
Modelnaam en clip-mask¶
Hieronder hebben we een optie om een model te maken van heel Nederland of een stukje Nederland:
model_name
: naam voor het Ribasim modelmask_poly
: GeoPackage met polygon
In [2]:
Copied!
model_name="ribasim_model_flevoland"
mask_gdf = gpd.read_file(DATA_DIR / "mask.gpkg")
mask_poly = mask_gdf.iloc[0].geometry
bbox = mask_poly.bounds
#bbox = None
model_name="ribasim_model_flevoland"
mask_gdf = gpd.read_file(DATA_DIR / "mask.gpkg")
mask_poly = mask_gdf.iloc[0].geometry
bbox = mask_poly.bounds
#bbox = None
1. Inlezen LHM-netwerk¶
In [3]:
Copied!
lhm_netwerk = Network.from_gpkg(DATA_DIR / "lhm-netwerk.gpkg", bbox=bbox)
lhm_netwerk.to_gpkg(DATA_DIR/ f"lhm_{model_name}.gpkg")
lhm_netwerk.explore()
lhm_netwerk = Network.from_gpkg(DATA_DIR / "lhm-netwerk.gpkg", bbox=bbox)
lhm_netwerk.to_gpkg(DATA_DIR/ f"lhm_{model_name}.gpkg")
lhm_netwerk.explore()
Out[3]:
Make this Notebook Trusted to load map: File -> Trust Notebook
2. Basin-knopen¶
We lezen de DM en LSW "knopen" in als Basins (1
) en maken de basin-profielen voor LSWS (2
) en DM (3
). We maken hiervoor een GeoDataFrame (4
) dat RIBASIM in kan.
2.1. Basin knopen¶
In [4]:
Copied!
# Inlezen knopen uit netwerk en simplied mozart SAQH tables
nodes_gdf = lhm_netwerk.nodes.copy()
input_mozart_ds = xr.open_dataset(DATA_DIR / r"ribasim_testmodel/simplified_SAQh.nc")
# knopen filteren op beschikbare data in mozart simplified
lsw_nodes_mask = nodes_gdf["index"].str.startswith("MZlsw_")
lsw_table_indices =[f"MZlsw_{i}" for i in input_mozart_ds["node"].values.astype(int)]
nodes_gdf = nodes_gdf[~lsw_nodes_mask | nodes_gdf["index"].isin(lsw_table_indices)]
nodes_gdf["lhm_index"] = nodes_gdf["index"]
nodes_gdf["type"] = "Basin"
nodes_gdf.reset_index(inplace=True, drop=True)
nodes_gdf.index = nodes_gdf.index + 1
nodes_gdf["index"] = nodes_gdf.index
# tonen aan de gebruiker
nodes_gdf
# Inlezen knopen uit netwerk en simplied mozart SAQH tables
nodes_gdf = lhm_netwerk.nodes.copy()
input_mozart_ds = xr.open_dataset(DATA_DIR / r"ribasim_testmodel/simplified_SAQh.nc")
# knopen filteren op beschikbare data in mozart simplified
lsw_nodes_mask = nodes_gdf["index"].str.startswith("MZlsw_")
lsw_table_indices =[f"MZlsw_{i}" for i in input_mozart_ds["node"].values.astype(int)]
nodes_gdf = nodes_gdf[~lsw_nodes_mask | nodes_gdf["index"].isin(lsw_table_indices)]
nodes_gdf["lhm_index"] = nodes_gdf["index"]
nodes_gdf["type"] = "Basin"
nodes_gdf.reset_index(inplace=True, drop=True)
nodes_gdf.index = nodes_gdf.index + 1
nodes_gdf["index"] = nodes_gdf.index
# tonen aan de gebruiker
nodes_gdf
Out[4]:
index | geometry | lhm_index | type | |
---|---|---|---|---|
1 | 1 | POINT (140573.890 476409.651) | MZlsw_21013 | Basin |
2 | 2 | POINT (142103.697 476901.516) | MZlsw_21011 | Basin |
3 | 3 | POINT (138645.749 477151.537) | MZlsw_21012 | Basin |
4 | 4 | POINT (144905.649 477910.142) | MZlsw_21010 | Basin |
5 | 5 | POINT (142435.955 478818.415) | MZlsw_21006 | Basin |
... | ... | ... | ... | ... |
134 | 134 | POINT (179324.333 504346.375) | MZlsw_260077 | Basin |
135 | 135 | POINT (177279.558 504319.119) | MZlsw_260078 | Basin |
136 | 136 | POINT (174454.178 507388.300) | MZlsw_260094 | Basin |
137 | 137 | POINT (176609.670 510253.423) | MZlsw_260093 | Basin |
138 | 138 | POINT (172140.739 512527.561) | MZlsw_260099 | Basin |
138 rows × 4 columns
2.2. LSW SAQ(h) profielen inlezen¶
In [5]:
Copied!
# tabellen toevoegen
da = input_mozart_ds.profile.transpose("node", "profile_col", "profile_row")
lsw_basin_profiles_df = pd.DataFrame(
[item.T for sublist in da.values for item in sublist.T],
columns = ["storage", "area","discharge", "level"]
)
lsw_basin_profiles_df["lhm_index"] = [f"MZlsw_{str(int(x))}" for x in da.node.values for _ in range(4)]
lsw_basin_profiles_df["remarks"] = "uit simplified_SAQh.nc"
lsw_basin_profiles_df = lsw_basin_profiles_df[lsw_basin_profiles_df["lhm_index"].isin(nodes_gdf["lhm_index"])]
lsw_basin_profiles_df
# tabellen toevoegen
da = input_mozart_ds.profile.transpose("node", "profile_col", "profile_row")
lsw_basin_profiles_df = pd.DataFrame(
[item.T for sublist in da.values for item in sublist.T],
columns = ["storage", "area","discharge", "level"]
)
lsw_basin_profiles_df["lhm_index"] = [f"MZlsw_{str(int(x))}" for x in da.node.values for _ in range(4)]
lsw_basin_profiles_df["remarks"] = "uit simplified_SAQh.nc"
lsw_basin_profiles_df = lsw_basin_profiles_df[lsw_basin_profiles_df["lhm_index"].isin(nodes_gdf["lhm_index"])]
lsw_basin_profiles_df
Out[5]:
storage | area | discharge | level | lhm_index | remarks | |
---|---|---|---|---|---|---|
2348 | 0.000000 | 68056.046875 | 0.000000 | -1.566615 | MZlsw_20024 | uit simplified_SAQh.nc |
2349 | 9365.777344 | 68056.046875 | 0.000000 | -1.366615 | MZlsw_20024 | uit simplified_SAQh.nc |
2350 | 18731.554688 | 68056.046875 | 0.024826 | -1.166615 | MZlsw_20024 | uit simplified_SAQh.nc |
2351 | 945943.500000 | 68056.046875 | 2.482551 | 18.633385 | MZlsw_20024 | uit simplified_SAQh.nc |
2352 | 0.000000 | 47411.921875 | 0.000000 | -0.968000 | MZlsw_20025 | uit simplified_SAQh.nc |
... | ... | ... | ... | ... | ... | ... |
31655 | 852750.562500 | 47374.812500 | 1.131316 | 14.410000 | MZlsw_260115 | uit simplified_SAQh.nc |
31656 | 0.000000 | 9373.941406 | 0.000000 | -4.593889 | MZlsw_260116 | uit simplified_SAQh.nc |
31657 | 1344.484253 | 9373.941406 | 0.000000 | -4.393889 | MZlsw_260116 | uit simplified_SAQh.nc |
31658 | 2688.968506 | 9373.941406 | 0.015075 | -4.193889 | MZlsw_260116 | uit simplified_SAQh.nc |
31659 | 135792.906250 | 9373.941406 | 1.507501 | 15.606111 | MZlsw_260116 | uit simplified_SAQh.nc |
524 rows × 6 columns
2.3. DM profielen inlezen¶
In [6]:
Copied!
# % uit DM nds file
def default_profile():
return pd.DataFrame(
data = [[-5, 0 ,0, None], [5, 1000000 , 10000000, None]],
columns= ["level", "area", "storage", "remarks"]
)
def lav(row):
dm_id = row.lhm_index[5:]
if dm_id not in dm_nds_df.index:
profile = default_profile()
profile["remarks"] = "default-profile"
else:
node_row = dm_nds_df.loc[dm_id]
if node_row.lav is None:
profile = default_profile()
remarks = []
if not (node_row.ar == 0.):
if not pd.isna(node_row.ar):
profile["area"] = node_row.ar
remarks += ["constant area" ]
if not (node_row.vo == 0.):
if not pd.isna(node_row.vo):
profile["storage"] = node_row.vo
remarks += ["constant volume"]
if remarks:
profile["remarks"] = ",".join(remarks)
else:
profile["remarks"] = "default-profile"
else:
profile = node_row.lav.copy(deep=True)
profile.index.name= "level"
profile.reset_index(inplace=True)
profile["remarks"] = "uit dm nds.txt"
profile["lhm_index"] = f"{row.lhm_index}"
return profile
dm_nds_df = read_dm_nds(LHM_DIR / r"dm/txtfiles_git/nds.txt").set_index("id")
dm_nodes_gdf = nodes_gdf[nodes_gdf["lhm_index"].str.startswith("DMnd_")]
dm_basin_profiles_df = pd.concat([lav(row) for row in dm_nodes_gdf.itertuples()]).reset_index(drop=True)
dm_basin_profiles_df
# % uit DM nds file
def default_profile():
return pd.DataFrame(
data = [[-5, 0 ,0, None], [5, 1000000 , 10000000, None]],
columns= ["level", "area", "storage", "remarks"]
)
def lav(row):
dm_id = row.lhm_index[5:]
if dm_id not in dm_nds_df.index:
profile = default_profile()
profile["remarks"] = "default-profile"
else:
node_row = dm_nds_df.loc[dm_id]
if node_row.lav is None:
profile = default_profile()
remarks = []
if not (node_row.ar == 0.):
if not pd.isna(node_row.ar):
profile["area"] = node_row.ar
remarks += ["constant area" ]
if not (node_row.vo == 0.):
if not pd.isna(node_row.vo):
profile["storage"] = node_row.vo
remarks += ["constant volume"]
if remarks:
profile["remarks"] = ",".join(remarks)
else:
profile["remarks"] = "default-profile"
else:
profile = node_row.lav.copy(deep=True)
profile.index.name= "level"
profile.reset_index(inplace=True)
profile["remarks"] = "uit dm nds.txt"
profile["lhm_index"] = f"{row.lhm_index}"
return profile
dm_nds_df = read_dm_nds(LHM_DIR / r"dm/txtfiles_git/nds.txt").set_index("id")
dm_nodes_gdf = nodes_gdf[nodes_gdf["lhm_index"].str.startswith("DMnd_")]
dm_basin_profiles_df = pd.concat([lav(row) for row in dm_nodes_gdf.itertuples()]).reset_index(drop=True)
dm_basin_profiles_df
Out[6]:
level | area | storage | remarks | lhm_index | |
---|---|---|---|---|---|
0 | -5 | 0.0 | 0.000000e+00 | default-profile | DMnd_591 |
1 | 5 | 1000000.0 | 1.000000e+07 | default-profile | DMnd_591 |
2 | -5 | 0.0 | 0.000000e+00 | default-profile | DMnd_592 |
3 | 5 | 1000000.0 | 1.000000e+07 | default-profile | DMnd_592 |
4 | -1.0 | 660600000.0 | 2.001700e+09 | uit dm nds.txt | DMnd_6058 |
5 | -0.6 | 660600000.0 | 2.266000e+09 | uit dm nds.txt | DMnd_6058 |
6 | -0.25 | 660600000.0 | 2.497200e+09 | uit dm nds.txt | DMnd_6058 |
7 | 0.25 | 660600000.0 | 2.827500e+09 | uit dm nds.txt | DMnd_6058 |
8 | 1.00 | 660600000.0 | 3.322900e+09 | uit dm nds.txt | DMnd_6058 |
9 | -5 | 0.0 | 0.000000e+00 | default-profile | DMnd_6058_1 |
10 | 5 | 1000000.0 | 1.000000e+07 | default-profile | DMnd_6058_1 |
11 | -1.0 | 73500000.0 | 2.227000e+08 | uit dm nds.txt | DMnd_60581 |
12 | -0.6 | 73500000.0 | 2.520000e+08 | uit dm nds.txt | DMnd_60581 |
13 | -0.25 | 73500000.0 | 2.778000e+08 | uit dm nds.txt | DMnd_60581 |
14 | 0.25 | 73500000.0 | 3.145000e+08 | uit dm nds.txt | DMnd_60581 |
15 | 1.00 | 73500000.0 | 3.697000e+08 | uit dm nds.txt | DMnd_60581 |
16 | -1.0 | 62700000.0 | 3.225000e+07 | uit dm nds.txt | DMnd_6059 |
17 | -0.6 | 62700000.0 | 5.733000e+07 | uit dm nds.txt | DMnd_6059 |
18 | -0.1 | 62700000.0 | 8.868000e+07 | uit dm nds.txt | DMnd_6059 |
19 | 0.40 | 62700000.0 | 1.200000e+08 | uit dm nds.txt | DMnd_6059 |
20 | 1.00 | 62700000.0 | 1.576000e+08 | uit dm nds.txt | DMnd_6059 |
21 | -5 | 0.0 | 0.000000e+00 | default-profile | DMnd_6059_1 |
22 | 5 | 1000000.0 | 1.000000e+07 | default-profile | DMnd_6059_1 |
2.4 Profielen combineren¶
In [7]:
Copied!
# Samenvoegen
basin_profiles_gdf = gpd.GeoDataFrame(
pd.concat(
[lsw_basin_profiles_df,
dm_basin_profiles_df],
ignore_index=True
),
geometry=gpd.GeoSeries(),
crs=28992
)
# DataType goed zetten
basin_profiles_gdf["level"] = basin_profiles_gdf["level"].astype(float)
basin_profiles_gdf["node_id"] = basin_profiles_gdf["lhm_index"].apply(lambda x: nodes_gdf[nodes_gdf["lhm_index"] == x].index.values[0])
# Tonen aan gebruiker
basin_profiles_gdf
# Samenvoegen
basin_profiles_gdf = gpd.GeoDataFrame(
pd.concat(
[lsw_basin_profiles_df,
dm_basin_profiles_df],
ignore_index=True
),
geometry=gpd.GeoSeries(),
crs=28992
)
# DataType goed zetten
basin_profiles_gdf["level"] = basin_profiles_gdf["level"].astype(float)
basin_profiles_gdf["node_id"] = basin_profiles_gdf["lhm_index"].apply(lambda x: nodes_gdf[nodes_gdf["lhm_index"] == x].index.values[0])
# Tonen aan gebruiker
basin_profiles_gdf
Out[7]:
storage | area | discharge | level | lhm_index | remarks | geometry | node_id | |
---|---|---|---|---|---|---|---|---|
0 | 0.000000e+00 | 6.805605e+04 | 0.000000 | -1.566615 | MZlsw_20024 | uit simplified_SAQh.nc | None | 27 |
1 | 9.365777e+03 | 6.805605e+04 | 0.000000 | -1.366615 | MZlsw_20024 | uit simplified_SAQh.nc | None | 27 |
2 | 1.873155e+04 | 6.805605e+04 | 0.024826 | -1.166615 | MZlsw_20024 | uit simplified_SAQh.nc | None | 27 |
3 | 9.459435e+05 | 6.805605e+04 | 2.482551 | 18.633385 | MZlsw_20024 | uit simplified_SAQh.nc | None | 27 |
4 | 0.000000e+00 | 4.741192e+04 | 0.000000 | -0.968000 | MZlsw_20025 | uit simplified_SAQh.nc | None | 7 |
... | ... | ... | ... | ... | ... | ... | ... | ... |
542 | 8.868000e+07 | 6.270000e+07 | NaN | -0.100000 | DMnd_6059 | uit dm nds.txt | None | 121 |
543 | 1.200000e+08 | 6.270000e+07 | NaN | 0.400000 | DMnd_6059 | uit dm nds.txt | None | 121 |
544 | 1.576000e+08 | 6.270000e+07 | NaN | 1.000000 | DMnd_6059 | uit dm nds.txt | None | 121 |
545 | 0.000000e+00 | 0.000000e+00 | NaN | -5.000000 | DMnd_6059_1 | default-profile | None | 122 |
546 | 1.000000e+07 | 1.000000e+06 | NaN | 5.000000 | DMnd_6059_1 | default-profile | None | 122 |
547 rows × 8 columns
3. Bouw netwerk¶
3.1. Tussen LHM links¶
In [8]:
Copied!
links_gdf = lhm_netwerk.links.copy()
links_gdf = links_gdf[links_gdf["node_from"].isin(nodes_gdf["lhm_index"]) & links_gdf["node_to"].isin(nodes_gdf["lhm_index"])]
passed_indexes = []
tabulated_profiles_df = pd.DataFrame(columns=["lhm_index", "level", "discharge", "remarks", "node_id"])
nodes = []
edges = []
manning_resistance = []
fractional_flow = []
lhm_nodes_gdf = nodes_gdf.reset_index().set_index("lhm_index")
lhm_indices = {v:k for k,v in nodes_gdf["lhm_index"].to_dict().items()}
def get_flow_point(geometries, node_type):
if len(geometries) > 1:
point = geometries.iloc[0].boundary.geoms[0]
point = Point(point.x, point.y - 100)
elif node_type == "TabulatedRatingCurve":
line = geometries.iloc[0]
point = line.interpolate(min(100, line.length / 2))
else:
line = geometries.iloc[0]
point = line.centroid
return point
nbr_groups = len(list(links_gdf.groupby("node_from")))
for idx, (lhm_id, df) in enumerate(links_gdf.groupby("node_from")):
report_progress(idx+1, nbr_groups, interval=5, print_at_interval=True)
if lhm_id in lhm_nodes_gdf.index:
#print(f"links vanaf {lhm_id}")
# inlezen van de eigenschappen van de start-knoop
node_from = lhm_nodes_gdf.loc[lhm_id]
# Als er een QH-relatie is, dan pakken we die
if lhm_id in lsw_table_indices:
# we make sure we get a new, available index.
tr_index = next_index(list(lhm_indices.values()))
lhm_indices[f"tr_{tr_index}"] = tr_index
# we add the node
node_type = "TabulatedRatingCurve"
tr_point = get_flow_point(df.geometry, node_type)
nodes += [(tr_index, lhm_id, node_type, tr_point)]
# we add the profile
prof_df = basin_profiles_gdf[basin_profiles_gdf["lhm_index"] == lhm_id][["level","discharge"]]
prof_df["lhm_id"] = lhm_id
prof_df["node_id"] = int(tr_index)
tabulated_profiles_df = pd.concat([tabulated_profiles_df, prof_df])
# we define a edge between the upstream node and
edges += [(node_from["index"], tr_index, LineString([node_from.geometry, tr_point]))]
else: # we don't need an extra node
node_type = "ManningResistance"
# add fraction or manning nodes
if len(df) > 1:
fraction = 1 / len(df)
for row in df.itertuples():
# we get, or make a node_to index
node_to = lhm_nodes_gdf.loc[row.node_to]
to_index = node_to["index"]
if node_type == "TabulatedRatingCurve":
if len(df) > 1:
frac_point = row.geometry.centroid
frac_index = next_index(list(lhm_indices.values()))
lhm_indices[f"frac_{frac_index}"] = frac_index
nodes += [(frac_index, lhm_id, "FractionalFlow", frac_point)]
fractional_flow += [(frac_index, fraction)]
edges += [
(tr_index, frac_index, LineString([tr_point, frac_point])),
(frac_index, to_index, LineString([frac_point, node_to.geometry]))
]
else:
edges += [(tr_index, to_index, LineString([tr_point, node_to.geometry]))]
else:
manning_point = row.geometry.centroid
manning_index = next_index(list(lhm_indices.values()))
lhm_indices[f"manning_{manning_index}"] = manning_index
nodes += [(manning_index, lhm_id, "ManningResistance", manning_point)]
manning_resistance += [(manning_index, 1000,0.04,50,2)]
edges += [
(node_from["index"], manning_index, LineString([node_from.geometry, manning_point])),
(manning_index, to_index, LineString([manning_point, node_to.geometry]))
]
tabulated_profiles_df["remarks"] = "uit simplified_SAQh.nc"
nodes_gdf = pd.concat(
[nodes_gdf,
gpd.GeoDataFrame(
nodes,
columns=["index","lhm_index","type","geometry"],
crs=28992
).set_index("index")]
)
edges_gdf = gpd.GeoDataFrame(edges, columns=["from_node_id","to_node_id","geometry"], crs=28992)
links_gdf = lhm_netwerk.links.copy()
links_gdf = links_gdf[links_gdf["node_from"].isin(nodes_gdf["lhm_index"]) & links_gdf["node_to"].isin(nodes_gdf["lhm_index"])]
passed_indexes = []
tabulated_profiles_df = pd.DataFrame(columns=["lhm_index", "level", "discharge", "remarks", "node_id"])
nodes = []
edges = []
manning_resistance = []
fractional_flow = []
lhm_nodes_gdf = nodes_gdf.reset_index().set_index("lhm_index")
lhm_indices = {v:k for k,v in nodes_gdf["lhm_index"].to_dict().items()}
def get_flow_point(geometries, node_type):
if len(geometries) > 1:
point = geometries.iloc[0].boundary.geoms[0]
point = Point(point.x, point.y - 100)
elif node_type == "TabulatedRatingCurve":
line = geometries.iloc[0]
point = line.interpolate(min(100, line.length / 2))
else:
line = geometries.iloc[0]
point = line.centroid
return point
nbr_groups = len(list(links_gdf.groupby("node_from")))
for idx, (lhm_id, df) in enumerate(links_gdf.groupby("node_from")):
report_progress(idx+1, nbr_groups, interval=5, print_at_interval=True)
if lhm_id in lhm_nodes_gdf.index:
#print(f"links vanaf {lhm_id}")
# inlezen van de eigenschappen van de start-knoop
node_from = lhm_nodes_gdf.loc[lhm_id]
# Als er een QH-relatie is, dan pakken we die
if lhm_id in lsw_table_indices:
# we make sure we get a new, available index.
tr_index = next_index(list(lhm_indices.values()))
lhm_indices[f"tr_{tr_index}"] = tr_index
# we add the node
node_type = "TabulatedRatingCurve"
tr_point = get_flow_point(df.geometry, node_type)
nodes += [(tr_index, lhm_id, node_type, tr_point)]
# we add the profile
prof_df = basin_profiles_gdf[basin_profiles_gdf["lhm_index"] == lhm_id][["level","discharge"]]
prof_df["lhm_id"] = lhm_id
prof_df["node_id"] = int(tr_index)
tabulated_profiles_df = pd.concat([tabulated_profiles_df, prof_df])
# we define a edge between the upstream node and
edges += [(node_from["index"], tr_index, LineString([node_from.geometry, tr_point]))]
else: # we don't need an extra node
node_type = "ManningResistance"
# add fraction or manning nodes
if len(df) > 1:
fraction = 1 / len(df)
for row in df.itertuples():
# we get, or make a node_to index
node_to = lhm_nodes_gdf.loc[row.node_to]
to_index = node_to["index"]
if node_type == "TabulatedRatingCurve":
if len(df) > 1:
frac_point = row.geometry.centroid
frac_index = next_index(list(lhm_indices.values()))
lhm_indices[f"frac_{frac_index}"] = frac_index
nodes += [(frac_index, lhm_id, "FractionalFlow", frac_point)]
fractional_flow += [(frac_index, fraction)]
edges += [
(tr_index, frac_index, LineString([tr_point, frac_point])),
(frac_index, to_index, LineString([frac_point, node_to.geometry]))
]
else:
edges += [(tr_index, to_index, LineString([tr_point, node_to.geometry]))]
else:
manning_point = row.geometry.centroid
manning_index = next_index(list(lhm_indices.values()))
lhm_indices[f"manning_{manning_index}"] = manning_index
nodes += [(manning_index, lhm_id, "ManningResistance", manning_point)]
manning_resistance += [(manning_index, 1000,0.04,50,2)]
edges += [
(node_from["index"], manning_index, LineString([node_from.geometry, manning_point])),
(manning_index, to_index, LineString([manning_point, node_to.geometry]))
]
tabulated_profiles_df["remarks"] = "uit simplified_SAQh.nc"
nodes_gdf = pd.concat(
[nodes_gdf,
gpd.GeoDataFrame(
nodes,
columns=["index","lhm_index","type","geometry"],
crs=28992
).set_index("index")]
)
edges_gdf = gpd.GeoDataFrame(edges, columns=["from_node_id","to_node_id","geometry"], crs=28992)
[####################] 100.0% completed
3.1. Toevoegen boundaries¶
In [9]:
Copied!
level_boundary = []
linear_resistance = []
nodes= []
edges = []
level_boundary_mask = nodes_gdf.index.isin(edges_gdf.to_node_id) & ~nodes_gdf.index.isin(edges_gdf.from_node_id)
for row in nodes_gdf.loc[level_boundary_mask].itertuples():
resist_point = Point(row.geometry.x, row.geometry.y + 50)
resist_index = next_index(list(lhm_indices.values()))
lhm_indices[f"resist_{resist_index}"] = resist_index
bnd_point = Point(row.geometry.x, row.geometry.y + 100)
bnd_index = next_index(list(lhm_indices.values()))
lhm_indices[f"bnd_{bnd_index}"] = bnd_index
nodes += [
(resist_index, lhm_id, "LinearResistance", resist_point),
(bnd_index, lhm_id, "LevelBoundary", bnd_point)
]
edges += [
(row.Index, resist_index, LineString([row.geometry, resist_point])),
(resist_index, bnd_index, LineString([resist_point, bnd_point]))
]
level_boundary += [(bnd_index, -0.01)]
linear_resistance += [(resist_index, 5000)]
nodes_gdf = pd.concat(
[nodes_gdf,
gpd.GeoDataFrame(
nodes,
columns=["index","origin","type","geometry"],
crs=28992
).set_index("index")]
)
edges_gdf = pd.concat(
[edges_gdf,
gpd.GeoDataFrame(
edges,
columns=["from_node_id","to_node_id","geometry"],
crs=28992)]
)
edges_gdf["edge_type"] = "flow"
level_boundary = []
linear_resistance = []
nodes= []
edges = []
level_boundary_mask = nodes_gdf.index.isin(edges_gdf.to_node_id) & ~nodes_gdf.index.isin(edges_gdf.from_node_id)
for row in nodes_gdf.loc[level_boundary_mask].itertuples():
resist_point = Point(row.geometry.x, row.geometry.y + 50)
resist_index = next_index(list(lhm_indices.values()))
lhm_indices[f"resist_{resist_index}"] = resist_index
bnd_point = Point(row.geometry.x, row.geometry.y + 100)
bnd_index = next_index(list(lhm_indices.values()))
lhm_indices[f"bnd_{bnd_index}"] = bnd_index
nodes += [
(resist_index, lhm_id, "LinearResistance", resist_point),
(bnd_index, lhm_id, "LevelBoundary", bnd_point)
]
edges += [
(row.Index, resist_index, LineString([row.geometry, resist_point])),
(resist_index, bnd_index, LineString([resist_point, bnd_point]))
]
level_boundary += [(bnd_index, -0.01)]
linear_resistance += [(resist_index, 5000)]
nodes_gdf = pd.concat(
[nodes_gdf,
gpd.GeoDataFrame(
nodes,
columns=["index","origin","type","geometry"],
crs=28992
).set_index("index")]
)
edges_gdf = pd.concat(
[edges_gdf,
gpd.GeoDataFrame(
edges,
columns=["from_node_id","to_node_id","geometry"],
crs=28992)]
)
edges_gdf["edge_type"] = "flow"
4. Wegschrijven model¶
In [10]:
Copied!
ribasim_node = ribasim.Node(static=nodes_gdf[['lhm_index', 'geometry', 'type']])
ribasim_edge = ribasim.Edge(static=edges_gdf)
# 2 mm/d precipitation, 1 mm/d evaporation
seconds_in_day = 24 * 3600
precipitation = 0.002 / seconds_in_day
evaporation = 0.001 / seconds_in_day
static_df = pd.DataFrame(nodes_gdf[nodes_gdf["type"] == "Basin"].reset_index()["index"].values, columns=["node_id"])
static_df["drainage"] = 0.0
static_df["potential_evaporation"] = evaporation
static_df["infiltration"] = 0.0
static_df["precipitation"] = precipitation
static_df["urban_runoff"] = 0.0
static_df["target_level"] = None
ribasim_basin = ribasim.Basin(
profile=basin_profiles_gdf,
static=static_df
)
ribasim_rating_curve = ribasim.TabulatedRatingCurve(static=tabulated_profiles_df)
ribasim_fractional_flow = ribasim.FractionalFlow(
static=pd.DataFrame(fractional_flow, columns=["node_id", "fraction"])
)
ribasim_manning_resistance = ribasim.ManningResistance(
static=pd.DataFrame(manning_resistance,columns=["node_id","length","manning_n","profile_width","profile_slope"])
)
ribasim_flow_boundary = None
ribasim_level_boundary = ribasim.LevelBoundary(
static=pd.DataFrame(
level_boundary,
columns=["node_id", "level"]
)
)
ribasim_linear_resistance = ribasim.LinearResistance(
static=pd.DataFrame(
linear_resistance,
columns=["node_id", "resistance"]
)
)
model = ribasim.Model(
modelname=model_name,
node=ribasim_node,
edge=ribasim_edge,
basin=ribasim_basin,
level_boundary=ribasim_level_boundary,
flow_boundary=ribasim_flow_boundary,
manning_resistance=ribasim_manning_resistance,
linear_resistance=ribasim_linear_resistance,
tabulated_rating_curve=ribasim_rating_curve,
fractional_flow=ribasim_fractional_flow,
starttime="2020-01-01 00:00:00",
endtime="2021-01-01 00:00:00",
)
model.write(MODEL_DIR / model_name)
ribasim_node = ribasim.Node(static=nodes_gdf[['lhm_index', 'geometry', 'type']])
ribasim_edge = ribasim.Edge(static=edges_gdf)
# 2 mm/d precipitation, 1 mm/d evaporation
seconds_in_day = 24 * 3600
precipitation = 0.002 / seconds_in_day
evaporation = 0.001 / seconds_in_day
static_df = pd.DataFrame(nodes_gdf[nodes_gdf["type"] == "Basin"].reset_index()["index"].values, columns=["node_id"])
static_df["drainage"] = 0.0
static_df["potential_evaporation"] = evaporation
static_df["infiltration"] = 0.0
static_df["precipitation"] = precipitation
static_df["urban_runoff"] = 0.0
static_df["target_level"] = None
ribasim_basin = ribasim.Basin(
profile=basin_profiles_gdf,
static=static_df
)
ribasim_rating_curve = ribasim.TabulatedRatingCurve(static=tabulated_profiles_df)
ribasim_fractional_flow = ribasim.FractionalFlow(
static=pd.DataFrame(fractional_flow, columns=["node_id", "fraction"])
)
ribasim_manning_resistance = ribasim.ManningResistance(
static=pd.DataFrame(manning_resistance,columns=["node_id","length","manning_n","profile_width","profile_slope"])
)
ribasim_flow_boundary = None
ribasim_level_boundary = ribasim.LevelBoundary(
static=pd.DataFrame(
level_boundary,
columns=["node_id", "level"]
)
)
ribasim_linear_resistance = ribasim.LinearResistance(
static=pd.DataFrame(
linear_resistance,
columns=["node_id", "resistance"]
)
)
model = ribasim.Model(
modelname=model_name,
node=ribasim_node,
edge=ribasim_edge,
basin=ribasim_basin,
level_boundary=ribasim_level_boundary,
flow_boundary=ribasim_flow_boundary,
manning_resistance=ribasim_manning_resistance,
linear_resistance=ribasim_linear_resistance,
tabulated_rating_curve=ribasim_rating_curve,
fractional_flow=ribasim_fractional_flow,
starttime="2020-01-01 00:00:00",
endtime="2021-01-01 00:00:00",
)
model.write(MODEL_DIR / model_name)