Skip to content

Commit

Permalink
Merge branch 'main' into docs_rws_regio
Browse files Browse the repository at this point in the history
  • Loading branch information
visr authored Oct 1, 2024
2 parents f5a24c6 + 71697aa commit 3da2de2
Show file tree
Hide file tree
Showing 18 changed files with 891 additions and 241 deletions.
10 changes: 3 additions & 7 deletions .gitignore
Original file line number Diff line number Diff line change
Expand Up @@ -20,17 +20,13 @@ dmypy.json
# pytest --basetemp
src/hydamo/tests/temp/
src/ribasim_nl/tests/temp/
src/peilbeheerst_model/tests/temp/

src/peilbeheerst_model/*.html
src/peilbeheerst_model/*.code-workspace
src/peilbeheerst_model/.vscode

src/peilbeheerst_model/*.jpeg
src/peilbeheerst_model/*.gpkg
src/peilbeheerst_model/tests_results
/src/peilbeheerst_model/Output_zdrive
/src/peilbeheerst_model/Rekenend_Model_Test
/src/peilbeheerst_model/vervallen
/src/peilbeheerst_model/Parametrize/ribasim
/src/peilbeheerst_model/Parametrize/fix_FF_HHSK.ipynb
/src/peilbeheerst_model/01_test_parse_crossings.ipynb

notebooks/rijkswaterstaat/plots/
202 changes: 21 additions & 181 deletions notebooks/de_dommel/05_add_berging.py
Original file line number Diff line number Diff line change
@@ -1,15 +1,11 @@
# %%

import geopandas as gpd
import numpy as np
import numpy.typing as npt
import pandas as pd
import rasterio
from rasterio.windows import from_bounds
from rasterstats import zonal_stats
from ribasim import Node
from ribasim.nodes import basin, tabulated_rating_curve
from ribasim.nodes import basin
from ribasim_nl import CloudStorage, Model
from ribasim_nl.berging import add_basin_statistics, get_basin_profile, get_rating_curve
from ribasim_nl.geometry import basin_to_point
from shapely.geometry import LineString

Expand All @@ -18,190 +14,19 @@
ribasim_toml = cloud.joinpath("DeDommel", "modellen", "DeDommel_parameterized", "model.toml")
model = Model.read(ribasim_toml)

banden = {
"maaiveld": 1,
"bodemhoogte_primair_winter": 2,
"bodemhoogte_primair_zomer": 3,
"bodemhoogte_secundair_winter": 4,
"bodemhoogte_secundair_zomer": 5,
"bodemhoogte_tertiair_winter": 6,
"bodemhoogte_tertiair_zomer": 7,
"ghg_2010-2019": 8,
"glg_2010-2019": 9,
"opp_primair": 10,
"opp_secundair": 11,
"opp_tertiair": 12,
}

basin_area_df = gpd.read_file(
cloud.joinpath("DeDommel", "verwerkt", "basin_area.gpkg"), engine="pyogrio", fid_as_index=True
)
basin_area_df.set_index("node_id", inplace=True)

lhm_rasters = cloud.joinpath("Basisgegevens", "LHM", "4.3", "input", "LHM_data.tif")
ma_raster = cloud.joinpath("Basisgegevens", "VanDerGaast_QH", "spafvoer1.tif")


def sample_raster(
raster_file,
df,
band=1,
fill_value: float | None = None,
all_touched=False,
stats="mean",
maaiveld_data: npt.ArrayLike | None = None,
):
with rasterio.open(raster_file) as raster_src:
# read band
data = raster_src.read(band)

if maaiveld_data is not None:
data = maaiveld_data - data

# fill nodata
if fill_value is not None:
data = np.where(data == raster_src.nodata, fill_value, data)

affine = raster_src.transform

return zonal_stats(df, data, affine=affine, stats=stats, nodata=raster_src.nodata, all_touched=all_touched)


def get_rating_curve(row, min_level, maaiveld: None | float = None):
flow_rate = np.round([0, row.ma * 0.2, row.ma * 0.33, row.ma / 2, row.ma * 2], decimals=2)
depth = np.round([row.glg + 1, row.glg, row.ghg, row.ghg / 2, 0], decimals=2)

# set GxG < 0 to 0
depth[depth < 0] = 0

# level relative to maaiveld
if maaiveld is not None:
level = maaiveld - depth
else:
level = row.maaiveld - depth

# make sure level >= min_level
level[level < min_level] = min_level

# flow_rate in m3/s
flow_rate = flow_rate / 1000 * row.geometry.area / 86400

df = pd.DataFrame({"level": np.round(level, decimals=2), "flow_rate": np.round(flow_rate, decimals=5)})
df.drop_duplicates("level", keep="first", inplace=True)
df.drop_duplicates("flow_rate", keep="last", inplace=True)

return tabulated_rating_curve.Static(level=df.level, flow_rate=df.flow_rate)


def get_basin_profile(basin_polygon, polygon, max_level, min_level):
with rasterio.open(lhm_rasters) as src:
level = np.array([], dtype=float)
area = np.array([], dtype=float)

# Get the window and its transform
window = from_bounds(*basin_polygon.bounds, transform=src.transform)

if (window.width) < 1 or (window.height < 1):
window = from_bounds(*basin_polygon.centroid.buffer(125).bounds, transform=src.transform)
window_transform = src.window_transform(window)

# Primary water bottom-level
window_data = src.read(3, window=window)

# We don't want hoofdwater / doorgaand water to be in profile
if (polygon is None) | (window_data.size == 0):
mask = ~np.isnan(window_data)
else:
mask = rasterio.features.geometry_mask(
[polygon], window_data.shape, window_transform, all_touched=True, invert=True
)
# Include nodata as False in mask
mask[np.isnan(window_data)] = False

# add levels
level = np.concat([level, window_data[mask].ravel()])

# add areas on same mask
window_data = src.read(10, window=window)
area = np.concat([area, window_data[mask].ravel()])

# Secondary water
window_data = src.read(5, window=window)
mask = ~np.isnan(window_data)
level = np.concat([level, window_data[mask].ravel()])

window_data = src.read(11, window=window)
area = np.concat([area, window_data[mask].ravel()])

# Tertiary water water
window_data = src.read(7, window=window)
mask = ~np.isnan(window_data)
level = np.concat([level, window_data[mask].ravel()])

window_data = src.read(12, window=window)
area = np.concat([area, window_data[mask].ravel()])

# Make sure area is never larger than polygon-area
area[area > basin_polygon.area] = basin_polygon.area

# If area is empty, we add min_level at 5% of polygon-area
if area.size == 0:
level = np.append(level, min_level)
area = np.append(area, basin_polygon.area * 0.05)

# Add extra row with max_level at basin_polygon.area
level = np.append(level, max_level)
area = np.append(area, basin_polygon.area)

# In pandas for magic
df = pd.DataFrame({"level": np.round(level, decimals=2), "area": np.round(area)})
df.sort_values(by="level", inplace=True)
df = df.set_index("level").cumsum().reset_index()
df.dropna(inplace=True)
df.drop_duplicates("level", keep="last", inplace=True)
lhm_raster_file = cloud.joinpath("Basisgegevens", "LHM", "4.3", "input", "LHM_data.tif")
ma_raster_file = cloud.joinpath("Basisgegevens", "VanDerGaast_QH", "spafvoer1.tif")

# Return profile
return basin.Profile(area=df.area, level=df.level)

basin_area_df = add_basin_statistics(df=basin_area_df, lhm_raster_file=lhm_raster_file, ma_raster_file=ma_raster_file)

# %% add columns

with rasterio.open(lhm_rasters) as raster_src:
# read band
maaiveld_data = raster_src.read(banden["maaiveld"])

# sample rasters
ghg = sample_raster(
raster_file=lhm_rasters,
df=basin_area_df,
band=banden["ghg_2010-2019"],
all_touched=True,
maaiveld_data=maaiveld_data,
)
basin_area_df.loc[:, ["ghg"]] = pd.Series(dtype=float)
basin_area_df.loc[:, ["ghg"]] = [i["mean"] for i in ghg]

glg = sample_raster(
raster_file=lhm_rasters,
df=basin_area_df,
band=banden["glg_2010-2019"],
all_touched=True,
maaiveld_data=maaiveld_data,
)
basin_area_df.loc[:, ["glg"]] = pd.Series(dtype=float)
basin_area_df.loc[:, ["glg"]] = [i["mean"] for i in glg]

ma = sample_raster(raster_file=ma_raster, df=basin_area_df, all_touched=True, fill_value=37) # 37mm/dag is
basin_area_df.loc[:, ["ma"]] = pd.Series(dtype=float)
basin_area_df.loc[:, ["ma"]] = [i["mean"] for i in ma]

maaiveld = sample_raster(
raster_file=lhm_rasters, df=basin_area_df, band=banden["maaiveld"], stats="mean min max", all_touched=True
)
basin_area_df.loc[:, ["maaiveld"]] = pd.Series(dtype=float)
basin_area_df.loc[:, ["maaiveld"]] = [i["mean"] if pd.isna(i["mean"]) else i["mean"] for i in maaiveld]
basin_area_df.loc[:, ["maaiveld_max"]] = [i["max"] if pd.isna(i["max"]) else i["max"] for i in maaiveld]
basin_area_df.loc[:, ["maaiveld_min"]] = [i["min"] if pd.isna(i["min"]) else i["min"] for i in maaiveld]
# %%update model
edge_id = model.edge.df.index.max() + 1
for row in model.basin.node.df.itertuples():
Expand Down Expand Up @@ -233,7 +58,13 @@ def get_basin_profile(basin_polygon, polygon, max_level, min_level):
min_level = max_level = basin_area_df.at[node_id, "maaiveld_min"]
if min_level == max_level:
min_level -= 0.1
basin_profile = get_basin_profile(basin_polygon, polygon, max_level=max_level, min_level=min_level)
basin_profile = get_basin_profile(
basin_polygon=basin_polygon,
polygon=polygon,
max_level=max_level,
min_level=min_level,
lhm_raster_file=lhm_raster_file,
)
data = [
basin_profile,
basin.State(level=[basin_profile.df.level.min() + 0.1]),
Expand Down Expand Up @@ -267,6 +98,15 @@ def get_basin_profile(basin_polygon, polygon, max_level, min_level):
else:
print(f"Geen basin-vlak voor {node_id}")

# %%
df = pd.DataFrame({"node_id": model.basin.node.df.index.to_list()})
df.index.name = "fid"
df.loc[:, "precipitation"] = 5.787037e-08
df.loc[:, "potential_evaporation"] = 1.157407e-08
df.loc[:, "drainage"] = 0
df.loc[:, "infiltration"] = 0
model.basin.static.df = df

# %%
ribasim_toml = cloud.joinpath("DeDommel", "modellen", "DeDommel_bergend", "model.toml")

Expand Down
Loading

0 comments on commit 3da2de2

Please sign in to comment.