Calibration of subcatchments defined by multiple gauges in a catchment

Jean-Michel Perraud 2020-01-28

Calibration of subcatchments defined by multiple gauges in a catchment

Use case

This vignette demonstrates how one can calibrate a catchment using multiple gauging points available within this catchment. Instead of setting up a whole-of-catchment calibration definition, it makes sense, at least in a system where subareas above a gauge points do not have a behavior dependent on other catchment processes (meaning mostly, no managed reservoirs). SWIFT offers capabilities to calibrate such subcatchments sequentially, feeding the input flow of upstream and already calibrated subcatchments to other subcatchments, thus cutting down on the complexity and runtime of the overall catchment calibration.

Data

The sample data that comes with the package contains a model definition for the South Esk catchment, including a short subset of the climate and flow record data.

library(swift)
library(DiagrammeR)

modelId <- 'GR4J'
siteId <- 'South_Esk'
simulation <- sampleCatchmentModel(siteId=siteId, configId='catchment')
# simulation <- swapModel(simulation, 'MuskingumNonLinear', 'channel_routing')
simulation <- swapModel(simulation, 'LagAndRoute', 'channel_routing')

A visual of the catchment structure (note: may not render yet through GitHub)

DiagrammeR(getCatchmentDotGraph(simulation))

seClimate <- sampleSeries(siteId=siteId, varName='climate')
seFlows <- sampleSeries(siteId=siteId, varName='flow')

The names of the climate series is already set to the climate input identifiers of the model simulation, so setting them as inputs is easy:

playInput(simulation, seClimate)
setSimulationSpan(simulation, start(seClimate), end(seClimate))
setSimulationTimeStep(simulation, 'hourly')

Moving on to define the parameters, free or fixed. We will use (for now - may change) the package calibragem, companion to SWIFT.

configureHourlyGr4j(simulation)

We define a function creating a realistic feasible parameter space. This is not the main object of this vignette, so we do not describe in details.

createMetaParameterizer <- function(simulation, refArea=250, timeSpan=3600L) {
  timeSpan <- as.integer(timeSpan)
  parameterizer <- defineGr4jScaledParameter(refArea, timeSpan)
  
  # Let's define _S0_ and _R0_ parameters such that for each GR4J model instance, _S = S0 * x1_ and _R = R0 * x3_
  pStates <- linearParameterizer(
                      paramName=c("S0","R0"), 
                      stateName=c("S","R"), 
                      scalingVarName=c("x1","x3"),
                      minPval=c(0.0,0.0), 
                      maxPval=c(1.0,1.0), 
                      value=c(0.9,0.9), 
                      selectorType='each subarea')
  
  initParameterizer <- makeStateInitParameterizer(pStates)
  parameterizer <- concatenateParameterizers(parameterizer, initParameterizer)
  
  lagAndRouteParameterizer <- function() {
    p <- data.frame(Name = c('alpha', 'inverse_velocity'),
        Value = c(1, 1),
        Min = c(1e-3, 1e-3),
        Max = c(1e2, 1e2),
        stringsAsFactors = FALSE)
    p <- createParameterizer('Generic links', p)
    return(p)
  }

  # Lag and route has several discrete storage type modes. One way to set up the modeP
  setupStorageType <- function(simulation) {
    p <- data.frame(Name = c('storage_type'), 
        Value = 1, 
        Min = 1,
        Max = 1,
        stringsAsFactors = FALSE)
    p <- createParameterizer('Generic links', p)
    applySysConfig(p, simulation)
  }
  setupStorageType(simulation)

  # transfer reach lengths to the Lag and route model
  linkIds <- getLinkIds(simulation)
  reachLengths <- getStateValue(simulation, paste0('link.', linkIds, '.Length'))
  setStateValue(simulation, paste0('link.', linkIds, '.reach_length'), reachLengths) 

  lnrp <- lagAndRouteParameterizer()
  parameterizer <- concatenateParameterizers(parameterizer, lnrp)  
  return(parameterizer)
}

parameterizer <- createMetaParameterizer(simulation)
parameterizerAsDataFrame(parameterizer)

##               Name       Min        Max     Value
## 1           log_x4  0.000000   2.380211 0.3054223
## 2           log_x1  0.000000   3.778151 0.5066903
## 3           log_x3  0.000000   3.000000 0.3154245
## 4         asinh_x2 -3.989327   3.989327 2.6377523
## 5               R0  0.000000   1.000000 0.9000000
## 6               S0  0.000000   1.000000 0.9000000
## 7            alpha  0.001000 100.000000 1.0000000
## 8 inverse_velocity  0.001000 100.000000 1.0000000

