Skip to content

Commit

Permalink
Squashed commit of the following:
Browse files Browse the repository at this point in the history 
commit e17fa6a
Author: Wietse van Gerwen <[email protected]>
Date:   Tue Mar 19 15:38:10 2024 +0100

    remove unnecessary import

commit 738e340
Merge: 96d4b9f 01af596
Author: Wietse van Gerwen <[email protected]>
Date:   Tue Mar 19 15:29:08 2024 +0100

    Merge branch 'main' into one_d_two_d-refactor

commit 96d4b9f
Author: Wietse van Gerwen <[email protected]>
Date:   Thu Feb 15 14:11:45 2024 +0100

    Update calculate_raster.py

commit 3b8b1b0
Author: Daniel <[email protected]>
Date:   Wed Feb 14 12:30:51 2024 +0100

    Delete one_d_two_d_exception.py

commit c025f5b
Author: Daniel <[email protected]>
Date:   Wed Feb 14 12:11:05 2024 +0100

    updated environment

commit ae38640
Author: Daniel <[email protected]>
Date:   Fri Feb 9 14:52:56 2024 +0100

    Create one_d_two_d_exception.py

commit 1b76dda
Author: wietsevangerwen <[email protected]>
Date:   Thu Jan 18 18:51:25 2024 +0100

    tests dont fail locally.

commit 166ba31
Author: wvangerwen <[email protected]>
Date:   Wed Dec 13 11:52:24 2023 +0100

    fix tests

commit 665a593
Author: wvangerwen <[email protected]>
Date:   Tue Dec 12 10:58:33 2023 +0100

    remove old code

commit 282346d
Author: wvangerwen <[email protected]>
Date:   Mon Dec 11 13:01:04 2023 +0100

    update grid and raster creation

commit f4955e0
Author: wvangerwen <[email protected]>
Date:   Mon Dec 11 12:58:33 2023 +0100

    change column order grid gpkg

commit c0783b2
Author: wvangerwen <[email protected]>
Date:   Thu Dec 7 12:51:53 2023 +0100

    remove unused code
  • Loading branch information
@wvangerwen
wvangerwen committed Mar 20, 2024
1 parent 3b93943 commit 8795f0b
Show file tree
Hide file tree
Showing 10 changed files with 577 additions and 313 deletions.
399 changes: 399 additions & 0 deletions deprecated/core/checks/one_d_two_d.py
View file

Large diffs are not rendered by default.

5 changes: 4 additions & 1 deletion envs/environment_test.yml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,7 +12,10 @@ dependencies:
- geopandas=0.12.2
- pandas=1.1.3
- scipy=1.6.2
- h5py=2.10.0
- h5py=3.8.0
- fiona=1.9.1
- shapely=2.0.1
- gdal=3.6.2

#User folder
- jupyterlab=3.6.3
Expand Down
8 changes: 2 additions & 6 deletions hhnk_threedi_tools/core/checks/grid_result_metadata.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,9 +5,7 @@


def calculate_rain_days(rain):
"""
Calculates days dry before and after rain
"""
"""Calculate days dry before and after rain"""
detected_rain = [i for i, e in enumerate(rain) if e > 1e-5]
# Collect indexes of items in rain where rain falls (every index represents an hour)
if detected_rain:
Expand All @@ -23,9 +21,7 @@ def calculate_rain_days(rain):


def get_rain_properties(results):
"""
Calculates the rain scenario used for this result
"""
"""Calculate the rain scenario used for this result"""
try:
# Calculates the mean of steps between timestamps (in seconds), then converts to minutes
dt = round(np.mean(np.diff(results.nodes.timestamps)) / 60, 0)
Expand Down
249 changes: 41 additions & 208 deletions hhnk_threedi_tools/core/checks/one_d_two_d.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,50 +3,40 @@
import hhnk_research_tools as hrt
import numpy as np
import pandas as pd
from hhnk_research_tools.variables import all_2d, t_end_rain_col, t_end_sum_col, t_start_rain_col
from shapely.geometry import LineString, box
from hhnk_research_tools.variables import t_end_rain_col, t_end_sum_col, t_start_rain_col
from shapely.geometry import LineString

# Local imports
import hhnk_threedi_tools.core.checks.grid_result_metadata as grid_result_metadata
from hhnk_threedi_tools.core.folders import Folders
from hhnk_threedi_tools.variables.database_aliases import df_geo_col
from hhnk_threedi_tools.core.result_rasters.calculate_raster import BaseCalculatorGPKG

