Fix subgrid tif

docs/changelog.rst

@@ -6,8 +6,25 @@ Distinction is made between new methods (Added), changes to existing methods (Ch
 The format is based on `Keep a Changelog`_, and this project adheres to
 `Semantic Versioning`_.
 
-v1.0.1 (3 August 2023)
+v1.0.2 (unreleased)
 ===================
+
+Added
+-----
+
+Changed
+-------
+
+Fixed
+-----
+- writing COG files in `SfincsModel.setup_subgrid` (the COG driver settings were wrong) PR #117
+
+Deprecated
+----------
+
+
+v1.0.1 (3 August 2023)
+======================
 This release contains several new features, such as burning in river bathymetry into the subgrid, setting up drainage structures and adding wind and pressure forcing.
 It also contains several bugfixes and improvements to existing methods. 
 It is recommended to use this release together with the latest version of the `SFINCS model <https://github.com/Deltares/SFINCS/releases/tag/v2.0.2_Blockhaus_release>`_.
@@ -22,7 +39,7 @@ Added
 
 Changed
 -------
-- `SfincsModel.setup_subgrid`` now supports the 'riv_datasets' to burn in river bathymetry into the subgrid. PR #84
+- `SfincsModel.setup_subgrid` now supports the 'riv_datasets' to burn in river bathymetry into the subgrid. PR #84
 - `SfincsModel.setup_mask_active` argument reset_mask default to True PR #94
 - `SfincsModel.read_config` allows to use a template input file from a directory different than the model root. PR #102
 - Added the option to use landuse/landcover data combined with a reclass table to `SfincsModel.setup_constant_infiltration`.  PR #103

hydromt_sfincs/__init__.py

@@ -3,7 +3,7 @@
 from os.path import dirname, join, abspath
 
 
-__version__ = "1.0.1"
+__version__ = "1.0.2.dev0"
 
 DATADIR = join(dirname(abspath(__file__)), "data")
 

hydromt_sfincs/subgrid.py

@@ -286,18 +286,21 @@ def build(
             1 / nr_subgrid_pixels
         )
 
+        # create COGs for topobathy/manning
         profile = dict(
-            driver="COG",
+            driver="GTiff",
             width=output_width,
             height=output_height,
             count=1,
             dtype=np.float32,
             crs=da_mask.raster.crs,
             transform=output_transform,
+            tiled=True,
             blockxsize=256,
             blockysize=256,
             compress="deflate",
             predictor=2,
+            profile="COG",
             nodata=np.nan,
             BIGTIFF="YES",  # Add the BIGTIFF option here
         )
@@ -389,7 +392,6 @@ def build(
                 if nactive == 0:  # not active cells in block
                     logger.debug("Skip block - No active cells")
                     continue
-                logger.debug(f"Processing block with {nactive} active cells..")
                 transform = da_mask_block.raster.transform
                 # add refi cells overlap in both dimensions for u and v in last row/col
                 reproj_kwargs = dict(
@@ -417,7 +419,9 @@ def build(
                 # raise warning if NaN values in active cells
                 if np.any(np.isnan(da_dep.values[da_mask_sbg > 0])) > 0:
                     npx = int(np.sum(np.isnan(da_dep.values[da_mask_sbg > 0])))
-                    logger.warning(f"Interpolate data at {npx} subgrid pixels")
+                    logger.warning(
+                        f"Interpolate elevation data at {npx} subgrid pixels"
+                    )
                 # always interpolate/extrapolate to avoid NaN values
                 da_dep = da_dep.raster.interpolate_na(
                     method="rio_idw", extrapolate=True
@@ -433,9 +437,9 @@ def build(
                     )
                     # raise warning if NaN values in active cells
                     if np.isnan(da_man.values[da_mask_sbg > 0]).any():
-                        logger.debug(
-                            "Missing values in manning roughness array, "
-                            "fill with default values"
+                        npx = int(np.sum(np.isnan(da_man.values[da_mask_sbg > 0])))
+                        logger.warning(
+                            f"Fill manning roughness data at {npx} subgrid pixels with default values"
                         )
                     # always fill based on land/sea elevation to avoid NaN values
                     da_man0 = xr.where(
@@ -489,6 +493,7 @@ def build(
                     yg = np.repeat(np.atleast_2d(yg), da_dep.raster.shape[0], axis=0)
 
                 # Now compute subgrid properties
+                logger.debug(f"Processing subgrid tables for {nactive} active cells..")
                 sn, sm = slice(bn0, bn1), slice(bm0, bm1)
                 (
                     self.z_zmin[sn, sm],