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:\DSDNov2014\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,False) AddToProject(hydroModel) flowModel = GetModelByName("water flow 1d") # 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_3_Measured.csv") # 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 = 20 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(1,len(timeSeriesToCompare)-1): # case 0 is data so don't use it. case = listCalibration[i-1] 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" lineResultsBestFit = CreateLineSeries(timeSeriesToCompare[indexBestFit]) 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"]) 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): Gui.DocumentViewsResolver.OpenViewForData# Show model on map view = _st.OpenView(model) view # Add satellite image layer satLayer = Gui.DocumentViews[Gui.DocumentViews.Count-1]_mf.CreateSatelliteImageLayer() satLayer.ExcludeFromMapExtent = True map = view.MapView.Map map.Layers.Add(satLayer) map.CoordinateSystem = _mf.Map.CoordinateSystemFactory.CreateFromEPSGCreateCoordinateSystem(3857) # webmercator # Add background (Satellite image) map layer = BingLayer() layer.Name = "Satellite Hybrid (Bing Maps)" map.Layers.Add(layer)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 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: timeMeasurement = timeValueList[0] valueMeasurement = timeValueList[1] for i in range(len(ts2Lookup))[minValueIndex:]: 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)) # average deviation |