Versions Compared
compared with
Key
- This line was added.
- This line was removed.
- Formatting was changed.
Here is the material corresponding to the SOBEK 3 scripting tutorial.
Data for the 2.1x SOBEK model.
Measurements of water level at gauging station.
Data to build the SOBEK 3 model step by step.
Script corresponding to this tutorial. The folder structure assumed in this script can be different from the one in your computer. So the script should be adapted before being run.
Code Block | |||
---|---|---|---|
|
...
|
...
| |
#region STEP 1: Import libraries and set data path |
...
from datetime import datetime, time |
...
from Libraries.ChartFunctions import *
from Libraries.Sobek2Functions import *
from Libraries.NetworkFunctions import BranchObjectType
from Libraries.SobekWaterFlowFunctions import *
from Libraries.StandardFunctions import *
import Libraries.MapFunctions as mf |
...
from Libraries.SOBEKTutorialLibrary import *
rootPath = r"D:\ |
...
Workshop\data" |
...
#endregion
|
...
#region STEP 2: Set up flow model |
...
#region STEP 2, OPTION A: Import SOBEK 2.13 model (Zwolle)
|
...
model213Path = rootPath + r"\SW_max.lit\2"
hydroModel = ImportSobek2Model(model213Path + r"\NETWORK.TP",False,True,False |
...
)
AddToProject(hydroModel)
|
...
flowModel = GetModelByName(" |
...
Flow1D") |
...
# Select location for which measurements are available
|
...
dataPoint = GetComputationGridLocationByName(flowModel, "3")
dataPoint.Name = "Gauging Station"
|
...
#endregion
#region STEP 2, OPTION B: Import model from multiple source files
|
...
#region STEP 2.B.1: Set up model network |
...
network = CreateNetworkFromBranchShapeFile(rootPath + r"\network_Branches.shp") |
...
AddCrossSectionsToNetwork(rootPath + r"\network_Cross Sections.shp", rootPath + r"\CrossSectionProfiles.csv", network)
AddBridgesToNetwork(rootPath + r"\network_Bridges.shp", rootPath + r"\BridgeProfiles.csv", network)
AddBranchObjectsFromShapeFile(BranchObjectType.LateralSource, rootPath + r"\network_Lateral Sources.shp", network)
AddBranchObjectsFromShapeFile(BranchObjectType.Weir , rootPath + r"\network_Weirs.shp", network)
|
...
flowModel = WaterFlowModel1D()
flowModel.Network = network
flowModel.Name = "Flow model"
AddToProject(flowModel)
|
...
# Show model on map
view = OpenView(flowModel)
# Add satellite image layer
satLayer = mf.CreateSatelliteImageLayer()
satLayer.ExcludeFromMapExtent = True
map = view.MapView.Map
map.Layers.Add(satLayer)
map.CoordinateSystem = mf.CreateCoordinateSystem(3857) # EPSG code => WGS 84 / Pseudo-Mercator
map.ZoomToExtents()
|
...
#endregion
|
...
#region STEP 2.B.2: Add Boundary and Initial Conditions
|
...
#region Import and assign BC and Laterals
|
...
ImportBoundaryConditions(rootPath + r"\boundaryConditions.csv", flowModel)
ImportLateralData(rootPath + r"\laterals.csv", flowModel)
#endregion
|
...
#region Roughness
SetDefaultRoughness(flowModel, "Main", RoughnessType.StricklerKs, 30) |
...
sectionFloodPlain = AddNewRoughnessSection(flowModel,"FloodPlain") |
...
crs1 = GetItemByName(network.CrossSections, "prof_SW2815-SW2870_Bo")
minY, maxY = GetMinYMaxYofCrossSectionProfile(crs1)
AddCrossSectionRoughness(flowModel, crs1, minY, maxY, sectionFloodPlain)
crs2 = GetItemByName(network.CrossSections, "prof_D20060515-DP-295")
minY, maxY = GetMinYMaxYofCrossSectionProfile(crs2)
AddCrossSectionRoughness(flowModel, crs2, minY, maxY, sectionFloodPlain) |
...
#endregion |
...
|
...
# set initial conditions
SetInitialConditionType(flowModel, InitialConditionType.Depth)
flowModel.DefaultInitialDepth = 3
|
...
#endregion
|
...
#region STEP 2.B.3: Space and time discretization |
...
# create computational grid with calculation points at every 0.5 m
CreateComputationalGrid(flowModel, gridAtFixedLength = True, fixedLength = 200)
|
...
# Select location for which measurements are available
dataPoint = GetComputationGridLocationByName(flowModel, "stadsgrachten Zwolle Midden_0.000")
dataPoint.Name = "Gauging Station"
|
...
SetModelTimes(flowModel, datetime(2007, 1, 15, 1, 0, 0), datetime(2007, 1, 17, 1, 0, 0), time(0,10,0,))
|
...
#endregion
|
...
#endregion
|
...
#endregion
|
...
#region STEP 3: Run model |
...
RunModel(flowModel, showDialog=True)
|
...
#endregion
|
...
#region STEP 4: Get Water Level measurements and model results at Gauging Station |
...
measuredTimeSeries = GetTimeSeriesFromCSVFile(rootPath + "\waterLevel_at_ |
...
GaugingStation.csv" |
...
, dateTimeFormat ="%d-%b-%y %H:%M:%S") # Get timeseries at calculation Gauging Station (begin Zwollen center channels) calculationPoint = GetComputationGridLocationByName(flowModel, "Gauging Station") waterlevelResults = GetTimeSeriesFromWaterFlowModel(flowModel, calculationPoint, "Water level") |
...
#endregion
|
...
#region STEP 5: Plot model results and compare with measurement data
|
...
dataSeries = CreatePointSeries(measuredTimeSeries)
lineResults = CreateLineSeries(waterlevelResults)
|
...
chart = CreateChart([dataSeries, lineResults])
chart.Name = "Water level at center Zwolle channels"
|
...
# Show the chart
view = OpenView(chart) |
...
#endregion |
...
#region STEP 6: Make chart nicer |
...
lineResults.Title = "Model results"
dataSeries.Title = "Measurements"
|
...
dataSeries.Color = Color.Red
dataSeries.Size = 4
dataSeries.LineVisible = True
dataSeries.LineColor = Color.Black
|
...
chart.BottomAxis.Title = "Time"
chart.LeftAxis.Title = "Waterlevel [m]"
chart.LeftAxis.Automatic = False
chart.LeftAxis.Maximum = 2
chart.LeftAxis.Minimum = 0
|
...
chart.Legend.Visible = True
chart.Legend.Alignment = LegendAlignment.Bottom
chart.Title = "Comparison measurements vs results"
chart.TitleVisible = True
|
...
#endregion
|
...
#region STEP 7: Main loop calibration
|
...
# Save measurements and initial model results in a list of time series to compare
timeSeriesToCompare = [measuredTimeSeries, waterlevelResults]
|
...
# Selection of parameters used for calibration and their range of variation
bcCalibration = GetBoundaryDataByName(flowModel, "Boxbergen")
rangeBC = [8, 10, 12]
roughnessMain = GetItemByName(flowModel.RoughnessSections, "Main").RoughnessNetworkCoverage
rangeRoughness = [24, 26, 28, 30, 32, 34]
|
...
# Deviation of current model results with respect to measurements
|
...
roughnessInitialValue = roughnessMain.DefaultValue
bcInitialValue = bcCalibration.Flow
warmUpTimeSteps = |
...
23 dev = GetAverageDeviation(measuredTimeSeries, waterlevelResults , startAt = warmUpTimeSteps) |
...
listCalibration = [[bcInitialValue, roughnessInitialValue, dev]] |
...
|
...
# Main calibration loop
|
...
for bcValue in rangeBC:
|
...
|
...
# Adjust boundary condition value |
...
bcCalibration.Flow = bcValue
|
...
|
...
for roughnessValue in rangeRoughness:
|
...
|
...
# Adjust default roughness value for roughness section Main
|
...
roughnessMain.DefaultValue = roughnessValue
|
...
|
...
RunModel(flowModel, showDialog=True) |
...
|
...
# Add deviation of results of current run with respect to measurements
|
...
calibResults = GetTimeSeriesFromWaterFlowModel(flowModel, calculationPoint, "Water level")
|
...
deviation = GetAverageDeviation(measuredTimeSeries, calibResults, startAt = warmUpTimeSteps) |
...
listCalibration.append([bcValue, roughnessValue, deviation]) |
...
|
...
# Add timeseries to timeSeriesToCompare list |
...
timeSeriesToCompare.append(calibResults) |
...
|
...
#endregion
|
...
#region STEP 8: show calibration results in chart
|
...
# create line series for timeseries |
...
lineFirstResults = CreateLineSeries(timeSeriesToCompare[1]) |
...
lineFirstResults.Title = "Initial model"
|
...
lineFirstResults.ShowInLegend = True
|
...
lineFirstResults.PointerVisible = False
|
...
lineFirstResults.Color = Color.LightGreen
|
...
series = [lineFirstResults]
|
...
for timeseries in timeSeriesToCompare[2:]: |
...
lineResults = CreateLineSeries(timeseries) |
...
lineResults.ShowInLegend = False
|
...
lineResults.Color = Color.Orange
|
...
lineResults.PointerVisible = False
|
...
series.append(lineResults)
|
...
lineResults.ShowInLegend = True |
...
lineResults.Title = "Calibration results" |
...
|
...
# create series for measurements
|
...
measurementsSeries = CreatePointSeries(measuredTimeSeries)
|
...
measurementsSeries.Title = "Measurements"
|
...
measurementsSeries.Color = Color.Red
|
...
measurementsSeries.Size = 4 |
...
measurementsSeries.LineVisible = True |
...
measurementsSeries.LineColor = Color.Black
|
...
series.append(measurementsSeries)
|
...
chartCalibration = CreateChart(series)
|
...
chartCalibration.Name = "Calibration using water level at center Zwolle channels" |
...
chartCalibration.Title = "Comparison measurements vs results" |
...
chartCalibration.TitleVisible = True
|
...
chartCalibration.Legend.Visible = True
|
...
chartCalibration.BottomAxis.Title = "Time"
|
...
chartCalibration.LeftAxis.Title = "Waterlevel [m]" |
...
chartCalibration.LeftAxis.Automatic = False |
...
chartCalibration.LeftAxis.Maximum = 2
|
...
chartCalibration.LeftAxis.Minimum = -0.5
|
...
# Show the chart
|
...
viewCalibration = OpenView(chartCalibration)
|
...
#endregion |
...
|
...
#region STEP 9: calculate goodness and select best fit |
...
indexBestFit = -1 |
...
bcBestFit = 1e100 |
...
roughnessBestFit = 1e100 |
...
devBestFit = 1e100
|
...
for i in range( |
...
len( |
...
listCalibration) |
...
): # case 0 is data so don't use it.
|
...
case = listCalibration[i |
...
]
|
...
if case[2] < devBestFit:
|
...
indexBestFit = i |
...
bcBestFit = case[0] |
...
roughnessBestFit = case[1] |
...
devBestFit = case[2] |
...
|
...
print "selected bc: " + str(bcBestFit) + "m3/s"
|
...
print "selected roughness " + str(roughnessBestFit) + "m^(1/3) s^-1" print indexBestFit print listCalibration |
...
lineResultsBestFit = CreateLineSeries(timeSeriesToCompare[indexBestFit+1]) |
...
lineResultsBestFit.Color = Color.Blue
|
...
lineResultsBestFit.Width = 2
|
...
lineResultsBestFit.Title = "Best fit"
|
...
lineResultsBestFit.ShowInLegend = True
|
...
lineResultsBestFit.PointerVisible = False
|
...
chartCalibration.Series.Add(lineResultsBestFit)
|
...
#endregion |
...
#region STEP 10: export results
|
...
ExportListToCsvFile(rootPath + r"\calibration.csv", listCalibration, ["BC value", "Roughness", "Average deviation"], delimiterChar = ",")
|
...
ExportListToCsvFile(rootPath + r"\bestFit.csv", timeSeriesToCompare[indexBestFit], ["Time", "Water level at Gauging station"])
|
...
ExportListToCsvFile(rootPath + r"\allTimeSeries.csv", timeSeriesToCompare[indexBestFit], ["Time", "Water level at Gauging station"]) values = [] for i in range(len(timeSeriesToCompare[1])): tempValues = [] for j in range(1,len(timeSeriesToCompare)): tempValues.append(timeSeriesToCompare[j][i][1]) values.append(tempValues) index = 0 for ts in timeSeriesToCompare: if index == 0: flag = " |
...
(initial set up) " elif index == indexBestFit: flag = " (best fit) " else: flag = " " fileName = rootPath + r"\Water Level at Gauging Station case " + str(index) + flag + ".csv" |
...
ExportListToCsvFile(fileName, ts, ["Time", "Water level at Gauging station"], ",") |
...
index += 1 |
...
|
...
chartCalibration.ExportAsImage(rootPath + r"\Calibration.png",1250,650) #endregion |
Script corresponding to the library for this tutorial. The folder structure assumed in this script can be different from the one in your computer. So the script should be adapted before being run.
Code Block | ||
---|---|---|
| ||
import csv as _csv from datetime import datetime import Libraries.ChartFunctions as _cf import Libraries.SobekWaterFlowFunctions as _wff import Libraries.StandardFunctions as _st import Libraries.NetworkFunctions as _nf import Libraries.MapFunctions as _mf from Libraries.NetworkFunctions import * from DelftTools.Hydro.Structures import BridgeFrictionType def CreateNetworkFromBranchShapeFile(branchShapeFile): """Creates a network using the provided shapefile (with line geometries).""" # create network network = _nf.HydroNetwork(Name = "Network") network.CoordinateSystem = _mf.GetShapeFileCoordinateSystem(branchShapeFile) number = 1 for feature in _mf.GetShapeFileFeatures(branchShapeFile): firstCoordinate = feature.Geometry.Coordinates[0] lastCoordinate = feature.Geometry.Coordinates[feature.Geometry.Coordinates.Length -1] # create nodes startNode = _nf.HydroNode(Name= feature.Attributes["Source"] ,Geometry = _mf.CreatePointGeometry(firstCoordinate.X, firstCoordinate.Y)) endNode = _nf.HydroNode(Name= feature.Attributes["Target"], Geometry = _mf.CreatePointGeometry(lastCoordinate.X, lastCoordinate.Y)) # create channel channel = _nf.Channel(startNode, endNode, Name = feature.Attributes["Name"]) channel.Geometry = feature.Geometry # add channel and nodes to network network.Nodes.AddRange([startNode, endNode]) network.Branches.Add(channel) number += 1 _nf.MergeNodesWithSameGeometry(network) return network def _CreateCsvLookup(profileCsvFile): lookup = {} with open(profileCsvFile) as csvfile: lines = _csv.reader(csvfile, delimiter=',') dataRows = None name = None for line in lines: if (len(line) == 0): # empty row continue if (len(line) == 1): # new definition if (name != None and dataRows != None): lookup[name] = dataRows # add definition to lookup name = line[0] dataRows = [] continue dataRows.append([float(s) for s in line]) # commit last row lookup[name] = dataRows return lookup def AddBranchObjectsFromShapeFile(branchObjectType, shapeFileName, network): for feature in _mf.GetShapeFileFeatures(shapeFileName): name = feature.Attributes ["Name"] branchName = feature.Attributes ["Branch"] chainage = feature.Attributes ["Chainage"] longName = feature.Attributes ["LongName"] branch = _st.GetItemByName(network.Branches, branchName) branchObject = _nf.CreateBranchObjectOnBranchUsingChainage(branchObjectType, name, branch, chainage) branchObject.LongName = longName if branchObjectType == _nf.BranchObjectType.Weir: branchObject.CrestWidth = feature.Attributes ["CrestWidth"] branchObject.CrestLevel = feature.Attributes ["CrestLevel"] def AddCrossSectionsToNetwork(csShapeFile, profileCsvFile, network): # get crossSection definitions as lookup (name, values) lookup = _CreateCsvLookup(profileCsvFile) for feature in _mf.GetShapeFileFeatures(csShapeFile): # get attributes name = feature.Attributes ["Name"] branchName = feature.Attributes ["Branch"] chainage = feature.Attributes ["Chainage"] crossSectionType = feature.Attributes ["CrossSectio"] thalweg = feature.Attributes ["Thalweg"] definitionName = feature.Attributes ["DefName"] # search branch branch = _st.GetItemByName(network.Branches, branchName) if (crossSectionType == "ZW"): type = _nf.BranchObjectType.CrossSectionZW elif (crossSectionType == "YZ"): type = _nf.BranchObjectType.CrossSectionYZ else : continue crossSection = _nf.CreateBranchObjectOnBranchUsingChainage(type, name, branch, chainage) _nf.SetCrossSectionProfile(crossSection, lookup[definitionName], thalweg) def AddBridgesToNetwork(csShapeFile, profileCsvFile, network): # get crossSection definitions as lookup (name, values) lookup = _CreateCsvLookup(profileCsvFile) for feature in _mf.GetShapeFileFeatures(csShapeFile): # get attributes name = feature.Attributes ["Name"] branchName = feature.Attributes ["Branch"] chainage = feature.Attributes ["Chainage"] length = feature.Attributes ["Length"] roughType = feature.Attributes ["RoughType"] roughnessValue = feature.Attributes ["Roughness"] inletLoss = feature.Attributes ["InLossCoef"] outletLoss = feature.Attributes ["OutLossCoef"] groundLayerRoughness = feature.Attributes ["GLRoughness"] pillarWidth= feature.Attributes ["PillarWidth"] shapeFactor= feature.Attributes ["ShapeFactor"] longName= feature.Attributes ["LongName"] # search branch branch = _st.GetItemByName(network.Branches, branchName) bridge = _nf.CreateBranchObjectOnBranchUsingChainage(_nf.BranchObjectType.Bridge, name, branch, chainage) bridge.Length = length bridge.InletLossCoefficient = inletLoss bridge.OutletLossCoefficient = outletLoss """brFrictionTypeDictionary = {"StricklerKs":BridgeFrictionType.StricklerKs, "Chezy":BridgeFrictionType.Chezy, "Manning":BridgeFrictionType.Manning, "StricklerKn":BridgeFrictionType.StricklerKn, "WhiteColebrook":BridgeFrictionType.WhiteColebrook}""" bridge.FrictionType = _st.ParseToEnum(roughType,BridgeFrictionType) bridge.Friction = roughnessValue bridge.GroundLayerRoughness = groundLayerRoughness bridge.PillarWidth = pillarWidth bridge.ShapeFactor = shapeFactor bridge.LongName = longName if (lookup.has_key(name)): bridge.BridgeType = _nf.BridgeType.Tabulated bridge.TabulatedCrossSectionDefinition.ZWDataTable.Clear() for item in lookup[name]: bridge.TabulatedCrossSectionDefinition.ZWDataTable.AddCrossSectionZWRow(item[0], item[1], item[2]) def ImportBoundaryConditions(path, flowModel): bcDictionary = {"H": _wff.BoundaryConditionType.WaterLevelConstant, "Q": _wff.BoundaryConditionType.FlowConstant, "None": _wff.BoundaryConditionType.NoBoundary} with open(path) as csvfile: lines = _csv.reader(csvfile,delimiter = ",") for line in lines: bcName = line[0] bcType = bcDictionary[line[1]] bcValue = float(line[2]) _wff.SetBoundaryCondition(flowModel, bcName, bcType, bcValue) def ImportLateralData(path, flowModel): latDictionary = {"FlowConstant": _wff.LateralDataType.FlowConstant} with open(path) as csvfile: lines = _csv.reader(csvfile,delimiter = ",") for line in lines: latName = line[0] latType = latDictionary[line[1]] latValue = float(line[2]) _wff.SetLateralData(flowModel, latName, latType, latValue) def GetMinYMaxYofCrossSectionProfile(crossSection): min = crossSection.Definition.YZDataTable[0].Yq max = crossSection.Definition.YZDataTable[crossSection.Definition.YZDataTable.Count-1].Yq return min, max def CreateOutputChart(flowModel, outputName, locationNames, colors): """Creates an output chart for the provided locations.""" chart = _cf.CreateChart([]) index = 0 for locationName in locationNames: location = _wff.GetComputationGridLocationByName(flowModel, locationName) timeSeries = _wff.GetTimeSeriesFromWaterFlowModel(flowModel, location, outputName) lineSeries = _cf.CreateLineSeries(timeSeries) lineSeries.Title = locationName lineSeries.Color = colors[index] lineSeries.Width = 3 lineSeries.PointerVisible = False lineSeries.Transparency = 50 chart.Series.Add(lineSeries) index += 1 chart.Name = outputName chart.TitleVisible = True chart.Title = "Compare " + outputName + " zwolle" chart.Legend.Visible = True chart.BottomAxis.Title = "Time" chart.LeftAxis.Title = outputName _st.OpenView(chart) return chart def OpenModelMapViewWithBackground(model): # Show model on map view = _st.OpenView(model) # Add satellite image layer satLayer = _mf.CreateSatelliteImageLayer() satLayer.ExcludeFromMapExtent = True map = view.MapView.Map map.Layers.Add(satLayer) map.CoordinateSystem = _mf.CreateCoordinateSystem(3857) # EPSG code => WGS 84 / Pseudo-Mercator map.ZoomToExtents() def _GetInterpolatedValue(timeLeft, valueLeft, timeRight, valueRight, timeToSearch): perc = (timeToSearch - timeLeft).total_seconds() / (timeRight - timeLeft).total_seconds() return ((valueRight - valueLeft) * perc) + valueLeft #region definition def GetAverageDeviation(ts1, ts2, startAt = 0): ts1Lookup = sorted(ts1, key=lambda tup: tup[0]) # measurements ts2Lookup = sorted(ts2, key=lambda tup: tup[0]) # calculated values ts1Min = ts1Lookup[0][0] ts1Max = ts1Lookup[-1][0] ts2Min = ts2Lookup[0][0] ts2Max = ts2Lookup[-1][0] if (ts1Min < ts2Min or ts1Max > ts2Max): raise Exception("Measurements out of range") totalSum = 0 minValueIndex = 0 for timeValueList in ts1Lookup[startAt:]: timeMeasurement = timeValueList[0] valueMeasurement = timeValueList[1] for i in range(minValueIndex,len(ts2Lookup)): timeCalc = ts2Lookup[i][0] valueCalc = ts2Lookup[i][1] if (timeCalc >= timeMeasurement): minValueIndex = i -1 timeleft = ts2Lookup[i-1][0] valueleft = ts2Lookup[i-1][1] totalSum += abs(valueMeasurement - _GetInterpolatedValue(timeleft, valueleft, timeCalc, valueCalc, timeMeasurement)) break return totalSum/ (len(ts1Lookup[startAt:])) # average deviation #endregion |