Compare 1D and 2D model output¶
This notebook show examples of how to use the Compare1D2D
utility to produce figures and tables.
Reading 1D and 2D output¶
Model output can be quite large, so we only want to read output once. In this notebook, we assume the utility has already initialised and 1D_H.csv
, 2D_H.csv
files are already generated. However, see the commented-out code on how to generate these csv files from netCDF model output.
from fm2prof import Project, utils
project = Project(fr'../../tests/test_data/compare1d2d/cases/case1/fm2prof.ini');
# To convert model output to csv files, initialise like this:
"""
plotter = utils.Compare1D2D(project=project,
path_1d='../path/to/dir/that/contains/dimr.xml',
path_2d='../path/to/his/netcdf/file' )
"""
# The minimal information you need to provide is an fm2prof config file.
plotter = utils.Compare1D2D(project=project)
╔═════╣ 2024-05-07 15:41 Starting new task: Loading configuration file ║ DEBUG T+ 0.01s Received ini file: ..\..\tests\test_data\compare1d2d\cases\case1\fm2prof.ini ║ ERROR T+ 0.03s Could not find input file: 2DMapOutput ║ ERROR T+ 0.06s Unexpected error reading input files. Check config file ║ WARNING T+ 0.11s exportmapfiles is not a known key ║ ERROR T+ 0.11s Unexpected error reading (debug) parameters. Check config file ╚═════╣ Task finished in 0.13sec ================================================================================ FM2PROF version 2.3.2 Documentation: https://deltares.github.io/Fm2Prof/ Authors: Koen Berends, Asako Fujisaki, Carles Soriano Perez, Ilia Awakimjan Contact: koen.berends@deltares.nl License: LPGL license. For more info see LICENSE.txt Copyright 2016-2020, University of Twente & Deltares ================================================================================ [input] 2DMapOutput = # Output file from FM2D model (.net file) CrossSectionLocationFile = # .csv or .txt file (, as delimiter) that contains four columns: X_coordinate,Y_coordinate,BranchName,Length,Chainage. RegionPolygonFile = # User defined polygons in json format SectionPolygonFile = # User defined polygons in json format [parameters] CaseName = # Subdirectory name under OutputDir (no space); if not specified, 'CaseXX' (01, 02..) are used MaximumPointsInProfile = 20 # Number of points which are used to generate cross-sections ConveyanceDetectionMethod = 1 # [0] method based on mean flow velocity, [1] method based on maximum flow velocity (default) AbsoluteVelocityThreshold = 0.01 # Absolute velocity threshold in m/s. Used to determine storage RelativeVelocityThreshold = 0.03 # Relative velocity threshold (percentage). Used to determine storage MinimumDepthThreshold = 0.02 # Minimum depth (m) for storage identification BedlevelCriterium = 0.05 # Ignore the lowest percentage of bed level points LakeTimesteps = 10 # Number of timesteps that are used for identifying lakes ExtrapolateStorage = True # Add storage to water level independent cross-section section SDCorrection = True # Use summerdike volume correction SDFloodplainBase = 0.5 # minimum distance between floodplain base level and crest level in meters SDTransitionHeight = 0.5 # Transition height at the summer dike (m) SDOptimisationMethod = 0 # [0] Optimise on Total volume, [1] optimise on flow volume, [2] optimise on both FrictionWeighingMethod = 0 # Options. [0] arithmetric mean, [1] Weighted average SkipMaps = 0 # number of timesteps to skip at begin of 2D output. Use to skip problems related to initial conditions. ClassificationMethod = 0 # How to classify 2D output using Section and Region polygons. Options: [0] Do not classify regions or sections. [1] Use DeltaShell MinimumTotalWidth = 0.5 # Minimum width in meters. Zero width may lead to numerical instability in 1D solvers [debug] ExportCSSData = False # If True, exports pickled (1) input data for cross-section object, (2) flow mask for conveyance/storage. ExportMapFiles = False # Export detailed map output. Only use for small models or limited number of cross-sections b/c output can be huge. CssSelection = # provide comma separated list of ints; e.g. 38, 39, 40. Leave empty to generate all cross-sections [output] OutputDirectory = output # Output files are saved in OutputDir/CaseName; if not specified, subdirectory is created in the current directory ║ DEBUG T+ 0.21s 1D netCDF file does not exist or is not provided. Input provided: None. ║ DEBUG T+ 0.23s 2D netCDF file does not exist or is not provided. Input provided: None. ║ INFO T+ 0.86s Using existing flow1d csv files ║ INFO T+ 0.87s using existing 1d-2d map ║ INFO T+ 1.23s Using existing flow2d csv files ║ INFO T+ 1.24s Using existing digitized file for 1d ║ INFO T+ 1.27s Using existing digitized file for 2d
Showing results along a route¶
The example data contains results for the Rijn river. A route along the Rhine we might inspect is the Bovenrijn (BR) to Waal (WL). Below, we set savefig=False
, so that the figure renders in the notebook. If you set savefig=True
, it will save as a png to the project
output folder.
plotter.figure_longitudinal(route=["BR", "WL"], savefig=False)
║ WARNING T+ 1.62s skipped labelling BR_863.9_R_LMW-H_Lobith-haven because too close to previous station
FigureOutput(fig=<Figure size 1200x1200 with 2 Axes>, axes=array([<Axes: title={'center': 'route: BR-WL'}, xlabel='Rivierkilometers', ylabel='Waterstand [m+NAP]'>, <Axes: xlabel='Rivierkilometers', ylabel='Verschil 1D-2D [m]'>], dtype=object), legend=<matplotlib.legend.Legend object at 0x000002A235E22B90>)
These results show that the first shown time is not really initialised very well. In fact, for this simulations it makes sense to start showing results only from January 1st. We can do this by giving the plotter a start time:
from datetime import datetime
plotter.start_time = datetime(year=2000, month=1, day=5)
plotter.figure_longitudinal(route=["BR", "WL"], savefig=False)
║ WARNING T+ 2.93s skipped labelling BR_863.9_R_LMW-H_Lobith-haven because too close to previous station
FigureOutput(fig=<Figure size 1200x1200 with 2 Axes>, axes=array([<Axes: title={'center': 'route: BR-WL'}, xlabel='Rivierkilometers', ylabel='Waterstand [m+NAP]'>, <Axes: xlabel='Rivierkilometers', ylabel='Verschil 1D-2D [m]'>], dtype=object), legend=<matplotlib.legend.Legend object at 0x000002A235D9A310>)
Combining multiple simulations in a single graph¶
To show the results of multiple simulations, it makes sense to plot a single line per simulation. Two available stats are implemented: max13
plots the average of the 13 highest values for each location, which makes sense for a Dynamic standard simulation. The last25
shows the average of the last 25 output, which makes most sense for a steady simulation.
# Loading output from two sources
project = Project(fr'../../tests/test_data/compare1d2d/cases/case1/fm2prof.ini')
project2 = Project(fr'../../tests/test_data/compare1d2d/cases/case2/fm2prof.ini')
plotter = utils.Compare1D2D(project=project,
start_time=datetime(year=2000, month=1, day=5))
plotter2 = utils.Compare1D2D(project=project2,
start_time=datetime(year=2000, month=1, day=5))
# Where setting savefig to False to get the Figure outptu
Figure = plotter.figure_longitudinal(route=['BR', "PK", "IJ"],
stat="last25",
label = "D8000",
savefig=False)
# The second potter takes the `FigureOutput` from the first plotter as argument for the `add_to_fig` parameter.
plotter2.figure_longitudinal(route=['BR', "PK", "IJ"],
stat="last25",
label = "D16000",
savefig=False,
add_to_fig=Figure)
║ INFO T+ 3.99s Starting new task: Loading configuration file ║ DEBUG T+ 4.01s Received ini file: ..\..\tests\test_data\compare1d2d\cases\case1\fm2prof.ini ║ ERROR T+ 4.03s Could not find input file: 2DMapOutput ║ ERROR T+ 4.04s Unexpected error reading input files. Check config file ║ WARNING T+ 4.11s exportmapfiles is not a known key ║ ERROR T+ 4.11s Unexpected error reading (debug) parameters. Check config file ║ INFO T+ 4.13s ║ INFO T+ 4.15s ================================================================================ ║ INFO T+ 4.16s FM2PROF version 2.3.2 ║ INFO T+ 4.18s Documentation: https://deltares.github.io/Fm2Prof/ ║ INFO T+ 4.20s Authors: Koen Berends, Asako Fujisaki, Carles Soriano Perez, Ilia Awakimjan ║ INFO T+ 4.21s Contact: koen.berends@deltares.nl ║ INFO T+ 4.23s License: LPGL license. For more info see LICENSE.txt ║ INFO T+ 4.23s Copyright 2016-2020, University of Twente & Deltares ║ INFO T+ 4.25s ================================================================================ ║ INFO T+ 4.26s ║ INFO T+ 4.27s [input] 2DMapOutput = # Output file from FM2D model (.net file) CrossSectionLocationFile = # .csv or .txt file (, as delimiter) that contains four columns: X_coordinate,Y_coordinate,BranchName,Length,Chainage. RegionPolygonFile = # User defined polygons in json format SectionPolygonFile = # User defined polygons in json format [parameters] CaseName = # Subdirectory name under OutputDir (no space); if not specified, 'CaseXX' (01, 02..) are used MaximumPointsInProfile = 20 # Number of points which are used to generate cross-sections ConveyanceDetectionMethod = 1 # [0] method based on mean flow velocity, [1] method based on maximum flow velocity (default) AbsoluteVelocityThreshold = 0.01 # Absolute velocity threshold in m/s. Used to determine storage RelativeVelocityThreshold = 0.03 # Relative velocity threshold (percentage). Used to determine storage MinimumDepthThreshold = 0.02 # Minimum depth (m) for storage identification BedlevelCriterium = 0.05 # Ignore the lowest percentage of bed level points LakeTimesteps = 10 # Number of timesteps that are used for identifying lakes ExtrapolateStorage = True # Add storage to water level independent cross-section section SDCorrection = True # Use summerdike volume correction SDFloodplainBase = 0.5 # minimum distance between floodplain base level and crest level in meters SDTransitionHeight = 0.5 # Transition height at the summer dike (m) SDOptimisationMethod = 0 # [0] Optimise on Total volume, [1] optimise on flow volume, [2] optimise on both FrictionWeighingMethod = 0 # Options. [0] arithmetric mean, [1] Weighted average SkipMaps = 0 # number of timesteps to skip at begin of 2D output. Use to skip problems related to initial conditions. ClassificationMethod = 0 # How to classify 2D output using Section and Region polygons. Options: [0] Do not classify regions or sections. [1] Use DeltaShell MinimumTotalWidth = 0.5 # Minimum width in meters. Zero width may lead to numerical instability in 1D solvers [debug] ExportCSSData = False # If True, exports pickled (1) input data for cross-section object, (2) flow mask for conveyance/storage. ExportMapFiles = False # Export detailed map output. Only use for small models or limited number of cross-sections b/c output can be huge. CssSelection = # provide comma separated list of ints; e.g. 38, 39, 40. Leave empty to generate all cross-sections [output] OutputDirectory = output # Output files are saved in OutputDir/CaseName; if not specified, subdirectory is created in the current directory ║ INFO T+ 4.28s Starting new task: Loading configuration file ║ DEBUG T+ 4.30s Received ini file: ..\..\tests\test_data\compare1d2d\cases\case2\fm2prof.ini ║ ERROR T+ 4.31s Could not find input file: 2DMapOutput ║ ERROR T+ 4.33s Unexpected error reading input files. Check config file ║ WARNING T+ 4.36s exportmapfiles is not a known key ║ ERROR T+ 4.38s Unexpected error reading (debug) parameters. Check config file ║ INFO T+ 4.38s ║ INFO T+ 4.39s ================================================================================ ║ INFO T+ 4.40s FM2PROF version 2.3.2 ║ INFO T+ 4.41s Documentation: https://deltares.github.io/Fm2Prof/ ║ INFO T+ 4.42s Authors: Koen Berends, Asako Fujisaki, Carles Soriano Perez, Ilia Awakimjan ║ INFO T+ 4.42s Contact: koen.berends@deltares.nl ║ INFO T+ 4.43s License: LPGL license. For more info see LICENSE.txt ║ INFO T+ 4.44s Copyright 2016-2020, University of Twente & Deltares ║ INFO T+ 4.45s ================================================================================ ║ INFO T+ 4.46s ║ INFO T+ 4.47s [input] 2DMapOutput = # Output file from FM2D model (.net file) CrossSectionLocationFile = # .csv or .txt file (, as delimiter) that contains four columns: X_coordinate,Y_coordinate,BranchName,Length,Chainage. RegionPolygonFile = # User defined polygons in json format SectionPolygonFile = # User defined polygons in json format [parameters] CaseName = # Subdirectory name under OutputDir (no space); if not specified, 'CaseXX' (01, 02..) are used MaximumPointsInProfile = 20 # Number of points which are used to generate cross-sections ConveyanceDetectionMethod = 1 # [0] method based on mean flow velocity, [1] method based on maximum flow velocity (default) AbsoluteVelocityThreshold = 0.01 # Absolute velocity threshold in m/s. Used to determine storage RelativeVelocityThreshold = 0.03 # Relative velocity threshold (percentage). Used to determine storage MinimumDepthThreshold = 0.02 # Minimum depth (m) for storage identification BedlevelCriterium = 0.05 # Ignore the lowest percentage of bed level points LakeTimesteps = 10 # Number of timesteps that are used for identifying lakes ExtrapolateStorage = True # Add storage to water level independent cross-section section SDCorrection = True # Use summerdike volume correction SDFloodplainBase = 0.5 # minimum distance between floodplain base level and crest level in meters SDTransitionHeight = 0.5 # Transition height at the summer dike (m) SDOptimisationMethod = 0 # [0] Optimise on Total volume, [1] optimise on flow volume, [2] optimise on both FrictionWeighingMethod = 0 # Options. [0] arithmetric mean, [1] Weighted average SkipMaps = 0 # number of timesteps to skip at begin of 2D output. Use to skip problems related to initial conditions. ClassificationMethod = 0 # How to classify 2D output using Section and Region polygons. Options: [0] Do not classify regions or sections. [1] Use DeltaShell MinimumTotalWidth = 0.5 # Minimum width in meters. Zero width may lead to numerical instability in 1D solvers [debug] ExportCSSData = False # If True, exports pickled (1) input data for cross-section object, (2) flow mask for conveyance/storage. ExportMapFiles = False # Export detailed map output. Only use for small models or limited number of cross-sections b/c output can be huge. CssSelection = # provide comma separated list of ints; e.g. 38, 39, 40. Leave empty to generate all cross-sections [output] OutputDirectory = output # Output files are saved in OutputDir/CaseName; if not specified, subdirectory is created in the current directory ║ DEBUG T+ 4.48s 1D netCDF file does not exist or is not provided. Input provided: None. ║ DEBUG T+ 4.49s 2D netCDF file does not exist or is not provided. Input provided: None. ║ INFO T+ 5.07s Using existing flow1d csv files ║ INFO T+ 5.08s using existing 1d-2d map ║ INFO T+ 5.40s Using existing flow2d csv files ║ INFO T+ 5.40s Using existing digitized file for 1d ║ INFO T+ 5.43s Using existing digitized file for 2d ║ DEBUG T+ 5.48s 1D netCDF file does not exist or is not provided. Input provided: None. ║ DEBUG T+ 5.48s 2D netCDF file does not exist or is not provided. Input provided: None. ║ INFO T+ 5.88s Using existing flow1d csv files ║ INFO T+ 5.88s using existing 1d-2d map ║ INFO T+ 6.17s Using existing flow2d csv files ║ INFO T+ 6.18s Using existing digitized file for 1d ║ INFO T+ 6.21s Using existing digitized file for 2d ║ WARNING T+ 6.65s skipped labelling BR_863.9_R_LMW-H_Lobith-haven because too close to previous station ║ WARNING T+ 6.65s skipped labelling PK_871.8_L_LMW-H_Pannerden because too close to previous station ║ WARNING T+ 6.66s skipped labelling PK_878.5_R_LMW-H_IJsselkop because too close to previous station ║ WARNING T+ 6.67s skipped labelling IJ_879.6_R_LMW-H_Hondsbroeksche-Pleij-IJssel because too close to previous station ║ WARNING T+ 6.67s skipped labelling IJ_881.1_R_LMW-H_Westervoort-IJsseldijkerwaard because too close to previous station ║ WARNING T+ 6.68s skipped labelling IJ_931.2_R_LMW-H_Eefde-beneden because too close to previous station
FigureOutput(fig=<Figure size 1200x1200 with 2 Axes>, axes=array([<Axes: title={'center': 'route: BR-PK-IJ'}, xlabel='Rivierkilometers', ylabel='Waterstand [m+NAP]'>, <Axes: xlabel='Rivierkilometers', ylabel='Verschil 1D-2D [m]'>], dtype=object), legend=<matplotlib.legend.Legend object at 0x000002A235EE17D0>)
Showing results for a single station¶
To compare results at a single station, use the following method
plotter.figure_at_station("WL_900.00", savefig=False)
If you are not sure what stations are available, use the stations
method to list all stations. The code below returns all stations that contain the abbreviation "WL"
[station for station in plotter.stations() if "WL" in station]
['WL_914.00', 'WL_915.00', 'WL_916.00', 'WL_917.00', 'WL_918.00', 'WL_919.00', 'WL_920.00', 'WL_921.00', 'WL_922.00', 'WL_923.00', 'WL_924.00', 'WL_925.00', 'WL_926.00', 'WL_927.00', 'WL_928.00', 'WL_929.00', 'WL_930.00', 'WL_931.00', 'WL_932.00', 'WL_933.00', 'WL_934.00', 'WL_935.00', 'WL_936.00', 'WL_937.00', 'WL_938.00', 'WL_939.00', 'WL_940.00', 'WL_941.00', 'WL_942.00', 'WL_943.00', 'WL_944.00', 'WL_945.00', 'WL_946.00', 'WL_947.00', 'WL_948.00', 'WL_949.00', 'WL_950.00', 'WL_951.00', 'WL_952.00', 'WL_926.1_L_LMW-H_Sint-Andries-Waal-g6', 'WL_934.8_L_LMW-H_Zaltbommel', 'WL_951.8_R_LMW-H_Vuren', 'WL_867.00', 'WL_913.3_R_LMW-H_Tiel-Waal', 'WL_868.00', 'WL_869.00', 'WL_870.00', 'WL_871.00', 'WL_872.00', 'WL_873.00', 'WL_874.00', 'WL_875.00', 'WL_876.00', 'WL_877.00', 'WL_878.00', 'WL_879.00', 'WL_880.00', 'WL_881.00', 'WL_882.00', 'WL_883.00', 'WL_884.00', 'WL_885.00', 'WL_886.00', 'WL_887.00', 'WL_888.00', 'WL_889.00', 'WL_890.00', 'WL_891.00', 'WL_892.00', 'WL_893.00', 'WL_894.00', 'WL_895.00', 'WL_896.00', 'WL_897.00', 'WL_898.00', 'WL_899.00', 'WL_900.00', 'WL_901.00', 'WL_902.00', 'WL_903.00', 'WL_904.00', 'WL_905.00', 'WL_906.00', 'WL_907.00', 'WL_908.00', 'WL_909.00', 'WL_910.00', 'WL_911.00', 'WL_912.00', 'WL_913.00', 'WL_884.9_L_LMW-H_Nijmegen-haven', 'WL_901.4_R_LMW-H_Dodewaard']
Compare discharge distribution¶
To plot discharge distribution between two stations, use the following code
plotter.figure_compare_discharge_at_stations(stations=['WL_869.00', 'PK_869.00'],
title="Pannerdensche Kop",
savefig=False)
FigureOutput(fig=<Figure size 1200x1000 with 3 Axes>, axes=array([<Axes: title={'center': 'tijdserie'}, ylabel='afvoer [m$^3$/s]'>, <Axes: title={'center': 'afvoerverdeling'}, xlabel='afvoer bovenstrooms [m$^3$/s]', ylabel='percentage t.o.v. totaal'>], dtype=object), legend=<matplotlib.legend.Legend object at 0x000002A235673410>)