# Local imports
from hhnk_threedi_tools.core.result_rasters.netcdf_to_gridgpkg import NetcdfToGPKG
from hhnk_threedi_tools.variables.default_variables import DEF_TRGT_CRS
from hhnk_threedi_tools.variables.one_d_two_d import (
content_type_col,
end_rain_sfx,
id_col,
kcu_col,
max_area_col,
max_sfx,
minimal_dem_col,
node_type_col,
one_d,
one_d_boundary_col,
one_d_two_d,
pump_capacity_m3_s_col,
q_m3_s_col,
spatialite_id_col,
start_rain_sfx,
storage_mm_col,
suffixes_list,
twelve_hr_after_rain_sfx,
two_d,
vel_m_s_col,
volume_m3_col,
wet_area_m2_col,
wtrlvl_col,
wtrlvl_m_col,
)


# TODO functies weer in class onderbrengen, class nu buiten gebruik.
class OneDTwoDTest:
TIMESTEPS = [1, 3, 15] # hours, 1=start rain, 3=end rain, 15=end calculation

def __init__(self, folder: Folders, revision=0, dem_path=None):
self.fenv = folder
self.folder = folder
self.revision = revision

self.result_fd = self.fenv.threedi_results.one_d_two_d[self.revision]
self.output_fd = self.fenv.output.one_d_two_d[self.revision]
self.result_fd = self.folder.threedi_results.one_d_two_d[self.revision]
self.output_fd = self.folder.output.one_d_two_d[self.revision]

self.grid_result = self.result_fd.grid

Expand All @@ -59,193 +49,46 @@ def __init__(self, folder: Folders, revision=0, dem_path=None):
self.timestep_df,
) = grid_result_metadata.construct_scenario(self.grid_result)

# if output_path:
# self.output_path = output_path
# self.layer_path = output_path + "/Layers"
# self.log_path = output_path + "/Logs"
# else:
# self.layer_path = str(folder.output.one_d_two_d[self.revision].layers)
# self.log_path = str(folder.output.one_d_two_d[self.revision].logs)

if dem_path:
self.dem = hrt.Raster(dem_path)
else:
self.dem = self.fenv.model.schema_base.rasters.dem

self.iresults = {}
self.dem = self.folder.model.schema_base.rasters.dem

@classmethod
def from_path(cls, path_to_polder, **kwargs):
return cls(Folders(path_to_polder), **kwargs)

@property
def results(self):
return self.iresults