Now, checking that a default parameter set works structurally on the simulation:

setParameterValue(parameterizer, 'asinh_x2', 0)
applySysConfig(parameterizer, simulation)
execSimulation(simulation)

We are now ready to enter the main topic of this vignette, subsetting the catchment into subcatchments for calibration purposes.

The sample gauge data flow contains identifiers that are of course distinct from the network node identifiers. We create a map between them (note - this information used to be in the NodeLink file in swiftv1), and we use these node as splitting points to derive subcatchments

gauges <- as.character(c( 92106, 592002, 18311, 93044,    25,   181))
names(gauges) <- paste0('node.', as.character( c(7,   12,   25,   30,   40,   43) ))   
splitElementIds <- names(gauges)
subCats <- splitToSubcatchments(simulation, splitElementIds)
str(subCats)

## List of 6
##  $ node.40:Formal class 'ExternalObjRef' [package "cinterop"] with 2 slots
##   .. ..@ obj :<externalptr> 
##   .. ..@ type: chr "MODEL_SIMULATION_PTR"
##  $ node.25:Formal class 'ExternalObjRef' [package "cinterop"] with 2 slots
##   .. ..@ obj :<externalptr> 
##   .. ..@ type: chr "MODEL_SIMULATION_PTR"
##  $ node.12:Formal class 'ExternalObjRef' [package "cinterop"] with 2 slots
##   .. ..@ obj :<externalptr> 
##   .. ..@ type: chr "MODEL_SIMULATION_PTR"
##  $ node.7 :Formal class 'ExternalObjRef' [package "cinterop"] with 2 slots
##   .. ..@ obj :<externalptr> 
##   .. ..@ type: chr "MODEL_SIMULATION_PTR"
##  $ node.30:Formal class 'ExternalObjRef' [package "cinterop"] with 2 slots
##   .. ..@ obj :<externalptr> 
##   .. ..@ type: chr "MODEL_SIMULATION_PTR"
##  $ node.43:Formal class 'ExternalObjRef' [package "cinterop"] with 2 slots
##   .. ..@ obj :<externalptr> 
##   .. ..@ type: chr "MODEL_SIMULATION_PTR"

The resulting list of subcatchment simulations is already ordered in an upstream to downstream order by SWIFT.

If we are to set up the first step of the sequential calibration:

elementId <- names(subCats)[1]

gaugeId <- gauges[elementId]
gaugeFlow <- seFlows[,gaugeId]
sc <- subCats[[elementId]]
applySysConfig(parameterizer,sc)
varId <- 'Catchment.StreamflowRate'
recordState(sc,varId)

DiagrammeR(getCatchmentDotGraph(sc))

Let’s view the default, uncalibrated output

obsVsCalc <- function(obs, calc, ylab="streamflow (m3/s)") {
    joki::plotTwoSeries(obs, calc, ylab=ylab, startTime = start(obs), endTime = end(obs))
}

execSimulation(sc)
obsVsCalc(gaugeFlow, getRecorded(sc, varId))

Now, setting up an objective (NSE) and optimizer:

objectiveId <- 'NSE'
objective <- createObjective(sc, varId, observation=gaugeFlow, objectiveId, start(seFlows), end(seFlows))
score <- getScore(objective,parameterizer)  

termination <- getMarginalTermination( tolerance = 1e-04, cutoffNoImprovement = 30, maxHours = 2/60) 
termination <- swift::CreateSceTerminationWila_Pkg_R('relative standard deviation', c('0.05','0.0167'))
sceParams <- getDefaultSceParameters()
params <- parameterizerAsDataFrame(parameterizer)
npars <- length(which(abs(params$Max-params$Min)>0))
sceParams <- SCEParameters(npars)
optimizer <- createSceOptimSwift(objective,terminationCriterion = termination, populationInitializer = parameterizer,SCEpars = sceParams)
calibLogger <- setCalibrationLogger(optimizer,"dummy")

optimStartTime <- lubridate::now();
calibResults <- executeOptimization(optimizer)
optimEndTime <- lubridate::now();
optimWallClock <- lubridate::as.duration(lubridate::interval(optimStartTime, optimEndTime))

optimWallClock

## [1] "11.4364669322968s"

And the resulting hydrograph follows. The NSE score is decent, but the magnitude of the peak is not well represented. We used a uniform value for the routing parameters; having a scaling based on link properties may be a line of enquiry.

sortedResults <- sortByScore(calibResults, 'NSE')
head(scoresAsDataFrame(sortedResults))

