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.
In [1]:
Copied!
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)
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-08-30 14:10 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.05s Unexpected error reading input files. Check config file ║ WARNING T+ 0.13s exportmapfiles is not a known key ║ ERROR T+ 0.14s Unexpected error reading (debug) parameters. Check config file ╚═════╣ Task finished in 0.15sec ================================================================================ FM2PROF version 2.3.3 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.23s 1D netCDF file does not exist or is not provided. Input provided: None. ║ DEBUG T+ 0.24s 2D netCDF file does not exist or is not provided. Input provided: None. ║ INFO T+ 1.07s Using existing flow1d csv files ║ INFO T+ 1.10s using existing 1d-2d map ║ INFO T+ 2.30s Using existing flow2d csv files ║ INFO T+ 2.31s Using existing digitized file for 1d ║ INFO T+ 2.37s 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.
In [2]:
Copied!
plotter.figure_longitudinal(route=["BR", "WL"], savefig=False)
plotter.figure_longitudinal(route=["BR", "WL"], savefig=False)
║ WARNING T+ 13.87s skipped labelling BR_863.9_R_LMW-H_Lobith-haven because too close to previous station
Out[2]:
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 0x00000247A18A7C50>)
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:
In [4]:
Copied!
from datetime import datetime
plotter.start_time = datetime(year=2000, month=1, day=5)
plotter.figure_longitudinal(route=["BR", "WL"], savefig=False)
from datetime import datetime
plotter.start_time = datetime(year=2000, month=1, day=5)
plotter.figure_longitudinal(route=["BR", "WL"], savefig=False)
║ WARNING T+ 38.98s skipped labelling BR_863.9_R_LMW-H_Lobith-haven because too close to previous station
Out[4]:
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 0x00000247AE31AE50>)
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.
In [5]:
Copied!
# 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="last3",
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="last3",
label = "D16000",
savefig=False,
add_to_fig=Figure)
# 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="last3",
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="last3",
label = "D16000",
savefig=False,
add_to_fig=Figure)
║ INFO T+ 44.24s Starting new task: Loading configuration file ║ DEBUG T+ 44.25s Received ini file: ..\..\tests\test_data\compare1d2d\cases\case1\fm2prof.ini ║ ERROR T+ 44.27s Could not find input file: 2DMapOutput ║ ERROR T+ 44.28s Unexpected error reading input files. Check config file ║ WARNING T+ 44.32s exportmapfiles is not a known key ║ ERROR T+ 44.35s Unexpected error reading (debug) parameters. Check config file ║ INFO T+ 44.37s ║ INFO T+ 44.38s ================================================================================ ║ INFO T+ 44.40s FM2PROF version 2.3.3 ║ INFO T+ 44.41s Documentation: https://deltares.github.io/Fm2Prof/ ║ INFO T+ 44.45s Authors: Koen Berends, Asako Fujisaki, Carles Soriano Perez, Ilia Awakimjan ║ INFO T+ 44.47s Contact: koen.berends@deltares.nl ║ INFO T+ 44.48s License: LPGL license. For more info see LICENSE.txt ║ INFO T+ 44.50s Copyright 2016-2020, University of Twente & Deltares ║ INFO T+ 44.51s ================================================================================ ║ INFO T+ 44.52s ║ INFO T+ 44.52s [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+ 44.55s Starting new task: Loading configuration file ║ DEBUG T+ 44.57s Received ini file: ..\..\tests\test_data\compare1d2d\cases\case2\fm2prof.ini ║ ERROR T+ 44.59s Could not find input file: 2DMapOutput ║ ERROR T+ 44.60s Unexpected error reading input files. Check config file ║ WARNING T+ 44.64s exportmapfiles is not a known key ║ ERROR T+ 44.68s Unexpected error reading (debug) parameters. Check config file ║ INFO T+ 44.69s ║ INFO T+ 44.71s ================================================================================ ║ INFO T+ 44.72s FM2PROF version 2.3.3 ║ INFO T+ 44.73s Documentation: https://deltares.github.io/Fm2Prof/ ║ INFO T+ 44.74s Authors: Koen Berends, Asako Fujisaki, Carles Soriano Perez, Ilia Awakimjan ║ INFO T+ 44.75s Contact: koen.berends@deltares.nl ║ INFO T+ 44.76s License: LPGL license. For more info see LICENSE.txt ║ INFO T+ 44.78s Copyright 2016-2020, University of Twente & Deltares ║ INFO T+ 44.80s ================================================================================ ║ INFO T+ 44.81s ║ INFO T+ 44.82s [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+ 44.83s 1D netCDF file does not exist or is not provided. Input provided: None. ║ DEBUG T+ 44.84s 2D netCDF file does not exist or is not provided. Input provided: None. ║ INFO T+ 46.07s Using existing flow1d csv files ║ INFO T+ 46.08s using existing 1d-2d map ║ INFO T+ 46.86s Using existing flow2d csv files ║ INFO T+ 46.86s Using existing digitized file for 1d ║ INFO T+ 46.89s Using existing digitized file for 2d ║ DEBUG T+ 46.93s 1D netCDF file does not exist or is not provided. Input provided: None. ║ DEBUG T+ 46.95s 2D netCDF file does not exist or is not provided. Input provided: None. ║ INFO T+ 48.12s Using existing flow1d csv files ║ INFO T+ 48.13s using existing 1d-2d map ║ INFO T+ 48.50s Using existing flow2d csv files ║ INFO T+ 48.51s Using existing digitized file for 1d ║ INFO T+ 48.55s Using existing digitized file for 2d ║ WARNING T+ 50.02s skipped labelling BR_863.9_R_LMW-H_Lobith-haven because too close to previous station ║ WARNING T+ 50.03s skipped labelling PK_871.8_L_LMW-H_Pannerden because too close to previous station ║ WARNING T+ 50.05s skipped labelling PK_878.5_R_LMW-H_IJsselkop because too close to previous station ║ WARNING T+ 50.06s skipped labelling IJ_879.6_R_LMW-H_Hondsbroeksche-Pleij-IJssel because too close to previous station ║ WARNING T+ 50.08s skipped labelling IJ_881.1_R_LMW-H_Westervoort-IJsseldijkerwaard because too close to previous station ║ WARNING T+ 50.10s skipped labelling IJ_931.2_R_LMW-H_Eefde-beneden because too close to previous station
Out[5]:
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 0x00000247AE770410>)
Showing results for a single station¶
To compare results at a single station, use the following method
In [5]:
Copied!
plotter.figure_at_station("WL_900.00", savefig=False)
plotter.figure_at_station("WL_900.00", savefig=False)
Out[5]:
FigureOutput(fig=<Figure size 1200x500 with 2 Axes>, axes=<Axes: title={'center': 'WL_900.00\nTijdreeks'}, ylabel='Waterstand [m+NAP]'>, legend=<matplotlib.legend.Legend object at 0x000002AC022D2110>)
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"
In [6]:
Copied!
[station for station in plotter.stations() if "WL" in station]
[station for station in plotter.stations() if "WL" in station]
Out[6]:
['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
In [6]:
Copied!
plotter.figure_compare_discharge_at_stations(stations=['WL_869.00', 'PK_869.00'],
title="Pannerdensche Kop",
savefig=False)
plotter.figure_compare_discharge_at_stations(stations=['WL_869.00', 'PK_869.00'],
title="Pannerdensche Kop",
savefig=False)
Out[6]:
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 0x00000247ADD02B10>)
Statistics¶
All statistics, like bias, mean average error and 'max13' can also be written to a csv file for analysis. The following method will immediately write t file. It can also take an optional argument file_path where you can change the name of the output
In [7]:
Copied!
plotter.statistics_to_file()
plotter.statistics_to_file()
║ INFO T+ 76.20s statistics written to ..\..\tests\test_data\compare1d2d\cases\case1\output\error_statistics.csv
In [ ]:
Copied!