def run_wlvl_depth_at_timesteps(self, overwrite=False):
"""Transform netcdf to grid gpkg and apply wlvl correction
Then create waterlevel and depth rasters at 3 timesteps:
1h : start rain
3h : end rain
15h : end calculation
"""
Deze functie bepaalt de waterstanden op de gegeven tijdstappen op basis van het 3di resultaat.
Vervolgens wordt op basis van de DEM en de waterstand per tijdstap de waterdiepte bepaald.
"""

def _create_depth_raster(self, windows, band_out, **kwargs):
"""hrt.Raster_calculator custom_run_window_function"""
self.dem = self.raster1
self.wlvl = self.raster2

block_dem = self.dem._read_array(window=windows["raster1"])
block_wlvl = self.wlvl._read_array(window=windows["raster2"])

# Calculate output
block_depth = np.subtract(block_wlvl, block_dem)

# Mask output
nodatamask = (block_dem == self.dem.nodata) | (block_dem == 10) | (block_depth < 0)
block_depth[nodatamask] = self.raster_out.nodata

# Get the window of the small raster
window_small = windows[[k for k, v in self.raster_mapping.items() if v == "small"][0]]

# Write to file
band_out.WriteArray(block_depth, xoff=window_small[0], yoff=window_small[1])

try:
timesteps_arr = [
self.timestep_df["t_start_rain"].value,
self.timestep_df["t_end_rain"].value,
self.timestep_df["t_end_sum"].value,
]
# hours since start of calculation
timestrings = [int(round(self.grid_result.nodes.timestamps[t] / 60 / 60, 0)) for t in timesteps_arr]

assert timestrings == [1, 3, 15]

# For each timestring calculate wlvl and depth raster.
for timestep, timestr in zip(timesteps_arr, timestrings):
wlvl_raster = getattr(self.output_fd, f"waterstand_T{timestr}")
depth_raster = getattr(self.output_fd, f"waterdiepte_T{timestr}")

# Calculate wlvl raster
nodes_2d_wlvl = self._read_2node_wlvl_at_timestep(timestep)
hrt.gdf_to_raster(
gdf=nodes_2d_wlvl,
value_field=wtrlvl_col,
raster_out=wlvl_raster,
nodata=self.dem.nodata,
metadata=self.dem.metadata,
read_array=False,
overwrite=overwrite,
netcdf_gpkg = NetcdfToGPKG.from_folder(folder=self.folder, threedi_result=self.result_fd)

# Convert netcdf to grid gpkg
netcdf_gpkg.run(
output_file=self.output_fd.grid_nodes_2d,
timesteps_seconds=[T * 3600 for T in self.TIMESTEPS],
overwrite=True,
)

# Create depth and wlvl rasters for each timestep.
grid_gdf = gpd.read_file(self.output_fd.grid_nodes_2d.path)
for T in self.TIMESTEPS:
with BaseCalculatorGPKG(
dem_path=self.folder.model.schema_base.rasters.dem,
grid_gdf=grid_gdf,
wlvl_column=f"wlvl_{T}h",
) as raster_calc:
raster_calc.run(
output_file=getattr(self.output_fd, f"waterdiepte_T{T}"), mode="MODE_WDEPTH", overwrite=overwrite
)

# Calculate depth raster
depth_calculator = hrt.RasterCalculator(
raster1=self.dem,
raster2=wlvl_raster,
raster_out=depth_raster,
custom_run_window_function=_create_depth_raster,
output_nodata=0,
verbose=False,
)

depth_calculator.run(overwrite=overwrite)
except Exception as e:
raise e from None

def _read_2node_wlvl_at_timestep(self, timestep):
"""timesteps is the index of the time in the timeseries you want to use
to calculate the wlvl and depth raster"""
nodes_2d = gpd.GeoDataFrame()
# * inputs every element from row as a new function argument.
nodes_2d[df_geo_col] = [box(*row) for row in self.grid_result.cells.subset(all_2d).cell_coords.T]
# waterstand
nodes_2d[wtrlvl_col] = self.grid_result.nodes.subset(all_2d).timeseries(indexes=[timestep]).s1[0]
return nodes_2d

def run_node_stats(self):
"""
Deze functie leest alle 2d nodes uit het 3di resultaat en berekent de volgende waarden:
* de minimale DEM waarde binnen het gebied van de betreffende node (geometrie is omgezet naar een vierkant)
* het totale oppervlak dat de node beslaat
Vervolgens wordt op drie tijdstappen (het begin van de regen het einde van de regen en het einde van de som
de volgende informatie berekend:
* de waterstand op de genoemde tijdstappen
* de hoeveelheid water (volume in m3) per tijdstap
* het natte oppervlak per tijdstap (in m2)
* opslag van regen in het gebied van de node (hoeveelheid water / totale oppervlak gebied)
"""
try:
nodes_wlvl = hrt.threedi.grid_nodes_to_gdf(self.grid_result)
# We need to keep a reference to the crs to restore it later
crs_orig = nodes_wlvl.crs
nodes_wlvl[id_col] = self.grid_result.nodes.id
nodes_wlvl[spatialite_id_col] = self.grid_result.nodes.content_pk
nodes_wlvl[node_type_col] = self.grid_result.nodes.node_type
# Replace numbers with human readable values
nodes_wlvl[node_type_col].replace([1, 3, 7], [two_d, one_d, one_d_boundary_col], inplace=True)

# totaal oppervlak
nodes_wlvl[max_area_col] = self.grid_result.nodes.sumax

# Load grid_result
# waterstand
wlvl = self.grid_result.nodes.timeseries(
indexes=[
self.timestep_df[t_start_rain_col].value,
self.timestep_df[t_end_rain_col].value,
self.timestep_df[t_end_sum_col].value,
]
).s1
volume = self.grid_result.nodes.timeseries(
indexes=[
self.timestep_df[t_start_rain_col].value,
self.timestep_df[t_end_rain_col].value,
self.timestep_df[t_end_sum_col].value,
]
).vol
# actueel nat oppervlak
wet_area = self.grid_result.nodes.timeseries(
indexes=[
self.timestep_df[t_start_rain_col].value,
self.timestep_df[t_end_rain_col].value,
self.timestep_df[t_end_sum_col].value,
]
).su

# Add grid_result to dataframe
args_lst = [start_rain_sfx, end_rain_sfx, twelve_hr_after_rain_sfx]
for index, time_str in enumerate(args_lst):
nodes_wlvl[wtrlvl_m_col + time_str] = np.round(wlvl[index], 2)
for index, time_str in enumerate(args_lst):
nodes_wlvl[wet_area_m2_col + time_str] = np.round(wet_area[index], 2)
for index, time_str in enumerate(args_lst):
nodes_wlvl[volume_m3_col + time_str] = np.round(volume[index], 2)
for index, time_str in enumerate(args_lst):
nodes_wlvl[storage_mm_col + time_str] = np.round(
nodes_wlvl[volume_m3_col + time_str] / nodes_wlvl[max_area_col], 2
raster_calc.run(
output_file=getattr(self.output_fd, f"waterstand_T{T}"), mode="MODE_WLVL", overwrite=overwrite
)

# select 2d nodes and create polygons for plotting.
nodes_2d = self._2d_nodes_to_grid(nodes_wlvl=nodes_wlvl)
orig_geom = nodes_2d[df_geo_col]
nodes_2d_gdf = gpd.GeoDataFrame(nodes_2d, geometry=orig_geom, crs=crs_orig)

# self.iresults["node_stats"] = nodes_2d_gdf
return nodes_2d_gdf
except Exception as e:
raise e from None

# TODO staat dit niet al in hhnk_research_tools
def _2d_nodes_to_grid(self, nodes_wlvl):
"""Transfer the nodes into polygons of the grid."""
try:
nodes_2d = nodes_wlvl[nodes_wlvl[node_type_col] == two_d].copy()
# replace geometry with polygons of the cells
nodes_2d.loc[:, df_geo_col] = [box(*row) for row in self.grid_result.cells.subset(all_2d).cell_coords.T]
nodes_2d.loc[:, minimal_dem_col] = self.grid_result.cells.subset(all_2d).z_coordinate
return nodes_2d
except Exception as e:
raise e from None

def run_flowline_stats(self):
"""
Deze functie leest alle stroom lijnen in uit het 3di resultaat. Vervolgens wordt gekeken naar het type van de lijn
Expand Down Expand Up @@ -379,20 +222,10 @@ def _read_pumpline_results(self):
if __name__ == "__main__":
from pathlib import Path

from hhnk_threedi_tools import Folders

TEST_MODEL = Path(__file__).parent.parent.parent.parent.full_path(r"tests/data/model_test/")
TEST_MODEL = Path(__file__).parents[3].joinpath(r"tests/data/model_test/")
folder = Folders(TEST_MODEL)
# %%
self = OneDTwoDTest.from_path(TEST_MODEL)

# def test_run_depth_at_timesteps_test(self):
"""test of de 0d1d test werkt"""
output = self.run_levels_depths_at_timesteps()

assert len(output) > 0
assert output[0] == 1
assert "waterdiepte_T15.tif" in self.test_1d2d.fenv.output.one_d_two_d[0].content
overwrite = True

# %%
folder.threedi_results.one_d_two_d.revisions
output = self.run_wlvl_depth_at_timesteps()
8 changes: 4 additions & 4 deletions hhnk_threedi_tools/core/climate_scenarios/klimaatsommen_prep.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -89,7 +89,7 @@ def netcdf_to_grid(
overwrite=False,
):
"""Transform netcdf to grid gpkg and apply wlvl correction
output will be stored in wlvl_max_corr column
output will be stored in wlvl_corr_max column
"""
# Select result
netcdf_gpkg = NetcdfToGPKG.from_folder(
Expand All @@ -110,7 +110,7 @@ def calculate_raster(
threedi_result: hrt.ThreediResult,
mode: str,
grid_filename: str = "grid_wlvl.gpkg",
wlvl_col_name: str = "wlvl_max_corr",
wlvl_col_name: str = "wlvl_corr_max",
overwrite=False,
):
"""Mode options are: 'MODE_WDEPTH', 'MODE_WLVL'"""
Expand All @@ -133,7 +133,7 @@ def calculate_depth(
scenario,
threedi_result: hrt.ThreediResult,
grid_filename: str,
wlvl_col_name="wlvl_max_corr",
wlvl_col_name="wlvl_corr_max",
overwrite=False,
):
scenario_raster = scenario.depth_max
Expand All @@ -151,7 +151,7 @@ def calculate_wlvl(
scenario,
threedi_result: hrt.ThreediResult,
grid_filename: str,
wlvl_col_name="wlvl_max_corr",
wlvl_col_name="wlvl_corr_max",
overwrite=False,
):
scenario_raster = scenario.wlvl_max
Expand Down
Loading

0 comments on commit 8795f0b

Please sign in to comment.