Rioxarray (#129)

* testing rioxarray

* more tests

* map_blocks working

* add hrt test

* test on different dem

* map_blocks speed testing

Generate blocks_gdf from a raster, specify which blocks we should do calculation, then skip the other blocks, returning nodata. Needs testing if it is effective.

* more blocktests

* tests

* speed tests home pc

Made some changes afterwards. rxr no blocks seems to be the prefered route.

* create overviews on raster

for faster rendering in gis. Will be added to options in hrt

* Update rxr_testing.py

* geocube gdf to raster

* hrt seems faster than geocube without interpolation.

* Create wlvl_writer.py

* deprecate hrt.create_meta_from_gdf

* formatting

* add DamageDem pytest

* DamageDem with rxr

* add oracledb

* landuse clip to rxr

* waterdelen rasterizen

* damagedem rxr reproject refactored

* GridToWaterDepth en GridToWaterLevel

* fix misc issues

* statistics to rxr

* Update grid_to_raster.py

* storage raster a la rastercalculator

* working storage raster with tests

* clean code

add OLD calculation. Dont merge this, used now to test speed differences, old method seems faster with the test dataset.

* set variables at init

* clean test

* correctly apply masks

* nodatamask works with nan and to func

* introduce RasterCalculatorRxr

* remove .pl

* refactor GridToWaterLevel test

* schema versie back to 217

* fix storage_raster test

* check wlvl raster on 1d2d test

Pytest not working yet. Result wlvl raster is quite different.

* Create to mkdir

* ModelbuilderRastersRxr opzet

* remove unused code

* add polder_tif to calculation rasters

* add waterdeel.tif to calculation_rasters

* fix tests waterdeel polder tif

* add ClipModelRasterCalc

Maakt clip van bronraster met de polder polygon. Optie om watervlakken ook een specifieke waarde te geven.

* clean schema rasters, working tests

* remove nodata_keys from RasterCalculatorRxr

* add landuse to clipmodelraster

* nodata_keys as input for more general behaviour

* improve readability

* remove create_landuse_polder_clip

Is now obsolete since it can be done with create_schematisation_rasters.ClipModelRasterCalc

* fix tests

* add todo wdepth klimaatsommen

* ruff to pyproject

* ruff command

* klimaatsommenprep uses new wlvl calc

Have to make a version that does the old one so we can compare.

* drooglegging sqlite to rxr

* add old wdepth creation as option to klimaatsommen

All tests working now.
Deprecate old wdepth calculation on later release

* add GridToRaster to init

* review changes

* No change to run tests

* fix tests

Something going on with PATH_NEW_FOLDER

---------

Co-authored-by: Daniel Tollenaar <[email protected]>

bin/run_ruff.bat

@@ -0,0 +1,2 @@
+REM force-exclude to make sure it uses the extend-exclude from pyproject.
+python -m ruff format ../hhnk_threedi_tools/**/*.py --force-exclude

envs/environment_test.yml

@@ -42,7 +42,7 @@ dependencies:
   - cached_property=1.5.2
   - rtree=1.0.1
   - pip=24.0
-  - oracledb==2.4.1
+  - oracledb=2.4.1
 
   - pip: # 3Di packages are not conda installable and are therefore installed by PIP
     - threedi_modelchecker==2.6.3 #threedi_results_analysis 3.9.0

hhnk_threedi_tools/__init__.py

@@ -16,10 +16,8 @@
 from hhnk_threedi_tools.core.checks.zero_d_one_d import ZeroDOneDTest
 from hhnk_threedi_tools.core.folders import Folders
 from hhnk_threedi_tools.core.folders_modelbuilder import FoldersModelbuilder
-from hhnk_threedi_tools.core.result_rasters.grid_to_raster import GridToRaster
-from hhnk_threedi_tools.core.modelbuilder.create_landuse_polder_clip import (
-    create_landuse_polder_clip,
-)
+from hhnk_threedi_tools.core.result_rasters.grid_to_raster import GridToWaterDepth, GridToWaterLevel
+from hhnk_threedi_tools.core.result_rasters.grid_to_raster_old import GridToRaster  # TODO deprecate
 from hhnk_threedi_tools.core.result_rasters.netcdf_to_gridgpkg import NetcdfToGPKG
 from hhnk_threedi_tools.core.schematisation import (
     migrate,

hhnk_threedi_tools/core/api/calculation.py

@@ -872,7 +872,6 @@ def download_sqlite(self, output_folder_sqlite=None):
         if self.model is None:
             return "define self.model_id first"
 
-
         output_path = Path(output_folder_sqlite).joinpath("tempfiles", f"model_{self.model_id}.zip")
         output_path.parent.mkdir(parents=True, exist_ok=True)  # Create parent folder
 
@@ -934,7 +933,7 @@ def simulation_info(self, str_type="text"):
             return f"Simulation: {sim.url}\
                     {newline}Scenario name: {sim.name}\
                     {newline}Organisation name: {sim.organisation_name}\
-                    {newline}Duration: {sim.duration}s ({sim.duration/3600}h)\
+                    {newline}Duration: {sim.duration}s ({sim.duration / 3600}h)\
                     {newline}Rain events: {self.data.rain}\
                     {newline}Control structures count: {len(self.data.structure_control)}\t(used={self.tracker.structure_control})\
                     {newline}Laterals count: {len(self.data.laterals)}\t(used={self.tracker.laterals})\
@@ -957,7 +956,7 @@ def simulation_info(self, str_type="text"):
                 f"Simulation id: <a href={sim.url}>{sim.id}</a>\
                     {newline}Scenario name: {sim.name}\
                     {newline}Organisation name: {sim.organisation_name}\
-                    {newline}Duration: {sim.duration}s ({sim.duration/3600}h)\
+                    {newline}Duration: {sim.duration}s ({sim.duration / 3600}h)\
                     {newline}Rain events: {self.data.rain}\
                     {newline}Control structures count: {len(self.data.structure_control)}\t(used={self.tracker.structure_control})\
                     {newline}Laterals count: {len(self.data.laterals)}\t(used={self.tracker.laterals})\

hhnk_threedi_tools/core/api/download_functions.py

@@ -1,4 +1,5 @@
 """Functions that interact with the results API"""
+
 import math
 import os
 import zipfile

hhnk_threedi_tools/core/api/download_gui_class.py

@@ -797,8 +797,8 @@ def download_batch(action):
             # batch_fd = Folders(batch_folder).
             # Create destination folder
 
-            self.vars.batch_fd.create()
-            self.vars.batch_fd.downloads.create()
+            self.vars.batch_fd.mkdir()
+            self.vars.batch_fd.downloads.mkdir()
 
             # Temporary disable download button
             self.download_batch.button.style.button_color = "orange"
@@ -879,7 +879,7 @@ def download_batch(action):
                     output_folder = getattr(self.vars.batch_fd.downloads, row["dl_name"]).netcdf
 
                     # Create destination folder
-                    output_folder.create()
+                    output_folder.mkdir()
                     output_folder = output_folder.path
 
                     # Start downloading of the files

hhnk_threedi_tools/core/checks/model_state.py

@@ -8,6 +8,7 @@
 used in folder structures
 
 """
+
 # system imports
 import os
 

hhnk_threedi_tools/core/checks/one_d_two_d.py

@@ -8,7 +8,7 @@
 
 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.core.result_rasters.grid_to_raster import GridToRaster
+from hhnk_threedi_tools.core.result_rasters.grid_to_raster import GridToWaterDepth, GridToWaterLevel
 
 # Local imports
 from hhnk_threedi_tools.core.result_rasters.netcdf_to_gridgpkg import NetcdfToGPKG
@@ -77,16 +77,22 @@ def run_wlvl_depth_at_timesteps(self, overwrite=False):
         # 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 GridToRaster(
+            with GridToWaterLevel(
                 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
+                level_raster = raster_calc.run(
+                    output_file=getattr(self.output_fd, f"waterstand_T{T}"),
+                    overwrite=overwrite,
                 )
-                raster_calc.run(
-                    output_file=getattr(self.output_fd, f"waterstand_T{T}"), mode="MODE_WLVL", overwrite=overwrite
+            with GridToWaterDepth(
+                dem_path=self.folder.model.schema_base.rasters.dem,
+                wlvl_path=level_raster,
+            ) as raster_calc:
+                _ = raster_calc.run(
+                    output_file=getattr(self.output_fd, f"waterdiepte_T{T}"),
+                    overwrite=overwrite,
                 )
 
     def run_flowline_stats(self):

hhnk_threedi_tools/core/checks/sqlite/sqlite_main.py

@@ -16,6 +16,7 @@
 import hhnk_research_tools as hrt
 import numpy as np
 import pandas as pd
+import xarray as xr
 from shapely import wkt
 from shapely.geometry import Point
 from threedigrid_builder import make_gridadmin
@@ -132,6 +133,7 @@
 ]
 
 
+# %%
 class SqliteCheck:
     def __init__(
         self,
@@ -194,72 +196,51 @@ def run_dem_max_value(self):
 
     def run_dewatering_depth(self, overwrite=False):
         """
-        Compares initial water level from fixed drainage level areas with
+        Compare initial water level from fixed drainage level areas with
         surface level in DEM of model. Initial water level should mostly
         be below surface level.
         """
+        drooglegging_raster = self.output_fd.drooglegging
+        create = hrt.check_create_new_file(output_file=drooglegging_raster.path, overwrite=overwrite)
+
+        if create:
+            # Rasterize fixeddrainage
+            wlvl_raster = self.output_fd.streefpeil
+            dtype = "float32"
+            nodata = hrt.variables.DEFAULT_NODATA_VALUES[dtype]
+
+            fixeddrainage_gdf = self.layer_fixeddrainage.load()
+            hrt.gdf_to_raster(
+                gdf=fixeddrainage_gdf,
+                value_field=COL_STREEFPEIL_BWN,
+                raster_out=wlvl_raster,
+                nodata=nodata,
+                metadata=self.dem.metadata,
+                read_array=False,
+                overwrite=overwrite,
+            )
 
-        def _create_drooglegging_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_dem, block_wlvl)
+            # Load data arrays
+            da_dem = self.dem.open_rxr()
+            da_wlvl = wlvl_raster.open_rxr()
+            crs = da_dem.rio.crs  # CRS is lost during xr.where; preserve it.
 
-            # Mask output
-            nodatamask = (block_dem == self.dem.nodata) | (block_dem == 10)
-            block_depth[nodatamask] = self.raster_out.nodata
+            da_out = da_dem - da_wlvl
 
-            # Get the window of the small raster
-            window_small = windows[[k for k, v in self.raster_mapping.items() if v == "small"][0]]
+            da_out = xr.where(da_dem == 10, nodata, da_out)
 
             # Write to file
-            band_out.WriteArray(block_depth, xoff=window_small[0], yoff=window_small[1])
-
-        try:
-            # Load layers
-            drooglegging_raster = self.output_fd.drooglegging
-
-            create = hrt.check_create_new_file(output_file=drooglegging_raster.path, overwrite=overwrite)
-
-            if create:
-                fixeddrainage_gdf = self.layer_fixeddrainage.load()
-                wlvl_raster = self.output_fd.streefpeil  # Rasterize fixeddrainage
-                hrt.gdf_to_raster(
-                    gdf=fixeddrainage_gdf,
-                    value_field=COL_STREEFPEIL_BWN,
-                    raster_out=wlvl_raster,
-                    nodata=self.dem.nodata,
-                    metadata=self.dem.metadata,
-                    read_array=False,
-                    overwrite=overwrite,
-                )
-
-                # Calculate drooglegging raster
-                drooglegging_calculator = hrt.RasterCalculator(
-                    raster1=self.dem,
-                    raster2=wlvl_raster,
-                    raster_out=drooglegging_raster,
-                    custom_run_window_function=_create_drooglegging_raster,
-                    output_nodata=self.dem.nodata,
-                    verbose=False,
-                )
-
-                drooglegging_calculator.run(overwrite=overwrite)
-
-                # remove temp files
-                wlvl_raster.unlink(missing_ok=True)
-                # self.results["dewatering_depth"] = output_file
-        except Exception as e:
-            raise e from None
+            da_out.rio.set_crs(crs)  # Reapply crs
+            drooglegging_raster = hrt.Raster.write(
+                drooglegging_raster,
+                result=da_out,
+                nodata=nodata,
+                dtype=dtype,
+                chunksize=None,
+            )
 
     def run_model_checks(self):
-        """
-        Collects all queries that are part of general model checks (see general_checks_queries file)
+        """Collect all queries that are part of general model checks (see general_checks_queries file)
         and executes them
         """
 
@@ -288,9 +269,10 @@ def run_geometry_checks(self):
         return result_db
 
     def run_imp_surface_area(self):
-        """
-        Calculates the impervious surface area (in the model), the area of the polder (based on the polder shapefile) and
-        the difference between the two.
+        """Calculate
+        1. the impervious surface area in the model
+        2. the area of the polder (based on the polder shapefile)
+        3. the difference between the two.
         """
         imp_surface_db = self.model.execute_sql_selection(impervious_surface_query)
         imp_surface_db.set_index("id", inplace=True)
@@ -352,9 +334,9 @@ def run_used_profiles(self):
         return channels_gdf
 
     def run_struct_channel_bed_level(self):
-        """
-        Checks whether the reference level of any of the adjacent cross section locations (channels) to a structure
-        is lower than the reference level for that structure (3di crashes if it is)
+        """Check whether the reference level of any of the adjacent cross section locations
+        (channels) to a structure is lower than the reference level for that structure
+        (3di crashes if it is)
         """
         datachecker_culvert_layer = self.fenv.source_data.datachecker.layers.culvert
         damo_duiker_sifon_layer = self.fenv.source_data.damo.layers.DuikerSifonHevel
@@ -453,9 +435,7 @@ def create_grid_from_sqlite(self, output_folder):
             hrt.gdf_write_to_geopackage(gdf, filepath=Path(output_folder) / f"{grid_type}.gpkg")
 
     def run_cross_section_duplicates(self, database):
-        """
-        Check for duplicate geometries in cross_section_locations.
-        """
+        """Check for duplicate geometries in cross_section_locations."""
         cross_section_point = database.execute_sql_selection(query=cross_section_location_query)
 
         # Make buffer of the points to identify if we have cross setion overlapping each other.
@@ -472,9 +452,7 @@ def run_cross_section_duplicates(self, database):
         return intersected_points
 
     def run_cross_section_no_vertex(self, database):
-        """
-        Check for cross_sections that are not located on the  vertex of a channel.
-        """
+        """Check for cross_sections that are not located on the  vertex of a channel."""
         # load cross locs and channels from sqlite
         cross_section_point = database.execute_sql_selection(query=cross_section_location_query)
         cross_section_point.rename({"geometry_point": "geometry"}, axis=1, inplace=True)
@@ -581,16 +559,14 @@ def calc_width_at_waterlevel(row):
 
 
 def split_round(item):
-    """
-    Split items in width and height columns by space, round all items in resulting list and converts to floats
+    """Split items in width and height columns by space,
+    round all items in resulting list and converts to floats
     """
     return [round(float(n), 2) for n in str(item).split(" ")]
 
 
 def get_max_depth(row):
-    """
-    calculates difference between initial waterlevel and reference level
-    """
+    """Calculate difference between initial waterlevel and reference level"""
     return round(float(row[initial_waterlevel_col]) - float(row[reference_level_col]), 2)
 
 
@@ -637,8 +613,7 @@ def add_datacheck_info(layer, gdf):
 
 
 def expand_multipolygon(df):
-    """
-    New version using explode, old version returned pandas dataframe not geopandas
+    """New version using explode, old version returned pandas dataframe not geopandas
     geodataframe (missing last line), I think it works now?
     """
     try: