#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