Page History

#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