Add reducing evaporation to the R-Instat water balance

[rdstern]

This has been srtongly requested by CIMH as an item in the original Instat that has yet to be added into R-Instat.

Here is the relevant page from the old Instat climatic guide. Notice the Simple and First in the old Instat, which is in 2 separate dialogues in R-Instat.

The Simple options from the old Instat (see above) are now all within the Climatic > Prepare > Transform dialogue.

I suggest we make the Options group box wider and put the Variable button on the same row as the Constant. Then underneath that we have a checkbox with label Reducing and default unchecked. If checked it opens an up-down with default value 0.5 . It has maximum 1 and minimum 0 and goes in steps of 0.01.

We still need the R calculation behind this, but we could start with these changes and then disable the Reducing checkbox and wait for the code.

The End of the Rains dialog is as follows:

That's a relief. We can have exactly the same changes here.

Then I need to enquire who can do the R part?

@N-thony you may wish to allocate?

[N-thony]

@MeSophie @NanyanziAlice can you provide an update on this?

[lilyclements]

@MeSophie are you available tomorrow morning for a meeting? If it is ready for the R-code part, I can then work on that tomorrow.

[MeSophie]

@MeSophie are you available tomorrow morning for a meeting? If it is ready for the R-code part, I can then work on that tomorrow.

@lilyclements I'm available tomorrow anytime. What I did soo far is just add controls in VB because I don't understand well how the dialog is working in R-Instat. I can+'t identified the package in the code so I don't even know which parameter this controls are using.

[lilyclements]

@MeSophie great, I can absolutely understand - this is a bit different to hte usual standard R code we use. It might be midday/early afternoon that I ring, we can organise on MatterMost :)

[MeSophie]

@MeSophie great, I can absolutely understand - this is a bit different to hte usual standard R code we use. It might be midday/early afternoon that I ring, we can organise on MatterMost :)

@lilylements well received. thank you

[lilyclements]

@rdstern before continuing, I want to ensure I have understood the code we want to implement here!

• Soil Holding Capacity = Capacity nud * Reducing nud
• What is Water Content at the start? I think this is 0 at the start, right?
• For each day, the Water Capacity is then the WB of the previous day + rain of today - evaporation of today.
• The value for evaporation is then what is calculated from this new "Reducing" checkbox:

I'm going to write this twice, once for you, once for me when I read this back!

• if Water Content is larger than the Soil Holding Capacity then the evaporation value is 5mm a day
• otherwise:
• if it is a rainy day and the rainfall exceeds the evaporation value given in the input box in R-Instat, then we run another if statement:
• let's set X = 5 / (water content) + rain that day. if X > 5mm, we lose 5mm (evaporation = 5), otherwise we lose X (evaporation = X)
• if it is not a rainy day or the rainfall does not exceed the evaporation value given in the input box in R-Instat, we lose 5/(water content)mm (evaporation = 5/water content).

Does this sound about right to you?

If this makes sense, then you can ignore the rest of this - it is all for my own mind when I read it back!

if Water Content > Soil Holding Capacity 
    evaporation = 5mm a day
otherwise
   if rainy day **and** the rainfall > ucrInputEvaporation
(        let's set `X = 5 / (Water Content)  + rain that day`.)
         if X > 5mm, evaporation = 5mm
         otherwise evaporation = X
   otherwise
          evaporation = (5/Water Content) mm

if this is the case then -

Reduce(f=function(x, y) pmin(pmax(x + y, 0), 100), x=tail(x=rain_max - evap, n=-1), init=0, accumulate=TRUE)

x = rain_today
y = WB_yesterday
100 = capacity
evap = evaporation amount in input

#TODO:
we add a column in for the changing evap value and that's what we read in instead of the value in the input 

[rdstern]

I am not sure if I have understood correctly. The dialog and formula is all there already, and seems to be working well. And it is working whether evaporation is a constant or is another variable.

I assume it starts at zero, but that should all be there already.

So there is nothing for you to do when the wb is full. The current formula is simply

wb = wb(before) + rf on day n, - evap on day n. if wb > capacity then wb = capacity if wb < 0 wb = 0.

You don't need to change the main formula. Just adjust the actual evap on day n in relation to the full evap on day n , before applying it.

The amendment adds fraction (default = 0.5) Perhaps make limits 0.8 and 0.2 so never need to worry about endpoints.

Then if wb (before) >= fraction capacity then evap on day n = full value If wb(before) < fraction capacity and rf on day n = 0 then evap on day n = evap normally wb(before)/ (fraction capacity).

Still have to worry about the adjustment to the evap if there is rain (and wb (n-1) is below the threshold * capacity:

a) If rain > full evap on day n, then evap is normal evap on that day. It isn't reduced. b) If rain > 0 and less than full evap on day n, then evap on day n is (say) reduced to evap (on dry day) + remaining evap * (rainfall/full evap).

Hope that's clear? e.g. Capacity 200. Threshold 0.6 - evap = 6 then evap is 6 if wb >= 120 If wb = 80 (so 2/3 of the threshold of 120, down to zero) then 1) evap = 4 on dry day (so 2/3 of the full evap as the wb is just 1/3 down from the 120 towards zero. 2) evap = 6 on rain day if rain >= 6 3) evap = 4.5 if rain = 1.5

[lilyclements]

@rdstern great, thanks. Just have one question - in your example, "evap (on dry day)" is 3. Where is this defined?

[rdstern]

@lilyclements I didn't see evap is 3 above. I have added a bit more explanation above. Hope it helps.

[lilyclements]

Ok, thanks. I'm attaching a snippet of data, does this help? It is the daily_niger data, we can check the values are the same between R-Isntat and Instat This is for: frac = 0.6 capacity = 200 evaporation_value = 6

# Code - for my own sake, just in case!
x <- daily_niger %>%
  filter(date > as.Date("1945-07-09")) %>%
  # create rain_min with NA as 0
  dplyr::mutate(rain_min = ifelse(is.na(rain), 0, rain)) %>%
  dplyr::mutate(wb_min = Reduce(f = function(x, y) pmin(pmax(y - WB_evaporation(x, frac, capacity, evaporation_value, y) + x, 0), capacity),
                                x = utils::tail(rain_min, n=-1), init = 0, accumulate=TRUE)) %>%
  dplyr::mutate(rain_max = ifelse(is.na(rain), capacity, rain)) %>%
  dplyr::mutate(wb_max = Reduce(f = function(x, y) pmin(pmax(y - WB_evaporation(x, frac, capacity, evaporation_value, y) + x, 0), capacity),
                                x = utils::tail(rain_max, n=-1), init = 0, accumulate=TRUE))%>%
  dplyr::mutate(wb = ifelse((wb_min != wb_max) | is.na(rain), NA, wb_min))

WB_evaporation <- function(water_balance, frac, capacity, evaporation_value, rain){
  if (water_balance >= frac*capacity){
    evaporation <- evaporation_value
  } else {
    if (rain == 0){
      evaporation <- evaporation_value * ((water_balance)/(frac*capacity))
    } else {
      if (water_balance < frac*capacity){
        if (rain > evaporation_value){
          evaporation <- evaporation_value
        } else {
          evaporation <- evaporation_value * ((water_balance)/(frac*capacity))
          evaporation <- evaporation + ((evaporation_value - evaporation)*(rain/evaporation_value))
        }
      } else {
        evaporation <- evaporation_value
        }
    }
  }
  return(evaporation)
}

[lilyclements]

If I calculate it manually, this is making sense to me. But there's always potential that I misunderstood!

[MeSophie]

@rdstern , @lilyclements it this code ready to implement or we need to compare the result with Instat first?

[lilyclements]

@MeSophie it's ready to implement. See #8404

The change to make is if the checkbox is checked, then let frac = the value in ucrNudWB

We change from

# unchecked
Reduce(f = function(x, y) pmin(pmax(y + x, 0), capacity),
                                x = utils::tail(rain_min - evaporation_value, n=-1), init = 0, accumulate=TRUE))

To

# if checked
Reduce(f = function(x, y) pmin(pmax(y - WB_evaporation(x, frac, capacity, evaporation_value, y) + x, 0), capacity),
                                x = utils::tail(rain_min, n=-1), init = 0, accumulate=TRUE))

Here: frac = the value in the ucrnudWB capacity = value in capacity nud evaporation_value = value in the evaporation input or variable receiver

This needs to be done for rain_min and rain_max. The example above is rain_min. For rain_max, it is exactly the same, but the variable rain_min is now rain_max.

[MeSophie]

@lilyclements Thank you.

[MeSophie]

@lilyclements Please can you explain to me clearly what remains to be done in the dialogue. I'm still confused. I saw you already added the controls what is great. Since the work is not in stand_alone_functions.R, the code in the dialog still confused me. Do I need to change the parameter "WB_evap_value" tofrac? According to Roger comment in PR #8404 the ucrNudWB and ucrChkWB are already working but need some improvment. So is there something to do else a part from Roger comment on PR #8404?

[lilyclements]

@MeSophie great.

So if you currently look at the dialog code, you'll see we run "Reduce" as an R function somewhere, with some parameters:

f x init accumulate

If this new checkbox is not checked, then we run the code as usual:

f = function(x, y) pmin(pmax(y + x, 0), capacity) x = utils::tail(rain_min - evaporation_value, n=-1) init = 0 accumulate = TRUE

If it is checked, then we edit the f and the x functions: f = function(x, y) pmin(pmax(y - WB_evaporation(x, frac, capacity, evaporation_value, y) + x, 0), capacity) x = utils::tail(rain_min, n=-1)

Currently For f, for example, we currently have this:

        clsPMinWBMinFunction.SetRCommand("function(x, y) pmin")
        clsPMinWBMinFunction.AddParameter("0", clsRFunctionParameter:=clsPMaxFunction, iPosition:=0, bIncludeArgumentName:=False)
        clsPMinWBMinFunction.AddParameter("1", iCapacityDefault, iPosition:=1, bIncludeArgumentName:=False)

...

        clsPMaxFunction.SetRCommand("pmax")
        clsPMaxFunction.AddParameter("0", "x + y", iPosition:=0, bIncludeArgumentName:=False)
        clsPMaxFunction.AddParameter("1", "0", bIncludeArgumentName:=False)

So line by line: Line 1: clsPMinWBMinFunction.SetRCommand("function(x, y) pmin") We set f to be the function function(x, y) pmin. We don't need to change this, because that's still what we want

Line 2: clsPMinWBMinFunction.AddParameter("0", clsRFunctionParameter:=clsPMaxFunction, iPosition:=0, bIncludeArgumentName:=False) We add a parameter 0 which equals clsPMaxFunction. This is given on line 4.

Line 3 clsPMinWBMinFunction.AddParameter("1", iCapacityDefault, iPosition:=1, bIncludeArgumentName:=False) We add another parameter 1 which equals iCapacityDefault. function(x, y) pmin(iCapacityDefault)

Line 4: clsPMaxFunction.SetRCommand("pmax") We have set a new R command as pmax. Hence, using lines 1 and 2 we now have function(x, y) pmin(pmax(), iCapacityDefault)

Line 5: clsPMaxFunction.AddParameter("0", "x + y", iPosition:=0, bIncludeArgumentName:=False) We add a parameter to pmax. Parameter is "0", parameter value is "x+y". We now have: function(x, y) pmin(pmax(x + y), iCapacityDefault)

Line 6: clsPMaxFunction.AddParameter("1", "0", bIncludeArgumentName:=False). We add a parameter to pmax. Parameter is "1", parameter value is "0". We now have: function(x, y) pmin(pmax(x + y, 0), iCapacityDefault)

And so this bit of code has made our "f = function(x, y) pmin(pmax(y + x, 0), capacity)"

Now, we want to have instead: f = function(x, y) pmin(pmax(y - WB_evaporation(x, frac, capacity, evaporation_value, y) + x, 0), capacity)

So, the bit that changes is in the pmax function. Unchecked: pmax(x + y) Checked: pmax(y - WB_evaporation(x, frac, capacity, evaporation_value, y) + x, 0) So we want something like:

clsPMaxFunction.AddParameter("0", << some string >>, iPosition:=0, bIncludeArgumentName:=False)

Where our string is x + y - clsWBEvaporation

And clsWBEvaporation is: WB_evaporation(x, frac, capacity, evaporation_value, y) E.g.

clsWBEvaporation.GetRCommand(WB_evaporation)
clsWBEvaporation.AddParameter("water_balance", "x", iPosition:=0)
clsWBEvaporation.AddParameter("frac", << value in WB nud >>, iPosition:=1)
clsWBEvaporation.AddParameter("capacity", iCapacityDefault(? - check what this is! capacity value), iPosition:=2)
clsWBEvaporation.AddParameter("evaporation_value", << value in the evaporation input or variable receiver >>, iPosition:=3)
clsWBEvaporation.AddParameter("rain", "rain_min", iPosition:=4)  # and rain_max later for WB max.

The bits in "<< >>" are just me trying to write words, they wouldn't be in the VB code.

Have a read through and follow the VB code along side it

[MeSophie]

@lilyclements Please when the Reducing WB is checked and also the rdovalue is checked do I need to remove the evaporation value in Pmax function? I means the code that I should have is function(x, y) pmin(pmax((x + y-WB_evaporation(x, 0.5, 60, 5, rain=rain_min)) - 5, 0), 60) or function(x, y) pmin(pmax((x + y-WB_evaporation(x, 0.5, 60, 5, rain=rain_min)) , 0), 60)

And also I have this error trying to run the code

This is how I defined my functions

[lilyclements]

@MeSophie you make a good point in your comment - the error says that "could not find function "WB_evaporation"." That is completely my mistake - WB_evaporation is a function I created, and then ... didn't add in! Can you add the following to the stand_alone_functions.R file on the branch you are working on?

WB_evaporation <- function(water_balance, frac, capacity, evaporation_value, rain){
  if (water_balance >= frac*capacity){
    evaporation <- evaporation_value
  } else {
    if (rain == 0){
      evaporation <- evaporation_value * ((water_balance)/(frac*capacity))
    } else {
      if (water_balance < frac*capacity){
        if (rain > evaporation_value){
          evaporation <- evaporation_value
        } else {
          evaporation <- evaporation_value * ((water_balance)/(frac*capacity))
          evaporation <- evaporation + ((evaporation_value - evaporation)*(rain/evaporation_value))
        }
      } else {
        evaporation <- evaporation_value
        }
    }
  }
  return(evaporation)
}

[lilyclements]

@MeSophie good point I suggest that when the variable rdo is checked, and the variable is run, lets let the variable receiver be evaporation_variable.

Then when the checkbox is checked and the variable rdo is run, we do:

evaporation_value <- ifelse(test=is.na(x= evaporation_variable), yes=5, no= evaporation_variable)
Reduce(f = function(x, y) pmin(pmax(y - WB_evaporation(x, frac, capacity, evaporation_value, y) + x, 0), capacity),
                                x = utils::tail(rain_min, n=-1), init = 0, accumulate=TRUE))

if the value rdo is run, we do:

Reduce(f = function(x, y) pmin(pmax(y - WB_evaporation(x, frac, capacity, evaporation_value, y) + x, 0), capacity),
                                x = utils::tail(rain_min, n=-1), init = 0, accumulate=TRUE))

where evaporation_value is the value in the nud.

I need to make some amendments to the WB_evaporation function to accomodate for this, so for the meantime I suggest to get it set up for evaporation_value, not evaporation_variable

[MeSophie]

@lilyclements I still have problem with wb_evaporation function. This is the code that I have following your comment.

# Code generated by the dialog, Transform

transform_calculation <- instat_calculation$new(type="calculation", function_exp="Reduce(function(x, y) pmin(pmax((x + y - WB_evaporation(x, 0.5, 60, 5, rain)), 0), 60), utils::tail(x=rain, n=-1), 0, accumulate=TRUE)", result_name="water", manipulations=list(), save=2, before=FALSE, calculated_from=list("dodoma"="rain"), adjacent_column="rain")
data_book$run_instat_calculation(calc=transform_calculation, display=FALSE)
rm(transform_calculation)

And this is the error

Please can you check it again?

[lilyclements]

@MeSophie what data are you using so that I can try to replicate this fully? I have noticed you have x=rain. I think you actually want to have x = min_rain or x = max_rain depending if you are doing the min or max bit.

I also noticed you have manipulations=list(). I think in here you would want something referring to the min/max rain calculation that should be done before this.

[MeSophie]

@lilyclements you can use dodoma data.

[lilyclements]

@MeSophie great - thanks. Does it get solved if the min_rain and max_rain are called in instead of rain? (and then the relevant bits in manipulations = list(). I believe that is what is run when the reducing checkbox are not checked, and that might explain the issues with missing values?

[MeSophie]

@lilyclements before change rain_minto rainI tried to run the code and obtain an error

[lilyclements]

Sure, so this error is because "rain_min" doesn't exist. Look at the code that runs when you don't click the checkbox - there would be a line that creates the "rain_min" variable.

[MeSophie]

@lilyclements Here is the code when you don't check the box. To be sure that the code is correct I run the code on my branch than switch to another branche and the is no line which produce rain_min. Or I,m the one to put it. If so please tell me how.

[lilyclements]

@MeSophie Gotcha - you're looking at the Transform dialog, and I'm looking at the End of Season dialog! Makes sense.

Thanks for flagging this - I'll fix up this problem in the R code for the Transform dialog for when there is no rainfall.

In the meantime, can you look at and edit the End of Season dialog? That one does use min_rain and max_rain.

[MeSophie]

@lilyclements Please any change about R code for the Transform dialog?

[lilyclements]

@MeSophie thank you for chasing me up on this. I haven't got around to it yet for the Transform dialog, but it should work for the End of Season dialog. Are you okay to work on that in the meantime? EDIT: Just saw your PR on this. I'll look at this first thing tomorrow

Thank you for being patient here - this is a fairly challenging (but fun! I hope!) task you have. I'll try to be more responsive - I've been on/off on leave.

[MeSophie]

@lilyclements Thank you for updating me. I saw your comment about End of rain Dialog. I'm on It.

[rdstern]

@MeSophie is this now finished? If so, then can you close it?

[MeSophie]

@MeSophie is this now finished? If so, then can you close it?

@rdstern I am not sure that this issue is completely resolved. I worked mainly on the end of season dialogue and very little on the transform dialogue. The codes I had were for End of Rain dialogue. So @lilyclements could you please check the code for transform Dialogue to identify what's missing? Thank you.

[lilyclements]

@rdstern am I correct in thinking that we want the "transform" bit on water balance to perform exactly the same as "end of seasons" on the "end rains" dialog?

[lilyclements]

@rdstern the issues reported by the CIMH team on the End of Season code come from them having missing values in their "evapotranspiration" variable. I have replied to them giving a solution for now on how to replace these missing values in the meantime, but we should have a solution for the next release in the R code. How do you suggest we handle these missing values?

Earlier in this issue, you say:

if wb (before) >= fraction capacity then evap on day n = full value If wb(before) < fraction capacity and rf on day n = 0 then evap on day n = evap normally wb(before)/ (fraction capacity).

Still have to worry about the adjustment to the evap if there is rain (and wb (n-1) is below the threshold * capacity:

a) If rain > full evap on day n, then evap is normal evap on that day. It isn't reduced.

So I just want to draw attention to the "a". The issue occurs when we check if "rain > full evap on day n". If "full evap on day n" is missing, then how do we want to check this?