From 3f76bc29f059214dbc1ce8689d8e52c53ac056f1 Mon Sep 17 00:00:00 2001 From: Timo Schaffhauser Date: Tue, 18 Aug 2026 12:30:58 +0200 Subject: [PATCH 1/7] groundwater soil wording correction in docs --- docs/5_annex_output-files/index.md | 2 +- implementation_guide_scale_offset_packing.md | 354 ++++++ lisflood_optimization_report.md | 284 +++++ lisflood_optimization_report_v2.md | 1042 ++++++++++++++++++ 4 files changed, 1681 insertions(+), 1 deletion(-) create mode 100644 implementation_guide_scale_offset_packing.md create mode 100644 lisflood_optimization_report.md create mode 100644 lisflood_optimization_report_v2.md diff --git a/docs/5_annex_output-files/index.md b/docs/5_annex_output-files/index.md index fd7cfbf0..5a0c8761 100644 --- a/docs/5_annex_output-files/index.md +++ b/docs/5_annex_output-files/index.md @@ -61,7 +61,7 @@ Output time series can be classified in the following categories: | infiltration | $\frac{mm}{timestep}$ | InfiltrationAvUpsTS | infiltrationUps.tss | | preferential (bypass) flow | $\frac{mm}{timestep}$ | PrefFlowAvUpsTS | prefFlowUps.tss | | percolation upper to lower soil layer | $\frac{mm}{timestep}$ | PercolationAvUpsTS | dTopToSubUps.tss | -| percolation lower soil layer to subsoil | $\frac{mm}{timestep}$ | SeepSubToGWAvUpsTS | dSubToUzUps.tss | +| percolation lower soil layer to Upper Zone | $\frac{mm}{timestep}$ | SeepSubToGWAvUpsTS | dSubToUzUps.tss | | surface runoff | $\frac{mm}{timestep}$ | SurfaceRunoffAvUpsTS | surfaceRunoffUps.tss | | outflow from upper zone | $\frac{mm}{timestep}$ | UZOutflowAvUpsTS | qUzUps.tss | | outflow from lower zone | $\frac{mm}{timestep}$ | LZOutflowAvUpsTS | qLzUps.tss | diff --git a/implementation_guide_scale_offset_packing.md b/implementation_guide_scale_offset_packing.md new file mode 100644 index 00000000..efb01cfa --- /dev/null +++ b/implementation_guide_scale_offset_packing.md @@ -0,0 +1,354 @@ +# Implementation Guide: CF Scale/Offset Packing for LISFLOOD Outputs + +## Overview + +This guide walks through adding int16 scale/offset packing to LISFLOOD output NetCDF files. +The feature stores floating-point output variables as packed int16 values using the CF-convention +`scale_factor` and `add_offset` attributes. All CF-compliant readers (xarray, CDO, NCO, QGIS) +automatically unpack on read — no downstream changes needed. + +**Result:** 75-90% smaller output files compared to float64 defaults. + +**Design principles:** +- Global on/off via a simple `True`/`False` setting (`OutputPacking`) +- Per-variable scale/offset stored directly in the `ReportedMap` namedtuple attributes +- State/end maps (for warm start) are never packed +- Range documentation kept in a separate `.md` reference file + +--- + +## Step 1: Add the setting to the reference settings XML + +**File:** `src/lisfloodSettings_reference.xml` + +Find the block where `OutputMapsDataType` is defined and add after it: + +```xml + + +The option "OutputPacking" enables CF-convention scale_factor/add_offset packing +of output maps into int16 (2 bytes per value instead of 4 or 8). + - "False" (default): write raw floating-point values (dtype from OutputMapsDataType) + - "True": pack into signed 16-bit integers using per-variable scale/offset + Readers automatically unpack using: value = packed * scale_factor + add_offset +Note: State/end maps used for warm starts are NEVER packed (always full precision). + + +``` + +Also add the pass-through in the bindings section (where `OutputMapsDataType` is passed): + +```xml + +``` + + +--- + +## Step 2: Extend the `ReportedMap` namedtuple + +**File:** `src/lisflood/global_modules/default_options.py` + +Change: + +```python +ReportedMap = namedtuple('ReportedMap', 'name, output_var, unit, end, steps, all, restrictoption, monthly, yearly') +``` + +To: + +```python +ReportedMap = namedtuple('ReportedMap', 'name, output_var, unit, end, steps, all, restrictoption, monthly, yearly, scale_factor, add_offset') +ReportedMap.__new__.__defaults__ = (None, None) # scale_factor and add_offset default to None (= no packing) +``` + +Setting `__new__.__defaults__` means all existing `ReportedMap(...)` entries remain valid +without modification — they'll get `scale_factor=None, add_offset=None` automatically. +Only variables you explicitly want to pack need the extra two fields. + +--- + +## Step 3: Add scale/offset values to output variables + +**File:** `src/lisflood/global_modules/default_options.py` + +For each variable you want to pack, add `scale_factor` and `add_offset` at the end of its +`ReportedMap` entry. Use the helper formula: + +``` +scale_factor = (physical_max - physical_min) / 65534 +add_offset = physical_min + scale_factor * 32767 +``` + +Examples: + +```python +'DischargeMaps': ReportedMap(name='DischargeMaps', output_var='ChanQAvg', + unit='m3/s', end=[], steps=[], + all=['repDischargeMaps'], restrictoption=[], + monthly=False, yearly=False, + scale_factor=3.052, # range 0-200000 m3/s + add_offset=100001.5), + +'WaterLevelMaps': ReportedMap(name='WaterLevelMaps', output_var='WaterLevel', + unit='m', end=[], steps=[], + all=['repWaterLevelMaps'], restrictoption=['nonInit'], + monthly=False, yearly=False, + scale_factor=6.104e-4, # range -10 to +30 m + add_offset=10.0), + +'SnowMaps': ReportedMap(name='SnowMaps', output_var='SnowCover', + unit='mm', end=[], steps=[], + all=['repSnowMaps'], restrictoption=['nonInit'], + monthly=False, yearly=False, + scale_factor=0.04578, # range 0-3000 mm + add_offset=1500.0), +``` + +Variables that should NOT be packed (state/end maps) simply omit the fields: + +```python +'ChanQEnd': ReportedMap(name='ChanQEnd', output_var='ChanQ', unit='m3/s', + end=['repEndMaps'], steps=[], all=[], + restrictoption=[], monthly=False, yearly=False), + # No scale_factor/add_offset → defaults to None → never packed +``` + + +--- + +## Step 4: Modify `write_netcdf_header()` to support packing + +**File:** `src/lisflood/global_modules/netcdf.py` + +### 4a. Change the function signature + +Add `map_value=None` as a parameter: + +```python +def write_netcdf_header(settings, var_name, netfile, DtDay, + value_standard_name, value_long_name, value_unit, + start_date, rep_steps, frequency, + map_value=None): # <-- NEW +``` + +### 4b. Replace the variable creation logic + +Find the block at the end of the function where the NetCDF variable is created. +Replace with: + +```python + # Determine if packing is active for this variable + 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 frequency is not None: # output file with "time" dimension + if packing_enabled and has_packing: + # CF scale/offset packing into int16 + value = nf1.createVariable(var_name, 'i2', ('time', dim_lat_y, dim_lon_x), + zlib=True, fill_value=np.int16(-32768), + chunksizes=(1, nrow, ncol)) + value.scale_factor = np.float64(map_value.scale_factor) + value.add_offset = np.float64(map_value.add_offset) + else: + # Standard float output (current behaviour) + value = nf1.createVariable(var_name, dtype, ('time', dim_lat_y, dim_lon_x), + zlib=True, fill_value=-9999, chunksizes=(1, nrow, ncol)) + else: + # End/state maps — NEVER pack (need full precision for warm start) + value = nf1.createVariable(var_name, dtype, (dim_lat_y, dim_lon_x), + zlib=True, fill_value=-9999) +``` + + +--- + +## Step 5: Pass `map_value` through the output call chain + +**File:** `src/lisflood/global_modules/output.py` + +In both `NetcdfWriter.write()` and `NetcdfStepsWriter.write()`, pass `self.map_value` to +`write_netcdf_header`: + +```python +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, + map_value=self.map_value) # <-- ADD THIS +``` + +--- + +## Step 6: Modify the write methods to pack data before storing + +**File:** `src/lisflood/global_modules/output.py` + +### 6a. In `NetcdfStepsWriter.write()`: + +Replace the data-writing loop: + +```python +for step, data in zip(self.step_range, self.data_steps): + map_np = uncompress_array(data) + + nc_var = nf1.variables[self.map_name] + if nc_var.dtype == np.int16: + # Pack float → int16 + scale = nc_var.scale_factor + offset = nc_var.add_offset + packed = np.round((map_np - offset) / scale).astype(np.float64) + packed = np.clip(packed, -32767, 32767) + packed[map_np == -9999] = -32768 # fill value + nc_var.set_auto_maskandscale(False) # CRITICAL: prevent double-packing + nc_var[step, :, :] = packed.astype(np.int16) + else: + nc_var[step, :, :] = map_np +``` + +**Critical:** `set_auto_maskandscale(False)` MUST be called before writing. Without it, the +netCDF4 library applies scale/offset again on write, resulting in double-packing (garbage values). + + +--- + +## Step 7: Testing + +### Quick verification script + +```python +import xarray as xr +import numpy as np + +# Open a packed output file — xarray auto-unpacks +ds = xr.open_dataset('path/to/output.nc') + +print(ds['DischargeMaps'].dtype) # float64 (auto-unpacked) +print(ds['DischargeMaps'].values[:5]) # physically reasonable values +print(ds['DischargeMaps'].encoding) # shows scale_factor, add_offset, dtype=int16 +``` + +### Comparison test + +Run the same simulation with `OutputPacking = False` and `OutputPacking = True`: + +```python +import xarray as xr +import numpy as np + +ref = xr.open_dataset('output_nopacking.nc') +packed = xr.open_dataset('output_packed.nc') + +var = 'DischargeMaps' +diff = np.abs(ref[var].values - packed[var].values) +print(f"Max absolute difference: {diff.max():.4f}") +print(f"Max relative difference: {(diff / (np.abs(ref[var].values) + 1e-10)).max():.6e}") + +# For discharge with range 0-200000: +# scale_factor ≈ 3.05 +# Max error ≈ ±1.5 m3/s (half a quantization step) +``` + +### File size comparison + +```bash +dir output_nopacking.nc output_packed.nc +# Expected: packed file ≈ 20-25% of original size +``` + + +--- + +## Summary of files to modify + +| # | File | What to do | +|---|------|------------| +| 1 | `src/lisfloodSettings_reference.xml` | Add `OutputPacking` textvar (True/False) + pass-through | +| 2 | `src/lisflood/global_modules/default_options.py` | Extend `ReportedMap` namedtuple with `scale_factor, add_offset` | +| 3 | `src/lisflood/global_modules/default_options.py` | Add scale/offset values to specific output variable entries | +| 4 | `src/lisflood/global_modules/netcdf.py` | Add `map_value` param to `write_netcdf_header()`; int16 creation logic | +| 5 | `src/lisflood/global_modules/output.py` | Pass `map_value` to `write_netcdf_header()` | +| 6 | `src/lisflood/global_modules/output.py` | Add packing logic in write methods | +| 7 | `docs/packing_ranges_reference.md` | Create documentation file with range/precision table | + +--- + +## How the flag works + +``` +User sets OutputPacking = "True" in settings XML + │ + ▼ +write_netcdf_header() reads binding['OutputPacking'] + │ + ▼ + ┌─────────────────────────────────┐ + │ packing_enabled = True │ + │ has_packing = map_value has │ + │ scale_factor and add_offset? │ + └──────────┬──────────────────────┘ + │ + ┌──────┴──────┐ + │ │ + has_packing no packing attrs + = True = False (None) + │ │ + ▼ ▼ + Create int16 Create float + variable with variable (normal + scale/offset behaviour) +``` + +Three conditions must ALL be true for packing to happen: +1. `OutputPacking = "True"` in settings (global toggle) +2. The variable's `ReportedMap` has `scale_factor` and `add_offset` defined (per-variable control) +3. The output has a time dimension (`frequency is not None`) — state/end maps are excluded + +--- + +## Gotchas + +1. **Double-packing**: Call `nc_var.set_auto_maskandscale(False)` before writing raw int16 values. Otherwise netCDF4 applies scale/offset again. + +2. **Fill value**: Use `-32768` (int16 minimum). The packed data range uses -32767 to +32767 (65534 levels). NetCDF4/xarray will mask cells with fill value as NaN on read. + +3. **Out-of-range clipping**: If the model produces values outside the defined range, they get clipped. Set ranges generously — better to waste quantization levels than clip real data. + +4. **State files**: The `frequency is None` guard ensures end maps are never packed. Don't add `scale_factor`/`add_offset` to End/State `ReportedMap` entries either (belt and braces). + +5. **Variable name matching**: The `map_value` object is passed directly from the output writer, so no name-to-range matching is needed — the metadata travels with the variable. + +--- + +## Helper: Computing scale/offset from a physical range + +```python +def compute_packing_params(vmin, vmax): + """Compute CF-convention scale_factor and add_offset for int16 packing. + + Parameters + ---------- + vmin : float - Minimum physical value + vmax : float - Maximum physical value + + Returns + ------- + scale_factor : float + add_offset : float + """ + n_levels = 65534.0 # int16 usable range: -32767 to +32767 + scale_factor = (vmax - vmin) / n_levels + add_offset = vmin + scale_factor * 32767.0 + return scale_factor, add_offset +``` + +Example outputs: +``` +compute_packing_params(0, 200000) → (3.052, 100001.5) discharge m3/s +compute_packing_params(0, 1.0) → (1.526e-5, 0.5) fraction +compute_packing_params(-10, 30) → (6.104e-4, 10.0) water level m +compute_packing_params(0, 3000) → (0.04578, 1500.0) snow mm +compute_packing_params(0, 500) → (0.00763, 250.0) precip mm/day +compute_packing_params(0, 5000) → (0.07630, 2500.0) groundwater mm +``` diff --git a/lisflood_optimization_report.md b/lisflood_optimization_report.md new file mode 100644 index 00000000..32502307 --- /dev/null +++ b/lisflood_optimization_report.md @@ -0,0 +1,284 @@ +# LISFLOOD I/O Optimization Report + +## Executive Summary + +The LISFLOOD hydrological model codebase has been analyzed for I/O performance and storage optimization opportunities. Key findings: + +1. **Output Data Type**: Default is `float64` (8 bytes) but can be set to `float32` (4 bytes) - potential **50% storage reduction** for output maps +2. **NetCDF Compression**: Already uses `zlib=True` for output files - good practice +3. **Input Reading**: Uses xarray with chunking for NetCDF inputs, but data is loaded into memory as float64 by default +4. **Caching**: Has `MapsCaching` option for static maps but not enabled by default +5. **Output Chunks**: Uses `OutputMapsChunks` setting (default 1) - could benefit from larger chunks + +## Codebase Overview + +### File Structure +- **Language**: Python (with Numba for performance-critical sections) +- **Main modules**: + - `src/lisflood/global_modules/netcdf.py` - NetCDF I/O handling (584 lines) + - `src/lisflood/global_modules/output.py` - Output writing (586 lines) + - `src/lisflood/global_modules/add1.py` - Map loading/compression (986 lines) + - `src/lisflood/Lisflood_dynamic.py` - Main time loop (269 lines) + - `src/lisflood/hydrological_modules/` - Hydrological process modules + +### Key I/O Files +| File | Purpose | +|------|---------| +| [`netcdf.py`](src/lisflood/global_modules/netcdf.py) | NetCDF reading/writing, xarray chunked readers | +| [`output.py`](src/lisflood/global_modules/output.py) | Output writers (NetCDF, PCRaster) | +| [`add1.py`](src/lisflood/global_modules/add1.py) | Map loading, compression, decompression | +| [`zusatz.py`](src/lisflood/global_modules/zusatz.py) | NetCDF file access utilities | + +--- + +## I/O Bottlenecks + +### Identified Issues + +| Location | Issue | Estimated Impact | Complexity | +|----------|-------|------------------|------------| +| [`netcdf.py:478`](src/lisflood/global_modules/netcdf.py:478) | Output dtype from binding is `float64` by default | High (50% storage) | Low | +| [`netcdf.py:572-574`](src/lisflood/global_modules/netcdf.py:572) | zlib compression with chunksizes=(1,nrow,ncol) - time dimension chunk=1 | Medium | Medium | +| [`add1.py:282`](src/lisflood/global_modules/add1.py:282) | `compressArray` converts to float64: `mapC.astype(float)` | Medium | Low | +| [`readmeteo.py:40`](src/lisflood/hydrological_modules/readmeteo.py:40) | xarray readers created per variable in __init__ | Low | Low | +| [`output.py:102`](src/lisflood/global_modules/output.py:102) | Single timestep writes: `nf1.variables[self.map_name][:, :] = map_np` | Medium | Medium | +| [`netcdf.py:265`](src/lisflood/global_modules/netcdf.py:265) | `chunk.load()` loads data synchronously | Low | Medium | + +### Detailed Analysis + +#### 1. Output Data Type (High Impact) +**Location**: [`lisfloodSettings_reference.xml:157`](src/lisfloodSettings_reference.xml:157) +```xml + +``` + +**Current behavior**: All output maps are written as float64 (8 bytes per value) + +**Recommendation**: Change default to `float32` - sufficient for most hydrological variables + +#### 2. Compression Chunk Size (Medium Impact) +**Location**: [`netcdf.py:572`](src/lisflood/global_modules/netcdf.py:572) +```python +value = nf1.createVariable(var_name, dtype, ('time', dim_lat_y, dim_lon_x), + zlib=True, fill_value=-9999, chunksizes=(1, nrow, ncol)) +``` + +**Current behavior**: Time chunk size = 1 (one timestep per chunk) + +**Recommendation**: Use larger time chunks (e.g., daily or monthly) for better compression ratio + +#### 3. Array Type Conversion (Medium Impact) +**Location**: [`add1.py:282`](src/lisflood/global_modules/add1.py:282) +```python +return mapC.astype(float) # This defaults to float64 +``` + +**Current behavior**: All compressed arrays become float64 regardless of input precision + +**Recommendation**: Preserve input dtype or allow explicit dtype specification + +--- + +## Data Type Optimization Opportunities + +### Variable Analysis + +| Variable | Current dtype | Suggested dtype | Range | Storage Saving | Risk | +|----------|---------------|-----------------|-------|----------------|------| +| Discharge (ChanQ) | float64 | float32 | 0-100000 m³/s | 50% | Low | +| Soil Moisture (W1, W2) | float64 | float32 | 0-1 (fraction) | 50% | Very Low | +| Groundwater (UZ, LZ) | float64 | float32 | 0-∞ mm | 50% | Low | +| Precipitation | float64 | float32 | 0-500 mm/day | 50% | Very Low | +| Temperature | float64 | float32 | -50 to 50 °C | 50% | Very Low | +| Snow Cover | float64 | float32 | 0-100% | 50% | Very Low | +| Channel Storage | float64 | float32 | 0-∞ m³ | 50% | Low | +| LDD (flow directions) | int8 | int8 | 0-9 | 0% | N/A | +| Lake/Reservoir IDs | int16 | int16 | 0-65535 | 0% | N/A | +| Mask maps | boolean | boolean | 0-1 | 0% | N/A | + +### Compression Encoding Opportunities + +For integer variables that could benefit from scale_factor/add_offset (CF convention): + +| Variable | Current | Suggested | Notes | +|----------|---------|-----------|-------| +| LDD | int8 | int8 (no change) | Already optimal | +| Land Use | int8/int16 | int16 + scale_factor | Could use packing | +| Lake IDs | int16 | int16 (no change) | Already optimal | + +--- + +## Compression & Encoding Recommendations + +### 1. Enable float32 as Default +**File**: [`lisfloodSettings_reference.xml`](src/lisfloodSettings_reference.xml:157) +```xml + +``` + +### 2. Optimize NetCDF Chunk Sizes +**File**: [`netcdf.py`](src/lisflood/global_modules/netcdf.py:572) + +Current: +```python +chunksizes=(1, nrow, ncol) +``` + +Recommended: +```python +# For daily output: chunk = 1 day +# For monthly output: chunk = 30 days +optimal_time_chunk = min(30, n_timesteps) # Use up to 30 days +chunksizes=(optimal_time_chunk, nrow, ncol) +``` + +### 3. Add Compression Level Option +**File**: [`netcdf.py:572`](src/lisflood/global_modules/netcdf.py:572) + +Current: +```python +zlib=True +``` + +Recommended: +```python +# Add setting for compression level (1-9, default 4 for balanced speed/ratio) +compression_level = int(binding.get('OutputCompressionLevel', 4)) +value = nf1.createVariable(..., zlib=True, complevel=compression_level) +``` + +### 4. Preserve Input Data Types +**File**: [`add1.py:282`](src/lisflood/global_modules/add1.py:282) + +Current: +```python +return mapC.astype(float) +``` + +Recommended: +```python +# Preserve dtype or allow specification +def compressArray(map, pcr=True, name=None, dtype=None): + # ... existing code ... + if dtype is None: + dtype = map.dtype if hasattr(map, 'dtype') else np.float64 + return mapC.astype(dtype) +``` + +--- + +## Prioritized Action Plan + +| Rank | Recommendation | Storage Impact | Runtime Impact | Complexity | Risk | Score | +|------|----------------|----------------|----------------|------------|------|-------| +| 1 | Change default OutputMapsDataType to float32 | 50% | ~5% | Low | No | 55 | +| 2 | Optimize NetCDF time chunk sizes | 10-30% | 5-15% | Medium | No | 20 | +| 3 | Add compression level setting | 5-15% | -5-10% | Low | No | 15 | +| 4 | Enable MapsCaching by default | 0% | 10-30% | Low | No | 10 | +| 5 | Preserve input dtype in compressArray | 0-20% | 0% | Low | No | 5 | +| 6 | Add async I/O for output writing | 0% | 10-20% | High | Maybe | 3 | + +*Score = (Storage Impact + Runtime Impact) / Complexity* + +--- + +## Code Snippets / Diff Examples + +### Recommendation 1: Change Default Output Data Type + +**File**: `src/lisfloodSettings_reference.xml` + +```diff +- ++ +``` + +**Impact**: 50% reduction in output file sizes for all NetCDF maps + +**Risk**: Very low - float32 provides ~7 significant digits, sufficient for all hydrological variables + +--- + +### Recommendation 2: Optimize NetCDF Time Chunk Sizes + +**File**: `src/lisflood/global_modules/netcdf.py` + +```diff +@@ -569,7 +569,12 @@ 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)) ++ # Optimize chunk size for better compression ++ n_timesteps = steps.size ++ # Use up to 30 days of data per chunk, or less if total timesteps < 30 ++ time_chunk = min(30, n_timesteps) if n_timesteps > 1 else 1 ++ chunks = (time_chunk, nrow, ncol) ++ value = nf1.createVariable(var_name, dtype, ('time', dim_lat_y, dim_lon_x), zlib=True, fill_value=-9999, chunksizes=chunks) +``` + +**Impact**: 10-30% better compression ratio, faster I/O for large datasets + +**Risk**: Low - chunking is internal to NetCDF, doesn't affect model results + +--- + +### Recommendation 3: Add Compression Level Setting + +**File**: `src/lisflood/global_modules/netcdf.py` + +```diff +@@ -571,7 +571,10 @@ 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)) ++ # Get compression level from settings (default 4 for balanced speed/ratio) ++ comp_level = int(binding.get('OutputCompressionLevel', 4)) ++ comp_level = max(1, min(9, comp_level)) # Clamp to valid range ++ value = nf1.createVariable(var_name, dtype, ('time', dim_lat_y, dim_lon_x), zlib=True, fill_value=-9999, chunksizes=chunks, complevel=comp_level) +``` + +**File**: `src/lisfloodSettings_reference.xml` (add new setting) + +```xml + +The option "OutputCompressionLevel" sets the zlib compression level (1-9). +Higher values give better compression but slower I/O. Default is 4. + + +``` + +**Impact**: 5-15% additional compression, tunable performance tradeoff + +**Risk**: Low - compression is lossless, doesn't affect model results + +--- + +## Additional Recommendations + +### Enable MapsCaching by Default +The `MapsCaching` option already exists but defaults to "False". Enabling it by default would cache static maps (DEM, land use, soil properties) in memory, avoiding repeated disk reads. + +**File**: `src/lisfloodSettings_reference.xml` +```xml + +``` + +**Impact**: 10-30% faster initialization for large models with many static maps + +**Risk**: Low - only affects static data, doesn't change model results + +### Consider Float32 for Internal Calculations +For memory-constrained systems, consider using float32 for internal state variables. This would require: +1. Testing to ensure numerical accuracy is maintained +2. Modifying array creation in initialization code +3. Potential impact on accumulation variables over long simulations + +--- + +## Summary + +The LISFLOOD codebase has a solid foundation for I/O operations with NetCDF compression already enabled. The primary optimization opportunity is changing the default output data type from float64 to float32, which would provide immediate 50% storage savings with no impact on model results. Secondary optimizations around chunk sizes and compression levels can provide additional 10-30% improvements. + +The implementation complexity for all recommended changes is low to medium, with minimal risk to model accuracy. The recommended changes are backward-compatible and can be implemented incrementally. \ No newline at end of file diff --git a/lisflood_optimization_report_v2.md b/lisflood_optimization_report_v2.md new file mode 100644 index 00000000..0a7236aa --- /dev/null +++ b/lisflood_optimization_report_v2.md @@ -0,0 +1,1042 @@ +# LISFLOOD Optimization Report v2 + +## Executive Summary + +This report identifies concrete opportunities to reduce wall-clock execution time and storage requirements across the LISFLOOD hydrological model codebase. It supersedes the earlier I/O-focused optimization report by covering the full computational pipeline including the newly added Muskingum-Cunge-Todini (MCT) routing module. + +Key findings: + +1. **MCT Routing Inner Loop** — The `mct_routing` kernel uses `prange` over pixels within each topological order, but orders are processed serially. The Newton-Raphson solver in `MCTRouting_single` (called twice per pixel) dominates runtime for MCT-enabled runs. +2. **Kinematic Wave Parallel Routing** — Already well-optimized with Numba `@njit(parallel=True)`, but the `numexpr` constant-term evaluation before calling the Numba kernel adds Python overhead every sub-step. +3. **Redundant Array Copies in Routing Loop** — Up to 10+ full-domain `.copy()` calls per routing sub-step when MCT + mass-balance reporting are active. +4. **I/O: Output Data Type** — Default `float64` output is still the factory setting; switching to `float32` halves output storage with no impact on results. +5. **I/O: NetCDF Chunk Size** — Time-dimension chunk of 1 limits compression ratio and read-back performance. +6. **Numba JIT Cold-Start** — First call to each `@njit` function triggers compilation; `cache=True` is already set but AOT pre-compilation could eliminate startup cost entirely. +7. **Memory: `compressArray` forces float64** — All compressed arrays are cast to `float64` regardless of input dtype. + + +--- + +## Part A — Runtime Performance + +--- + +### A1. MCT Routing: Serial Order Loop with Parallel Inner Loop + +| Attribute | Detail | +|-----------|--------| +| **File** | `src/lisflood/hydrological_modules/mct.py` | +| **Function** | `mct_routing()` (lines ~113-175) | +| **Impact** | **HIGH** | + +**Why it is a bottleneck:** +The MCT routing kernel is decorated `@njit(parallel=True)` but only the inner `prange(first, last)` loop over pixels within a single topological order is parallelised. The outer `for order in range(num_orders)` loop is inherently serial (data dependency: downstream pixels need upstream results). For domains with many orders but few pixels per order (long river stems), most iterations run with minimal parallelism. + +Inside the inner loop, `MCTRouting_single` is called per pixel. It contains: +- A Newton-Raphson iterative solver (`hoq`) with up to 1000 iterations per call +- Two full iterations of the MCT parameter calibration (`for i in range(2)`) +- Multiple calls to `qoh` (Manning-based Q-h relationships) + +For EFAS (~7 million pixels, ~5000 MCT pixels, ~200 topological orders), the MCT kernel can account for 30-50% of the routing sub-step time. + +**Recommendations:** + +1. **Reduce Newton-Raphson iterations in `hoq`**: The convergence tolerance is `1e-6` and max iterations is 1000. Profile to check actual average iteration count. If typically < 20, the overhead is acceptable. If convergence is slow for certain geometries, consider providing a better initial guess from the previous timestep's water depth (store `y_prev` as state). + +2. **Vectorise `MCTRouting_single` for batches**: Instead of calling a scalar function per pixel inside `prange`, restructure to pass arrays of pixel data for each order-batch and process them with vectorised NumPy/Numba array operations. This would enable SIMD and reduce function-call overhead. + +3. **Pre-compute static derived quantities**: `np.arctan(1 / ChanSdXdY[kinpix])` is computed every timestep for every pixel. Store `ANalv` as a pre-computed array in `MCTWave.__init__`. + +4. **Consider adaptive sub-stepping**: When Courant number < 0.5 for most pixels, the MCT solution is over-resolved. Allow the routing module to skip MCT computation for pixels where flow conditions have not materially changed (delta-Q threshold). + + +--- + +### A2. MCT Routing: Redundant `arctan` Computation Per Pixel Per Timestep + +| Attribute | Detail | +|-----------|--------| +| **File** | `src/lisflood/hydrological_modules/mct.py` | +| **Function** | `mct_routing()`, inner loop body | +| **Impact** | **MEDIUM** | + +**Current code (inside prange loop):** +```python +ANalv = np.arctan(1 / ChanSdXdY[kinpix]) +``` + +**Why it matters:** +`arctan` is a transcendental function computed per pixel, per routing sub-step, per model timestep. For EFAS with 4 routing sub-steps and 5000 MCT pixels, this is 20000 `arctan` calls per model timestep — all with static input. The result never changes. + +**Recommendation:** +Pre-compute `ANalv` once during `MCTWave.__init__()` and pass it as an array to `mct_routing`: + +```python +# In MCTWave.__init__: +self.ANalv = np.arctan(1.0 / ChanSdXdY) + +# In mct_routing signature: add ANalv parameter +# In the inner loop: replace np.arctan(...) with ANalv[kinpix] +``` + +**Estimated saving:** ~2-5% of MCT kernel time (removes transcendental from hot loop). + + +--- + +### A3. Kinematic Wave: `numexpr` Overhead Before Numba Kernel + +| Attribute | Detail | +|-----------|--------| +| **File** | `src/lisflood/hydrological_modules/kinematic_wave_parallel.py` | +| **Function** | `kinematicWave.kinematicWaveRouting()` | +| **Impact** | **MEDIUM** | + +**Current code:** +```python +lateral_inflow = nx.evaluate("q * dx", local_dict={"q": specific_lateral_inflow, "dx": self.space_delta}) +constant = nx.evaluate("a_dx_div_dt * Qold ** b + lateral_inflow", local_dict={...}) +``` + +**Why it matters:** +These `numexpr.evaluate` calls are executed every routing sub-step (typically 2-8 times per model timestep). Each call: +- Creates temporary arrays (full domain size) +- Involves Python-level dict construction and string parsing +- Runs a multi-threaded expression evaluator that competes with Numba's own thread pool + +For the kinematic routing (which applies to ALL channel pixels including those later overwritten by MCT), this adds measurable Python overhead on every sub-step. + +**Recommendation:** +Move the constant-term computation into the Numba `kinematicRouting` kernel itself. The expressions are simple element-wise operations that Numba can fuse into the main loop without allocating intermediates: + +```python +# Inside kinematicRouting (already @njit parallel): +# Replace 'constant' parameter with raw inputs +# Compute constant[pix] = a_dx_div_dt[pix] * discharge[pix]**beta + lateral_inflow[pix] +# directly at point of use +``` + +Alternatively, replace `numexpr` with a small `@njit` helper that computes `constant` in-place: + +```python +@njit(parallel=True, cache=True) +def compute_constant(constant, a_dx_div_dt, discharge, beta, lateral_inflow, space_delta, specific_lateral_inflow): + for pix in prange(constant.size): + constant[pix] = a_dx_div_dt[pix] * discharge[pix]**beta + specific_lateral_inflow[pix] * space_delta[pix] +``` + +**Estimated saving:** 5-10% of kinematic routing time (eliminates temporary array allocations and Python-level overhead per sub-step). + + +--- + +### A4. Routing Loop: Excessive `.copy()` Calls Per Sub-Step + +| Attribute | Detail | +|-----------|--------| +| **File** | `src/lisflood/hydrological_modules/routing.py` | +| **Function** | `routing.dynamic()` | +| **Impact** | **MEDIUM** | + +**Current code pattern (inside the `NoRoutSteps` loop):** +```python +SideflowChanM3 = self.var.ToChanM3RunoffDt.copy() +ChanQ_0 = self.var.ChanQ.copy() +ChanM3_0 = self.var.ChanM3.copy() +ChanQ = self.var.ChanQKin.copy() +ChanM3 = self.var.ChanM3Kin.copy() +ChanQAvgDt = self.var.ChanQKinAvgDt.copy() +ChanQAvgDt_old = self.var.ChanQAvgDt.copy() +# + additional copies inside repMBTs block +``` + +**Why it matters:** +Each `.copy()` allocates a new array of size `num_channel_pixels` (e.g. ~1.8M for EFAS). With MCT + mass-balance reporting active, there are ~10-12 full-domain copies per sub-step. With 4 sub-steps, that is ~40-48 array allocations per model timestep — pure memory allocation + memcpy overhead. + +**Recommendation:** +1. Pre-allocate scratch buffers once in `routing.initial()` and reuse them: +```python +# In initial(): +self._buf_ChanQ_0 = np.empty_like(self.var.ChanQ) +self._buf_ChanM3_0 = np.empty_like(self.var.ChanM3) + +# In dynamic(): +np.copyto(self._buf_ChanQ_0, self.var.ChanQ) # reuse buffer, no allocation +``` + +2. For `SideflowChanM3`: it starts from `ToChanM3RunoffDt` and then has values added/subtracted. Use an in-place pattern: +```python +SideflowChanM3 = self._buf_sideflow +np.copyto(SideflowChanM3, self.var.ToChanM3RunoffDt) +if option['openwaterevapo']: + SideflowChanM3 -= self.var.EvaAddM3Dt +# ... +``` + +3. For the `repMBTs` block: many `.copy()` calls exist solely to avoid modifying the original. Use indexing with `np.where` or masked assignment instead. + +**Estimated saving:** 3-8% of total routing time (reduces GC pressure and memcpy for large domains). + + +--- + +### A5. Numba JIT Cold-Start Compilation Overhead + +| Attribute | Detail | +|-----------|--------| +| **Files** | `mct.py`, `kinematic_wave_parallel_tools.py`, `soilloop.py` | +| **Functions** | All `@njit` decorated functions | +| **Impact** | **MEDIUM** (one-time cost, significant for short runs / calibration) | + +**Why it matters:** +Although `cache=True` is set on all Numba-compiled functions, the cache is invalidated whenever: +- The source file changes (even a comment) +- Numba or NumPy is upgraded +- The function signature changes due to different input dtypes + +During calibration workflows (1000s of short runs), the first run in each new environment pays 10-30 seconds of JIT compilation. For the MCT module alone, there are 6 `@njit` functions. + +**Recommendations:** + +1. **Ahead-of-Time (AOT) compilation**: Use `numba.pycc` to pre-compile the performance-critical kernels into a shared library. This eliminates JIT overhead entirely: +```python +from numba.pycc import CC +cc = CC('lisflood_routing_compiled') + +@cc.export('mct_routing', '...') +def mct_routing(...): ... + +cc.compile() +``` + +2. **Warm-up script**: Provide a lightweight `warmup_numba.py` that imports all JIT functions and calls them once with tiny dummy arrays. Run this as part of container/environment setup. + +3. **Pin Numba + NumPy versions** in Docker/conda environments to avoid cache invalidation across runs. + +**Estimated saving:** 10-30 seconds per cold start; near-zero for warm cache. + + +--- + +### A6. Surface Routing: Three Separate Kinematic Wave Calls + +| Attribute | Detail | +|-----------|--------| +| **File** | `src/lisflood/hydrological_modules/surface_routing.py` | +| **Function** | `surface_routing.dynamic()` | +| **Impact** | **MEDIUM** | + +**Current code:** +```python +self.direct_surface_router.kinematicWaveRouting(self.var.OFQDirectAvg, self.var.OFQDirect, SideflowDirect) +self.other_surface_router.kinematicWaveRouting(self.var.OFQOtherAvg, self.var.OFQOther, SideflowOther) +self.forest_surface_router.kinematicWaveRouting(self.var.OFQForestAvg, self.var.OFQForest, SideflowForest) +``` + +**Why it matters:** +Three independent kinematic wave routing calls are made sequentially, each with its own Numba kernel invocation. Each call processes the full overland-flow pixel domain. Since these three calls are independent (different land-use fractions with no inter-dependency), they could be batched or run concurrently. + +**Recommendations:** + +1. **Batch into a single kernel call**: Modify the kinematic routing kernel to accept a "batch" dimension (3 land-use types). Process all three in a single Numba `prange` call, tripling the available parallelism per order: +```python +# Single call handles all 3 land-use fractions +kinematicRoutingBatched(discharge_avg_batch, discharge_batch, lateral_inflow_batch, ...) +``` + +2. **Alternatively, use Python threading**: Since the three calls release the GIL (Numba `nogil=True` is implicit in `parallel=True`), they can be dispatched to a `ThreadPoolExecutor` with 3 workers. Each call then runs on its own Numba thread pool subset. + +**Estimated saving:** 10-25% of surface routing time (better core utilisation). + + +--- + +### A7. Soil Loop: Numba Thread Contention with Numexpr + +| Attribute | Detail | +|-----------|--------| +| **File** | `src/lisflood/hydrological_modules/soilloop.py` | +| **Functions** | `soilColumnsWaterBalance`, `interception_water_balance` | +| **Impact** | **LOW-MEDIUM** | + +**Why it matters:** +The soil module uses both `numexpr` (via `nx.evaluate`) and Numba `@njit(parallel=True)` with `prange`. Both libraries spawn their own thread pools: +- Numba uses the TBB or OpenMP backend +- numexpr uses its own thread pool (default: number of cores) + +When both are active in the same process, they compete for CPU cores. This over-subscription can cause context-switching overhead, especially on HPC nodes with many cores. + +**Recommendations:** + +1. **Limit numexpr threads**: Set `numexpr.set_num_threads(1)` when Numba parallel functions are the primary workload. Or coordinate thread counts: `NUMEXPR_MAX_THREADS=4` and `NUMBA_NUM_THREADS=N-4`. + +2. **Replace remaining numexpr calls with Numba**: The soil module already has extensive Numba coverage. The few remaining `nx.evaluate` calls (in `kinematic_wave_parallel.py`) can be absorbed into Numba kernels, eliminating the second thread pool entirely. + +3. **Set `NUMBA_THREADING_LAYER=tbb`** explicitly in the environment to ensure deterministic thread management. + +**Estimated saving:** 2-5% overall on many-core systems (16+ cores). + + +--- + +### A8. I/O: Synchronous Chunk Loading in XarrayChunked + +| Attribute | Detail | +|-----------|--------| +| **File** | `src/lisflood/global_modules/netcdf.py` | +| **Function** | `XarrayChunked.load_next_chunk()` and `XarrayChunked.__getitem__()` | +| **Impact** | **LOW-MEDIUM** | + +**Current code:** +```python +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() # blocks until data is in memory +``` + +**Why it matters:** +When a new temporal chunk boundary is crossed, the model blocks while xarray loads the next chunk from disk. For 5 forcing variables (Precip, Tavg, ET0, ES0, E0), this happens synchronously in `readmeteo.dynamic()`. Each chunk load involves NetCDF decompression (zlib) which is CPU-bound. + +**Recommendations:** + +1. **Prefetch next chunk asynchronously**: Use a background thread to load the next chunk before it is needed: +```python +import threading + +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() + +def prefetch_next_chunk(self): + if self.ichunk + 1 < len(self.chunk_indexes) - 1: + self._prefetch_thread = threading.Thread(target=self._prefetch) + self._prefetch_thread.start() +``` + +2. **Increase chunk size**: The setting `NetCDFTimeChunks` controls chunk size. For daily forcing over a year, set to 365 (load entire year at once). Memory cost is modest: 5 variables * 1.8M pixels * 365 days * 4 bytes = ~13 GB for float32 EFAS forcing loaded fully in memory. + +3. **Use `MapsCaching=True`** for forcing data when memory permits — loads the entire time series at initialization. + +**Estimated saving:** 1-5% for chunked reads; up to 10% if chunk boundaries align with expensive computation. + + +--- + +### A9. Output Writing: Per-Timestep NetCDF Writes + +| Attribute | Detail | +|-----------|--------| +| **File** | `src/lisflood/global_modules/output.py` | +| **Functions** | `NetcdfStepsWriter.stage()`, `NetcdfStepsWriter.write()` | +| **Impact** | **LOW-MEDIUM** | + +**Current behaviour:** +The `OutputMapsChunks` setting (default 1) controls how many timesteps are buffered before writing. With default=1, every reporting timestep triggers: +1. `uncompress_array()` — reconstruct 2D map from compressed 1D array +2. File open (or keep-open via `iterOpenNetcdf`) +3. Single-slice write to NetCDF variable +4. File close (if chunk boundary) + +**Why it matters:** +For runs reporting every timestep (e.g. hourly discharge maps for a year = 8760 writes), the overhead of repeated file I/O becomes significant. Each write involves: +- Python-level overhead of netCDF4 library calls +- OS-level file metadata updates +- zlib compression of each 2D slice + +**Recommendations:** + +1. **Increase `OutputMapsChunks`**: Set to 30-365 depending on output frequency. This buffers multiple timesteps in memory and writes them in a single batch, improving compression ratio and reducing file I/O overhead. + +2. **Use the existing `OutputMapsFactoryThreads` class**: The codebase already contains an async-write implementation using `ThreadPool`. It is documented as "NOT FULLY TESTED" but the approach is sound. Validate and enable it as an option: +```xml + +``` + +3. **Batch `uncompress_array` calls**: When writing multiple output variables at the same timestep, the decompression mask operation is repeated for each variable. Cache the mask indexing. + +**Estimated saving:** 5-15% for I/O-heavy configurations (many output maps, frequent reporting). + + +--- + +### A10. Dynamic Loop: Python-Level Overhead Per Timestep + +| Attribute | Detail | +|-----------|--------| +| **File** | `src/lisflood/Lisflood_dynamic.py` | +| **Function** | `LisfloodModel_dyn.dynamic()` | +| **Impact** | **LOW** | + +**Why it matters:** +Each model timestep involves: +- `datetime.timedelta` computation and `strftime` formatting +- Multiple `if option[...]` checks (dict lookups) +- `sys.stdout.write` and `sys.stdout.flush` for progress reporting + +For sub-hourly simulations (e.g. 15-minute timesteps over multiple years), there are 100K+ timesteps. Python-level overhead accumulates. + +**Recommendations:** + +1. **Cache option flags as local booleans** at the start of `dynamic()`: +```python +opt_MCT = option['MCTRouting'] +opt_split = option['SplitRouting'] +opt_repMBTs = option['repMBTs'] +# use local variables in if-statements +``` + +2. **Reduce progress output frequency**: Print progress every N steps rather than every step: +```python +if not flags['quiet'] and i % 100 == 0: + sys.stdout.write(...) +``` + +3. **Avoid `uuid.uuid4()` on first step**: The CDFFlags initialization creates a UUID on timestep 1. This is a system call that can be slow on some platforms. Move to `initial()`. + +**Estimated saving:** 1-2% for very long simulations with many timesteps. + + +--- + +### A11. MCT Routing: Calibration Points Check Inside Hot Loop + +| Attribute | Detail | +|-----------|--------| +| **File** | `src/lisflood/hydrological_modules/mct.py` | +| **Function** | `mct_routing()`, inner loop | +| **Impact** | **LOW-MEDIUM** | + +**Current code:** +```python +for ups_ix in range(num_upstream_pixels[kinpix]): + ups_pix = upstream_pixels[ups_ix] + if np.any(CalibPointsIds == ups_pix): + ql += ChanQAvgDt[ups_pix] + else: + q00 += ChanQ_0[ups_pix] + q0m += ChanQAvgDt[ups_pix] + q01 += ChanQ[ups_pix] +``` + +**Why it matters:** +`np.any(CalibPointsIds == ups_pix)` performs a linear scan of the `CalibPointsIds` array for every upstream pixel of every MCT pixel, on every routing sub-step. If there are K calibration points and U upstream connections per pixel, this is O(K * U * num_MCT_pixels * num_substeps) comparisons. + +For EFAS with ~50 calibration points and 5000 MCT pixels, this check dominates the inner loop when K is non-trivial. + +**Recommendation:** +Replace the linear scan with a pre-computed boolean lookup array: +```python +# In MCTWave.__init__: +is_calib_point = np.zeros(num_all_pixels, dtype=np.bool_) +is_calib_point[CalibPointsIds] = True + +# In mct_routing inner loop: +if is_calib_point[ups_pix]: + ... +``` + +This changes the check from O(K) to O(1) per upstream pixel. + +**Estimated saving:** 5-15% of MCT kernel time when calibration points are active. + + +--- + +### A12. Water Balance Module: Repeated `np.bincount` + `np.take` Pattern + +| Attribute | Detail | +|-----------|--------| +| **File** | `src/lisflood/hydrological_modules/waterbalance.py` | +| **Function** | `waterbalance.dynamic()` | +| **Impact** | **LOW** | + +**Current code pattern (repeated ~10 times):** +```python +WaterIn += np.take(np.bincount(self.var.Catchments, weights=some_array), self.var.Catchments) +``` + +**Why it matters:** +Each `np.bincount(..., weights=...)` creates a temporary array of size `max(Catchments)+1`, then `np.take` expands it back to full domain size. With 10+ such calls per timestep, this creates substantial temporary memory allocation and cache pressure. + +**Recommendations:** + +1. **Accumulate weights first, then do a single bincount**: Instead of calling bincount separately for each variable, sum the weighted arrays first: +```python +total_weights = self.var.TotalPrecipitationWB * self.var.MMtoM3 +if option['inflow']: + total_weights += self.var.sumInWB +WaterIn = np.take(np.bincount(self.var.Catchments, weights=total_weights), self.var.Catchments) +``` + +2. **Pre-allocate the bincount result buffer**: Reuse the same output array across calls. + +3. **Consider skipping mass balance for production runs**: The `repMBTs` option enables extensive per-catchment accounting. For operational forecasting (where mass balance is not the focus), disabling this saves all associated computation. + +**Estimated saving:** 1-3% of total timestep time. + + +--- + +## Part B — Storage & Memory Optimization + +--- + +### B1. Output Data Type: Default float64 + +| Attribute | Detail | +|-----------|--------| +| **File** | `src/lisfloodSettings_reference.xml` (line ~206) | +| **Setting** | `OutputMapsDataType` | +| **Impact** | **HIGH** (50% output storage reduction) | + +**Current default:** +```xml + +``` + +**Why it matters:** +All output NetCDF maps are written as float64 (8 bytes per value). For typical hydrological variables (discharge 0-100000 m³/s, soil moisture 0-1, temperature -50 to 50°C), float32 provides ~7 significant digits — more than sufficient. + +For an EFAS daily discharge output: 1.8M pixels * 8 bytes * 365 days = ~4.8 GB/year in float64 vs 2.4 GB/year in float32. + +**Recommendation:** +Change the default to `float32`: +```xml + +``` + +**Risk:** Very low. Float32 precision (7 significant digits) exceeds the physical accuracy of any hydrological variable LISFLOOD produces. No model results are affected — this is output-only. + +**Estimated saving:** 50% reduction in all output map file sizes. + + +--- + +### B2. NetCDF Time Chunk Size + +| Attribute | Detail | +|-----------|--------| +| **File** | `src/lisflood/global_modules/netcdf.py` (line ~572) | +| **Function** | `write_netcdf_header()` | +| **Impact** | **MEDIUM** (10-30% better compression) | + +**Current code:** +```python +value = nf1.createVariable(var_name, dtype, ('time', dim_lat_y, dim_lon_x), + zlib=True, fill_value=-9999, chunksizes=(1, nrow, ncol)) +``` + +**Why it matters:** +A time-chunk of 1 means each timestep is compressed independently. Hydrological fields exhibit strong temporal autocorrelation — consecutive timesteps are similar. Compressing multiple timesteps together exploits this redundancy for better compression ratios. + +Additionally, reading back data for post-processing (e.g. extracting a time series at one pixel) requires decompressing full spatial slices one at a time when time-chunk=1. + +**Recommendation:** +Use a time chunk that matches typical access patterns: +```python +# For daily output files: chunk 30 days together +time_chunk = min(30, n_timesteps) if n_timesteps > 1 else 1 +chunksizes = (time_chunk, nrow, ncol) +``` + +For hourly outputs (sub-daily), consider `time_chunk = min(24, n_timesteps)` (one day of hours). + +Make this configurable via a new setting: +```xml + +``` + +**Risk:** Low — chunking is internal to NetCDF format; does not affect data values. + +**Estimated saving:** 10-30% better compression ratio; faster time-series extraction in post-processing. + + +--- + +### B3. Compression Level Setting + +| Attribute | Detail | +|-----------|--------| +| **File** | `src/lisflood/global_modules/netcdf.py` | +| **Function** | `write_netcdf_header()` | +| **Impact** | **LOW-MEDIUM** | + +**Current code:** +```python +zlib=True # uses default complevel (typically 4 in netCDF4 library) +``` + +**Why it matters:** +The `complevel` parameter is not explicitly set, defaulting to whatever the netCDF4 library uses (usually 4). For operational runs where output I/O time is a bottleneck, a lower compression level (1-2) can significantly speed up writes at the cost of slightly larger files. For archival runs, a higher level (6-9) produces smaller files at the cost of write speed. + +**Recommendation:** +Add an explicit compression level setting: +```python +comp_level = int(binding.get('OutputCompressionLevel', 4)) +value = nf1.createVariable(..., zlib=True, complevel=comp_level, ...) +``` + +```xml + +``` + +**Trade-off guidance:** +| Level | Write speed | File size | Use case | +|-------|-------------|-----------|----------| +| 1 | Fast | Larger (+20%) | Real-time forecasting | +| 4 | Balanced | Baseline | General use | +| 6-9 | Slow | Smaller (-10-20%) | Long-term archival | + +**Estimated saving:** Configurable — up to 20% smaller files (level 6+) or 30% faster writes (level 1). + + +--- + +### B4. `compressArray` Forces float64 on All Arrays + +| Attribute | Detail | +|-----------|--------| +| **File** | `src/lisflood/global_modules/add1.py` (line ~291) | +| **Function** | `compressArray()` | +| **Impact** | **MEDIUM** (memory footprint) | + +**Current code:** +```python +return mapC.astype(float) # float = float64 in NumPy +``` + +**Why it matters:** +Every map loaded through `compressArray` is cast to float64 regardless of its original dtype. This includes: +- Boolean masks (could be bool/int8 → 1 byte becomes 8 bytes) +- Integer maps like LDD, land-use classes (int8/int16 → float64) +- Maps read as float32 from NetCDF files (4 bytes → 8 bytes) + +For EFAS with ~1.8M pixels, each unnecessary float64 promotion costs 14.4 MB (from float32) or 12.6 MB (from int8). With dozens of static maps loaded at initialization, this adds up to 200-500 MB of wasted memory. + +**Recommendation:** +Preserve the input dtype or allow explicit specification: +```python +def compressArray(map, pcr=True, name=None, force_load_with_nans=False, dtype=None): + # ... existing logic ... + if dtype is not None: + return mapC.astype(dtype) + elif hasattr(mapC, 'dtype') and not np.issubdtype(mapC.dtype, np.floating): + # Keep integer/boolean maps in their native type + return mapC + else: + return mapC.astype(np.float64) # backward compatible default for float maps +``` + +For a less invasive change, at minimum avoid promoting boolean and integer maps: +```python +if np.issubdtype(mapC.dtype, np.integer) or np.issubdtype(mapC.dtype, np.bool_): + return mapC +return mapC.astype(float) +``` + +**Risk:** Low if done carefully. Some downstream code may assume float64; a `grep` for `.astype(float)` patterns in dependent code should be checked. + +**Estimated saving:** 200-500 MB RAM reduction for large EFAS/GloFAS domains. + + +--- + +### B5. MCT State Variables: Extra Arrays Per Pixel + +| Attribute | Detail | +|-----------|--------| +| **File** | `src/lisflood/hydrological_modules/routing.py` | +| **Variables** | `PrevCm0`, `PrevDm0`, `ChanQ_0`, `ChanM3_0` (MCT-specific state) | +| **Impact** | **LOW** | + +**Why it matters:** +MCT routing introduces 4 additional full-domain float64 arrays as state variables: +- `PrevCm0` (Courant number) — only meaningful at MCT pixels +- `PrevDm0` (Reynolds number) — only meaningful at MCT pixels +- `ChanQ_0` (previous discharge copy) — full domain +- `ChanM3_0` (previous storage copy) — full domain + +For EFAS (~1.8M pixels): 4 arrays * 1.8M * 8 bytes = ~57 MB. However, only ~5000 pixels are MCT pixels (0.3% of domain). + +**Recommendation:** +Store `PrevCm0` and `PrevDm0` as MCT-only compressed arrays (5000 pixels) rather than full-domain arrays: +```python +# In MCTWave.__init__: +self.PrevCm0_mct = np.ones(num_mct_pixels) +self.PrevDm0_mct = np.zeros(num_mct_pixels) +``` + +Map between MCT-local and full-domain indices using the existing `mapping_mct` array. This reduces memory by ~25 MB for these two arrays alone. + +For `ChanQ_0` and `ChanM3_0`: these are full-domain copies needed as "previous state" inputs. They could be eliminated by reorganizing the MCT call to read directly from `self.var.ChanQ` before kinematic routing overwrites it (restructure call order). + +**Estimated saving:** 25-57 MB RAM. + + +--- + +### B6. Enable MapsCaching by Default for Static Maps + +| Attribute | Detail | +|-----------|--------| +| **File** | `src/lisfloodSettings_reference.xml` | +| **Setting** | `MapsCaching` | +| **Impact** | **MEDIUM** (runtime), **LOW** (memory trade-off) | + +**Current default:** +```xml + +``` + +**Why it matters:** +When `MapsCaching` is `False`, the `XarrayChunked` reader is used for forcing data. When `True`, the `XarrayCached` class (which uses `@Cache` decorator) stores the full dataset in memory, avoiding repeated disk reads. + +For static maps loaded via `loadmap()`, caching avoids re-reading the same NetCDF file if it is referenced multiple times (e.g. soil parameters used in both initialization and dynamic sections). + +**Recommendation:** +Enable caching by default: +```xml + +``` + +For memory-constrained environments, document the memory cost and allow users to disable it. + +**Memory cost:** For EFAS forcing data (5 variables, 1.8M pixels, 365 days, float32): ~13 GB. This is acceptable on modern HPC nodes (128-512 GB) but may be prohibitive on smaller machines. + +**Estimated saving:** 10-30% faster initialization; eliminates repeated disk reads for static data during the run. + + +--- + +### B7. Output Variable Pruning: Map vs TSS Trade-offs + +| Attribute | Detail | +|-----------|--------| +| **File** | Settings XML (user configuration) | +| **Impact** | **HIGH** (operational storage) | + +**Why it matters:** +LISFLOOD can output dozens of 2D map variables at every timestep. For EFAS operational runs, the output volume is often the dominant storage cost. Many variables are only needed at specific gauge points (time series) rather than as full 2D maps. + +**Recommendations:** + +1. **Prefer TSS over Maps for point-based diagnostics**: Variables like discharge, water level, and lake/reservoir levels are typically only needed at gauge locations. Use TSS (Time Series) output instead of full maps: + - Storage for 1 variable, 1000 gauges, 1 year daily: ~3 MB (TSS) vs ~2.4 GB (map, float32) + - Factor 800x reduction + +2. **Use monthly/yearly frequency for slow-changing variables**: Soil moisture, groundwater, snow cover change slowly. Report them monthly rather than daily: + - 12x reduction in output volume for these variables + +3. **Document output configurations**: Provide template settings files for common use cases: + - `settings_minimal_output.xml` — discharge TSS only + - `settings_operational.xml` — discharge maps + key state variables + - `settings_research.xml` — full variable set + +4. **Consider lossy compression for diagnostic maps**: For variables like soil moisture or temperature that are only used for visualization, NetCDF's `least_significant_digit` option can dramatically improve compression: +```python +value = nf1.createVariable(..., zlib=True, least_significant_digit=3) +# Keeps 3 significant digits, enables much better compression +``` + +**Estimated saving:** 50-95% output volume reduction depending on configuration. + + +--- + +### B8. Internal Computation Precision: Float64 Everywhere + +| Attribute | Detail | +|-----------|--------| +| **Files** | All hydrological modules | +| **Impact** | **LOW** (memory), **LOW** (runtime on modern CPUs) | + +**Current state:** +All internal state variables and computations use float64. This is the safe default for numerical stability, especially for: +- Accumulation variables over long simulations (mass balance) +- Channel routing where small flows interact with large storages +- The MCT solver which is sensitive to precision in Courant/Reynolds numbers + +**Recommendation:** +Do NOT change internal precision to float32 for state variables. The risk of numerical drift over multi-year simulations outweighs the modest memory savings. Float64 arithmetic is nearly as fast as float32 on modern x86-64 CPUs (same SIMD width for scalar operations; only vectorised code benefits from float32's 2x throughput). + +**Exception:** Intermediate temporary arrays that are not accumulated (e.g. `SideflowChan`, `SideflowDirect`) could safely use float32 if memory pressure is extreme, but the savings (~14 MB per temporary for EFAS) rarely justify the added complexity. + +**Conclusion:** Keep float64 for internal computation. Focus memory optimization on B4 (dtype preservation for non-float maps) and B1 (output dtype). + + +--- + +## Part C — Prioritized Action Plan + +| Rank | ID | Recommendation | Runtime Impact | Storage Impact | Complexity | Risk | +|------|----|----------------|----------------|----------------|------------|------| +| 1 | B9 | CF scale/offset packing (int16 output) | ~5% (I/O) | **75-90%** | Medium | Low | +| 2 | B1 | Change default OutputMapsDataType to float32 | ~5% (I/O) | **50%** | Low | None | +| 2 | A11 | MCT: Replace linear CalibPoints scan with boolean lookup | 5-15% MCT | — | Low | None | +| 3 | A2 | MCT: Pre-compute ANalv (arctan) | 2-5% MCT | — | Low | None | +| 4 | A3 | Kinematic wave: Replace numexpr with Numba helper | 5-10% routing | — | Medium | Low | +| 5 | B2 | Optimize NetCDF output time chunk size | 5% (I/O) | **10-30%** | Low | None | +| 6 | A4 | Routing: Pre-allocate scratch buffers, reduce .copy() | 3-8% routing | — | Medium | Low | +| 7 | B4 | compressArray: Preserve input dtype | — | **200-500 MB** | Low | Low | +| 8 | A6 | Surface routing: Batch 3 kinematic calls | 10-25% surface | — | Medium | Low | +| 9 | B7 | Output variable pruning (TSS vs Maps) | 5-15% (I/O) | **50-95%** | Config only | None | +| 10 | A1 | MCT: Better initial guess / vectorise batches | 10-20% MCT | — | High | Medium | +| 11 | B3 | Add configurable compression level | ±10% (I/O) | **±20%** | Low | None | +| 12 | A5 | Numba AOT / warm-up script | 10-30s startup | — | Medium | Low | +| 13 | B6 | Enable MapsCaching by default | 10-30% init | +memory | Low | Low | +| 14 | A8 | Async chunk prefetch for forcing I/O | 1-5% | — | Medium | Low | +| 15 | A9 | Increase OutputMapsChunks / enable async writes | 5-15% (I/O) | — | Low | Low | +| 16 | A7 | Resolve Numba/numexpr thread contention | 2-5% | — | Low | None | +| 17 | A12 | Water balance: reduce bincount calls | 1-3% | — | Low | None | +| 18 | A10 | Dynamic loop: cache options, reduce stdout | 1-2% | — | Low | None | +| 19 | B5 | MCT: Compress PrevCm0/PrevDm0 to MCT-only arrays | — | **25-57 MB** | Medium | Low | + + +--- + +## Part D — Quick Wins (Implementable in < 1 day each) + +### D1. Change default output dtype (B1) +**File:** `src/lisfloodSettings_reference.xml` +```diff +- ++ +``` + +### D2. Pre-compute ANalv in MCT (A2) +**File:** `src/lisflood/hydrological_modules/mct.py` + +In `MCTWave.__init__()`, add: +```python +self.ANalv = np.arctan(1.0 / ChanSdXdY) +``` + +Add `ANalv` as a parameter to `mct_routing()` and in the inner loop replace: +```python +ANalv = np.arctan(1 / ChanSdXdY[kinpix]) +``` +with: +```python +ANalv = precomputed_ANalv[kinpix] +``` + +### D3. Replace CalibPoints linear scan with boolean array (A11) +**File:** `src/lisflood/hydrological_modules/mct.py` + +In `MCTWave.__init__()` or `routing.initialMCT()`: +```python +num_pixels = len(ChanLength) +is_calib_point = np.zeros(num_pixels, dtype=np.bool_) +if len(CalibPointsIds) > 0: + is_calib_point[CalibPointsIds] = True +``` + +Pass `is_calib_point` to `mct_routing` and replace: +```python +if np.any(CalibPointsIds == ups_pix): +``` +with: +```python +if is_calib_point[ups_pix]: +``` + +### D4. Add compression level setting (B3) +**File:** `src/lisflood/global_modules/netcdf.py`, in `write_netcdf_header()`: +```python +comp_level = int(binding.get('OutputCompressionLevel', 4)) +value = nf1.createVariable(var_name, dtype, dims, zlib=True, complevel=comp_level, fill_value=-9999, chunksizes=chunks) +``` + +### D5. Increase default OutputMapsChunks (A9) +**File:** `src/lisfloodSettings_reference.xml` +```diff +- ++ +``` + + +--- + +## Part E — Changes from Previous Report + +| Previous Recommendation | Status | Notes | +|------------------------|--------|-------| +| Change OutputMapsDataType to float32 | **Still valid** | Retained as #1 priority | +| Optimize NetCDF time chunk sizes | **Still valid** | Code unchanged; chunksizes=(1,nrow,ncol) persists | +| Add compression level setting | **Still valid** | Not yet implemented | +| Enable MapsCaching by default | **Still valid** | Still defaults to False | +| Preserve input dtype in compressArray | **Still valid** | `mapC.astype(float)` unchanged | +| Add async I/O for output writing | **Partially implemented** | `OutputMapsFactoryThreads` exists but is not enabled by default | +| Consider float32 for internal calculations | **Superseded** | Now explicitly recommended AGAINST (see B8) — risk outweighs benefit | + +### New sections added (not in original report): + +| Section | Reason | +|---------|--------| +| A1, A2, A11 | MCT routing module did not exist | +| A3 | numexpr usage in kinematic_wave_parallel.py is new | +| A4 | .copy() pattern in routing loop amplified by MCT additions | +| A5 | Numba cold-start now affects more functions (MCT adds 6 new @njit) | +| A6 | Surface routing now uses parallel kinematic wave (was PCRaster before) | +| A7 | Thread contention newly relevant with parallel soil + parallel routing | +| A10 | Dynamic loop analysis refreshed against current code | +| A12 | Water balance module rewritten with new bincount patterns | +| B5 | MCT-specific state variable overhead is new | +| B7 | Output pruning guidance (operational vs research) | +| B8 | Explicit recommendation against float32 internals | + + +--- + +## Appendix: Profiling Guidance + +To validate the estimates in this report, the following profiling approach is recommended: + +### Wall-clock timing per module +Add timing instrumentation in `Lisflood_dynamic.py`: +```python +import time +t0 = time.perf_counter() +self.readmeteo_module.dynamic() +t_meteo = time.perf_counter() - t0 + +t0 = time.perf_counter() +self.soilloop_module.dynamic_soil() +t_soil = time.perf_counter() - t0 + +# ... etc for each module call +``` + +### Numba kernel timing +Use `numba.core.config.DEVELOPER_MODE = 1` and Numba's built-in timing, or wrap calls: +```python +t0 = time.perf_counter() +kwpt.kinematicRouting(...) +t_kin = time.perf_counter() - t0 +``` + +### Memory profiling +Use `tracemalloc` or `memory_profiler` to identify peak memory usage: +```python +import tracemalloc +tracemalloc.start() +# ... run model ... +snapshot = tracemalloc.take_snapshot() +top_stats = snapshot.statistics('lineno') +``` + +### I/O profiling +Monitor disk I/O with system tools (`iostat`, `iotop`) or Python's `cProfile` focused on netCDF4 calls. + +--- + +*Report generated: July 2026* +*Codebase version: current HEAD (post-MCT integration)* +*Previous report: `lisflood_optimization_report.md` (I/O-focused, pre-MCT)* + +--- + +### B9. CF-Convention Scale/Offset Packing (int16 storage of float data) + +| Attribute | Detail | +|-----------|--------| +| **File** | `src/lisflood/global_modules/netcdf.py` | +| **Function** | `write_netcdf_header()` | +| **Impact** | **HIGH** (75% output storage reduction vs float64; 50% vs float32) | + +**Current state:** +Output variables are stored as raw floating-point values (float64 or float32). No packing is applied. + +**What is scale/offset packing?** +The CF conventions define two variable attributes — `scale_factor` and `add_offset` — that allow floating-point data to be stored as smaller integer types (typically int16 or uint16). On read, the client library automatically reconstructs the original value: + +``` +unpacked_value = packed_value * scale_factor + add_offset +``` + +A 16-bit integer provides 65536 distinct values. By choosing `scale_factor` and `add_offset` to cover the physical range of the variable, you get fixed-precision storage at **2 bytes per value** instead of 4 (float32) or 8 (float64). + +**Why it matters:** +- float64 → int16: **75% reduction** (8 bytes → 2 bytes per value) +- float32 → int16: **50% reduction** (4 bytes → 2 bytes per value) +- Combined with zlib compression, the actual on-disk savings can exceed 80% because integer data with limited range compresses far better than floating-point data. +- This is standard practice in climate/meteorological data (CMIP6, ERA5, EFAS forcing data already use this). + +**Variable-specific packing parameters:** + +| Variable | Physical Range | scale_factor | add_offset | Precision | Suitable? | +|----------|---------------|--------------|------------|-----------|-----------| +| Discharge (m³/s) | 0 – 100,000 | 1.53 | 50,000 | ±0.8 m³/s | Yes (large rivers) | +| Discharge (m³/s) | 0 – 10,000 | 0.153 | 5,000 | ±0.08 m³/s | Yes (medium rivers) | +| Soil Moisture (fraction) | 0 – 1 | 1.53e-5 | 0.5 | ±8e-6 | Yes | +| Snow Water Equiv. (mm) | 0 – 2000 | 0.031 | 1000 | ±0.015 mm | Yes | +| Temperature (°C) | -50 – +50 | 0.00153 | 0 | ±0.001 °C | Yes | +| ET (mm/day) | 0 – 20 | 3.05e-4 | 10 | ±1.5e-4 mm | Yes | +| Water Level (m) | -5 – +20 | 3.82e-4 | 7.5 | ±0.2 mm | Yes | +| Groundwater (mm) | 0 – 5000 | 0.076 | 2500 | ±0.04 mm | Yes | + +**Implementation:** + +In `write_netcdf_header()`, add packing when a packing configuration is provided: + +```python +# Determine packing based on variable metadata or user settings +packing = binding.get('OutputPacking', 'none') # 'none', 'int16', 'auto' + +if packing == 'int16' and frequency is not None: + # Use variable-specific range or a safe default + vmin, vmax = get_variable_range(var_name) # lookup table + scale = (vmax - vmin) / 65534.0 # leave room for fill_value + offset = vmin + scale * 32767.0 + + value = nf1.createVariable(var_name, 'i2', ('time', dim_lat_y, dim_lon_x), + zlib=True, fill_value=-32767, + chunksizes=(time_chunk, nrow, ncol)) + value.scale_factor = scale + value.add_offset = offset +else: + # Current behaviour (float32 or float64) + value = nf1.createVariable(var_name, dtype, ('time', dim_lat_y, dim_lon_x), + zlib=True, fill_value=-9999, + chunksizes=(1, nrow, ncol)) +``` + +The write step must also pack the data before storing: +```python +# In NetcdfStepsWriter.write() or equivalent: +if hasattr(nf1.variables[self.map_name], 'scale_factor'): + scale = nf1.variables[self.map_name].scale_factor + offset = nf1.variables[self.map_name].add_offset + packed = np.round((map_np - offset) / scale).astype(np.int16) + packed[map_np == -9999] = -32767 # fill value + nf1.variables[self.map_name][step, :, :] = packed +else: + nf1.variables[self.map_name][step, :, :] = map_np +``` + +**Settings XML:** +```xml + + +``` + +**Compatibility:** +- All modern NetCDF readers (xarray, CDO, NCO, QGIS, Python netCDF4) automatically apply `scale_factor` and `add_offset` on read — no user action needed. +- CF-compliant: follows CF-1.6+ conventions exactly. +- Backward compatible: if `OutputPacking=none`, behaviour is unchanged. + +**Risks:** +- **Precision loss**: int16 provides ~4.8 significant digits (vs ~7 for float32, ~15 for float64). For most hydrological variables this is adequate. For variables with very large dynamic range (e.g. discharge spanning 0.001 to 100,000 m³/s), consider using per-timestep adaptive scaling or splitting into sub-ranges. +- **State/restart files should NOT be packed**: Warm-start state variables (used to restart the model) must retain full float64 precision to avoid drift. Only apply packing to reporting/output maps. + +**Estimated saving:** +- vs current float64 default: **75% output file size reduction** +- vs float32: **50% additional reduction on top of B1** +- Combined with zlib (already enabled): effective on-disk ratios of 85-90% reduction are achievable + From 06dfe39ca8dec59cf8af0ced578510cebf16b124 Mon Sep 17 00:00:00 2001 From: Timo Schaffhauser Date: Fri, 28 Aug 2026 15:08:15 +0200 Subject: [PATCH 2/7] 2 typos --- docs/4_Static-Maps_reservoirs-lakes/index.md | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/docs/4_Static-Maps_reservoirs-lakes/index.md b/docs/4_Static-Maps_reservoirs-lakes/index.md index bdde1f94..483d8df6 100644 --- a/docs/4_Static-Maps_reservoirs-lakes/index.md +++ b/docs/4_Static-Maps_reservoirs-lakes/index.md @@ -54,7 +54,7 @@ The reservoirs map shows the outflow location of each reservoir: each outflow po |Reservoir normal outflow| res_normal_outflow.txt;
2 columms: ID VALUE; 1 row for each reservoir|Units: m3/s|Normal outflow| |Reseervoir minimum outflow| res_min_outflow.txt;
2 columms: ID VALUE; 1 row for each reservoir|Units: m3/s|Minimum outflow| -The well-known Global Reservoir and Dam Database[GDW](https://www.globaldamwatch.org/grand) now superseeded by the Global Dam Watch [GDW](https://www.globaldamwatch.org/database) is a relevant example of source of data for lakes map and tables. +The well-known Global Reservoir and Dam Database [GDW](https://www.globaldamwatch.org/grand) now superseeded by the Global Dam Watch [GDW](https://www.globaldamwatch.org/database) is a relevant example of source of data for lakes map and tables. ### Methodology As a first step, it is recommended to create a file including all reservoir information required by OS LISFLOOD and some relevant metadata that can help with model analysis and results description. @@ -137,7 +137,7 @@ The following paragraphs provide guidelines for the generation of the lake map a Lake unique identifier (1) and coordinates of the outlet mapped on the OS LISFLOOD local drainage direction map ([ldd](../4_Static-Maps_topography/index.md)) (3) are required to generate the lake map. Geographic coordinates of the lake outlet (2) and OS LISFLOOD local drainage direction map ([ldd](../4_Static-Maps_topography/index.md)) are essential to generate (3). Adequate model representation requires the agreement between lake catchment area (7) and OS LISFLOOD [upstream area map](../4_Static-Maps_topography/index.md). -Lake surface area can be retrieved from local datasets or global datasets such as HydroLAKES](https://www.hydrosheds.org/products/hydrolakes), [GLWD](https://www.hydrosheds.org/products/glwd), [GRAND](https://www.globaldamwatch.org/grand). +Lake surface area can be retrieved from local datasets or global datasets such as [HydroLAKES](https://www.hydrosheds.org/products/hydrolakes), [GLWD](https://www.hydrosheds.org/products/glwd), [GRAND](https://www.globaldamwatch.org/grand). Where lake outlet width cannot be retrieved from external dataset, it can be measured with GIS tools. From b54fcb61fbdbe52bab2b6bf3c0027432b6916832 Mon Sep 17 00:00:00 2001 From: Timo Schaffhauser Date: Mon, 7 Sep 2026 17:44:06 +0200 Subject: [PATCH 3/7] Update scale offset documentation --- docs/3_step6_model-output/index.md | 92 ++++++++++++++++++++++++++++++ 1 file changed, 92 insertions(+) diff --git a/docs/3_step6_model-output/index.md b/docs/3_step6_model-output/index.md index 08e6982d..15faf7cd 100644 --- a/docs/3_step6_model-output/index.md +++ b/docs/3_step6_model-output/index.md @@ -44,3 +44,95 @@ Time series file have *.tss* extension and they can be opened with any text edit It is important to be aware of the spatial domain for which each time series is computed. For instance, all *rate variables* are reported as pixel-average values. Soil moisture and groundwater storage are reported for the permeable fraction of each pixel only. The reported snow cover is the average of the snow depths in snow zones A, B and C. This [**Annex**](../5_annex_output-files/index.md) summarises most of the options to report additional output maps. + +## Output data compression (int16 packing) + +By default, output maps are written as floating-point values (`float32` or `float64`, controlled by the `OutputMapsDataType` setting). For large simulations (long time series, global domains, many variables) this can produce very large files. + +LISFLOOD offers an optional lossy compression that stores output maps as signed 16-bit integers (`int16`) using the CF-convention `scale_factor`/`add_offset` packing scheme. This roughly halves the file size compared to `float32` (2 bytes per value instead of 4) and reduces it fourfold compared to `float64`. + +### How to enable it + +Add the following to the settings file: + +```xml + +``` + +- `"False"` (default): output maps are written as raw floating-point values. +- `"True"`: output maps are packed into `int16` using per-variable `scale_factor` and `add_offset`. + +No further user action is required. The `scale_factor` and `add_offset` for each variable are defined internally (see `default_options.py`). Standard netCDF readers (xarray, CDO, NCO, ncview, Panoply) automatically unpack the values on read, so from a user perspective the data appears as normal floating-point. + +### Packing and unpacking formulas + +When writing, each floating-point value is encoded into an integer: + +``` +packed_int = round( (value - add_offset) / scale_factor ) +``` + +When reading, the CF-compliant reader restores the physical value: + +``` +value = packed_int * scale_factor + add_offset +``` + +The packed integer is stored in the range `[-32766, +32767]`. The value `-32767` is reserved as the fill value (`_FillValue`) for masked / outside-domain pixels, so it is never used for real data. + +### Precision and value range + +The `scale_factor` directly determines both the resolution and the representable range of a packed variable: + +- **Resolution (quantization step):** equal to `scale_factor`. A value read back can differ from the original by at most `scale_factor / 2`. +- **Representable range:** + - minimum = `add_offset + scale_factor * (-32766)` + - maximum = `add_offset + scale_factor * (+32767)` + +For a target physical range `[min, max]`, the parameters are chosen as: + +``` +scale_factor = (max - min) / 65533 +add_offset = (min + max) / 2 +``` + +(65533 is the number of usable integer levels, i.e. `32767 - (-32766)`.) + +**Example** — actual soil evaporation (`ESActMaps`) uses `scale_factor = 3.1e-4`, `add_offset = 10.0`: +- minimum = `10.0 + 3.1e-4 * (-32766)` ≈ `0.0 mm/day` +- maximum = `10.0 + 3.1e-4 * (+32767)` ≈ `+20.2 mm/day` +- resolution ≈ `3.1e-4 mm/day` (values accurate to about `±1.5e-4 mm/day`) + +The `add_offset` is placed at the centre of the target range (here 10.0, the midpoint of roughly 0–20 mm/day), so the available integer levels are distributed symmetrically around it. + +Values falling outside the representable range are clipped to the nearest bound, and a warning is issued once per variable during the run. Because packing is lossy, a tiny quantization error is expected: for instance an input of exactly `0.0` may be read back as a value on the order of `1e-4`. This is far below any hydrologically meaningful threshold and is normal for packed datasets (the same behaviour applies to ERA5, CMIP6, and other CF-packed products). + +### Variables that are never packed + +- **State and "end" maps** (used to initialise warm-start runs) are always written at full floating-point precision, regardless of the `OutputPacking` setting. This preserves the exact model state needed for a bit-reproducible restart. +- Variables whose value range is too wide to pack meaningfully (for example **discharge**, which spans several orders of magnitude globally) are intentionally left unpacked to avoid an unacceptable loss of precision. + +## Temporal aggregation of output maps + +Instead of writing a value at every model time step, LISFLOOD can aggregate selected output variables to monthly or yearly means or sums. This dramatically reduces output volume for long simulations while retaining the climatologically relevant signal. + +Four settings control temporal aggregation, each taking a semicolon-separated list of output variable names: + +```xml + + + + +``` + +- **Mean** aggregation writes the average of all time steps within each period. +- **Sum** aggregation writes the accumulated total over each period. +- A value is written only at the end of a completed period (month-end or year-end). Partial periods at the end of a simulation are not written. +- When a variable is listed for aggregation, its normal per-time-step output is suppressed to avoid duplication. + +The aggregation boundaries are detected directly from the calendar date and therefore work for daily as well as sub-daily time steps. + +**Interaction with packing:** + +- **Mean** aggregates stay within the same value range as the daily values, so `int16` packing (if `OutputPacking = True`) is applied as usual. +- **Sum** aggregates can greatly exceed the daily value range (a monthly sum can be up to ~30× a daily value). To prevent overflow, packing is automatically disabled for sum-aggregated outputs, which are always written as floating-point. From e1e423453b33c5948450e7520cb035ce369f72a8 Mon Sep 17 00:00:00 2001 From: Carlo Russo <51730707+doc78@users.noreply.github.com> Date: Tue, 8 Sep 2026 10:26:32 +0200 Subject: [PATCH 4/7] Add temporary markdown files to .gitignore --- .gitignore | 3 +++ 1 file changed, 3 insertions(+) diff --git a/.gitignore b/.gitignore index ef5f354a..252b2998 100644 --- a/.gitignore +++ b/.gitignore @@ -27,3 +27,6 @@ lisflood_model.egg-info .vscode/ *.ipynb .ipynb_checkpoints/ +implementation_guide_scale_offset_packing.md +lisflood_optimization_report.md +lisflood_optimization_report_v2.md From 9d9103cb372b5bc0b065412fab0fa104664fabcc Mon Sep 17 00:00:00 2001 From: Carlo Russo <51730707+doc78@users.noreply.github.com> Date: Tue, 8 Sep 2026 10:28:41 +0200 Subject: [PATCH 5/7] Delete implementation_guide_scale_offset_packing.md --- implementation_guide_scale_offset_packing.md | 354 ------------------- 1 file changed, 354 deletions(-) delete mode 100644 implementation_guide_scale_offset_packing.md diff --git a/implementation_guide_scale_offset_packing.md b/implementation_guide_scale_offset_packing.md deleted file mode 100644 index efb01cfa..00000000 --- a/implementation_guide_scale_offset_packing.md +++ /dev/null @@ -1,354 +0,0 @@ -# Implementation Guide: CF Scale/Offset Packing for LISFLOOD Outputs - -## Overview - -This guide walks through adding int16 scale/offset packing to LISFLOOD output NetCDF files. -The feature stores floating-point output variables as packed int16 values using the CF-convention -`scale_factor` and `add_offset` attributes. All CF-compliant readers (xarray, CDO, NCO, QGIS) -automatically unpack on read — no downstream changes needed. - -**Result:** 75-90% smaller output files compared to float64 defaults. - -**Design principles:** -- Global on/off via a simple `True`/`False` setting (`OutputPacking`) -- Per-variable scale/offset stored directly in the `ReportedMap` namedtuple attributes -- State/end maps (for warm start) are never packed -- Range documentation kept in a separate `.md` reference file - ---- - -## Step 1: Add the setting to the reference settings XML - -**File:** `src/lisfloodSettings_reference.xml` - -Find the block where `OutputMapsDataType` is defined and add after it: - -```xml - - -The option "OutputPacking" enables CF-convention scale_factor/add_offset packing -of output maps into int16 (2 bytes per value instead of 4 or 8). - - "False" (default): write raw floating-point values (dtype from OutputMapsDataType) - - "True": pack into signed 16-bit integers using per-variable scale/offset - Readers automatically unpack using: value = packed * scale_factor + add_offset -Note: State/end maps used for warm starts are NEVER packed (always full precision). - - -``` - -Also add the pass-through in the bindings section (where `OutputMapsDataType` is passed): - -```xml - -``` - - ---- - -## Step 2: Extend the `ReportedMap` namedtuple - -**File:** `src/lisflood/global_modules/default_options.py` - -Change: - -```python -ReportedMap = namedtuple('ReportedMap', 'name, output_var, unit, end, steps, all, restrictoption, monthly, yearly') -``` - -To: - -```python -ReportedMap = namedtuple('ReportedMap', 'name, output_var, unit, end, steps, all, restrictoption, monthly, yearly, scale_factor, add_offset') -ReportedMap.__new__.__defaults__ = (None, None) # scale_factor and add_offset default to None (= no packing) -``` - -Setting `__new__.__defaults__` means all existing `ReportedMap(...)` entries remain valid -without modification — they'll get `scale_factor=None, add_offset=None` automatically. -Only variables you explicitly want to pack need the extra two fields. - ---- - -## Step 3: Add scale/offset values to output variables - -**File:** `src/lisflood/global_modules/default_options.py` - -For each variable you want to pack, add `scale_factor` and `add_offset` at the end of its -`ReportedMap` entry. Use the helper formula: - -``` -scale_factor = (physical_max - physical_min) / 65534 -add_offset = physical_min + scale_factor * 32767 -``` - -Examples: - -```python -'DischargeMaps': ReportedMap(name='DischargeMaps', output_var='ChanQAvg', - unit='m3/s', end=[], steps=[], - all=['repDischargeMaps'], restrictoption=[], - monthly=False, yearly=False, - scale_factor=3.052, # range 0-200000 m3/s - add_offset=100001.5), - -'WaterLevelMaps': ReportedMap(name='WaterLevelMaps', output_var='WaterLevel', - unit='m', end=[], steps=[], - all=['repWaterLevelMaps'], restrictoption=['nonInit'], - monthly=False, yearly=False, - scale_factor=6.104e-4, # range -10 to +30 m - add_offset=10.0), - -'SnowMaps': ReportedMap(name='SnowMaps', output_var='SnowCover', - unit='mm', end=[], steps=[], - all=['repSnowMaps'], restrictoption=['nonInit'], - monthly=False, yearly=False, - scale_factor=0.04578, # range 0-3000 mm - add_offset=1500.0), -``` - -Variables that should NOT be packed (state/end maps) simply omit the fields: - -```python -'ChanQEnd': ReportedMap(name='ChanQEnd', output_var='ChanQ', unit='m3/s', - end=['repEndMaps'], steps=[], all=[], - restrictoption=[], monthly=False, yearly=False), - # No scale_factor/add_offset → defaults to None → never packed -``` - - ---- - -## Step 4: Modify `write_netcdf_header()` to support packing - -**File:** `src/lisflood/global_modules/netcdf.py` - -### 4a. Change the function signature - -Add `map_value=None` as a parameter: - -```python -def write_netcdf_header(settings, var_name, netfile, DtDay, - value_standard_name, value_long_name, value_unit, - start_date, rep_steps, frequency, - map_value=None): # <-- NEW -``` - -### 4b. Replace the variable creation logic - -Find the block at the end of the function where the NetCDF variable is created. -Replace with: - -```python - # Determine if packing is active for this variable - 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 frequency is not None: # output file with "time" dimension - if packing_enabled and has_packing: - # CF scale/offset packing into int16 - value = nf1.createVariable(var_name, 'i2', ('time', dim_lat_y, dim_lon_x), - zlib=True, fill_value=np.int16(-32768), - chunksizes=(1, nrow, ncol)) - value.scale_factor = np.float64(map_value.scale_factor) - value.add_offset = np.float64(map_value.add_offset) - else: - # Standard float output (current behaviour) - value = nf1.createVariable(var_name, dtype, ('time', dim_lat_y, dim_lon_x), - zlib=True, fill_value=-9999, chunksizes=(1, nrow, ncol)) - else: - # End/state maps — NEVER pack (need full precision for warm start) - value = nf1.createVariable(var_name, dtype, (dim_lat_y, dim_lon_x), - zlib=True, fill_value=-9999) -``` - - ---- - -## Step 5: Pass `map_value` through the output call chain - -**File:** `src/lisflood/global_modules/output.py` - -In both `NetcdfWriter.write()` and `NetcdfStepsWriter.write()`, pass `self.map_value` to -`write_netcdf_header`: - -```python -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, - map_value=self.map_value) # <-- ADD THIS -``` - ---- - -## Step 6: Modify the write methods to pack data before storing - -**File:** `src/lisflood/global_modules/output.py` - -### 6a. In `NetcdfStepsWriter.write()`: - -Replace the data-writing loop: - -```python -for step, data in zip(self.step_range, self.data_steps): - map_np = uncompress_array(data) - - nc_var = nf1.variables[self.map_name] - if nc_var.dtype == np.int16: - # Pack float → int16 - scale = nc_var.scale_factor - offset = nc_var.add_offset - packed = np.round((map_np - offset) / scale).astype(np.float64) - packed = np.clip(packed, -32767, 32767) - packed[map_np == -9999] = -32768 # fill value - nc_var.set_auto_maskandscale(False) # CRITICAL: prevent double-packing - nc_var[step, :, :] = packed.astype(np.int16) - else: - nc_var[step, :, :] = map_np -``` - -**Critical:** `set_auto_maskandscale(False)` MUST be called before writing. Without it, the -netCDF4 library applies scale/offset again on write, resulting in double-packing (garbage values). - - ---- - -## Step 7: Testing - -### Quick verification script - -```python -import xarray as xr -import numpy as np - -# Open a packed output file — xarray auto-unpacks -ds = xr.open_dataset('path/to/output.nc') - -print(ds['DischargeMaps'].dtype) # float64 (auto-unpacked) -print(ds['DischargeMaps'].values[:5]) # physically reasonable values -print(ds['DischargeMaps'].encoding) # shows scale_factor, add_offset, dtype=int16 -``` - -### Comparison test - -Run the same simulation with `OutputPacking = False` and `OutputPacking = True`: - -```python -import xarray as xr -import numpy as np - -ref = xr.open_dataset('output_nopacking.nc') -packed = xr.open_dataset('output_packed.nc') - -var = 'DischargeMaps' -diff = np.abs(ref[var].values - packed[var].values) -print(f"Max absolute difference: {diff.max():.4f}") -print(f"Max relative difference: {(diff / (np.abs(ref[var].values) + 1e-10)).max():.6e}") - -# For discharge with range 0-200000: -# scale_factor ≈ 3.05 -# Max error ≈ ±1.5 m3/s (half a quantization step) -``` - -### File size comparison - -```bash -dir output_nopacking.nc output_packed.nc -# Expected: packed file ≈ 20-25% of original size -``` - - ---- - -## Summary of files to modify - -| # | File | What to do | -|---|------|------------| -| 1 | `src/lisfloodSettings_reference.xml` | Add `OutputPacking` textvar (True/False) + pass-through | -| 2 | `src/lisflood/global_modules/default_options.py` | Extend `ReportedMap` namedtuple with `scale_factor, add_offset` | -| 3 | `src/lisflood/global_modules/default_options.py` | Add scale/offset values to specific output variable entries | -| 4 | `src/lisflood/global_modules/netcdf.py` | Add `map_value` param to `write_netcdf_header()`; int16 creation logic | -| 5 | `src/lisflood/global_modules/output.py` | Pass `map_value` to `write_netcdf_header()` | -| 6 | `src/lisflood/global_modules/output.py` | Add packing logic in write methods | -| 7 | `docs/packing_ranges_reference.md` | Create documentation file with range/precision table | - ---- - -## How the flag works - -``` -User sets OutputPacking = "True" in settings XML - │ - ▼ -write_netcdf_header() reads binding['OutputPacking'] - │ - ▼ - ┌─────────────────────────────────┐ - │ packing_enabled = True │ - │ has_packing = map_value has │ - │ scale_factor and add_offset? │ - └──────────┬──────────────────────┘ - │ - ┌──────┴──────┐ - │ │ - has_packing no packing attrs - = True = False (None) - │ │ - ▼ ▼ - Create int16 Create float - variable with variable (normal - scale/offset behaviour) -``` - -Three conditions must ALL be true for packing to happen: -1. `OutputPacking = "True"` in settings (global toggle) -2. The variable's `ReportedMap` has `scale_factor` and `add_offset` defined (per-variable control) -3. The output has a time dimension (`frequency is not None`) — state/end maps are excluded - ---- - -## Gotchas - -1. **Double-packing**: Call `nc_var.set_auto_maskandscale(False)` before writing raw int16 values. Otherwise netCDF4 applies scale/offset again. - -2. **Fill value**: Use `-32768` (int16 minimum). The packed data range uses -32767 to +32767 (65534 levels). NetCDF4/xarray will mask cells with fill value as NaN on read. - -3. **Out-of-range clipping**: If the model produces values outside the defined range, they get clipped. Set ranges generously — better to waste quantization levels than clip real data. - -4. **State files**: The `frequency is None` guard ensures end maps are never packed. Don't add `scale_factor`/`add_offset` to End/State `ReportedMap` entries either (belt and braces). - -5. **Variable name matching**: The `map_value` object is passed directly from the output writer, so no name-to-range matching is needed — the metadata travels with the variable. - ---- - -## Helper: Computing scale/offset from a physical range - -```python -def compute_packing_params(vmin, vmax): - """Compute CF-convention scale_factor and add_offset for int16 packing. - - Parameters - ---------- - vmin : float - Minimum physical value - vmax : float - Maximum physical value - - Returns - ------- - scale_factor : float - add_offset : float - """ - n_levels = 65534.0 # int16 usable range: -32767 to +32767 - scale_factor = (vmax - vmin) / n_levels - add_offset = vmin + scale_factor * 32767.0 - return scale_factor, add_offset -``` - -Example outputs: -``` -compute_packing_params(0, 200000) → (3.052, 100001.5) discharge m3/s -compute_packing_params(0, 1.0) → (1.526e-5, 0.5) fraction -compute_packing_params(-10, 30) → (6.104e-4, 10.0) water level m -compute_packing_params(0, 3000) → (0.04578, 1500.0) snow mm -compute_packing_params(0, 500) → (0.00763, 250.0) precip mm/day -compute_packing_params(0, 5000) → (0.07630, 2500.0) groundwater mm -``` From 59aef519502a546ca81512527035a0108c13f136 Mon Sep 17 00:00:00 2001 From: Carlo Russo <51730707+doc78@users.noreply.github.com> Date: Tue, 8 Sep 2026 10:28:54 +0200 Subject: [PATCH 6/7] Delete lisflood_optimization_report.md --- lisflood_optimization_report.md | 284 -------------------------------- 1 file changed, 284 deletions(-) delete mode 100644 lisflood_optimization_report.md diff --git a/lisflood_optimization_report.md b/lisflood_optimization_report.md deleted file mode 100644 index 32502307..00000000 --- a/lisflood_optimization_report.md +++ /dev/null @@ -1,284 +0,0 @@ -# LISFLOOD I/O Optimization Report - -## Executive Summary - -The LISFLOOD hydrological model codebase has been analyzed for I/O performance and storage optimization opportunities. Key findings: - -1. **Output Data Type**: Default is `float64` (8 bytes) but can be set to `float32` (4 bytes) - potential **50% storage reduction** for output maps -2. **NetCDF Compression**: Already uses `zlib=True` for output files - good practice -3. **Input Reading**: Uses xarray with chunking for NetCDF inputs, but data is loaded into memory as float64 by default -4. **Caching**: Has `MapsCaching` option for static maps but not enabled by default -5. **Output Chunks**: Uses `OutputMapsChunks` setting (default 1) - could benefit from larger chunks - -## Codebase Overview - -### File Structure -- **Language**: Python (with Numba for performance-critical sections) -- **Main modules**: - - `src/lisflood/global_modules/netcdf.py` - NetCDF I/O handling (584 lines) - - `src/lisflood/global_modules/output.py` - Output writing (586 lines) - - `src/lisflood/global_modules/add1.py` - Map loading/compression (986 lines) - - `src/lisflood/Lisflood_dynamic.py` - Main time loop (269 lines) - - `src/lisflood/hydrological_modules/` - Hydrological process modules - -### Key I/O Files -| File | Purpose | -|------|---------| -| [`netcdf.py`](src/lisflood/global_modules/netcdf.py) | NetCDF reading/writing, xarray chunked readers | -| [`output.py`](src/lisflood/global_modules/output.py) | Output writers (NetCDF, PCRaster) | -| [`add1.py`](src/lisflood/global_modules/add1.py) | Map loading, compression, decompression | -| [`zusatz.py`](src/lisflood/global_modules/zusatz.py) | NetCDF file access utilities | - ---- - -## I/O Bottlenecks - -### Identified Issues - -| Location | Issue | Estimated Impact | Complexity | -|----------|-------|------------------|------------| -| [`netcdf.py:478`](src/lisflood/global_modules/netcdf.py:478) | Output dtype from binding is `float64` by default | High (50% storage) | Low | -| [`netcdf.py:572-574`](src/lisflood/global_modules/netcdf.py:572) | zlib compression with chunksizes=(1,nrow,ncol) - time dimension chunk=1 | Medium | Medium | -| [`add1.py:282`](src/lisflood/global_modules/add1.py:282) | `compressArray` converts to float64: `mapC.astype(float)` | Medium | Low | -| [`readmeteo.py:40`](src/lisflood/hydrological_modules/readmeteo.py:40) | xarray readers created per variable in __init__ | Low | Low | -| [`output.py:102`](src/lisflood/global_modules/output.py:102) | Single timestep writes: `nf1.variables[self.map_name][:, :] = map_np` | Medium | Medium | -| [`netcdf.py:265`](src/lisflood/global_modules/netcdf.py:265) | `chunk.load()` loads data synchronously | Low | Medium | - -### Detailed Analysis - -#### 1. Output Data Type (High Impact) -**Location**: [`lisfloodSettings_reference.xml:157`](src/lisfloodSettings_reference.xml:157) -```xml - -``` - -**Current behavior**: All output maps are written as float64 (8 bytes per value) - -**Recommendation**: Change default to `float32` - sufficient for most hydrological variables - -#### 2. Compression Chunk Size (Medium Impact) -**Location**: [`netcdf.py:572`](src/lisflood/global_modules/netcdf.py:572) -```python -value = nf1.createVariable(var_name, dtype, ('time', dim_lat_y, dim_lon_x), - zlib=True, fill_value=-9999, chunksizes=(1, nrow, ncol)) -``` - -**Current behavior**: Time chunk size = 1 (one timestep per chunk) - -**Recommendation**: Use larger time chunks (e.g., daily or monthly) for better compression ratio - -#### 3. Array Type Conversion (Medium Impact) -**Location**: [`add1.py:282`](src/lisflood/global_modules/add1.py:282) -```python -return mapC.astype(float) # This defaults to float64 -``` - -**Current behavior**: All compressed arrays become float64 regardless of input precision - -**Recommendation**: Preserve input dtype or allow explicit dtype specification - ---- - -## Data Type Optimization Opportunities - -### Variable Analysis - -| Variable | Current dtype | Suggested dtype | Range | Storage Saving | Risk | -|----------|---------------|-----------------|-------|----------------|------| -| Discharge (ChanQ) | float64 | float32 | 0-100000 m³/s | 50% | Low | -| Soil Moisture (W1, W2) | float64 | float32 | 0-1 (fraction) | 50% | Very Low | -| Groundwater (UZ, LZ) | float64 | float32 | 0-∞ mm | 50% | Low | -| Precipitation | float64 | float32 | 0-500 mm/day | 50% | Very Low | -| Temperature | float64 | float32 | -50 to 50 °C | 50% | Very Low | -| Snow Cover | float64 | float32 | 0-100% | 50% | Very Low | -| Channel Storage | float64 | float32 | 0-∞ m³ | 50% | Low | -| LDD (flow directions) | int8 | int8 | 0-9 | 0% | N/A | -| Lake/Reservoir IDs | int16 | int16 | 0-65535 | 0% | N/A | -| Mask maps | boolean | boolean | 0-1 | 0% | N/A | - -### Compression Encoding Opportunities - -For integer variables that could benefit from scale_factor/add_offset (CF convention): - -| Variable | Current | Suggested | Notes | -|----------|---------|-----------|-------| -| LDD | int8 | int8 (no change) | Already optimal | -| Land Use | int8/int16 | int16 + scale_factor | Could use packing | -| Lake IDs | int16 | int16 (no change) | Already optimal | - ---- - -## Compression & Encoding Recommendations - -### 1. Enable float32 as Default -**File**: [`lisfloodSettings_reference.xml`](src/lisfloodSettings_reference.xml:157) -```xml - -``` - -### 2. Optimize NetCDF Chunk Sizes -**File**: [`netcdf.py`](src/lisflood/global_modules/netcdf.py:572) - -Current: -```python -chunksizes=(1, nrow, ncol) -``` - -Recommended: -```python -# For daily output: chunk = 1 day -# For monthly output: chunk = 30 days -optimal_time_chunk = min(30, n_timesteps) # Use up to 30 days -chunksizes=(optimal_time_chunk, nrow, ncol) -``` - -### 3. Add Compression Level Option -**File**: [`netcdf.py:572`](src/lisflood/global_modules/netcdf.py:572) - -Current: -```python -zlib=True -``` - -Recommended: -```python -# Add setting for compression level (1-9, default 4 for balanced speed/ratio) -compression_level = int(binding.get('OutputCompressionLevel', 4)) -value = nf1.createVariable(..., zlib=True, complevel=compression_level) -``` - -### 4. Preserve Input Data Types -**File**: [`add1.py:282`](src/lisflood/global_modules/add1.py:282) - -Current: -```python -return mapC.astype(float) -``` - -Recommended: -```python -# Preserve dtype or allow specification -def compressArray(map, pcr=True, name=None, dtype=None): - # ... existing code ... - if dtype is None: - dtype = map.dtype if hasattr(map, 'dtype') else np.float64 - return mapC.astype(dtype) -``` - ---- - -## Prioritized Action Plan - -| Rank | Recommendation | Storage Impact | Runtime Impact | Complexity | Risk | Score | -|------|----------------|----------------|----------------|------------|------|-------| -| 1 | Change default OutputMapsDataType to float32 | 50% | ~5% | Low | No | 55 | -| 2 | Optimize NetCDF time chunk sizes | 10-30% | 5-15% | Medium | No | 20 | -| 3 | Add compression level setting | 5-15% | -5-10% | Low | No | 15 | -| 4 | Enable MapsCaching by default | 0% | 10-30% | Low | No | 10 | -| 5 | Preserve input dtype in compressArray | 0-20% | 0% | Low | No | 5 | -| 6 | Add async I/O for output writing | 0% | 10-20% | High | Maybe | 3 | - -*Score = (Storage Impact + Runtime Impact) / Complexity* - ---- - -## Code Snippets / Diff Examples - -### Recommendation 1: Change Default Output Data Type - -**File**: `src/lisfloodSettings_reference.xml` - -```diff -- -+ -``` - -**Impact**: 50% reduction in output file sizes for all NetCDF maps - -**Risk**: Very low - float32 provides ~7 significant digits, sufficient for all hydrological variables - ---- - -### Recommendation 2: Optimize NetCDF Time Chunk Sizes - -**File**: `src/lisflood/global_modules/netcdf.py` - -```diff -@@ -569,7 +569,12 @@ 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)) -+ # Optimize chunk size for better compression -+ n_timesteps = steps.size -+ # Use up to 30 days of data per chunk, or less if total timesteps < 30 -+ time_chunk = min(30, n_timesteps) if n_timesteps > 1 else 1 -+ chunks = (time_chunk, nrow, ncol) -+ value = nf1.createVariable(var_name, dtype, ('time', dim_lat_y, dim_lon_x), zlib=True, fill_value=-9999, chunksizes=chunks) -``` - -**Impact**: 10-30% better compression ratio, faster I/O for large datasets - -**Risk**: Low - chunking is internal to NetCDF, doesn't affect model results - ---- - -### Recommendation 3: Add Compression Level Setting - -**File**: `src/lisflood/global_modules/netcdf.py` - -```diff -@@ -571,7 +571,10 @@ 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)) -+ # Get compression level from settings (default 4 for balanced speed/ratio) -+ comp_level = int(binding.get('OutputCompressionLevel', 4)) -+ comp_level = max(1, min(9, comp_level)) # Clamp to valid range -+ value = nf1.createVariable(var_name, dtype, ('time', dim_lat_y, dim_lon_x), zlib=True, fill_value=-9999, chunksizes=chunks, complevel=comp_level) -``` - -**File**: `src/lisfloodSettings_reference.xml` (add new setting) - -```xml - -The option "OutputCompressionLevel" sets the zlib compression level (1-9). -Higher values give better compression but slower I/O. Default is 4. - - -``` - -**Impact**: 5-15% additional compression, tunable performance tradeoff - -**Risk**: Low - compression is lossless, doesn't affect model results - ---- - -## Additional Recommendations - -### Enable MapsCaching by Default -The `MapsCaching` option already exists but defaults to "False". Enabling it by default would cache static maps (DEM, land use, soil properties) in memory, avoiding repeated disk reads. - -**File**: `src/lisfloodSettings_reference.xml` -```xml - -``` - -**Impact**: 10-30% faster initialization for large models with many static maps - -**Risk**: Low - only affects static data, doesn't change model results - -### Consider Float32 for Internal Calculations -For memory-constrained systems, consider using float32 for internal state variables. This would require: -1. Testing to ensure numerical accuracy is maintained -2. Modifying array creation in initialization code -3. Potential impact on accumulation variables over long simulations - ---- - -## Summary - -The LISFLOOD codebase has a solid foundation for I/O operations with NetCDF compression already enabled. The primary optimization opportunity is changing the default output data type from float64 to float32, which would provide immediate 50% storage savings with no impact on model results. Secondary optimizations around chunk sizes and compression levels can provide additional 10-30% improvements. - -The implementation complexity for all recommended changes is low to medium, with minimal risk to model accuracy. The recommended changes are backward-compatible and can be implemented incrementally. \ No newline at end of file From c7a5575b1fc7f5323de64f96a77b43e8fd88d498 Mon Sep 17 00:00:00 2001 From: Carlo Russo <51730707+doc78@users.noreply.github.com> Date: Tue, 8 Sep 2026 10:29:04 +0200 Subject: [PATCH 7/7] Delete lisflood_optimization_report_v2.md --- lisflood_optimization_report_v2.md | 1042 ---------------------------- 1 file changed, 1042 deletions(-) delete mode 100644 lisflood_optimization_report_v2.md diff --git a/lisflood_optimization_report_v2.md b/lisflood_optimization_report_v2.md deleted file mode 100644 index 0a7236aa..00000000 --- a/lisflood_optimization_report_v2.md +++ /dev/null @@ -1,1042 +0,0 @@ -# LISFLOOD Optimization Report v2 - -## Executive Summary - -This report identifies concrete opportunities to reduce wall-clock execution time and storage requirements across the LISFLOOD hydrological model codebase. It supersedes the earlier I/O-focused optimization report by covering the full computational pipeline including the newly added Muskingum-Cunge-Todini (MCT) routing module. - -Key findings: - -1. **MCT Routing Inner Loop** — The `mct_routing` kernel uses `prange` over pixels within each topological order, but orders are processed serially. The Newton-Raphson solver in `MCTRouting_single` (called twice per pixel) dominates runtime for MCT-enabled runs. -2. **Kinematic Wave Parallel Routing** — Already well-optimized with Numba `@njit(parallel=True)`, but the `numexpr` constant-term evaluation before calling the Numba kernel adds Python overhead every sub-step. -3. **Redundant Array Copies in Routing Loop** — Up to 10+ full-domain `.copy()` calls per routing sub-step when MCT + mass-balance reporting are active. -4. **I/O: Output Data Type** — Default `float64` output is still the factory setting; switching to `float32` halves output storage with no impact on results. -5. **I/O: NetCDF Chunk Size** — Time-dimension chunk of 1 limits compression ratio and read-back performance. -6. **Numba JIT Cold-Start** — First call to each `@njit` function triggers compilation; `cache=True` is already set but AOT pre-compilation could eliminate startup cost entirely. -7. **Memory: `compressArray` forces float64** — All compressed arrays are cast to `float64` regardless of input dtype. - - ---- - -## Part A — Runtime Performance - ---- - -### A1. MCT Routing: Serial Order Loop with Parallel Inner Loop - -| Attribute | Detail | -|-----------|--------| -| **File** | `src/lisflood/hydrological_modules/mct.py` | -| **Function** | `mct_routing()` (lines ~113-175) | -| **Impact** | **HIGH** | - -**Why it is a bottleneck:** -The MCT routing kernel is decorated `@njit(parallel=True)` but only the inner `prange(first, last)` loop over pixels within a single topological order is parallelised. The outer `for order in range(num_orders)` loop is inherently serial (data dependency: downstream pixels need upstream results). For domains with many orders but few pixels per order (long river stems), most iterations run with minimal parallelism. - -Inside the inner loop, `MCTRouting_single` is called per pixel. It contains: -- A Newton-Raphson iterative solver (`hoq`) with up to 1000 iterations per call -- Two full iterations of the MCT parameter calibration (`for i in range(2)`) -- Multiple calls to `qoh` (Manning-based Q-h relationships) - -For EFAS (~7 million pixels, ~5000 MCT pixels, ~200 topological orders), the MCT kernel can account for 30-50% of the routing sub-step time. - -**Recommendations:** - -1. **Reduce Newton-Raphson iterations in `hoq`**: The convergence tolerance is `1e-6` and max iterations is 1000. Profile to check actual average iteration count. If typically < 20, the overhead is acceptable. If convergence is slow for certain geometries, consider providing a better initial guess from the previous timestep's water depth (store `y_prev` as state). - -2. **Vectorise `MCTRouting_single` for batches**: Instead of calling a scalar function per pixel inside `prange`, restructure to pass arrays of pixel data for each order-batch and process them with vectorised NumPy/Numba array operations. This would enable SIMD and reduce function-call overhead. - -3. **Pre-compute static derived quantities**: `np.arctan(1 / ChanSdXdY[kinpix])` is computed every timestep for every pixel. Store `ANalv` as a pre-computed array in `MCTWave.__init__`. - -4. **Consider adaptive sub-stepping**: When Courant number < 0.5 for most pixels, the MCT solution is over-resolved. Allow the routing module to skip MCT computation for pixels where flow conditions have not materially changed (delta-Q threshold). - - ---- - -### A2. MCT Routing: Redundant `arctan` Computation Per Pixel Per Timestep - -| Attribute | Detail | -|-----------|--------| -| **File** | `src/lisflood/hydrological_modules/mct.py` | -| **Function** | `mct_routing()`, inner loop body | -| **Impact** | **MEDIUM** | - -**Current code (inside prange loop):** -```python -ANalv = np.arctan(1 / ChanSdXdY[kinpix]) -``` - -**Why it matters:** -`arctan` is a transcendental function computed per pixel, per routing sub-step, per model timestep. For EFAS with 4 routing sub-steps and 5000 MCT pixels, this is 20000 `arctan` calls per model timestep — all with static input. The result never changes. - -**Recommendation:** -Pre-compute `ANalv` once during `MCTWave.__init__()` and pass it as an array to `mct_routing`: - -```python -# In MCTWave.__init__: -self.ANalv = np.arctan(1.0 / ChanSdXdY) - -# In mct_routing signature: add ANalv parameter -# In the inner loop: replace np.arctan(...) with ANalv[kinpix] -``` - -**Estimated saving:** ~2-5% of MCT kernel time (removes transcendental from hot loop). - - ---- - -### A3. Kinematic Wave: `numexpr` Overhead Before Numba Kernel - -| Attribute | Detail | -|-----------|--------| -| **File** | `src/lisflood/hydrological_modules/kinematic_wave_parallel.py` | -| **Function** | `kinematicWave.kinematicWaveRouting()` | -| **Impact** | **MEDIUM** | - -**Current code:** -```python -lateral_inflow = nx.evaluate("q * dx", local_dict={"q": specific_lateral_inflow, "dx": self.space_delta}) -constant = nx.evaluate("a_dx_div_dt * Qold ** b + lateral_inflow", local_dict={...}) -``` - -**Why it matters:** -These `numexpr.evaluate` calls are executed every routing sub-step (typically 2-8 times per model timestep). Each call: -- Creates temporary arrays (full domain size) -- Involves Python-level dict construction and string parsing -- Runs a multi-threaded expression evaluator that competes with Numba's own thread pool - -For the kinematic routing (which applies to ALL channel pixels including those later overwritten by MCT), this adds measurable Python overhead on every sub-step. - -**Recommendation:** -Move the constant-term computation into the Numba `kinematicRouting` kernel itself. The expressions are simple element-wise operations that Numba can fuse into the main loop without allocating intermediates: - -```python -# Inside kinematicRouting (already @njit parallel): -# Replace 'constant' parameter with raw inputs -# Compute constant[pix] = a_dx_div_dt[pix] * discharge[pix]**beta + lateral_inflow[pix] -# directly at point of use -``` - -Alternatively, replace `numexpr` with a small `@njit` helper that computes `constant` in-place: - -```python -@njit(parallel=True, cache=True) -def compute_constant(constant, a_dx_div_dt, discharge, beta, lateral_inflow, space_delta, specific_lateral_inflow): - for pix in prange(constant.size): - constant[pix] = a_dx_div_dt[pix] * discharge[pix]**beta + specific_lateral_inflow[pix] * space_delta[pix] -``` - -**Estimated saving:** 5-10% of kinematic routing time (eliminates temporary array allocations and Python-level overhead per sub-step). - - ---- - -### A4. Routing Loop: Excessive `.copy()` Calls Per Sub-Step - -| Attribute | Detail | -|-----------|--------| -| **File** | `src/lisflood/hydrological_modules/routing.py` | -| **Function** | `routing.dynamic()` | -| **Impact** | **MEDIUM** | - -**Current code pattern (inside the `NoRoutSteps` loop):** -```python -SideflowChanM3 = self.var.ToChanM3RunoffDt.copy() -ChanQ_0 = self.var.ChanQ.copy() -ChanM3_0 = self.var.ChanM3.copy() -ChanQ = self.var.ChanQKin.copy() -ChanM3 = self.var.ChanM3Kin.copy() -ChanQAvgDt = self.var.ChanQKinAvgDt.copy() -ChanQAvgDt_old = self.var.ChanQAvgDt.copy() -# + additional copies inside repMBTs block -``` - -**Why it matters:** -Each `.copy()` allocates a new array of size `num_channel_pixels` (e.g. ~1.8M for EFAS). With MCT + mass-balance reporting active, there are ~10-12 full-domain copies per sub-step. With 4 sub-steps, that is ~40-48 array allocations per model timestep — pure memory allocation + memcpy overhead. - -**Recommendation:** -1. Pre-allocate scratch buffers once in `routing.initial()` and reuse them: -```python -# In initial(): -self._buf_ChanQ_0 = np.empty_like(self.var.ChanQ) -self._buf_ChanM3_0 = np.empty_like(self.var.ChanM3) - -# In dynamic(): -np.copyto(self._buf_ChanQ_0, self.var.ChanQ) # reuse buffer, no allocation -``` - -2. For `SideflowChanM3`: it starts from `ToChanM3RunoffDt` and then has values added/subtracted. Use an in-place pattern: -```python -SideflowChanM3 = self._buf_sideflow -np.copyto(SideflowChanM3, self.var.ToChanM3RunoffDt) -if option['openwaterevapo']: - SideflowChanM3 -= self.var.EvaAddM3Dt -# ... -``` - -3. For the `repMBTs` block: many `.copy()` calls exist solely to avoid modifying the original. Use indexing with `np.where` or masked assignment instead. - -**Estimated saving:** 3-8% of total routing time (reduces GC pressure and memcpy for large domains). - - ---- - -### A5. Numba JIT Cold-Start Compilation Overhead - -| Attribute | Detail | -|-----------|--------| -| **Files** | `mct.py`, `kinematic_wave_parallel_tools.py`, `soilloop.py` | -| **Functions** | All `@njit` decorated functions | -| **Impact** | **MEDIUM** (one-time cost, significant for short runs / calibration) | - -**Why it matters:** -Although `cache=True` is set on all Numba-compiled functions, the cache is invalidated whenever: -- The source file changes (even a comment) -- Numba or NumPy is upgraded -- The function signature changes due to different input dtypes - -During calibration workflows (1000s of short runs), the first run in each new environment pays 10-30 seconds of JIT compilation. For the MCT module alone, there are 6 `@njit` functions. - -**Recommendations:** - -1. **Ahead-of-Time (AOT) compilation**: Use `numba.pycc` to pre-compile the performance-critical kernels into a shared library. This eliminates JIT overhead entirely: -```python -from numba.pycc import CC -cc = CC('lisflood_routing_compiled') - -@cc.export('mct_routing', '...') -def mct_routing(...): ... - -cc.compile() -``` - -2. **Warm-up script**: Provide a lightweight `warmup_numba.py` that imports all JIT functions and calls them once with tiny dummy arrays. Run this as part of container/environment setup. - -3. **Pin Numba + NumPy versions** in Docker/conda environments to avoid cache invalidation across runs. - -**Estimated saving:** 10-30 seconds per cold start; near-zero for warm cache. - - ---- - -### A6. Surface Routing: Three Separate Kinematic Wave Calls - -| Attribute | Detail | -|-----------|--------| -| **File** | `src/lisflood/hydrological_modules/surface_routing.py` | -| **Function** | `surface_routing.dynamic()` | -| **Impact** | **MEDIUM** | - -**Current code:** -```python -self.direct_surface_router.kinematicWaveRouting(self.var.OFQDirectAvg, self.var.OFQDirect, SideflowDirect) -self.other_surface_router.kinematicWaveRouting(self.var.OFQOtherAvg, self.var.OFQOther, SideflowOther) -self.forest_surface_router.kinematicWaveRouting(self.var.OFQForestAvg, self.var.OFQForest, SideflowForest) -``` - -**Why it matters:** -Three independent kinematic wave routing calls are made sequentially, each with its own Numba kernel invocation. Each call processes the full overland-flow pixel domain. Since these three calls are independent (different land-use fractions with no inter-dependency), they could be batched or run concurrently. - -**Recommendations:** - -1. **Batch into a single kernel call**: Modify the kinematic routing kernel to accept a "batch" dimension (3 land-use types). Process all three in a single Numba `prange` call, tripling the available parallelism per order: -```python -# Single call handles all 3 land-use fractions -kinematicRoutingBatched(discharge_avg_batch, discharge_batch, lateral_inflow_batch, ...) -``` - -2. **Alternatively, use Python threading**: Since the three calls release the GIL (Numba `nogil=True` is implicit in `parallel=True`), they can be dispatched to a `ThreadPoolExecutor` with 3 workers. Each call then runs on its own Numba thread pool subset. - -**Estimated saving:** 10-25% of surface routing time (better core utilisation). - - ---- - -### A7. Soil Loop: Numba Thread Contention with Numexpr - -| Attribute | Detail | -|-----------|--------| -| **File** | `src/lisflood/hydrological_modules/soilloop.py` | -| **Functions** | `soilColumnsWaterBalance`, `interception_water_balance` | -| **Impact** | **LOW-MEDIUM** | - -**Why it matters:** -The soil module uses both `numexpr` (via `nx.evaluate`) and Numba `@njit(parallel=True)` with `prange`. Both libraries spawn their own thread pools: -- Numba uses the TBB or OpenMP backend -- numexpr uses its own thread pool (default: number of cores) - -When both are active in the same process, they compete for CPU cores. This over-subscription can cause context-switching overhead, especially on HPC nodes with many cores. - -**Recommendations:** - -1. **Limit numexpr threads**: Set `numexpr.set_num_threads(1)` when Numba parallel functions are the primary workload. Or coordinate thread counts: `NUMEXPR_MAX_THREADS=4` and `NUMBA_NUM_THREADS=N-4`. - -2. **Replace remaining numexpr calls with Numba**: The soil module already has extensive Numba coverage. The few remaining `nx.evaluate` calls (in `kinematic_wave_parallel.py`) can be absorbed into Numba kernels, eliminating the second thread pool entirely. - -3. **Set `NUMBA_THREADING_LAYER=tbb`** explicitly in the environment to ensure deterministic thread management. - -**Estimated saving:** 2-5% overall on many-core systems (16+ cores). - - ---- - -### A8. I/O: Synchronous Chunk Loading in XarrayChunked - -| Attribute | Detail | -|-----------|--------| -| **File** | `src/lisflood/global_modules/netcdf.py` | -| **Function** | `XarrayChunked.load_next_chunk()` and `XarrayChunked.__getitem__()` | -| **Impact** | **LOW-MEDIUM** | - -**Current code:** -```python -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() # blocks until data is in memory -``` - -**Why it matters:** -When a new temporal chunk boundary is crossed, the model blocks while xarray loads the next chunk from disk. For 5 forcing variables (Precip, Tavg, ET0, ES0, E0), this happens synchronously in `readmeteo.dynamic()`. Each chunk load involves NetCDF decompression (zlib) which is CPU-bound. - -**Recommendations:** - -1. **Prefetch next chunk asynchronously**: Use a background thread to load the next chunk before it is needed: -```python -import threading - -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() - -def prefetch_next_chunk(self): - if self.ichunk + 1 < len(self.chunk_indexes) - 1: - self._prefetch_thread = threading.Thread(target=self._prefetch) - self._prefetch_thread.start() -``` - -2. **Increase chunk size**: The setting `NetCDFTimeChunks` controls chunk size. For daily forcing over a year, set to 365 (load entire year at once). Memory cost is modest: 5 variables * 1.8M pixels * 365 days * 4 bytes = ~13 GB for float32 EFAS forcing loaded fully in memory. - -3. **Use `MapsCaching=True`** for forcing data when memory permits — loads the entire time series at initialization. - -**Estimated saving:** 1-5% for chunked reads; up to 10% if chunk boundaries align with expensive computation. - - ---- - -### A9. Output Writing: Per-Timestep NetCDF Writes - -| Attribute | Detail | -|-----------|--------| -| **File** | `src/lisflood/global_modules/output.py` | -| **Functions** | `NetcdfStepsWriter.stage()`, `NetcdfStepsWriter.write()` | -| **Impact** | **LOW-MEDIUM** | - -**Current behaviour:** -The `OutputMapsChunks` setting (default 1) controls how many timesteps are buffered before writing. With default=1, every reporting timestep triggers: -1. `uncompress_array()` — reconstruct 2D map from compressed 1D array -2. File open (or keep-open via `iterOpenNetcdf`) -3. Single-slice write to NetCDF variable -4. File close (if chunk boundary) - -**Why it matters:** -For runs reporting every timestep (e.g. hourly discharge maps for a year = 8760 writes), the overhead of repeated file I/O becomes significant. Each write involves: -- Python-level overhead of netCDF4 library calls -- OS-level file metadata updates -- zlib compression of each 2D slice - -**Recommendations:** - -1. **Increase `OutputMapsChunks`**: Set to 30-365 depending on output frequency. This buffers multiple timesteps in memory and writes them in a single batch, improving compression ratio and reducing file I/O overhead. - -2. **Use the existing `OutputMapsFactoryThreads` class**: The codebase already contains an async-write implementation using `ThreadPool`. It is documented as "NOT FULLY TESTED" but the approach is sound. Validate and enable it as an option: -```xml - -``` - -3. **Batch `uncompress_array` calls**: When writing multiple output variables at the same timestep, the decompression mask operation is repeated for each variable. Cache the mask indexing. - -**Estimated saving:** 5-15% for I/O-heavy configurations (many output maps, frequent reporting). - - ---- - -### A10. Dynamic Loop: Python-Level Overhead Per Timestep - -| Attribute | Detail | -|-----------|--------| -| **File** | `src/lisflood/Lisflood_dynamic.py` | -| **Function** | `LisfloodModel_dyn.dynamic()` | -| **Impact** | **LOW** | - -**Why it matters:** -Each model timestep involves: -- `datetime.timedelta` computation and `strftime` formatting -- Multiple `if option[...]` checks (dict lookups) -- `sys.stdout.write` and `sys.stdout.flush` for progress reporting - -For sub-hourly simulations (e.g. 15-minute timesteps over multiple years), there are 100K+ timesteps. Python-level overhead accumulates. - -**Recommendations:** - -1. **Cache option flags as local booleans** at the start of `dynamic()`: -```python -opt_MCT = option['MCTRouting'] -opt_split = option['SplitRouting'] -opt_repMBTs = option['repMBTs'] -# use local variables in if-statements -``` - -2. **Reduce progress output frequency**: Print progress every N steps rather than every step: -```python -if not flags['quiet'] and i % 100 == 0: - sys.stdout.write(...) -``` - -3. **Avoid `uuid.uuid4()` on first step**: The CDFFlags initialization creates a UUID on timestep 1. This is a system call that can be slow on some platforms. Move to `initial()`. - -**Estimated saving:** 1-2% for very long simulations with many timesteps. - - ---- - -### A11. MCT Routing: Calibration Points Check Inside Hot Loop - -| Attribute | Detail | -|-----------|--------| -| **File** | `src/lisflood/hydrological_modules/mct.py` | -| **Function** | `mct_routing()`, inner loop | -| **Impact** | **LOW-MEDIUM** | - -**Current code:** -```python -for ups_ix in range(num_upstream_pixels[kinpix]): - ups_pix = upstream_pixels[ups_ix] - if np.any(CalibPointsIds == ups_pix): - ql += ChanQAvgDt[ups_pix] - else: - q00 += ChanQ_0[ups_pix] - q0m += ChanQAvgDt[ups_pix] - q01 += ChanQ[ups_pix] -``` - -**Why it matters:** -`np.any(CalibPointsIds == ups_pix)` performs a linear scan of the `CalibPointsIds` array for every upstream pixel of every MCT pixel, on every routing sub-step. If there are K calibration points and U upstream connections per pixel, this is O(K * U * num_MCT_pixels * num_substeps) comparisons. - -For EFAS with ~50 calibration points and 5000 MCT pixels, this check dominates the inner loop when K is non-trivial. - -**Recommendation:** -Replace the linear scan with a pre-computed boolean lookup array: -```python -# In MCTWave.__init__: -is_calib_point = np.zeros(num_all_pixels, dtype=np.bool_) -is_calib_point[CalibPointsIds] = True - -# In mct_routing inner loop: -if is_calib_point[ups_pix]: - ... -``` - -This changes the check from O(K) to O(1) per upstream pixel. - -**Estimated saving:** 5-15% of MCT kernel time when calibration points are active. - - ---- - -### A12. Water Balance Module: Repeated `np.bincount` + `np.take` Pattern - -| Attribute | Detail | -|-----------|--------| -| **File** | `src/lisflood/hydrological_modules/waterbalance.py` | -| **Function** | `waterbalance.dynamic()` | -| **Impact** | **LOW** | - -**Current code pattern (repeated ~10 times):** -```python -WaterIn += np.take(np.bincount(self.var.Catchments, weights=some_array), self.var.Catchments) -``` - -**Why it matters:** -Each `np.bincount(..., weights=...)` creates a temporary array of size `max(Catchments)+1`, then `np.take` expands it back to full domain size. With 10+ such calls per timestep, this creates substantial temporary memory allocation and cache pressure. - -**Recommendations:** - -1. **Accumulate weights first, then do a single bincount**: Instead of calling bincount separately for each variable, sum the weighted arrays first: -```python -total_weights = self.var.TotalPrecipitationWB * self.var.MMtoM3 -if option['inflow']: - total_weights += self.var.sumInWB -WaterIn = np.take(np.bincount(self.var.Catchments, weights=total_weights), self.var.Catchments) -``` - -2. **Pre-allocate the bincount result buffer**: Reuse the same output array across calls. - -3. **Consider skipping mass balance for production runs**: The `repMBTs` option enables extensive per-catchment accounting. For operational forecasting (where mass balance is not the focus), disabling this saves all associated computation. - -**Estimated saving:** 1-3% of total timestep time. - - ---- - -## Part B — Storage & Memory Optimization - ---- - -### B1. Output Data Type: Default float64 - -| Attribute | Detail | -|-----------|--------| -| **File** | `src/lisfloodSettings_reference.xml` (line ~206) | -| **Setting** | `OutputMapsDataType` | -| **Impact** | **HIGH** (50% output storage reduction) | - -**Current default:** -```xml - -``` - -**Why it matters:** -All output NetCDF maps are written as float64 (8 bytes per value). For typical hydrological variables (discharge 0-100000 m³/s, soil moisture 0-1, temperature -50 to 50°C), float32 provides ~7 significant digits — more than sufficient. - -For an EFAS daily discharge output: 1.8M pixels * 8 bytes * 365 days = ~4.8 GB/year in float64 vs 2.4 GB/year in float32. - -**Recommendation:** -Change the default to `float32`: -```xml - -``` - -**Risk:** Very low. Float32 precision (7 significant digits) exceeds the physical accuracy of any hydrological variable LISFLOOD produces. No model results are affected — this is output-only. - -**Estimated saving:** 50% reduction in all output map file sizes. - - ---- - -### B2. NetCDF Time Chunk Size - -| Attribute | Detail | -|-----------|--------| -| **File** | `src/lisflood/global_modules/netcdf.py` (line ~572) | -| **Function** | `write_netcdf_header()` | -| **Impact** | **MEDIUM** (10-30% better compression) | - -**Current code:** -```python -value = nf1.createVariable(var_name, dtype, ('time', dim_lat_y, dim_lon_x), - zlib=True, fill_value=-9999, chunksizes=(1, nrow, ncol)) -``` - -**Why it matters:** -A time-chunk of 1 means each timestep is compressed independently. Hydrological fields exhibit strong temporal autocorrelation — consecutive timesteps are similar. Compressing multiple timesteps together exploits this redundancy for better compression ratios. - -Additionally, reading back data for post-processing (e.g. extracting a time series at one pixel) requires decompressing full spatial slices one at a time when time-chunk=1. - -**Recommendation:** -Use a time chunk that matches typical access patterns: -```python -# For daily output files: chunk 30 days together -time_chunk = min(30, n_timesteps) if n_timesteps > 1 else 1 -chunksizes = (time_chunk, nrow, ncol) -``` - -For hourly outputs (sub-daily), consider `time_chunk = min(24, n_timesteps)` (one day of hours). - -Make this configurable via a new setting: -```xml - -``` - -**Risk:** Low — chunking is internal to NetCDF format; does not affect data values. - -**Estimated saving:** 10-30% better compression ratio; faster time-series extraction in post-processing. - - ---- - -### B3. Compression Level Setting - -| Attribute | Detail | -|-----------|--------| -| **File** | `src/lisflood/global_modules/netcdf.py` | -| **Function** | `write_netcdf_header()` | -| **Impact** | **LOW-MEDIUM** | - -**Current code:** -```python -zlib=True # uses default complevel (typically 4 in netCDF4 library) -``` - -**Why it matters:** -The `complevel` parameter is not explicitly set, defaulting to whatever the netCDF4 library uses (usually 4). For operational runs where output I/O time is a bottleneck, a lower compression level (1-2) can significantly speed up writes at the cost of slightly larger files. For archival runs, a higher level (6-9) produces smaller files at the cost of write speed. - -**Recommendation:** -Add an explicit compression level setting: -```python -comp_level = int(binding.get('OutputCompressionLevel', 4)) -value = nf1.createVariable(..., zlib=True, complevel=comp_level, ...) -``` - -```xml - -``` - -**Trade-off guidance:** -| Level | Write speed | File size | Use case | -|-------|-------------|-----------|----------| -| 1 | Fast | Larger (+20%) | Real-time forecasting | -| 4 | Balanced | Baseline | General use | -| 6-9 | Slow | Smaller (-10-20%) | Long-term archival | - -**Estimated saving:** Configurable — up to 20% smaller files (level 6+) or 30% faster writes (level 1). - - ---- - -### B4. `compressArray` Forces float64 on All Arrays - -| Attribute | Detail | -|-----------|--------| -| **File** | `src/lisflood/global_modules/add1.py` (line ~291) | -| **Function** | `compressArray()` | -| **Impact** | **MEDIUM** (memory footprint) | - -**Current code:** -```python -return mapC.astype(float) # float = float64 in NumPy -``` - -**Why it matters:** -Every map loaded through `compressArray` is cast to float64 regardless of its original dtype. This includes: -- Boolean masks (could be bool/int8 → 1 byte becomes 8 bytes) -- Integer maps like LDD, land-use classes (int8/int16 → float64) -- Maps read as float32 from NetCDF files (4 bytes → 8 bytes) - -For EFAS with ~1.8M pixels, each unnecessary float64 promotion costs 14.4 MB (from float32) or 12.6 MB (from int8). With dozens of static maps loaded at initialization, this adds up to 200-500 MB of wasted memory. - -**Recommendation:** -Preserve the input dtype or allow explicit specification: -```python -def compressArray(map, pcr=True, name=None, force_load_with_nans=False, dtype=None): - # ... existing logic ... - if dtype is not None: - return mapC.astype(dtype) - elif hasattr(mapC, 'dtype') and not np.issubdtype(mapC.dtype, np.floating): - # Keep integer/boolean maps in their native type - return mapC - else: - return mapC.astype(np.float64) # backward compatible default for float maps -``` - -For a less invasive change, at minimum avoid promoting boolean and integer maps: -```python -if np.issubdtype(mapC.dtype, np.integer) or np.issubdtype(mapC.dtype, np.bool_): - return mapC -return mapC.astype(float) -``` - -**Risk:** Low if done carefully. Some downstream code may assume float64; a `grep` for `.astype(float)` patterns in dependent code should be checked. - -**Estimated saving:** 200-500 MB RAM reduction for large EFAS/GloFAS domains. - - ---- - -### B5. MCT State Variables: Extra Arrays Per Pixel - -| Attribute | Detail | -|-----------|--------| -| **File** | `src/lisflood/hydrological_modules/routing.py` | -| **Variables** | `PrevCm0`, `PrevDm0`, `ChanQ_0`, `ChanM3_0` (MCT-specific state) | -| **Impact** | **LOW** | - -**Why it matters:** -MCT routing introduces 4 additional full-domain float64 arrays as state variables: -- `PrevCm0` (Courant number) — only meaningful at MCT pixels -- `PrevDm0` (Reynolds number) — only meaningful at MCT pixels -- `ChanQ_0` (previous discharge copy) — full domain -- `ChanM3_0` (previous storage copy) — full domain - -For EFAS (~1.8M pixels): 4 arrays * 1.8M * 8 bytes = ~57 MB. However, only ~5000 pixels are MCT pixels (0.3% of domain). - -**Recommendation:** -Store `PrevCm0` and `PrevDm0` as MCT-only compressed arrays (5000 pixels) rather than full-domain arrays: -```python -# In MCTWave.__init__: -self.PrevCm0_mct = np.ones(num_mct_pixels) -self.PrevDm0_mct = np.zeros(num_mct_pixels) -``` - -Map between MCT-local and full-domain indices using the existing `mapping_mct` array. This reduces memory by ~25 MB for these two arrays alone. - -For `ChanQ_0` and `ChanM3_0`: these are full-domain copies needed as "previous state" inputs. They could be eliminated by reorganizing the MCT call to read directly from `self.var.ChanQ` before kinematic routing overwrites it (restructure call order). - -**Estimated saving:** 25-57 MB RAM. - - ---- - -### B6. Enable MapsCaching by Default for Static Maps - -| Attribute | Detail | -|-----------|--------| -| **File** | `src/lisfloodSettings_reference.xml` | -| **Setting** | `MapsCaching` | -| **Impact** | **MEDIUM** (runtime), **LOW** (memory trade-off) | - -**Current default:** -```xml - -``` - -**Why it matters:** -When `MapsCaching` is `False`, the `XarrayChunked` reader is used for forcing data. When `True`, the `XarrayCached` class (which uses `@Cache` decorator) stores the full dataset in memory, avoiding repeated disk reads. - -For static maps loaded via `loadmap()`, caching avoids re-reading the same NetCDF file if it is referenced multiple times (e.g. soil parameters used in both initialization and dynamic sections). - -**Recommendation:** -Enable caching by default: -```xml - -``` - -For memory-constrained environments, document the memory cost and allow users to disable it. - -**Memory cost:** For EFAS forcing data (5 variables, 1.8M pixels, 365 days, float32): ~13 GB. This is acceptable on modern HPC nodes (128-512 GB) but may be prohibitive on smaller machines. - -**Estimated saving:** 10-30% faster initialization; eliminates repeated disk reads for static data during the run. - - ---- - -### B7. Output Variable Pruning: Map vs TSS Trade-offs - -| Attribute | Detail | -|-----------|--------| -| **File** | Settings XML (user configuration) | -| **Impact** | **HIGH** (operational storage) | - -**Why it matters:** -LISFLOOD can output dozens of 2D map variables at every timestep. For EFAS operational runs, the output volume is often the dominant storage cost. Many variables are only needed at specific gauge points (time series) rather than as full 2D maps. - -**Recommendations:** - -1. **Prefer TSS over Maps for point-based diagnostics**: Variables like discharge, water level, and lake/reservoir levels are typically only needed at gauge locations. Use TSS (Time Series) output instead of full maps: - - Storage for 1 variable, 1000 gauges, 1 year daily: ~3 MB (TSS) vs ~2.4 GB (map, float32) - - Factor 800x reduction - -2. **Use monthly/yearly frequency for slow-changing variables**: Soil moisture, groundwater, snow cover change slowly. Report them monthly rather than daily: - - 12x reduction in output volume for these variables - -3. **Document output configurations**: Provide template settings files for common use cases: - - `settings_minimal_output.xml` — discharge TSS only - - `settings_operational.xml` — discharge maps + key state variables - - `settings_research.xml` — full variable set - -4. **Consider lossy compression for diagnostic maps**: For variables like soil moisture or temperature that are only used for visualization, NetCDF's `least_significant_digit` option can dramatically improve compression: -```python -value = nf1.createVariable(..., zlib=True, least_significant_digit=3) -# Keeps 3 significant digits, enables much better compression -``` - -**Estimated saving:** 50-95% output volume reduction depending on configuration. - - ---- - -### B8. Internal Computation Precision: Float64 Everywhere - -| Attribute | Detail | -|-----------|--------| -| **Files** | All hydrological modules | -| **Impact** | **LOW** (memory), **LOW** (runtime on modern CPUs) | - -**Current state:** -All internal state variables and computations use float64. This is the safe default for numerical stability, especially for: -- Accumulation variables over long simulations (mass balance) -- Channel routing where small flows interact with large storages -- The MCT solver which is sensitive to precision in Courant/Reynolds numbers - -**Recommendation:** -Do NOT change internal precision to float32 for state variables. The risk of numerical drift over multi-year simulations outweighs the modest memory savings. Float64 arithmetic is nearly as fast as float32 on modern x86-64 CPUs (same SIMD width for scalar operations; only vectorised code benefits from float32's 2x throughput). - -**Exception:** Intermediate temporary arrays that are not accumulated (e.g. `SideflowChan`, `SideflowDirect`) could safely use float32 if memory pressure is extreme, but the savings (~14 MB per temporary for EFAS) rarely justify the added complexity. - -**Conclusion:** Keep float64 for internal computation. Focus memory optimization on B4 (dtype preservation for non-float maps) and B1 (output dtype). - - ---- - -## Part C — Prioritized Action Plan - -| Rank | ID | Recommendation | Runtime Impact | Storage Impact | Complexity | Risk | -|------|----|----------------|----------------|----------------|------------|------| -| 1 | B9 | CF scale/offset packing (int16 output) | ~5% (I/O) | **75-90%** | Medium | Low | -| 2 | B1 | Change default OutputMapsDataType to float32 | ~5% (I/O) | **50%** | Low | None | -| 2 | A11 | MCT: Replace linear CalibPoints scan with boolean lookup | 5-15% MCT | — | Low | None | -| 3 | A2 | MCT: Pre-compute ANalv (arctan) | 2-5% MCT | — | Low | None | -| 4 | A3 | Kinematic wave: Replace numexpr with Numba helper | 5-10% routing | — | Medium | Low | -| 5 | B2 | Optimize NetCDF output time chunk size | 5% (I/O) | **10-30%** | Low | None | -| 6 | A4 | Routing: Pre-allocate scratch buffers, reduce .copy() | 3-8% routing | — | Medium | Low | -| 7 | B4 | compressArray: Preserve input dtype | — | **200-500 MB** | Low | Low | -| 8 | A6 | Surface routing: Batch 3 kinematic calls | 10-25% surface | — | Medium | Low | -| 9 | B7 | Output variable pruning (TSS vs Maps) | 5-15% (I/O) | **50-95%** | Config only | None | -| 10 | A1 | MCT: Better initial guess / vectorise batches | 10-20% MCT | — | High | Medium | -| 11 | B3 | Add configurable compression level | ±10% (I/O) | **±20%** | Low | None | -| 12 | A5 | Numba AOT / warm-up script | 10-30s startup | — | Medium | Low | -| 13 | B6 | Enable MapsCaching by default | 10-30% init | +memory | Low | Low | -| 14 | A8 | Async chunk prefetch for forcing I/O | 1-5% | — | Medium | Low | -| 15 | A9 | Increase OutputMapsChunks / enable async writes | 5-15% (I/O) | — | Low | Low | -| 16 | A7 | Resolve Numba/numexpr thread contention | 2-5% | — | Low | None | -| 17 | A12 | Water balance: reduce bincount calls | 1-3% | — | Low | None | -| 18 | A10 | Dynamic loop: cache options, reduce stdout | 1-2% | — | Low | None | -| 19 | B5 | MCT: Compress PrevCm0/PrevDm0 to MCT-only arrays | — | **25-57 MB** | Medium | Low | - - ---- - -## Part D — Quick Wins (Implementable in < 1 day each) - -### D1. Change default output dtype (B1) -**File:** `src/lisfloodSettings_reference.xml` -```diff -- -+ -``` - -### D2. Pre-compute ANalv in MCT (A2) -**File:** `src/lisflood/hydrological_modules/mct.py` - -In `MCTWave.__init__()`, add: -```python -self.ANalv = np.arctan(1.0 / ChanSdXdY) -``` - -Add `ANalv` as a parameter to `mct_routing()` and in the inner loop replace: -```python -ANalv = np.arctan(1 / ChanSdXdY[kinpix]) -``` -with: -```python -ANalv = precomputed_ANalv[kinpix] -``` - -### D3. Replace CalibPoints linear scan with boolean array (A11) -**File:** `src/lisflood/hydrological_modules/mct.py` - -In `MCTWave.__init__()` or `routing.initialMCT()`: -```python -num_pixels = len(ChanLength) -is_calib_point = np.zeros(num_pixels, dtype=np.bool_) -if len(CalibPointsIds) > 0: - is_calib_point[CalibPointsIds] = True -``` - -Pass `is_calib_point` to `mct_routing` and replace: -```python -if np.any(CalibPointsIds == ups_pix): -``` -with: -```python -if is_calib_point[ups_pix]: -``` - -### D4. Add compression level setting (B3) -**File:** `src/lisflood/global_modules/netcdf.py`, in `write_netcdf_header()`: -```python -comp_level = int(binding.get('OutputCompressionLevel', 4)) -value = nf1.createVariable(var_name, dtype, dims, zlib=True, complevel=comp_level, fill_value=-9999, chunksizes=chunks) -``` - -### D5. Increase default OutputMapsChunks (A9) -**File:** `src/lisfloodSettings_reference.xml` -```diff -- -+ -``` - - ---- - -## Part E — Changes from Previous Report - -| Previous Recommendation | Status | Notes | -|------------------------|--------|-------| -| Change OutputMapsDataType to float32 | **Still valid** | Retained as #1 priority | -| Optimize NetCDF time chunk sizes | **Still valid** | Code unchanged; chunksizes=(1,nrow,ncol) persists | -| Add compression level setting | **Still valid** | Not yet implemented | -| Enable MapsCaching by default | **Still valid** | Still defaults to False | -| Preserve input dtype in compressArray | **Still valid** | `mapC.astype(float)` unchanged | -| Add async I/O for output writing | **Partially implemented** | `OutputMapsFactoryThreads` exists but is not enabled by default | -| Consider float32 for internal calculations | **Superseded** | Now explicitly recommended AGAINST (see B8) — risk outweighs benefit | - -### New sections added (not in original report): - -| Section | Reason | -|---------|--------| -| A1, A2, A11 | MCT routing module did not exist | -| A3 | numexpr usage in kinematic_wave_parallel.py is new | -| A4 | .copy() pattern in routing loop amplified by MCT additions | -| A5 | Numba cold-start now affects more functions (MCT adds 6 new @njit) | -| A6 | Surface routing now uses parallel kinematic wave (was PCRaster before) | -| A7 | Thread contention newly relevant with parallel soil + parallel routing | -| A10 | Dynamic loop analysis refreshed against current code | -| A12 | Water balance module rewritten with new bincount patterns | -| B5 | MCT-specific state variable overhead is new | -| B7 | Output pruning guidance (operational vs research) | -| B8 | Explicit recommendation against float32 internals | - - ---- - -## Appendix: Profiling Guidance - -To validate the estimates in this report, the following profiling approach is recommended: - -### Wall-clock timing per module -Add timing instrumentation in `Lisflood_dynamic.py`: -```python -import time -t0 = time.perf_counter() -self.readmeteo_module.dynamic() -t_meteo = time.perf_counter() - t0 - -t0 = time.perf_counter() -self.soilloop_module.dynamic_soil() -t_soil = time.perf_counter() - t0 - -# ... etc for each module call -``` - -### Numba kernel timing -Use `numba.core.config.DEVELOPER_MODE = 1` and Numba's built-in timing, or wrap calls: -```python -t0 = time.perf_counter() -kwpt.kinematicRouting(...) -t_kin = time.perf_counter() - t0 -``` - -### Memory profiling -Use `tracemalloc` or `memory_profiler` to identify peak memory usage: -```python -import tracemalloc -tracemalloc.start() -# ... run model ... -snapshot = tracemalloc.take_snapshot() -top_stats = snapshot.statistics('lineno') -``` - -### I/O profiling -Monitor disk I/O with system tools (`iostat`, `iotop`) or Python's `cProfile` focused on netCDF4 calls. - ---- - -*Report generated: July 2026* -*Codebase version: current HEAD (post-MCT integration)* -*Previous report: `lisflood_optimization_report.md` (I/O-focused, pre-MCT)* - ---- - -### B9. CF-Convention Scale/Offset Packing (int16 storage of float data) - -| Attribute | Detail | -|-----------|--------| -| **File** | `src/lisflood/global_modules/netcdf.py` | -| **Function** | `write_netcdf_header()` | -| **Impact** | **HIGH** (75% output storage reduction vs float64; 50% vs float32) | - -**Current state:** -Output variables are stored as raw floating-point values (float64 or float32). No packing is applied. - -**What is scale/offset packing?** -The CF conventions define two variable attributes — `scale_factor` and `add_offset` — that allow floating-point data to be stored as smaller integer types (typically int16 or uint16). On read, the client library automatically reconstructs the original value: - -``` -unpacked_value = packed_value * scale_factor + add_offset -``` - -A 16-bit integer provides 65536 distinct values. By choosing `scale_factor` and `add_offset` to cover the physical range of the variable, you get fixed-precision storage at **2 bytes per value** instead of 4 (float32) or 8 (float64). - -**Why it matters:** -- float64 → int16: **75% reduction** (8 bytes → 2 bytes per value) -- float32 → int16: **50% reduction** (4 bytes → 2 bytes per value) -- Combined with zlib compression, the actual on-disk savings can exceed 80% because integer data with limited range compresses far better than floating-point data. -- This is standard practice in climate/meteorological data (CMIP6, ERA5, EFAS forcing data already use this). - -**Variable-specific packing parameters:** - -| Variable | Physical Range | scale_factor | add_offset | Precision | Suitable? | -|----------|---------------|--------------|------------|-----------|-----------| -| Discharge (m³/s) | 0 – 100,000 | 1.53 | 50,000 | ±0.8 m³/s | Yes (large rivers) | -| Discharge (m³/s) | 0 – 10,000 | 0.153 | 5,000 | ±0.08 m³/s | Yes (medium rivers) | -| Soil Moisture (fraction) | 0 – 1 | 1.53e-5 | 0.5 | ±8e-6 | Yes | -| Snow Water Equiv. (mm) | 0 – 2000 | 0.031 | 1000 | ±0.015 mm | Yes | -| Temperature (°C) | -50 – +50 | 0.00153 | 0 | ±0.001 °C | Yes | -| ET (mm/day) | 0 – 20 | 3.05e-4 | 10 | ±1.5e-4 mm | Yes | -| Water Level (m) | -5 – +20 | 3.82e-4 | 7.5 | ±0.2 mm | Yes | -| Groundwater (mm) | 0 – 5000 | 0.076 | 2500 | ±0.04 mm | Yes | - -**Implementation:** - -In `write_netcdf_header()`, add packing when a packing configuration is provided: - -```python -# Determine packing based on variable metadata or user settings -packing = binding.get('OutputPacking', 'none') # 'none', 'int16', 'auto' - -if packing == 'int16' and frequency is not None: - # Use variable-specific range or a safe default - vmin, vmax = get_variable_range(var_name) # lookup table - scale = (vmax - vmin) / 65534.0 # leave room for fill_value - offset = vmin + scale * 32767.0 - - value = nf1.createVariable(var_name, 'i2', ('time', dim_lat_y, dim_lon_x), - zlib=True, fill_value=-32767, - chunksizes=(time_chunk, nrow, ncol)) - value.scale_factor = scale - value.add_offset = offset -else: - # Current behaviour (float32 or float64) - value = nf1.createVariable(var_name, dtype, ('time', dim_lat_y, dim_lon_x), - zlib=True, fill_value=-9999, - chunksizes=(1, nrow, ncol)) -``` - -The write step must also pack the data before storing: -```python -# In NetcdfStepsWriter.write() or equivalent: -if hasattr(nf1.variables[self.map_name], 'scale_factor'): - scale = nf1.variables[self.map_name].scale_factor - offset = nf1.variables[self.map_name].add_offset - packed = np.round((map_np - offset) / scale).astype(np.int16) - packed[map_np == -9999] = -32767 # fill value - nf1.variables[self.map_name][step, :, :] = packed -else: - nf1.variables[self.map_name][step, :, :] = map_np -``` - -**Settings XML:** -```xml - - -``` - -**Compatibility:** -- All modern NetCDF readers (xarray, CDO, NCO, QGIS, Python netCDF4) automatically apply `scale_factor` and `add_offset` on read — no user action needed. -- CF-compliant: follows CF-1.6+ conventions exactly. -- Backward compatible: if `OutputPacking=none`, behaviour is unchanged. - -**Risks:** -- **Precision loss**: int16 provides ~4.8 significant digits (vs ~7 for float32, ~15 for float64). For most hydrological variables this is adequate. For variables with very large dynamic range (e.g. discharge spanning 0.001 to 100,000 m³/s), consider using per-timestep adaptive scaling or splitting into sub-ranges. -- **State/restart files should NOT be packed**: Warm-start state variables (used to restart the model) must retain full float64 precision to avoid drift. Only apply packing to reporting/output maps. - -**Estimated saving:** -- vs current float64 default: **75% output file size reduction** -- vs float32: **50% additional reduction on top of B1** -- Combined with zlib (already enabled): effective on-disk ratios of 85-90% reduction are achievable -