Feature/offset scale outputs

.gitignore

@@ -1,6 +1,7 @@
 *.gem
 *.sublime-project
 *.sublime-workspace
+.kiro
 .bundle
 .DS_Store
 .jekyll-metadata
@@ -27,3 +28,6 @@ lisflood_model.egg-info
 .vscode/
 *.ipynb
 .ipynb_checkpoints/
+implementation_guide_scale_offset_packing.md
+lisflood_optimization_report.md
+lisflood_optimization_report_v2.md

docs/5_annex_output-files/index.md

@@ -157,7 +157,7 @@ To speed up the pre-run and to prevent that results are taken from the pre-run,
 | actual transpiration                    | repTaMaps              | $\frac{mm}{timestep}$      | TaMaps                                             | tact            |
 | rainfall interception                   | repInterceptionMaps    | $\frac{mm}{timestep}$      | InterceptionMaps <br> InterceptionForestMaps                                   | int <br> intF            |
 | evaporation of intercepted water        | repEWIntMaps           | $\frac{mm}{timestep}$      | EWIntMaps                                          | ewint           |
-| leaf drainage                           | repLeafDrainageMaps    | $\frac{mm}{timestep}$      | LeafDrainageMaps <br> LeafDrainageForestMaps                                   | ldra <br> draF            |
+| leaf drainage                           | repLeafDrainageMaps    | $\frac{mm}{timestep}$      | LeafDrainageMaps <br> LeafDrainageForestMaps                     no              | ldra <br> draF            |
 | infiltration                            | repInfiltrationMaps    | $\frac{mm}{timestep}$      | InfiltrationMaps <br> InfiltrationForestMaps                                   | inf <br> infF             |
 | preferential (bypass) flow              | repPrefFlowMaps        | $\frac{mm}{timestep}$      | PrefFlowMaps <br> PrefFlowtherMaps <br> PrefFlowForestMaps <br> PrefFlowIrrigationMaps                                       | pflowpixel <br> pflow <br> pflowF <br> pflowi           |
 | percolation upper to lower soil layer   | repPercolationMaps     | $\frac{mm}{timestep}$      | Percolation1ato1bOtherMaps <br> Percolation1to1bForestMaps <br> Percolation1ato1bIrrigationMaps <br> Percolation1bto2OtherMaps <br> Percolation1bto2ForestMaps <br> Percolation1bto2IrrigationMaps                                    | Percolation1ato1bOther <br> Percolation1ato1bForest <br> Percolation1to2Irrigation <br> Percolation1bto2Other <br> Percolation1bto2Forest <br> Percolation1bto2Irrigation           |

src/lisflood/global_modules/netcdf.py

@@ -2,6 +2,7 @@
 import glob
 import warnings
 import xarray as xr
+import dask
 import numpy as np
 import datetime
 import pcraster
@@ -260,9 +261,20 @@ def load_next_chunk(self):
         self.ichunk += 1
         begin = self.chunk_indexes[self.ichunk]
         end = self.chunk_indexes[self.ichunk+1]
-        
+
         chunk = self.dataset.isel(time=range(begin, end))
-        self.dataset_chunk = chunk.load()  # triggers xarray computation
+        # Force the single-threaded (synchronous) dask scheduler for this read.
+        # The netCDF4/HDF5 libraries shipped by conda-forge are built without
+        # thread-safety (nompi). Under dask's default threaded scheduler, several
+        # worker threads read from the same HDF5 file concurrently while, at the
+        # same time, garbage collection in another thread may close a netCDF file
+        # handle (NetCDF4DataStore.__del__ -> close). Concurrent read/close on a
+        # non-threadsafe HDF5 corrupts its global state and aborts the process
+        # (Fatal Python error: Aborted / SIGABRT). Reads were already serialized
+        # by xarray's HDF5 lock, so running synchronously removes the crash with
+        # negligible performance impact.
+        with dask.config.set(scheduler='synchronous'):
+            self.dataset_chunk = chunk.load()  # triggers xarray computation
 
     def __getitem__(self, step):
 
