diff --git a/.gitignore b/.gitignore index ef5f354a..3ec87b21 100644 --- a/.gitignore +++ b/.gitignore @@ -1,6 +1,7 @@ *.gem *.sublime-project *.sublime-workspace +.kiro .bundle .DS_Store .jekyll-metadata @@ -27,3 +28,6 @@ lisflood_model.egg-info .vscode/ *.ipynb .ipynb_checkpoints/ +implementation_guide_scale_offset_packing.md +lisflood_optimization_report.md +lisflood_optimization_report_v2.md diff --git a/Dockerfile b/Dockerfile index 92ad8fb4..eac44765 100755 --- a/Dockerfile +++ b/Dockerfile @@ -2,7 +2,7 @@ # docker push jrce1/lisflood FROM continuumio/miniconda3 -LABEL maintainer="Stefania Grimaldi, Cinzia Mazzetti, Carlo Russo, Valerio Lorini, Ad de Roo" +LABEL maintainer="Stefania Grimaldi, Timo Schaffhauser, Carlo Russo, Cinzia Mazzetti, Corentin Carton De Wiart" ENV DEBIAN_FRONTEND=noninteractive diff --git a/README.md b/README.md index a8e9f3d4..538db81f 100644 --- a/README.md +++ b/README.md @@ -1,3 +1,5 @@ +[![DOI](https://img.shields.io/badge/DOI-10.5281%2Fzenodo.21107672-blue.svg)](https://doi.org/10.5281/zenodo.21107672) + # Lisflood OS This repository hosts source code of LISFLOOD model. diff --git a/docs/5_annex_output-files/index.md b/docs/5_annex_output-files/index.md index a03575f9..a1ef8598 100644 --- a/docs/5_annex_output-files/index.md +++ b/docs/5_annex_output-files/index.md @@ -157,7 +157,7 @@ To speed up the pre-run and to prevent that results are taken from the pre-run, | actual transpiration | repTaMaps | $\frac{mm}{timestep}$ | TaMaps | tact | | rainfall interception | repInterceptionMaps | $\frac{mm}{timestep}$ | InterceptionMaps
InterceptionForestMaps | int
intF | | evaporation of intercepted water | repEWIntMaps | $\frac{mm}{timestep}$ | EWIntMaps | ewint | -| leaf drainage | repLeafDrainageMaps | $\frac{mm}{timestep}$ | LeafDrainageMaps
LeafDrainageForestMaps | ldra
draF | +| leaf drainage | repLeafDrainageMaps | $\frac{mm}{timestep}$ | LeafDrainageMaps
LeafDrainageForestMaps no | ldra
draF | | infiltration | repInfiltrationMaps | $\frac{mm}{timestep}$ | InfiltrationMaps
InfiltrationForestMaps | inf
infF | | preferential (bypass) flow | repPrefFlowMaps | $\frac{mm}{timestep}$ | PrefFlowMaps
PrefFlowtherMaps
PrefFlowForestMaps
PrefFlowIrrigationMaps | pflowpixel
pflow
pflowF
pflowi | | percolation upper to lower soil layer | repPercolationMaps | $\frac{mm}{timestep}$ | Percolation1ato1bOtherMaps
Percolation1to1bForestMaps
Percolation1ato1bIrrigationMaps
Percolation1bto2OtherMaps
Percolation1bto2ForestMaps
Percolation1bto2IrrigationMaps | Percolation1ato1bOther
Percolation1ato1bForest
Percolation1to2Irrigation
Percolation1bto2Other
Percolation1bto2Forest
Percolation1bto2Irrigation | diff --git a/setup.py b/setup.py index eb672278..ba2c8c77 100755 --- a/setup.py +++ b/setup.py @@ -18,7 +18,7 @@ IMPORTANT Note: To test pip installation: python setup.py testpypi -pip install --index-url https://test.pypi.org/simple/ lisflood-model==5.0.0 +pip install --index-url https://test.pypi.org/simple/ --extra-index-url https://pypi.org/simple lisflood-model==5.0.0 To publish on PyPi: diff --git a/src/lisflood/global_modules/default_options.py b/src/lisflood/global_modules/default_options.py index 32de212c..d3a955c9 100644 --- a/src/lisflood/global_modules/default_options.py +++ b/src/lisflood/global_modules/default_options.py @@ -1,7 +1,8 @@ from collections import namedtuple TimeSeries = namedtuple('TimeSeries', 'name, output_var, where, repoption, restrictoption, operation') -ReportedMap = namedtuple('ReportedMap', 'name, output_var, unit, end, steps, all, restrictoption, monthly, yearly') +ReportedMap = namedtuple('ReportedMap', 'name, output_var, unit, end, steps, all, restrictoption, monthly, yearly, scale_factor, add_offset') +ReportedMap.__new__.__defaults__ = (None, None) # scale_factor, add_offset default_options = {'EnKF': False, 'InitLisflood': False, @@ -351,7 +352,7 @@ output_var='DirectRunoff', unit='mm', end=[], steps=['repSurfaceRunoffMaps'], all=[''], restrictoption=['nonInit'], monthly=False, - yearly=False), + yearly=False, scale_factor=0.0062, add_offset=200.0), 'DischargeEnd': ReportedMap(name='DischargeEnd', output_var='ChanQAvg', unit='m3/s', end=[], steps=[], all=[], restrictoption=[], monthly=False, @@ -368,38 +369,46 @@ 'ESActForestMaps': ReportedMap(name='ESActForestMaps', output_var='ESAct[1]', unit='mm', end=[], steps=[''], all=['repESActMaps'], restrictoption=['nonInit'], - monthly=False, yearly=False), + monthly=False, yearly=False,scale_factor=0.00031, add_offset=10.), 'ESActMaps': ReportedMap(name='ESActMaps', output_var='ESActPixel', unit='mm', end=[], steps=[], all=['repESActMaps','repE2O1'], - restrictoption=['nonInit'], monthly=False, yearly=False), + restrictoption=['nonInit'], monthly=False, yearly=False, + scale_factor=0.00031, add_offset=10.0), 'ESRefMapsOut': ReportedMap(name='ESRefMapsOut', output_var='ESRef', unit='mm', end=[], steps=[], all=['repESRefMaps'], - restrictoption=[], monthly=False, yearly=False), + restrictoption=[], monthly=False, yearly=False, + scale_factor=0.00031, add_offset=10.0), 'ETActMaps': ReportedMap(name='ETActMaps', output_var='ESActPixel+self.var.TaPixel+self.var.TaInterceptionAll+self.var.EvaAddM3*self.var.M3toMM', unit='mm', end=[], steps=[], all=['repETActMaps','repE2O1'], - restrictoption=['nonInit'], monthly=False, yearly=False), + restrictoption=['nonInit'], monthly=False, yearly=False, + scale_factor=0.016, add_offset=495.0), 'ETActBudykoMaps': ReportedMap(name='ETActBudykoMaps', output_var='ESActNoIrrPixel+self.var.TaNoIrrPixel+self.var.TaInterceptionNoIrrPixel', unit='mm', end=[], steps=[], all=['repETActMaps','repE2O1'], - restrictoption=['nonInit'], monthly=False, yearly=False), + restrictoption=['nonInit'], monthly=False, yearly=False, + scale_factor=0.016, add_offset=495.0), 'ETRefMapsOut': ReportedMap(name='ETRefMapsOut', output_var='ETRef', unit='mm', end=[], steps=[], all=['repETRefMaps'], - restrictoption=[], monthly=False, yearly=False), + restrictoption=[], monthly=False, yearly=False, + scale_factor=0.00031, add_offset=10.0), 'EWIntForestMaps': ReportedMap(name='EWIntForestMaps', output_var='TaInterception[1]', unit='mm', end=[], steps=[], all=['repEWIntMaps'], restrictoption=['nonInit'], monthly=False, - yearly=False), + yearly=False, scale_factor=0.00031, add_offset=10.0), 'EWIntMaps': ReportedMap(name='EWIntMaps', output_var='TaInterceptionAll', unit='mm', end=[], steps=[], all=['repEWIntMaps','repE2O1'], - restrictoption=['nonInit'], monthly=False, yearly=False), + restrictoption=['nonInit'], monthly=False, yearly=False, + scale_factor=0.00031, add_offset=10.0), 'EWRefMapsOut': ReportedMap(name='EWRefMapsOut', output_var='EWRef', unit='mm', end=[], steps=[], all=['repEWRefMaps'], - restrictoption=[], monthly=False, yearly=False), + restrictoption=[], monthly=False, yearly=False, + scale_factor=0.00031, add_offset=10.0), 'EWater': ReportedMap(name='EWater', output_var='EvaAddM3*self.var.M3toMM', unit='mm', end=[], steps=[], all=['repE2O1'], - restrictoption=['nonInit'], monthly=False, yearly=False), + restrictoption=['nonInit'], monthly=False, yearly=False, + scale_factor=0.00031, add_offset=10.0), 'FalkenmarkM3Capita1': ReportedMap(name='FalkenmarkM3Capita1', output_var='FalkenmarkM3Capita1', unit='M3Capita', end=[], @@ -426,7 +435,7 @@ output_var='SurfaceRunoff+self.var.UZOutflowPixel', unit='mm', end=[], steps=[], all=['repFastRunoffMaps','repE2O1'], restrictoption=['nonInit'], monthly=False, - yearly=False), + yearly=False, scale_factor=0.0077, add_offset=250), 'FossilGroundwaterUsedM3': ReportedMap(name='FossilGroundwaterUsedM3', output_var='FossilGroundwaterUsedM3', unit='m3', end=[], @@ -460,49 +469,49 @@ end=[], steps=['repGwPercUZLZMaps'], all=['repGwPercUZLZMaps'], restrictoption=['nonInit'], monthly=False, - yearly=False), + yearly=False, scale_factor=4.6e-05, add_offset=1.5), 'GwPercUZLZIrrigationMaps': ReportedMap(name='GwPercUZLZForestMaps', output_var='GwPercUZLZ[2]', unit='mm', end=[], steps=['repGwPercUZLZMaps'], all=['repGwPercUZLZMaps'], restrictoption=['nonInit'], monthly=False, - yearly=False), + yearly=False, scale_factor=4.6e-05, add_offset=1.5), 'GwPercUZLZOtherMaps': ReportedMap(name='GwPercUZLZForestMaps', output_var='GwPercUZLZ[0]', unit='mm', end=[], steps=[], all=['repGwPercUZLZMaps'], restrictoption=['nonInit'], monthly=False, - yearly=False), + yearly=False, scale_factor=4.6e-05, add_offset=1.5), 'GwPercUZLZMaps': ReportedMap(name='GwPercUZLZMaps', output_var='GwPercUZLZPixel', unit='mm', end=[], steps=[], all=['repGwPercUZLZMaps','repE2O2'], restrictoption=['nonInit'], monthly=False, - yearly=False), + yearly=False, scale_factor=4.6e-05, add_offset=1.5), 'InfiltrationForestMaps': ReportedMap(name='InfiltrationForestMaps', output_var='Infiltration[1]', unit='mm', end=[], steps=[], all=['repInfiltrationMaps'], restrictoption=['nonInit'], monthly=False, - yearly=False), + yearly=False, scale_factor=0.0077, add_offset=250.), 'InfiltrationMaps': ReportedMap(name='InfiltrationMaps', output_var='Infiltration[0]', unit='mm', end=[], steps=[], all=['repInfiltrationMaps'], restrictoption=['nonInit'], monthly=False, - yearly=False), + yearly=False, scale_factor=0.0077, add_offset=250.), 'InterceptionForestMaps': ReportedMap(name='InterceptionForestMaps', output_var='Interception[1]', unit='mm', end=[], steps=[], all=['repInterceptionMaps'], restrictoption=['nonInit'], monthly=False, - yearly=False), + yearly=False, scale_factor=0.00023, add_offset=7.5), 'InterceptionMaps': ReportedMap(name='InterceptionMaps', output_var='Interception[0]', unit='mm', end=[], steps=[], all=['repInterceptionMaps'], restrictoption=['nonInit'], monthly=False, - yearly=False), + yearly=False, scale_factor=0.00023, add_offset=7.5), 'LocalFreshwaterM3': ReportedMap(name='LocalFreshwaterM3', output_var='LocalFreshwaterM3', unit='m3', end=[], @@ -530,15 +539,16 @@ restrictoption=[], monthly=False, yearly=False), 'LZMaps': ReportedMap(name='LZMaps', output_var='LZ', unit='mm', end=[], steps=[], all=['repLZMaps', 'repE2O2'], restrictoption=['nonInit'], - monthly=False, yearly=False), + monthly=False, yearly=False,scale_factor=0.23, add_offset=7500.0), 'LZOutflowMaps': ReportedMap(name='LZOutflowMaps', output_var='LZOutflowToChannelPixel', unit='mm', end=[], steps=[], all=['repLZOutflowMaps','repE2O1'], restrictoption=['nonInit'], monthly=False, - yearly=False), + yearly=False, scale_factor=4.6e-05, add_offset=1.5), 'LZState': ReportedMap(name='LZState', output_var='LZ', unit='mm', end=[], steps=['repStateMaps'], all=[], - restrictoption=['nonInit'], monthly=False, yearly=False), + restrictoption=['nonInit'], monthly=False, yearly=False, + scale_factor=0.092, add_offset=-1000.0), 'LakeAbstractionM3': ReportedMap(name='LakeAbstractionM3', output_var='LakeAbstractionM3', unit='m3', end=[], steps=['repTotalAbs'], all=[], @@ -697,34 +707,34 @@ output_var='SeepTopToSubA[0]', unit='mm', end=[], steps=[], all=['repPercolationMaps'], restrictoption=['nonInit'], monthly=False, - yearly=False), + yearly=False, scale_factor=4.6e-05, add_offset=1.5), 'Percolation1ato1bForestMaps': ReportedMap(name='Percolation1ato1bForestMaps', output_var='SeepTopToSubA[1]', unit='mm', end=[], steps=[], all=['repPercolationMaps'], restrictoption=['nonInit'], monthly=False, - yearly=False), + yearly=False, scale_factor=4.6e-05, add_offset=1.5), 'Percolation1ato1bIrrigationMaps': ReportedMap(name='Percolation1ato1bIrrigationMaps', output_var='SeepTopToSubA[2]', unit='mm', end=[], steps=[], all=['repPercolationMaps'], restrictoption=['nonInit'], monthly=False, - yearly=False), + yearly=False, scale_factor=4.6e-05, add_offset=1.5), 'Percolation1bto2OtherMaps': ReportedMap(name='Percolation1bto2OtherMaps', output_var='SeepTopToSubB[0]', unit='mm', end=[], steps=[], all=['repPercolationMaps'], restrictoption=['nonInit'], monthly=False, - yearly=False), + yearly=False, scale_factor=4.6e-05, add_offset=1.5), 'Percolation1bto2ForestMaps': ReportedMap(name='Percolation1bto2ForestMaps', output_var='SeepTopToSubB[1]', unit='mm', end=[], steps=[], all=['repPercolationMaps'], restrictoption=['nonInit'], monthly=False, - yearly=False), + yearly=False, scale_factor=4.6e-05, add_offset=1.5), 'Percolation1bto2IrrigationMaps': ReportedMap(name='Percolation1bto2IrrigationMaps', output_var='SeepTopToSubB[2]', unit='mm', end=[], steps=[], all=['repPercolationMaps'], restrictoption=['nonInit'], monthly=False, - yearly=False), + yearly=False, scale_factor=4.6e-05, add_offset=1.5), 'PrevCmMCTEnd': ReportedMap(name='PrevCmMCTEnd', output_var='PrevCm0', unit='-', end=['repEndMaps'], steps=[], all=[], restrictoption=['nonInit','MCTRouting'], monthly=False, yearly=False), @@ -759,7 +769,7 @@ unit='mm/timestep', end=[], steps=[], all=['repPrecipitationMaps'], restrictoption=[], monthly=False, - yearly=False), + yearly=False, scale_factor=0.016, add_offset=500.), 'PrefFlowOtherMaps': ReportedMap(name='PrefFlowOtherMaps', output_var='PrefFlow[0]', unit='mm', end=[], steps=[], all=['repPrefFlowMaps'], @@ -781,7 +791,8 @@ yearly=False), 'RainMaps': ReportedMap(name='RainMaps', output_var='Rain', unit='mm', end=[], steps=[], all=['repRainMaps','repE2O1'], - restrictoption=['nonInit'], monthly=False, yearly=False), + restrictoption=['nonInit'], monthly=False, yearly=False, + scale_factor=0.016, add_offset=500.), 'RegionMonthExternalInflowM3': ReportedMap(name='RegionMonthExternalInflowM3', output_var='RegionMonthExternalInflowM3', unit='M3', end=[], @@ -890,39 +901,39 @@ end=[], steps=[], all=['repSeepSubToGWMaps'], restrictoption=['nonInit'], monthly=False, - yearly=False), + yearly=False, scale_factor=0.00092, add_offset=30.), 'SeepSubToGWforestMaps': ReportedMap(name='SeepSubToGWForestMaps', output_var='SeepSubToGW[1]', unit='mm', end=[], steps=[], all=['repSeepSubToGWMaps'], restrictoption=['nonInit'], monthly=False, - yearly=False), + yearly=False, scale_factor=0.00092, add_offset=30.), 'SeepSubToGWirrigationMaps': ReportedMap(name='SeepSubToGWIrrigationMaps', output_var='SeepSubToGW[2]', unit='mm', end=[], steps=[], all=['repSeepSubToGWMaps'], restrictoption=['nonInit'], monthly=False, - yearly=False), + yearly=False, scale_factor=0.00092, add_offset=30.), 'SeepSubToGWMaps': ReportedMap(name='SeepSubToGWMaps', output_var='SeepSubToGWPixel', unit='mm', end=[], steps=[], all=['repSeepSubToGWMaps'], restrictoption=['nonInit'], monthly=False, - yearly=False), + yearly=False, scale_factor=0.00092, add_offset=30.), 'SeepTopToSubBAverageOtherMap': ReportedMap(name='SeepTopToSubBAverageOtherMap', output_var='SeepTopToSubBAv[0]', unit='mm/timestep', end=['InitLisflood'], steps=[], all=[], restrictoption=[], monthly=False, - yearly=False), + yearly=False, scale_factor=0.00092, add_offset=30.), 'SeepTopToSubBAverageForestMap': ReportedMap(name='SeepTopToSubBAverageForestMap', output_var='SeepTopToSubBAv[1]', unit='mm/timestep', end=['InitLisflood'], steps=[], all=[], restrictoption=[], monthly=False, - yearly=False), + yearly=False, scale_factor=0.00092, add_offset=30.), 'SeepTopToSubBAverageIrrigationMap': ReportedMap(name='SeepTopToSubBAverageIrrigationMap', output_var='SeepTopToSubBAv[2]', unit='mm/timestep', end=['InitLisflood'], steps=[], all=[], restrictoption=[], monthly=False, - yearly=False), + yearly=False, scale_factor=0.00092, add_offset=30.), 'SnowCoverAEnd': ReportedMap(name='SnowCoverAEnd', output_var='SnowCoverS[0]', unit='mm', end=['repEndMaps'], steps=[], all=[], restrictoption=[], monthly=False, @@ -931,7 +942,7 @@ output_var='SnowCoverS[0]', unit='mm', end=[], steps=['repStateMaps'], all=[], restrictoption=['nonInit'], monthly=False, - yearly=False), + yearly=False, scale_factor=0.16, add_offset=5000.), 'SnowCoverBEnd': ReportedMap(name='SnowCoverBEnd', output_var='SnowCoverS[1]', unit='mm', end=['repEndMaps'], steps=[], all=[], restrictoption=[], monthly=False, @@ -940,7 +951,7 @@ output_var='SnowCoverS[1]', unit='mm', end=[], steps=['repStateMaps'], all=[], restrictoption=['nonInit'], monthly=False, - yearly=False), + yearly=False, scale_factor=0.16, add_offset=5000.), 'SnowCoverCEnd': ReportedMap(name='SnowCoverCEnd', output_var='SnowCoverS[2]', unit='mm', end=['repEndMaps'], steps=[], all=[], restrictoption=[], monthly=False, @@ -949,39 +960,40 @@ output_var='SnowCoverS[2]', unit='mm', end=[], steps=['repStateMaps'], all=[], restrictoption=['nonInit'], monthly=False, - yearly=False), + yearly=False, scale_factor=0.16, add_offset=5000.), 'SnowCoverMaps': ReportedMap(name='SnowCoverMaps', output_var='SnowCover', unit='mm', end=[], steps=[], all=['repSnowCoverMaps','repE2O1'], restrictoption=['nonInit'], monthly=False, - yearly=False), + yearly=False, scale_factor=0.16, add_offset=5000.), 'SnowMaps': ReportedMap(name='SnowMaps', output_var='Snow', unit='mm', end=[], steps=[], all=['repSnowMaps'], - restrictoption=['nonInit'], monthly=False, yearly=False), + restrictoption=['nonInit'], monthly=False, yearly=False, + scale_factor=0.0031, add_offset=100.), 'SnowMeltMaps': ReportedMap(name='SnowMeltMaps', output_var='SnowMelt', unit='mm', end=[], steps=[], all=['repSnowMeltMaps'], restrictoption=['nonInit'], monthly=False, - yearly=False), + yearly=False, scale_factor=0.0039, add_offset=125.), 'SurfaceRunoffMaps': ReportedMap(name='SurfaceRunoffMaps', output_var='SurfaceRunoff', unit='mm', end=[], steps=[], all=['repSurfaceRunoffMaps'], restrictoption=['nonInit'], monthly=False, - yearly=False), + yearly=False, scale_factor=0.0077, add_offset=250.), 'TaOtherMaps': ReportedMap(name='TaOtherMaps', output_var='Ta[0]', unit='mm', end=[], steps=[], all=['repTaMaps'], restrictoption=['nonInit'], monthly=False, - yearly=False), + yearly=False, scale_factor=3.4e-4, add_offset=9.0), 'TaForestMaps': ReportedMap(name='TaForestMaps', output_var='Ta[1]', unit='mm', end=[], steps=[], all=['repTaMaps'], restrictoption=['nonInit'], monthly=False, - yearly=False), + yearly=False, scale_factor=3.4e-4, add_offset=9.0), 'TaIrrigationMaps': ReportedMap(name='TaIrrigationMaps', output_var='Ta[2]', unit='mm', end=[], steps=[], all=['repTaMaps'], restrictoption=['nonInit'], monthly=False, - yearly=False), + yearly=False, scale_factor=3.4e-4, add_offset=9.0), 'TaMaps': ReportedMap(name='TaMaps', output_var='TaPixel', unit='mm', end=[], steps=[], all=['repE2O1','repTaMaps'], restrictoption=['nonInit'], - monthly=False, yearly=False), + monthly=False, yearly=False, scale_factor=3.4e-4, add_offset=9.0), 'TavgMapsOut': ReportedMap(name='TavgMapsOut', output_var='Tavg', unit='degree', end=[], steps=[], all=['repTavgMaps'], restrictoption=[], monthly=False, yearly=False), @@ -996,12 +1008,12 @@ unit='', end=[], steps=[], all=['repThetaForestMaps','repThetaMaps'], restrictoption=['nonInit'], monthly=False, - yearly=False), + yearly=False, scale_factor=1.6e-05, add_offset=0.5), 'Theta1ForestState': ReportedMap(name='Theta1ForestState', output_var='Theta1a[1]', unit='', end=[], steps=['repStateMaps'], all=[], restrictoption=['nonInit'], monthly=False, - yearly=False), + yearly=False, scale_factor=1.6e-05, add_offset=0.5), 'Theta1IrrigationEnd': ReportedMap(name='Theta1IrrigationEnd', output_var='Theta1a[2]', unit='', end=['repEndMaps','InitLisflood'], steps=[], all=[], @@ -1012,20 +1024,20 @@ end=[], steps=[], all=['repThetaIrrigationMaps','repThetaMaps'], restrictoption=['nonInit'], monthly=False, - yearly=False), + yearly=False, scale_factor=1.6e-05, add_offset=0.5), 'Theta1IrrigationState': ReportedMap(name='Theta1IrrigationState', output_var='Theta1a[2]', unit='', end=[], steps=['repStateMaps'], all=[], restrictoption=['nonInit'], monthly=False, - yearly=False), + yearly=False, scale_factor=1.6e-05, add_offset=0.5), 'Theta1Maps': ReportedMap(name='Theta1Maps', output_var='Theta1a[0]', unit='', end=[], steps=[], all=['repThetaMaps','repE2O2'], restrictoption=['nonInit'], monthly=False, - yearly=False), + yearly=False, scale_factor=1.6e-05, add_offset=0.5), 'Theta1State': ReportedMap(name='Theta1State', output_var='Theta1a[0]', unit='', end=[], steps=['repStateMaps'], all=[], restrictoption=['nonInit'], monthly=False, - yearly=False), + yearly=False, scale_factor=1.6e-05, add_offset=0.5), 'Theta2End': ReportedMap(name='Theta2End', output_var='Theta1b[0]', unit='', end=['repEndMaps','InitLisflood'], steps=[], all=[], restrictoption=[], monthly=False, yearly=False), @@ -1037,12 +1049,12 @@ unit='', end=[], steps=[], all=['repThetaForestMaps','repThetaMaps'], restrictoption=['nonInit'], monthly=False, - yearly=False), + yearly=False, scale_factor=1.6e-05, add_offset=0.5), 'Theta2ForestState': ReportedMap(name='Theta2ForestState', output_var='Theta1b[1]', unit='', end=[], steps=['repStateMaps'], all=[], restrictoption=['nonInit'], monthly=False, - yearly=False), + yearly=False, scale_factor=1.6e-05, add_offset=0.5), 'Theta2IrrigationEnd': ReportedMap(name='Theta2IrrigationEnd', output_var='Theta1b[2]', unit='', end=['repEndMaps','InitLisflood'], steps=[], all=[], @@ -1053,20 +1065,20 @@ end=[], steps=[], all=['repThetaIrrigationMaps','repThetaMaps'], restrictoption=['nonInit'], monthly=False, - yearly=False), + yearly=False, scale_factor=1.6e-05, add_offset=0.5), 'Theta2IrrigationState': ReportedMap(name='Theta2IrrigationState', output_var='Theta1b[2]', unit='', end=[], steps=['repStateMaps'], all=[], restrictoption=['nonInit'], monthly=False, - yearly=False), + yearly=False, scale_factor=1.6e-05, add_offset=0.5), 'Theta2Maps': ReportedMap(name='Theta2Maps', output_var='Theta1b[0]', unit='', end=[], steps=[], all=['repThetaMaps','repE2O2'], restrictoption=['nonInit'], monthly=False, - yearly=False), + yearly=False, scale_factor=1.6e-05, add_offset=0.5), 'Theta2State': ReportedMap(name='Theta2State', output_var='Theta1b[0]', unit='', end=[], steps=['repStateMaps'], all=[], restrictoption=['nonInit'], monthly=False, - yearly=False), + yearly=False, scale_factor=1.6e-05, add_offset=0.5), 'Theta3End': ReportedMap(name='Theta3End', output_var='Theta2[0]', unit='', end=['repEndMaps','InitLisflood'], steps=[], all=[], restrictoption=[], monthly=False, yearly=False), @@ -1078,12 +1090,12 @@ unit='', end=[], steps=[], all=['repThetaForestMaps','repThetaMaps'], restrictoption=['nonInit'], monthly=False, - yearly=False), + yearly=False, scale_factor=1.6e-05, add_offset=0.5), 'Theta3ForestState': ReportedMap(name='Theta3ForestState', output_var='Theta2[1]', unit='', end=[], steps=['repStateMaps'], all=[], restrictoption=['nonInit'], monthly=False, - yearly=False), + yearly=False, scale_factor=1.6e-05, add_offset=0.5), 'Theta3IrrigationEnd': ReportedMap(name='Theta3IrrigationEnd', output_var='Theta2[2]', unit='', end=['repEndMaps','InitLisflood'], steps=[], all=[], @@ -1093,20 +1105,20 @@ output_var='Theta2[2]', unit='', end=[], steps=[], all=['repThetaIrrigationMaps','repThetaMaps'], restrictoption=['nonInit'], monthly=False, - yearly=False), + yearly=False, scale_factor=1.6e-05, add_offset=0.5), 'Theta3IrrigationState': ReportedMap(name='Theta3IrrigationState', output_var='Theta2[2]', unit='', end=[], steps=['repStateMaps'], all=[], restrictoption=['nonInit'], monthly=False, - yearly=False), + yearly=False, scale_factor=1.6e-05, add_offset=0.5), 'Theta3Maps': ReportedMap(name='Theta3Maps', output_var='Theta2[0]', unit='', end=[], steps=[], all=['repThetaMaps','repE2O2'], restrictoption=['nonInit'], monthly=False, - yearly=False), + yearly=False, scale_factor=1.6e-05, add_offset=0.5), 'Theta3State': ReportedMap(name='Theta3State', output_var='Theta2[0]', unit='', end=[], steps=['repStateMaps'], all=[], restrictoption=['nonInit'], monthly=False, - yearly=False), + yearly=False, scale_factor=1.6e-05, add_offset=0.5), 'TotalAbsGroundwater': ReportedMap(name='TotalAbsGroundwater', output_var='abstraction_GW_actual_M3*self.var.M3toMM', unit='mm', end=[], steps=['repTotalAbs'], @@ -1151,17 +1163,17 @@ output_var='TotalRunoff', unit='mm', end=[], steps=[], all=['repTotalRunoffMaps'], restrictoption=['nonInit'], - monthly=False, yearly=False), + monthly=False, yearly=False, scale_factor=0.0077, add_offset=250.), 'TotalToChanMaps': ReportedMap(name='TotalToChanMaps', output_var='ToChanM3Runoff*self.var.M3toMM', unit='mm', end=[], steps=[], all=['repTotalRunoffMaps'], restrictoption=['nonInit'], monthly=False, - yearly=False), + yearly=False, scale_factor=0.0054, add_offset=175.), 'TransLossMaps': ReportedMap(name='TransLossMaps', output_var='TransLossM3*self.var.M3toMM', unit='mm', end=[], steps=[], all=['repGwLossMaps','repE2O2'], restrictoption=['nonInit','TransLoss'], monthly=False, - yearly=False), + yearly=False, scale_factor=0.0092, add_offset=300.), 'UpstreamInflowUsedM3': ReportedMap(name='UpstreamInflowUsedM3', output_var='UpstreamInflowUsedM3', unit='M3', end=[], steps=['repWaterUse'], @@ -1182,52 +1194,53 @@ 'UZForestMaps': ReportedMap(name='UZForestMaps', output_var='UZ[1]', unit='mm', end=[], steps=[], all=['repUZMaps'], restrictoption=['nonInit'], monthly=False, - yearly=False), + yearly=False, scale_factor=0.054, add_offset=1750.), 'UZForestState': ReportedMap(name='UZForestState', output_var='UZ[1]', unit='mm', end=[], steps=['repStateMaps'], all=[], restrictoption=['nonInit'], monthly=False, - yearly=False), + yearly=False, scale_factor=0.054, add_offset=1750.), 'UZIrrigationEnd': ReportedMap(name='UZIrrigationEnd', output_var='UZ[2]', unit='mm', end=['repEndMaps','InitLisflood'], steps=[], all=[], restrictoption=[], - monthly=False, yearly=False), + monthly=False, yearly=False,), 'UZIrrigationMaps': ReportedMap(name='UZIrrigationMaps', output_var='UZ[2]', unit='mm', end=[], steps=[], all=['repUZMaps'], restrictoption=['nonInit'], - monthly=False, yearly=False), + monthly=False, yearly=False, scale_factor=0.054, add_offset=1750.), 'UZIrrigationState': ReportedMap(name='UZIrrigationState', output_var='UZ[2]', unit='mm', end=[], steps=['repStateMaps'], all=[], restrictoption=['nonInit'], - monthly=False, yearly=False), + monthly=False, yearly=False, scale_factor=0.054, add_offset=1750.), 'UZMaps': ReportedMap(name='UZMaps', output_var='UZ[0]', unit='mm', end=[], steps=[], all=['repUZMaps'], restrictoption=['nonInit'], - monthly=False, yearly=False), + monthly=False, yearly=False, scale_factor=0.054, add_offset=1750.), 'UZOutflowForestMaps': ReportedMap(name='UZOutflowForestMaps', output_var='UZOutflow[1]', unit='mm', end=[], steps=[], all=['repUZOutflowMaps'], restrictoption=['nonInit'], monthly=False, - yearly=False), + yearly=False, scale_factor=0.0046, add_offset=150.), 'UZOutflowIrrigationMaps': ReportedMap(name='UZOutflowIrrigationMaps', output_var='UZOutflow[2]', unit='mm', end=[], steps=[], all=['repUZOutflowMaps'], restrictoption=['nonInit'], monthly=False, - yearly=False), + yearly=False, scale_factor=0.0046, add_offset=150.), 'UZOutflowOtherMaps': ReportedMap(name='UZOutflowOtherMaps', output_var='UZOutflow[0]', unit='mm', end=[], steps=[], all=['repUZOutflowMaps'], restrictoption=['nonInit'], monthly=False, - yearly=False), + yearly=False, scale_factor=0.0046, add_offset=150.), 'UZOutflowMaps': ReportedMap(name='UZOutflowMaps', output_var='UZOutflowPixel', unit='mm', end=[], steps=[], all=['repUZOutflowMaps'], restrictoption=['nonInit'], monthly=False, - yearly=False), + yearly=False, scale_factor=0.0046, add_offset=150.), 'UZState': ReportedMap(name='UZState', output_var='UZ[0]', unit='mm', end=[], steps=['repStateMaps'], all=[], - restrictoption=['nonInit'], monthly=False, yearly=False), + restrictoption=['nonInit'], monthly=False, yearly=False, + scale_factor=0.054, add_offset=1750.), 'WEI_Abs': ReportedMap(name='WEI_Abs', output_var='WEI_Abs', unit='-', end=[], steps=['repWIndex'], all=[], restrictoption=['nonInit', 'wateruse', 'indicator'], @@ -1258,17 +1271,17 @@ unit='m', end=[], steps=[], all=['repWaterDepthMaps'], restrictoption=['nonInit'], monthly=False, - yearly=False), + yearly=False, scale_factor=0.0039, add_offset=125.), 'WaterDepthState': ReportedMap(name='WaterDepthState', output_var='WaterDepth', unit='m', end=[], steps=[], all=[], restrictoption=['nonInit'], - monthly=False, yearly=False), + monthly=False, yearly=False, scale_factor=0.0039, add_offset=125.), 'TotalWaterStorageMaps': ReportedMap(name='TotalWaterStorageMaps', output_var='TotalWaterStorageMM', unit='mm', end=[], steps=[], all=['repTotalWaterStorageMaps'], restrictoption=['nonInit'], - monthly=False, yearly=False), + monthly=False, yearly=False, scale_factor=5.4e+02, add_offset=17497500.0), 'WaterSecurityIndex': ReportedMap(name='WaterSecurityIndex', output_var='WaterSecurityIndex', unit='-', end=[], steps=['repWIndex'], all=[], diff --git a/src/lisflood/global_modules/netcdf.py b/src/lisflood/global_modules/netcdf.py index 5d08ed19..e890bd35 100644 --- a/src/lisflood/global_modules/netcdf.py +++ b/src/lisflood/global_modules/netcdf.py @@ -439,7 +439,7 @@ def write_netcdf_header(settings, start_date, rep_steps, frequency, - ): + map_value=None): """ Writes a netcdf header without the data inside @@ -465,7 +465,11 @@ def write_netcdf_header(settings, list of reporting steps frequency: output frequency (all, monthly or yearly) - + map_value: ReportedMap namedtuple or None, optional + Output variable metadata from default_options.py. When provided and OutputPacking is + enabled, its scale_factor and add_offset attributes are used for int16 packing. + If None (default), no packing is applied. + Returns ------- object @@ -569,7 +573,25 @@ def write_netcdf_header(settings, time.units = 'minutes since %s' % start_date.strftime("%Y-%m-%d %H:%M:%S.0") nf1.variables["time"][:] = date2num(time_stamps, time.units, time.calendar) - value = nf1.createVariable(var_name, dtype, ('time', dim_lat_y, dim_lon_x), zlib=True, fill_value=-9999, chunksizes=(1, nrow, ncol)) + # value = nf1.createVariable(var_name, dtype, ('time', dim_lat_y, dim_lon_x), zlib=True, fill_value=-9999, chunksizes=(1, nrow, ncol)) + # Packing: use int16 with CF scale/offset if enabled and variable has packing metadata + packing_enabled = binding.get('OutputPacking', 'False') == 'True' + has_packing = (map_value is not None + and getattr(map_value, 'scale_factor', None) is not None + and getattr(map_value, 'add_offset', None) is not None) + if packing_enabled and has_packing: + var_dtype = 'i2' + var_fill = default_fillvals['i2'] + else: + var_dtype = dtype + var_fill = -9999 + + value = nf1.createVariable(var_name, var_dtype, ('time', dim_lat_y, dim_lon_x), + zlib=True, fill_value=var_fill, chunksizes=(1, nrow, ncol)) + if packing_enabled and has_packing: + value.scale_factor = np.float64(map_value.scale_factor) + value.add_offset = np.float64(map_value.add_offset) + else: value = nf1.createVariable(var_name, dtype, (dim_lat_y, dim_lon_x), zlib=True, fill_value=-9999) diff --git a/src/lisflood/global_modules/output.py b/src/lisflood/global_modules/output.py index e238f21e..47640371 100644 --- a/src/lisflood/global_modules/output.py +++ b/src/lisflood/global_modules/output.py @@ -15,17 +15,27 @@ """ import os +import datetime import numpy as np from pcraster import ifthen, catchmenttotal, mapmaximum import sys +import warnings from .zusatz import TimeoutputTimeseries from .add1 import decompress, valuecell, loadmap, compressArray from .netcdf import write_netcdf_header, iterOpenNetcdf, nanCheckMap, uncompress_array from .errors import LisfloodFileError, LisfloodWarning -from .settings import inttodate, CDFFlags, LisSettings +from .settings import inttodate, CDFFlags, LisSettings, MaskInfo +from netCDF4 import default_fillvals +# ------------------------------------------------------------------------ +# Packing constants for int16 CF scale/offset encoding +# ------------------------------------------------------------------------ +PACK_FILL = np.int16(default_fillvals['i2']) # -32767 +PACK_MIN = PACK_FILL + 1 # -32766 +PACK_MAX = np.iinfo(np.int16).max # 32767 + # ------------------------------------------------------------------------ # Writer classes # ------------------------------------------------------------------------ @@ -95,7 +105,7 @@ def write(self, start_date, rep_steps): if self.data is not None: nf1 = write_netcdf_header(self.settings, self.map_name, self.map_path, self.var.DtDay, self.map_key, self.map_value.output_var, self.map_value.unit, - start_date, rep_steps, self.frequency) + start_date, rep_steps, self.frequency, map_value=self.map_value) map_np = uncompress_array(self.data) @@ -151,16 +161,38 @@ def write(self, start_date, rep_steps): if self.step_range[0] == 0: nf1 = write_netcdf_header(self.settings, self.map_name, self.map_path, self.var.DtDay, self.map_key, self.map_value.output_var, self.map_value.unit, - start_date, rep_steps, self.frequency) + start_date, rep_steps, self.frequency, map_value=self.map_value) else: nf1 = iterOpenNetcdf(self.map_path, "", 'a', format='NETCDF4') + nc_var = nf1.variables[self.map_name] + is_packed = nc_var.dtype == np.int16 + if is_packed: + nc_var.set_auto_maskandscale(False) + scale = nc_var.scale_factor + offset = nc_var.add_offset + nodata_mask = MaskInfo.instance().info.mask + for step, data in zip(self.step_range, self.data_steps): - nf1.variables[self.map_name][step, :, :] = uncompress_array(data) + map_np = uncompress_array(data) + if is_packed: + packed = np.round((map_np - offset) / scale).astype(np.float64) + clipped = ((packed < -32767) | (packed > 32767)) & (map_np != -9999) + if clipped.any(): + vmin = offset + scale * (-32767) + vmax = offset + scale * 32767 + warnings.warn(LisfloodWarning( + f"OutputPacking: {clipped.sum()} values in '{self.map_name}' outside " + f"packing range [{vmin:.4g}, {vmax:.4g}] and will be clipped." + )) + packed = np.clip(packed, PACK_MIN, PACK_MAX) + packed[nodata_mask] = PACK_FILL + nc_var[step, :, :] = packed.astype(np.int16) + else: + nc_var[step, :, :] = map_np nf1.close() - # clear lists for next chunk self.step_range.clear() self.data_steps.clear() else: @@ -381,6 +413,87 @@ def _start_date(self): def _rep_steps(self): return self._rep_steps_val + +class MapOutputAggregated(MapOutput): + """Handles temporal aggregation (monthly/yearly mean/sum) for a variable.""" + + def __init__(self, var, map_key, map_value, frequency, operation): + out_type = 'all' # accumulates every timestep + settings = LisSettings.instance() + binding = settings.binding + self._start_date_val = var.CalendarDayStart + self._rep_steps_val = range(binding['StepStartInt'], binding['StepEndInt'] + 1) + + self._operation = operation # 'mean' or 'sum' + self._accum_buffer = None + self._accum_count = 0 + self._write_step = 0 # own step counter for NetCDF time dimension + + # Disable int16 packing for sum aggregates — monthly/yearly sums can exceed + # the int16 range calibrated for daily values. Mean aggregates stay within + # the same value range as daily output, so packing remains valid. + if operation == 'sum': + map_value_no_pack = map_value._replace(scale_factor=None, add_offset=None) + else: + map_value_no_pack = map_value + + super().__init__(var, out_type, frequency, map_key, map_value_no_pack) + + # Force immediate write for aggregated outputs (one slice per period) + if hasattr(self, 'writer') and hasattr(self.writer, 'chunks'): + self.writer.chunks = 1 + + def _output_checkpoint(self): + """Always True — we accumulate every timestep.""" + return True + + @property + def _start_date(self): + return self._start_date_val + + @property + def _rep_steps(self): + return self._rep_steps_val + + def stage(self): + """Accumulate instead of storing instantaneous values.""" + self.step = self.var.currentTimeStep() + map_np = self.writer._extract_map() + + if self._accum_buffer is None: + self._accum_buffer = np.zeros_like(map_np) + + self._accum_buffer += map_np + self._accum_count += 1 + + def write(self): + """Write only at period boundary (month-end or year-end).""" + current_date = self.var.CalendarDate + next_date = current_date + datetime.timedelta(days=self.var.DtDay) + if self.frequency == 'monthly': + is_boundary = current_date.month != next_date.month + elif self.frequency == 'yearly': + is_boundary = current_date.year != next_date.year + else: + is_boundary = True + + if is_boundary and self._accum_buffer is not None: + # Finalize + if self._operation == 'mean': + result = self._accum_buffer / self._accum_count + else: # sum + result = self._accum_buffer + + # Stage the aggregated result into the writer + self.writer.data_steps.append(result) + self.writer.step_range.append(self._write_step) + self.writer.write(self._start_date, self._rep_steps) + + # Increment own step counter and reset accumulator + self._write_step += 1 + self._accum_buffer = None + self._accum_count = 0 + # ------------------------------------------------------------------------ # Output factory # ------------------------------------------------------------------------ @@ -428,11 +541,38 @@ def __init__(self, var): if out.is_valid(): outputs.append(out) - check_duplicates = [] + # --- Temporal aggregation outputs --- + binding = settings.binding + aggregation_configs = { + 'OutputMonthlyMean': ('monthly', 'mean'), + 'OutputMonthlySum': ('monthly', 'sum'), + 'OutputYearlyMean': ('yearly', 'mean'), + 'OutputYearlySum': ('yearly', 'sum'), + } + aggregated_vars = set() # track which vars are aggregated + + reportedmaps = settings.options['reportedmaps'] + for setting_key, (frequency, operation) in aggregation_configs.items(): + var_list = binding.get(setting_key, '').split(';') + for var_name in var_list: + var_name = var_name.strip() + if var_name and var_name in reportedmaps: + map_value = reportedmaps[var_name] + out = MapOutputAggregated(var, var_name, map_value, frequency, operation) + if out.is_valid(): + outputs.append(out) + aggregated_vars.add(var_name) + + # Remove normal outputs for variables that are now aggregated outputs_clean = [] + check_duplicates = [] for out in outputs: + # Skip normal Maps/All output if variable is aggregated + if hasattr(out, 'map_key') and out.map_key in aggregated_vars: + if not isinstance(out, (MapOutputEnd, MapOutputAggregated)): + continue if out.map_path in check_duplicates: - print(f'Warning! Output map {out.map_path} is duplicated, check list of outputs') + print(f'Warning! Output map {out.map_path} is duplicated') else: check_duplicates.append(out.map_path) outputs_clean.append(out) diff --git a/src/lisflood/global_modules/settings.py b/src/lisflood/global_modules/settings.py index d60a93bd..8285407f 100755 --- a/src/lisflood/global_modules/settings.py +++ b/src/lisflood/global_modules/settings.py @@ -841,6 +841,14 @@ def __init__(self, model): msg += "\t[X] The simulation output as specified in the settings file can be found in {}\n".format(out_dir) msg += "\t[X] Activated modules: {}\n".format(activated_options) msg += "\t[X] Report options: {}\n".format(activated_rep) + # Packing and aggregation info + binding = settings.binding + if binding.get('OutputPacking', 'False') == 'True': + msg += "\t[X] Output Packing: int16 scale/offset enabled\n" + for agg_key in ['OutputMonthlyMean', 'OutputMonthlySum', 'OutputYearlyMean', 'OutputYearlySum']: + agg_val = binding.get(agg_key, '').strip() + if agg_val: + msg += "\t[X] {}: {}\n".format(agg_key, agg_val) self._msg = '{}{}'.format(header, msg) def __str__(self): diff --git a/src/lisfloodSettings_reference.xml b/src/lisfloodSettings_reference.xml index 57097349..0746f45e 100644 --- a/src/lisfloodSettings_reference.xml +++ b/src/lisfloodSettings_reference.xml @@ -47,6 +47,7 @@ You can use builtin path variables in this template and reference to other paths + @@ -210,6 +211,17 @@ You can use builtin path variables in this template and reference to other paths + + + 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). + + + ************************************************************** PARALLELISATION WITH NUMBA (USED IN ROUTING AND SOILLOOP) @@ -343,6 +355,33 @@ You can use builtin path variables in this template and reference to other paths + + + Semicolon-separated list of variable names to output as monthly mean. + Variables listed here are excluded from normal ReportSteps output. + Example: "DischargeMaps;Theta1Maps;UZMaps" + Leave empty to disable. + + + + + + See OutputMonthlyMean + + + + + + See OutputMonthlyMean + + + + + + See OutputMonthlyMean + + + @@ -1723,6 +1762,13 @@ You can use builtin path variables in this template and reference to other paths + + + location of calibration points + OPTIONAL: nominal map with locations of calibration points + + + OPTIONAL: observed or simulated input hydrographs as time series [cu m / s] @@ -1831,6 +1877,7 @@ You can use builtin path variables in this template and reference to other paths + @@ -1906,7 +1953,34 @@ You can use builtin path variables in this template and reference to other paths Number of days used for internal spin-up (fluxes computations during prerun) - + + + + Semicolon-separated list of variable names to output as monthly mean. + Variables listed here are excluded from normal ReportSteps output. + Example: "DischargeMaps;Theta1Maps;UZMaps" + Leave empty to disable. + + + + + + See OutputMonthlyMean + + + + + + See OutputMonthlyMean + + + + + + See OutputMonthlyMean + + + @@ -5152,8 +5226,8 @@ You can use builtin path variables in this template and reference to other paths ************************************************************** TABLES WITH TOPSOIL SOIL PHYSICAL PARAMETERS (HYPRES) - Each parameter is defined for upper (1a and 1b) and - lower (2) soil layers + Each parameter is defined for upper (1 and 2) and + lower (3) soil layers ************************************************************** @@ -5257,7 +5331,7 @@ You can use builtin path variables in this template and reference to other paths ************************************************************** TABLES WITH TOPSOIL SOIL PHYSICAL PARAMETERS (HYPRES) for Forest - Each parameter is defined for 1a and 1b + Each parameter is defined for 1 and 2 Normal parameter is taken for the lower soil layer ************************************************************** @@ -5651,6 +5725,13 @@ You can use builtin path variables in this template and reference to other paths + + + location of calibration points + OPTIONAL: nominal map with locations of calibration points + + + Observed or simulated input hydrographs as diff --git a/tests/test_scale_offset_aggregation.py b/tests/test_scale_offset_aggregation.py new file mode 100644 index 00000000..898b87c0 --- /dev/null +++ b/tests/test_scale_offset_aggregation.py @@ -0,0 +1,574 @@ +""" +Tests for the scale/offset packing and temporal aggregation features. + +These tests verify that: +1. Output maps produced with int16 scale/offset packing are equivalent to + float outputs within the expected quantization tolerance. +2. Aggregated (monthly mean/sum) outputs match the result of aggregating + the daily time-step outputs after the fact. +""" +from __future__ import absolute_import + +import os +import shutil +import datetime + +import numpy as np +import pytest +from netCDF4 import Dataset + +from lisflood.main import lisfloodexe +from lisflood.global_modules.default_options import default_options +from lisflood.global_modules.output import PACK_MIN, PACK_MAX + + +# --------------------------------------------------------------------------- +# Helpers +# --------------------------------------------------------------------------- + +def _find_data_var(ds): + """Find the main data variable in a netCDF4 Dataset (the 3D variable).""" + for name in ds.variables: + var = ds.variables[name] + if len(var.dimensions) == 3: + return name + # fallback: return first variable that is not a dimension or projection + dims = set(ds.dimensions.keys()) + skip = dims | {'laea', 'lambert_azimuthal_equal_area'} + for name in ds.variables: + if name not in skip: + return name + raise ValueError(f"No data variable found in {ds.filepath()}") + +def setoptions_with_new_vars(settings_file, opts_to_set=None, opts_to_unset=None, + vars_to_set=None, new_vars=None): + """Extended version of setoptions that can inject NEW textvar entries. + + Parameters + ---------- + settings_file : str + Path to XML settings template. + opts_to_set / opts_to_unset : list of str + Options to enable/disable. + vars_to_set : dict + Existing textvars to modify (name -> value). + new_vars : dict + New textvars to inject (name -> value). They are added to the first + section if not already present. + """ + import uuid + from bs4 import BeautifulSoup + from lxml import etree + from lisflood.global_modules.settings import LisSettings, Singleton + from lisflood.global_modules.errors import LisfloodError + + if isinstance(opts_to_set, str): + opts_to_set = [opts_to_set] + if isinstance(opts_to_unset, str): + opts_to_unset = [opts_to_unset] + + opts_to_set = [] if opts_to_set is None else opts_to_set + opts_to_unset = [] if opts_to_unset is None else opts_to_unset + vars_to_set = {} if vars_to_set is None else vars_to_set + new_vars = {} if new_vars is None else new_vars + + with open(settings_file) as filetocheck: + etree.parse(filetocheck) + + with open(settings_file) as tpl: + soup = BeautifulSoup(tpl, 'lxml-xml') + + # Set/unset options + lfoptions = soup.find('lfoptions') + for opt in opts_to_set: + tag = soup.find("setoption", {'name': opt}) + if tag: + tag['choice'] = '1' + elif lfoptions: + # Option not in XML — inject it + new_tag = soup.new_tag("setoption", attrs={'name': opt, 'choice': '1'}) + lfoptions.append(new_tag) + for opt in opts_to_unset: + tag = soup.find("setoption", {'name': opt}) + if tag: + tag['choice'] = '0' + elif lfoptions: + new_tag = soup.new_tag("setoption", attrs={'name': opt, 'choice': '0'}) + lfoptions.append(new_tag) + + # Modify existing textvars + for textvar, value in vars_to_set.items(): + for tag in soup.find_all("textvar", {'name': textvar}): + tag['value'] = value + break + + # Inject new textvars (add to both and sections) + if new_vars: + lfuser = soup.find('lfuser') + lfbinding = soup.find('lfbinding') + for name, value in new_vars.items(): + # Add/update in lfuser + if lfuser: + existing = lfuser.find("textvar", {'name': name}) + if existing: + existing['value'] = value + else: + new_tag = soup.new_tag("textvar", attrs={'name': name, 'value': value}) + lfuser.append(new_tag) + # Add/update in lfbinding (this is where the model reads bindings from) + if lfbinding: + existing = lfbinding.find("textvar", {'name': name}) + if existing: + existing['value'] = value + else: + new_tag = soup.new_tag("textvar", attrs={'name': name, 'value': value}) + lfbinding.append(new_tag) + + # Write temporary settings file + uid = uuid.uuid4() + filename = os.path.join(os.path.dirname(settings_file), + './{}_{}.xml'.format(os.path.basename(settings_file), uid)) + with open(filename, 'w') as dest: + dest.write(soup.prettify()) + try: + Singleton._instances = {} + Singleton._current = {} + settings = LisSettings(filename) + options = settings.options + for opt in opts_to_set: + options[opt] = True + for opt in opts_to_unset: + options[opt] = False + except LisfloodError as e: + raise e + finally: + os.unlink(filename) + return settings + + +def read_nc_variable(nc_path, var_name=None): + """Read the data variable from a NetCDF file. + + Returns + ------- + data : np.ndarray + The data array (auto-unscaled by netCDF4 library). Masked/fill values become NaN. + var_name : str + Name of the variable read. + """ + with Dataset(nc_path, 'r') as ds: + if var_name is None: + var_name = _find_data_var(ds) + var = ds.variables[var_name] + data = var[:] + # Convert masked arrays to float with NaN for masked/fill values + if hasattr(data, 'filled'): + data = data.filled(np.nan).astype(np.float64) + else: + data = np.array(data, dtype=np.float64) + # Also treat legacy -9999 fill values as NaN + data[data == -9999.0] = np.nan + return data, var_name + + +def get_nc_packing_info(nc_path, var_name=None): + """Get packing metadata from a NetCDF file. + + Returns + ------- + dict with keys: dtype, scale_factor, add_offset (or None if not packed) + """ + with Dataset(nc_path, 'r') as ds: + ds.set_auto_maskandscale(False) + if var_name is None: + var_name = _find_data_var(ds) + var = ds.variables[var_name] + ncattrs = var.ncattrs() + info = { + 'dtype': var.dtype, + 'scale_factor': var.getncattr('scale_factor') if 'scale_factor' in ncattrs else None, + 'add_offset': var.getncattr('add_offset') if 'add_offset' in ncattrs else None, + } + return info + + +# --------------------------------------------------------------------------- +# Test class +# --------------------------------------------------------------------------- + +@pytest.mark.slow +class TestScaleOffsetPacking: + """Test that int16 scale/offset packing produces outputs equivalent to float.""" + + case_dir = os.path.join(os.path.dirname(__file__), 'data', 'LF_ETRS89_UseCase') + settings_file = os.path.join(case_dir, 'settings', 'full.xml') + + # Variables to test: these have scale_factor/add_offset defined + # and use simple report options available in the test catchment. + # Theta1Maps: scale_factor=1.526e-5, add_offset=0.5 (range 0-1) + # Using repThetaMaps which produces Theta1Maps + + def _run_lisflood(self, output_dir, packing=False, dt_sec='86400', + step_start='30/07/2016 06:00', step_end='05/08/2016 06:00', + opts_to_set=None, opts_to_unset=None, new_vars=None): + """Run lisflood with specified settings and return output path.""" + full_out_dir = os.path.join(self.case_dir, 'out', output_dir) + if os.path.exists(full_out_dir): + shutil.rmtree(full_out_dir, ignore_errors=True) + os.makedirs(full_out_dir, exist_ok=True) + + base_opts_set = list(opts_to_set or []) + base_opts_unset = list(opts_to_unset or []) + base_new_vars = dict(new_vars or {}) + + # Configure packing + if packing: + base_new_vars['OutputPacking'] = 'True' + else: + base_new_vars['OutputPacking'] = 'False' + + settings = setoptions_with_new_vars( + self.settings_file, + opts_to_set=base_opts_set, + opts_to_unset=base_opts_unset, + vars_to_set={ + 'StepStart': step_start, + 'StepEnd': step_end, + 'DtSec': dt_sec, + 'PathOut': full_out_dir, + }, + new_vars=base_new_vars, + ) + lisfloodexe(settings) + return full_out_dir + + def test_packing_all_variables(self): + """All variables with scale_factor/add_offset should match float output within tolerance. + + Dynamically discovers all packed variables from default_options and checks + every one that the model produces in this configuration. + """ + # Enable as many report options as possible to produce packed variables + opts = [ + 'repThetaMaps', 'repThetaForestMaps', 'repThetaIrrigationMaps', + 'repE2O2', 'repUZMaps', 'repGwPercUZLZMaps', + 'repSnowMaps', 'repSnowCoverMaps', 'repSnowMeltMaps', + 'repDischargeMaps', 'repSurfaceRunoffMaps', + 'repFastRunoffMaps', 'repInfiltrationMaps', 'repInterceptionMaps', + 'repTotalRunoffMaps', 'repWaterDepthMaps', + 'repESActMaps', 'repETActMaps', 'repETRefMaps', 'repEWIntMaps', + 'repTaMaps', 'repRainMaps', 'repPrefFlowMaps', + 'repSeepSubToGWMaps', 'repUZOutflowMaps', + ] + + # Run without packing (float64 output) + out_float = self._run_lisflood( + 'test_packing_float', + packing=False, + opts_to_set=opts, + ) + + # Run with packing (int16 output) + out_packed = self._run_lisflood( + 'test_packing_int16', + packing=True, + opts_to_set=opts, + ) + + # Build lookup: all variables that have scale_factor defined + reportedmaps = default_options['reportedmaps'] + packed_vars = { + name: rm for name, rm in reportedmaps.items() + if getattr(rm, 'scale_factor', None) is not None + } + + checked_count = 0 + + for filename in sorted(os.listdir(out_packed)): + if not filename.endswith('.nc'): + continue + + nc_float = os.path.join(out_float, filename) + nc_packed = os.path.join(out_packed, filename) + + if not os.path.exists(nc_float): + continue + + info_packed = get_nc_packing_info(nc_packed) + if info_packed['dtype'] != np.dtype('int16'): + continue # not packed, skip + + # Get tolerance from the file's own scale_factor/add_offset attributes + scale_factor = info_packed['scale_factor'] + add_offset = info_packed['add_offset'] + + data_float, _ = read_nc_variable(nc_float) + data_packed, _ = read_nc_variable(nc_packed) + + valid = ~np.isnan(data_float) & ~np.isnan(data_packed) + if not valid.any(): + continue + + # Assert no valid float values are outside the packing range + pack_min = add_offset + scale_factor * PACK_MIN + pack_max = add_offset + scale_factor * PACK_MAX + out_of_range = valid & ((data_float < pack_min) | (data_float > pack_max)) + assert not out_of_range.any(), \ + f"{filename}: {out_of_range.sum()} values outside packing range " \ + f"[{pack_min:.4g}, {pack_max:.4g}]. Scale/offset parameters need adjustment." + + # Compare packed vs float within quantization tolerance + tolerance = scale_factor # full step tolerance + max_diff = np.max(np.abs(data_float[valid] - data_packed[valid])) + assert max_diff <= tolerance, \ + f"Max difference {max_diff} exceeds tolerance {tolerance} for {filename}" + checked_count += 1 + + assert checked_count > 0, "No packed variables were checked" + + def teardown_method(self): + """Clean up output directories.""" + out_base = os.path.join(self.case_dir, 'out') + if os.path.exists(out_base): + for d in os.listdir(out_base): + if d.startswith('test_packing_'): + shutil.rmtree(os.path.join(out_base, d), ignore_errors=True) + + +@pytest.mark.slow +class TestTemporalAggregation: + """Test that temporal aggregation (monthly mean/sum) matches manual aggregation of daily outputs.""" + + case_dir = os.path.join(os.path.dirname(__file__), 'data', 'LF_ETRS89_UseCase') + settings_file = os.path.join(case_dir, 'settings', 'full.xml') + + def _run_lisflood(self, output_dir, dt_sec='86400', + step_start='02/07/2016 06:00', step_end='01/09/2016 06:00', + opts_to_set=None, opts_to_unset=None, new_vars=None): + """Run lisflood with specified settings.""" + full_out_dir = os.path.join(self.case_dir, 'out', output_dir) + if os.path.exists(full_out_dir): + shutil.rmtree(full_out_dir, ignore_errors=True) + os.makedirs(full_out_dir, exist_ok=True) + + base_opts_set = list(opts_to_set or []) + base_opts_unset = list(opts_to_unset or []) + base_new_vars = dict(new_vars or {}) + + settings = setoptions_with_new_vars( + self.settings_file, + opts_to_set=base_opts_set, + opts_to_unset=base_opts_unset, + vars_to_set={ + 'StepStart': step_start, + 'StepEnd': step_end, + 'DtSec': dt_sec, + 'PathOut': full_out_dir, + }, + new_vars=base_new_vars, + ) + lisfloodexe(settings) + return full_out_dir + + @pytest.mark.parametrize("frequency,operation", [ + ("monthly", "mean"), + ("monthly", "sum"), + ("yearly", "mean"), + ("yearly", "sum"), + ]) + def test_temporal_aggregation(self, frequency, operation): + """Temporal aggregation should match manual aggregation of daily outputs. + + Dynamically discovers all variables with scale_factor from default_options, + enables their report options, and verifies that the aggregated output + matches the manually computed aggregation from daily outputs. + + The aggregation writes at period boundaries (month-end or year-end), + so only complete periods produce output. + """ + # Choose time range: monthly needs ~2 months, yearly needs >1 year + if frequency == 'monthly': + step_start = '02/07/2016 06:00' + step_end = '01/09/2016 06:00' + else: # yearly + step_start = '02/01/2016 06:00' + step_end = '02/01/2017 06:00' + # Build list of all report options that produce packed variables + reportedmaps = default_options['reportedmaps'] + packed_vars = { + name: rm for name, rm in reportedmaps.items() + if getattr(rm, "scale_factor", None) is not None and rm.all and not rm.steps + } + # Collect all report options needed + all_report_opts = set() + for rm in packed_vars.values(): + all_report_opts.update(rm.all) + all_report_opts.discard('') + opts = list(all_report_opts) + + # Build semicolon-separated list of all packed variable names for aggregation + agg_var_list = ';'.join(packed_vars.keys()) + + # Select the aggregation setting key + agg_setting = f'Output{frequency.capitalize()}{operation.capitalize()}' + + # Run with daily all-steps output (no aggregation) + out_daily = self._run_lisflood( + f'test_agg_{frequency}_{operation}_daily', + step_start=step_start, + step_end=step_end, + opts_to_set=opts, + new_vars={ + 'OutputPacking': 'False', + 'OutputMonthlyMean': '', + 'OutputMonthlySum': '', + 'OutputYearlyMean': '', + 'OutputYearlySum': '', + }, + ) + + # Run with temporal aggregation for all packed variables + out_agg = self._run_lisflood( + f'test_agg_{frequency}_{operation}_agg', + step_start=step_start, + step_end=step_end, + opts_to_set=opts, + new_vars={ + 'OutputPacking': 'False', + 'OutputMonthlyMean': agg_var_list if agg_setting == 'OutputMonthlyMean' else '', + 'OutputMonthlySum': agg_var_list if agg_setting == 'OutputMonthlySum' else '', + 'OutputYearlyMean': agg_var_list if agg_setting == 'OutputYearlyMean' else '', + 'OutputYearlySum': agg_var_list if agg_setting == 'OutputYearlySum' else '', + }, + ) + + # Aggregation function for manual computation + agg_func = np.nanmean if operation == 'mean' else np.nansum + + checked_count = 0 + for filename in sorted(os.listdir(out_agg)): + if not filename.endswith('.nc'): + continue + + nc_daily_path = os.path.join(out_daily, filename) + nc_agg_path = os.path.join(out_agg, filename) + + if not os.path.exists(nc_daily_path): + continue + + # Load daily data + with Dataset(nc_daily_path, 'r') as ds: + time_var = ds.variables['time'] + times = time_var[:] + time_units = time_var.units + calendar = time_var.calendar + + from netCDF4 import num2date + dates = num2date(times, time_units, calendar) + + var_name = _find_data_var(ds) + daily_data = ds.variables[var_name][:] + + # Convert to float with NaN for masked values + if hasattr(daily_data, 'filled'): + daily_data = daily_data.filled(np.nan).astype(np.float64) + daily_data[daily_data == -9999.0] = np.nan + + # Group daily data by period (month or year) + dt_day = 1 + period_groups = {} + for i, d in enumerate(dates): + if frequency == 'monthly': + period_key = (d.year, d.month) + else: + period_key = d.year + if period_key not in period_groups: + period_groups[period_key] = [] + period_groups[period_key].append(daily_data[i]) + + # Only keep completed periods (boundary reached) + completed_periods = set() + for i, d in enumerate(dates): + next_d = d + datetime.timedelta(days=dt_day) + if frequency == 'monthly': + if d.month != next_d.month: + completed_periods.add((d.year, d.month)) + else: + if d.year != next_d.year: + completed_periods.add(d.year) + + # Compute aggregation for each completed period + manual_agg = [] + for key in sorted(completed_periods): + stack = np.array(period_groups[key]) + manual_agg.append(agg_func(stack, axis=0)) + manual_agg = np.array(manual_agg) + + # Read aggregated output + data_agg, _ = read_nc_variable(nc_agg_path) + + # Only compare if the aggregated output has fewer time steps than daily + # (confirms aggregation actually happened for this variable) + if data_agg.shape[0] >= daily_data.shape[0]: + continue + + if data_agg.shape != manual_agg.shape: + continue + + valid = ~np.isnan(manual_agg) & ~np.isnan(data_agg) + if not valid.any(): + continue + + # Compare: absolute tolerance for mean, relative for sum + if operation == 'mean': + max_diff = np.max(np.abs(data_agg[valid] - manual_agg[valid])) + assert max_diff < 1e-6, \ + f"{filename}: {frequency} {operation} max difference {max_diff} exceeds tolerance 1e-6" + else: + max_rel_diff = np.max(np.abs( + (data_agg[valid] - manual_agg[valid]) / + np.where(manual_agg[valid] != 0, manual_agg[valid], 1.0) + )) + assert max_rel_diff < 1e-5, \ + f"{filename}: {frequency} {operation} max relative difference {max_rel_diff} exceeds tolerance 1e-5" + checked_count += 1 + + assert checked_count > 0, f"No aggregated variables were checked for {frequency} {operation}" + + def test_aggregation_sum_disables_packing(self): + """Sum-aggregated outputs should NOT use int16 packing even when OutputPacking=True. + + Monthly/yearly sums can exceed the int16 range calibrated for daily values, + so packing is disabled for sum aggregates. + """ + out_dir = self._run_lisflood( + 'test_agg_no_packing', + step_start='02/07/2016 06:00', + step_end='01/09/2016 06:00', + opts_to_set=['repThetaMaps', 'repE2O2'], + new_vars={ + 'OutputPacking': 'True', + 'OutputMonthlyMean': '', + 'OutputMonthlySum': 'Theta1Maps', + 'OutputYearlyMean': '', + 'OutputYearlySum': '', + }, + ) + + nc_path = os.path.join(out_dir, 'tha.nc') + assert os.path.exists(nc_path), f"Output not found: {nc_path}" + + info = get_nc_packing_info(nc_path) + # Sum-aggregated outputs should remain float (not packed) + assert info['dtype'] != np.dtype('int16'), \ + "Sum-aggregated output should NOT be int16-packed" + assert info['scale_factor'] is None, \ + "Sum-aggregated output should not have scale_factor attribute" + + def teardown_method(self): + """Clean up output directories.""" + out_base = os.path.join(self.case_dir, 'out') + if os.path.exists(out_base): + for d in os.listdir(out_base): + if d.startswith('test_agg_'): + shutil.rmtree(os.path.join(out_base, d), ignore_errors=True)