##         NSE   log_x4    log_x1    log_x3  asinh_x2        R0        S0
## 1 0.4949697 2.082811 0.8537391 0.5311710 0.2138180 0.2610210 0.3526996
## 2 0.4949311 2.089526 0.8630253 0.5096557 0.2259685 0.1997594 0.3133773
## 3 0.4947418 2.080529 0.8402662 0.5453678 0.2241863 0.2453242 0.2806822
## 4 0.4947411 2.087740 0.8549421 0.4918872 0.2107718 0.2001515 0.3279347
## 5 0.4947155 2.082976 0.8540680 0.5326692 0.2326300 0.3102163 0.3848358
## 6 0.4947014 2.099333 0.8598832 0.4789433 0.1990987 0.2159813 0.3402485
##      alpha inverse_velocity
## 1 57.43339         44.65305
## 2 50.28133         52.15916
## 3 54.01892         47.18598
## 4 63.34453         52.55784
## 5 63.39566         46.43670
## 6 57.53490         50.23801

p <- swift::getScoreAtIndex(sortedResults, 1)
p <- GetSystemConfigurationWila_R(p)
parameterizerAsDataFrame(p)

##               Name       Min        Max      Value
## 1           log_x4  0.000000   2.380211  2.0828106
## 2           log_x1  0.000000   3.778151  0.8537391
## 3           log_x3  0.000000   3.000000  0.5311710
## 4         asinh_x2 -3.989327   3.989327  0.2138180
## 5               R0  0.000000   1.000000  0.2610210
## 6               S0  0.000000   1.000000  0.3526996
## 7            alpha  0.001000 100.000000 57.4333924
## 8 inverse_velocity  0.001000 100.000000 44.6530483

applySysConfig(p, sc)
execSimulation(sc)
obsVsCalc(gaugeFlow, getRecorded(sc, varId))

We can create a subcatchment parameterizer, such that when applied to the whole of the South Esk, only the states of the subareas, links and nodes of the subcatchment are potentially affected.

sp <- subcatchmentParameterizer(p, sc)
applySysConfig(sp, simulation)
getStateValue(simulation, paste0('subarea.', 34:40, '.x2'))

## subarea.34.x2 subarea.35.x2 subarea.36.x2 subarea.37.x2 subarea.38.x2 
##     0.0000000     0.0000000     0.0000000     0.1443521     0.1443521 
## subarea.39.x2 subarea.40.x2 
##     0.1443521     0.0000000

# saIds <- getSubareaIds(simulation)

spFile <- tempfile()
SaveParameterizer_R(sp, spFile)
# as of May 2016, cannot reload this on Linux. There is an issue leading to a segfault, possibly related to having added support nan, inf as numeric values.
# as of July 2016, this works on Windows however. Try and see.
sp2 <- LoadParameterizer_R(spFile)

if(file.exists(spFile)) { file.remove(spFile) }

## [1] TRUE

Whole of catchment calibration combining point gauges

gauges <- as.character(c( 92106, 592002, 18311, 93044,    25,   181))
names(gauges) <- paste0('node.', as.character( c(7,   12,   25,   30,   40,   43) ))
calibNodes <- paste0('node.', as.character( c(7,   12) ))

elementId <- names(subCats)[1]

gaugeId <- gauges[calibNodes]
gaugeFlow <- seFlows[,gaugeId]

varId <- paste0(calibNodes, '.OutflowRate')
recordState(simulation,varId)

objectiveId <- 'NSE'
objective <- createObjective(simulation, varId[1], observation=gaugeFlow[,1], objectiveId, start(seFlows), end(seFlows))

addToCompositeObjective <- function(compositeObj, objective, weight, name) {
  obj <- objective
  if (cinterop::isExternalObjRef(objective, 'OBJECTIVE_EVALUATOR_WILA_PTR')) {
    obj <- UnwrapObjectiveEvaluatorWila_R(objective)
  }
  AddSingleObservationObjectiveEvaluator_R(compositeObj, obj, weight, name)
}

co <- CreateEmptyCompositeObjectiveEvaluator_R()
addToCompositeObjective(co, objective, 1.0, varId[1])
objective <- createObjective(simulation, varId[2], observation=gaugeFlow[,2], objectiveId, start(seFlows), end(seFlows))
addToCompositeObjective(co, objective, 1.0, varId[2])

co <- swift::WrapObjectiveEvaluatorWila_R(co, clone=TRUE)

score <- getScore(co,parameterizer) 
# scoresAsDataFrame(score)