@@ -439,7 +451,7 @@ def write_netcdf_header(settings,
                         start_date,
                         rep_steps,
                         frequency,
-                        ):
+                        map_value=None):
 
     """ Writes a netcdf header without the data inside
 
@@ -465,7 +477,11 @@ def write_netcdf_header(settings,
         list of reporting steps
     frequency:
         output frequency (all, monthly or yearly)
-
+    map_value: ReportedMap namedtuple or None, optional
+        Output variable metadata from default_options.py. When provided and OutputPacking is
+        enabled, its scale_factor and add_offset attributes are used for int16 packing.
+        If None (default), no packing is applied.
+
     Returns
     -------
     object
@@ -569,7 +585,25 @@ def write_netcdf_header(settings,
             time.units = 'minutes since %s' % start_date.strftime("%Y-%m-%d %H:%M:%S.0")
         nf1.variables["time"][:] = date2num(time_stamps, time.units, time.calendar)
 
-        value = nf1.createVariable(var_name, dtype, ('time', dim_lat_y, dim_lon_x), zlib=True, fill_value=-9999, chunksizes=(1, nrow, ncol))
+        # value = nf1.createVariable(var_name, dtype, ('time', dim_lat_y, dim_lon_x), zlib=True, fill_value=-9999, chunksizes=(1, nrow, ncol))
+        # Packing: use int16 with CF scale/offset if enabled and variable has packing metadata
+        packing_enabled = binding.get('OutputPacking', 'False') == 'True'
+        has_packing = (map_value is not None
+                       and getattr(map_value, 'scale_factor', None) is not None
+                       and getattr(map_value, 'add_offset', None) is not None)
+        if packing_enabled and has_packing:
+            var_dtype = 'i2'
+            var_fill = default_fillvals['i2']
+        else:
+            var_dtype = dtype
+            var_fill = -9999
+
+        value = nf1.createVariable(var_name, var_dtype, ('time', dim_lat_y, dim_lon_x),
+                                   zlib=True, fill_value=var_fill, chunksizes=(1, nrow, ncol))
+        if packing_enabled and has_packing:
+            value.scale_factor = np.float64(map_value.scale_factor)
+            value.add_offset = np.float64(map_value.add_offset)
+
     else:
         value = nf1.createVariable(var_name, dtype, (dim_lat_y, dim_lon_x), zlib=True, fill_value=-9999)
 

src/lisflood/global_modules/output.py

@@ -15,17 +15,27 @@
 
 """
 import os
+import datetime
 import numpy as np
 from pcraster import ifthen, catchmenttotal, mapmaximum
 import sys
+import warnings
 
 from .zusatz import TimeoutputTimeseries
 from .add1 import decompress, valuecell, loadmap, compressArray
 from .netcdf import write_netcdf_header, iterOpenNetcdf, nanCheckMap, uncompress_array
 from .errors import LisfloodFileError, LisfloodWarning
-from .settings import inttodate, CDFFlags, LisSettings
+from .settings import inttodate, CDFFlags, LisSettings, MaskInfo
+from netCDF4 import default_fillvals
 
 
+# ------------------------------------------------------------------------
+# Packing constants for int16 CF scale/offset encoding
+# ------------------------------------------------------------------------
+PACK_FILL = np.int16(default_fillvals['i2'])      # -32767
+PACK_MIN = PACK_FILL + 1                          # -32766
+PACK_MAX = np.iinfo(np.int16).max                 #  32767
+
 # ------------------------------------------------------------------------
 # Writer classes
 # ------------------------------------------------------------------------
@@ -95,7 +105,7 @@ def write(self, start_date, rep_steps):
         if self.data is not None:
             nf1 = write_netcdf_header(self.settings, self.map_name, self.map_path, self.var.DtDay,
                                     self.map_key, self.map_value.output_var, self.map_value.unit,
-                                    start_date, rep_steps, self.frequency)
+                                    start_date, rep_steps, self.frequency, map_value=self.map_value)
 
             map_np = uncompress_array(self.data)
 
@@ -151,16 +161,41 @@ def write(self, start_date, rep_steps):
                 if self.step_range[0] == 0:
                     nf1 = write_netcdf_header(self.settings, self.map_name, self.map_path, self.var.DtDay,
                                             self.map_key, self.map_value.output_var, self.map_value.unit,
-                                            start_date, rep_steps, self.frequency)
+                                            start_date, rep_steps, self.frequency, map_value=self.map_value)
                 else:
                     nf1 = iterOpenNetcdf(self.map_path, "", 'a', format='NETCDF4')
 
+                nc_var = nf1.variables[self.map_name]
+                nc_var.set_auto_maskandscale(False)
+                is_packed = nc_var.dtype == np.int16
+                if is_packed:
+                    scale = nc_var.scale_factor
+                    offset = nc_var.add_offset
+                    nodata_mask = MaskInfo.instance().info.mask
+
                 for step, data in zip(self.step_range, self.data_steps):
-                    nf1.variables[self.map_name][step, :, :] = uncompress_array(data)
+                    map_np = uncompress_array(data)
+                    if is_packed:
+                        packed = np.round((map_np - offset) / scale).astype(np.float64)
+                        clipped = ((packed < -32767) | (packed > 32767)) & (map_np != -9999)
+                        if clipped.any():
+                            vmin = offset + scale * (-32767)
+                            vmax = offset + scale * 32767
+                            warnings.warn(LisfloodWarning(
+                                f"OutputPacking: {clipped.sum()} values in '{self.map_name}' outside "
+                                f"packing range [{vmin:.4g}, {vmax:.4g}] and will be clipped."
+                            ))
+                        packed = np.clip(packed, PACK_MIN, PACK_MAX)
+                        packed[nodata_mask] = PACK_FILL
+                        nc_var[step, :, :] = packed.astype(np.int16)
+                    else:
+                        # For non-packed: convert masked array to plain array with fill value
+                        if hasattr(map_np, 'filled'):
+                            map_np = map_np.filled(-9999)
+                        nc_var[step, :, :] = map_np
 
                 nf1.close()
 
-                # clear lists for next chunk
                 self.step_range.clear()
                 self.data_steps.clear()
             else:
@@ -381,6 +416,87 @@ def _start_date(self):
     def _rep_steps(self):
         return self._rep_steps_val
 
+
+class MapOutputAggregated(MapOutput):
+    """Handles temporal aggregation (monthly/yearly mean/sum) for a variable."""
+
+    def __init__(self, var, map_key, map_value, frequency, operation):
+        out_type = 'all'  # accumulates every timestep
+        settings = LisSettings.instance()
+        binding = settings.binding
+        self._start_date_val = var.CalendarDayStart
+        self._rep_steps_val = range(binding['StepStartInt'], binding['StepEndInt'] + 1)
+
+        self._operation = operation  # 'mean' or 'sum'
+        self._accum_buffer = None
+        self._accum_count = 0
+        self._write_step = 0  # own step counter for NetCDF time dimension
+
+        # Disable int16 packing for sum aggregates — monthly/yearly sums can exceed
+        # the int16 range calibrated for daily values. Mean aggregates stay within
+        # the same value range as daily output, so packing remains valid.
+        if operation == 'sum':
+            map_value_no_pack = map_value._replace(scale_factor=None, add_offset=None)
+        else:
+            map_value_no_pack = map_value
+
+        super().__init__(var, out_type, frequency, map_key, map_value_no_pack)
+
+        # Force immediate write for aggregated outputs (one slice per period)
+        if hasattr(self, 'writer') and hasattr(self.writer, 'chunks'):
+            self.writer.chunks = 1
+
+    def _output_checkpoint(self):
+        """Always True — we accumulate every timestep."""
+        return True
+
+    @property
+    def _start_date(self):
+        return self._start_date_val
+
+    @property
+    def _rep_steps(self):
+        return self._rep_steps_val
+
+    def stage(self):
+        """Accumulate instead of storing instantaneous values."""
+        self.step = self.var.currentTimeStep()
+        map_np = self.writer._extract_map()
+
+        if self._accum_buffer is None:
+            self._accum_buffer = np.zeros_like(map_np)
+
+        self._accum_buffer += map_np
+        self._accum_count += 1
+
+    def write(self):
+        """Write only at period boundary (month-end or year-end)."""
+        current_date = self.var.CalendarDate
+        next_date = current_date + datetime.timedelta(days=self.var.DtDay)
+        if self.frequency == 'monthly':
+            is_boundary = current_date.month != next_date.month
+        elif self.frequency == 'yearly':
+            is_boundary = current_date.year != next_date.year
+        else:
+            is_boundary = True
+
+        if is_boundary and self._accum_buffer is not None:
+            # Finalize
+            if self._operation == 'mean':
+                result = self._accum_buffer / self._accum_count
+            else:  # sum
+                result = self._accum_buffer
+
+            # Stage the aggregated result into the writer
+            self.writer.data_steps.append(result)
+            self.writer.step_range.append(self._write_step)
+            self.writer.write(self._start_date, self._rep_steps)
+
+            # Increment own step counter and reset accumulator
+            self._write_step += 1
+            self._accum_buffer = None
+            self._accum_count = 0
+
 # ------------------------------------------------------------------------
 # Output factory
 # ------------------------------------------------------------------------
@@ -428,11 +544,38 @@ def __init__(self, var):
             if out.is_valid():
                 outputs.append(out)
 
-        check_duplicates = []
+        # --- Temporal aggregation outputs ---
+        binding = settings.binding
+        aggregation_configs = {
+            'OutputMonthlyMean': ('monthly', 'mean'),
+            'OutputMonthlySum': ('monthly', 'sum'),
+            'OutputYearlyMean': ('yearly', 'mean'),
+            'OutputYearlySum': ('yearly', 'sum'),
+        }
+        aggregated_vars = set()  # track which vars are aggregated
+
+        reportedmaps = settings.options['reportedmaps']
+        for setting_key, (frequency, operation) in aggregation_configs.items():
+            var_list = binding.get(setting_key, '').split(';')
+            for var_name in var_list:
+                var_name = var_name.strip()
+                if var_name and var_name in reportedmaps:
+                    map_value = reportedmaps[var_name]
+                    out = MapOutputAggregated(var, var_name, map_value, frequency, operation)
+                    if out.is_valid():
+                        outputs.append(out)
+                        aggregated_vars.add(var_name)
+
+        # Remove normal outputs for variables that are now aggregated
         outputs_clean = []
+        check_duplicates = []
         for out in outputs:
+            # Skip normal Maps/All output if variable is aggregated
+            if hasattr(out, 'map_key') and out.map_key in aggregated_vars:
+                if not isinstance(out, (MapOutputEnd, MapOutputAggregated)):
+                    continue
             if out.map_path in check_duplicates:
-                print(f'Warning! Output map {out.map_path} is duplicated, check list of outputs')
+                print(f'Warning! Output map {out.map_path} is duplicated')
             else:
                 check_duplicates.append(out.map_path)
                 outputs_clean.append(out)

src/lisflood/global_modules/settings.py

@@ -841,6 +841,14 @@ def __init__(self, model):
         msg += "\t[X] The simulation output as specified in the settings file can be found in {}\n".format(out_dir)
         msg += "\t[X] Activated modules: {}\n".format(activated_options)
         msg += "\t[X] Report options: {}\n".format(activated_rep)
+        # Packing and aggregation info
+        binding = settings.binding
+        if binding.get('OutputPacking', 'False') == 'True':
+            msg += "\t[X] Output Packing: int16 scale/offset enabled\n"
+        for agg_key in ['OutputMonthlyMean', 'OutputMonthlySum', 'OutputYearlyMean', 'OutputYearlySum']:
+            agg_val = binding.get(agg_key, '').strip()
+            if agg_val:
+                msg += "\t[X] {}: {}\n".format(agg_key, agg_val)
         self._msg = '{}{}'.format(header, msg)
 
     def __str__(self):