Fixed Reducing Water Issue in Transform Climatic Dialogue

[MeSophie]

Fixes #8891 @rdstern , @lilyclements I made some change on Transform> Water Balance dioalogue. Please have a look. @lilyclements Could you please check the code of the windows script? I didn't introduce the .x and .int parameters because they weren't in the initial reduce function. I don't know if this is the script to get because the summary of the result is very large compared to the End of Saison Dialogue.

[rdstern]

@MeSophie good you are back to fixing this. a) Could you include the default in this water option of including the rain variable in the dialog if it finds one. b) I first tried the default of fixed 5. Seems ok. c) Then I put 5 in a variable called evap and tried. 1) It starts with -10 (there was a bug before so it started with -5. It should never be negative of course.
2) And you seem to be subtracting double the value in the column, so 10 always? d) Then constant again and it gives an error:

[MeSophie]

@rdstern Could you please test this dialogue again? @rdstern or @lilyclements I'm wondering why for a value below 17 for Perman data and 20 I think for dodoma and reducing box checked, the water values are all the same. At first I thought the dialogue wasn't working. Then I had the idea of doing the annual summary of the results and I noticed a change, not very visible with dodoma but visible for Perman. On the other hand, when the variable is used, the variation in reducing is very quickly noticeable.

[rdstern]

@MeSophie it is getting closer, though I can't yet check the reducing part, see c) below: a) While adding this component, I note that we now make the default capacity 100 in the End of season dialog. For consistency could you change the 60 in this dialog to 100. b) Second small issue. I check by adding a constant 5 column and calling it evap. Then, using this variable should give the same results as when we use the value 5. It does except for the very first day, when it gives -5. I think this may always have been there. Of course, with the formula we use, it should never be negative.
c) And the main one is that the reducing option doesn't work. I always get:

[MeSophie]

@rdstern The error message is due to missing values in the rain variable. And if I have understood correctly, part c) of your first comment ties is the same as part b) of your second comment. These negative values have been there since before the Reducing box was added. This may come from the calculation, I don't know if @lilyclements can help us here.

[rdstern]

@lilyclements and @MeSophie As you say, the problem of the -5 on day 1, when you allow for an evaporation variable should be examined, but there is a much larger problem when there are missing values in the rainfall variable. I am ok, if missing values in evaporation are not allowed, but we need to be able to "recover" from missing values in rainfall. Currently, if I have 1 missing day in 1935 for Dodoma, then thwe whole of the water variable is missing from then on. I am sure we specified something more general earlier, and hope we can be a bit cleverer in the water code now. I specify it mainly for the non-reducing situation, so we (at least get that working correctly. Lily, I specify in words and hope it is easy to translate into R code. a) If rain is missing then water is automatically missing.
b) Once the rain is not missing again we work out the water balance on two different scenarios. First if the water balance is 0 on the day before, and second if it is full on the day before. We continue to make water balance missing when these 2 are different. When they are the same that's the value of the water balance and we continue like that. I assume we may want to make those 2 variables for the whole record and then use them at the end to specify the missing values in the water balance? I suggest with the reducing option we ignore that aspect when working out the missing values for the water balance.

[lilyclements]

@rdstern I just want to confirm:

1. Currently, if the rainfall is missing, the water balance is missing for that day.
2. If the rainfall value of yesterday were missing, replace the missing rain value with 0mm, and calculate the WB. Otherwise, just calculate the WB on the rainfall given.
3. If the rainfall value of yesterday were missing, replace the missing rain value with the maximum (e.g., 100mm), and calculate the WB. Otherwise, just calculate the WB on the rainfall given.
4. We make water balance missing when these 2 are different. When they are the same that's the value of the water balance and we continue like that.

E.g.

wb_function <- function(rainfall){
  rain_min<- ifelse(test=is.na(x=rainfall), yes=0, no=rainfall)
  wb_min <- purrr::accumulate(.f= ~ pmin(pmax(..1 + ..2, 0), 60), .x=tail(x=rain_min - 5, n=-1), .init=0)
  rain_max <- ifelse(test=is.na(x=rainfall), yes=100, no=rainfall)
  wb_max <- purrr::accumulate(.f= ~ pmin(pmax(..1 + ..2, 0), 60), .x=tail(x=rain_max - 5, n=-1), .init=0)
  wb <- ifelse(test=(wb_min != wb_max) | is.na(x=rainfall), yes=NA, no=wb_min)
  return(wb)
}

rainfall <- c(10, NA, 20, 24)
wb_function(rainfall)
# gives 0, NA, NA, NA because days 3 and 4 did not have the same value for wb_min as for wb_max

rainfall <- c(60, NA, 60, 60)
wb_function(rainfall)
# gives 0, NA, NA, 60 because day 3 did not have the same value for wb_min as for wb_max,
# but day 4 did have the same value

This is currently what is done

[lilyclements]

@rdstern I'm not getting the same issue as you are for dodoma. If one day is NA, I don't get that every day is NA. This R script should work in R (if it doesn't then replace 0.7.16 with your R-Instat version)

# Code generated by the dialog, Import Dataset
new_RDS <- readRDS(file="C:/Program Files/R-Instat/0.7.16/static/Library/Climatic/Tanzania/dodoma.RDS")
dodoma <- new_RDS$get_data_frame("dodoma")
dodoma <- dodoma %>% filter(year == 1935)
dodoma$rain[10] <- NA

dodoma %>%
  mutate(rain_min = ifelse(test=is.na(x=rain), yes=0, no=rain),
         wb_min = purrr::accumulate(.f= ~ pmin(pmax(..1 + ..2, 0), 60), .x=tail(x=rain_min - 5, n=-1), .init=0),
         rain_max = ifelse(test=is.na(x=rain), yes=100, no=rain),
         wb_max = purrr::accumulate(.f= ~ pmin(pmax(..1 + ..2, 0), 60), .x=tail(x=rain_max - 5, n=-1), .init=0),
         wb = ifelse(test=(wb_min != wb_max) | is.na(x=rain), yes=NA, no=wb_min))

In that code I have set day 10 (randomly) to be NA. If you run that, you can see that rain_min replaces the NAs with a 0, rain_max replaces an NA with the maximum. wb_min, wb_max, and wb variables are then calculated, with wb as missing when wb_min != wb_max, or when rainfall is missing.

[rdstern]

@lilyclements that's a relief - I thought that was originally implemented. So let me check again on @MeSophie latest version.

[rdstern]

@lilyclements and @MeSophie I am using the latest version now. The R-code for the Dodoma data is here:

# Dialog: Transform
transform_calculation <- instat_calculation$new(type="calculation", function_exp="purrr::accumulate(.f=~ pmin(pmax(..1 + ..2 - 5, 0), 100), rain1, accumulate=TRUE)", result_name="water1", manipulations=list(), save=2, before=FALSE, calculated_from=list("dodoma"="rain1"), adjacent_column="rain1")
data_book$run_instat_calculation(calc=transform_calculation, display=FALSE)

I added NA in January 1935 here:

And it remains NA for ever. I don't see in the code above, how it could do otherwise?

I then tried with a column, called evap, and got the same result. Then with reducing and it said missing isn't yet allowed.

More generally I would like to have a new version of R-Instat as soon as this reducing evaporation issue, is ok, and evapotranspiration is confirmed to be working now. Then we can tell CIMH they could test. @jkmusyoka is evapotranspiration now ok? Perhaps that could be your starting point on the climatic guide, and the R-Instat help?

[lilyclements]

And it remains NA for ever. I don't see in the code above, how it could do otherwise?

@rdstern thanks for sharing. My mistake - I was looking at the "End of Seasons", not the "Transform" dialog. So, as you pointed out, we want to be running the bits where we replace the NAs with 0/max and check if the value is the same etc etc.

@MeSophie the code that should run here is the same as the End of Seasons dialog - that is,

rain_min <- instat_calculation$new(type="calculation", function_exp="ifelse(test=is.na(x=rainfall), yes=0, no=rainfall)", result_name="rain_min", calculated_from=list("ghana"="rainfall"))
wb_min <- instat_calculation$new(type="calculation", function_exp="purrr::accumulate(.f= ~ pmin(pmax(..1 + ..2, 0), 100), .x=tail(x=rain_min - 5, n=-1), .init=0)", result_name="wb_min", sub_calculations=list(rain_min))
rain_max <- instat_calculation$new(type="calculation", function_exp="ifelse(yes=100, test=is.na(x=rainfall), no=rainfall)", result_name="rain_max", calculated_from=list("ghana"="rainfall"))
wb_max <- instat_calculation$new(type="calculation", function_exp="purrr::accumulate(.f= ~ pmin(pmax(..1 + ..2, 0), 100), .x=tail(x=rain_max - 5, n=-1), .init=0)", result_name="wb_max", sub_calculations=list(rain_max))
wb <- instat_calculation$new(type="calculation", function_exp="ifelse(test=(wb_min != wb_max) | is.na(x=rainfall), yes=NA, no=wb_min)", result_name="wb", sub_calculations=list(wb_min, wb_max))
conditions_filter <- instat_calculation$new(type="filter", function_exp="(wb <= 0.5) | is.na(x=rainfall)", sub_calculations=list(wb))
grouping_by_station_year <- instat_calculation$new(type="by", calculated_from=list("ghana"="station","ghana"="year"))
end_season <- instat_calculation$new(type="summary", function_exp="ifelse(test=is.na(x=dplyr::first(x=wb)), yes=NA, no=dplyr::first(x=doy))", result_name="end_season", calculated_from=list("ghana"="doy"), save=2)
end_of_season_combined <- instat_calculation$new(type="combination", manipulations=list(conditions_filter, grouping_by_station_year), sub_calculations=list(end_season))
data_book$run_instat_calculation(calc=end_of_season_combined, display=FALSE)

[MeSophie]

And it remains NA for ever. I don't see in the code above, how it could do otherwise?

@rdstern thanks for sharing. My mistake - I was looking at the "End of Seasons", not the "Transform" dialog. So, as you pointed out, we want to be running the bits where we replace the NAs with 0/max and check if the value is the same etc etc.

@MeSophie the code that should run here is the same as the End of Seasons dialog - that is,

rain_min <- instat_calculation$new(type="calculation", function_exp="ifelse(test=is.na(x=rainfall), yes=0, no=rainfall)", result_name="rain_min", calculated_from=list("ghana"="rainfall"))
wb_min <- instat_calculation$new(type="calculation", function_exp="purrr::accumulate(.f= ~ pmin(pmax(..1 + ..2, 0), 100), .x=tail(x=rain_min - 5, n=-1), .init=0)", result_name="wb_min", sub_calculations=list(rain_min))
rain_max <- instat_calculation$new(type="calculation", function_exp="ifelse(yes=100, test=is.na(x=rainfall), no=rainfall)", result_name="rain_max", calculated_from=list("ghana"="rainfall"))
wb_max <- instat_calculation$new(type="calculation", function_exp="purrr::accumulate(.f= ~ pmin(pmax(..1 + ..2, 0), 100), .x=tail(x=rain_max - 5, n=-1), .init=0)", result_name="wb_max", sub_calculations=list(rain_max))
wb <- instat_calculation$new(type="calculation", function_exp="ifelse(test=(wb_min != wb_max) | is.na(x=rainfall), yes=NA, no=wb_min)", result_name="wb", sub_calculations=list(wb_min, wb_max))
conditions_filter <- instat_calculation$new(type="filter", function_exp="(wb <= 0.5) | is.na(x=rainfall)", sub_calculations=list(wb))
grouping_by_station_year <- instat_calculation$new(type="by", calculated_from=list("ghana"="station","ghana"="year"))
end_season <- instat_calculation$new(type="summary", function_exp="ifelse(test=is.na(x=dplyr::first(x=wb)), yes=NA, no=dplyr::first(x=doy))", result_name="end_season", calculated_from=list("ghana"="doy"), save=2)
end_of_season_combined <- instat_calculation$new(type="combination", manipulations=list(conditions_filter, grouping_by_station_year), sub_calculations=list(end_season))
data_book$run_instat_calculation(calc=end_of_season_combined, display=FALSE)

@lilyclements Okay thank you for this clarification.

[MeSophie]

@lilyclements Please does this mean that I have to add the control water balance to the transform dialogue or should I just assign it the default value 0.50 in the code?

[lilyclements]

@MeSophie I am pretty sure you should add the water control button to the dialog. @rdstern can you confirm?

also just to check - @rdstern should the output of this be on the daily data, or do we expect it to give it on the summary (by year) data frame?

[rdstern]

@MeSophie this dialog starts with daily data and is just a transformation of it. So the result is also daily data. It should give the same daily values that are then also produced and summarised in the end of season. So the controls are the same as there (except, or course you don't need the threshold value. The idea is that, if you want to understand the results given by the end of season, then you produce the transformation here, and can then see when (each year) it drops below any threshold that you have specified.

[lilyclements]

@rdstern thank you for explaining.

If I am understanding correctly then @MeSophie, in this case the code we would like is:

rain_min <- instat_calculation$new(type="calculation", function_exp="ifelse(test=is.na(x=rainfall), yes=0, no=rainfall)", result_name="rain_min", calculated_from=list("ghana"="rainfall"))
wb_min <- instat_calculation$new(type="calculation", function_exp="purrr::accumulate(.f= ~ pmin(pmax(..1 + ..2, 0), 100), .x=tail(x=rain_min - 5, n=-1), .init=0)", result_name="wb_min", sub_calculations=list(rain_min))
rain_max <- instat_calculation$new(type="calculation", function_exp="ifelse(yes=100, test=is.na(x=rainfall), no=rainfall)", result_name="rain_max", calculated_from=list("ghana"="rainfall"))
wb_max <- instat_calculation$new(type="calculation", function_exp="purrr::accumulate(.f= ~ pmin(pmax(..1 + ..2, 0), 100), .x=tail(x=rain_max - 5, n=-1), .init=0)", result_name="wb_max", sub_calculations=list(rain_max))
wb <- instat_calculation$new(type="calculation", function_exp="ifelse(test=(wb_min != wb_max) | is.na(x=rainfall), yes=NA, no=wb_min)", save = 2, adjacent_column = "rainfall", result_name="wb", sub_calculations=list(wb_min, wb_max))
data_book$run_instat_calculation(calc=wb, display=FALSE)

Here,

• rainfall = variable in the element column
• 100 = capacity (default is 60)
• 5 = value in the "Value" input (or similarly the code for variable that we do in end of seasons)

Note that this is the same code as in end of seasons for rain_min, wb_min, rain_max, wb_max.

The code is pretty much the same as in end of seasons for wb,except that here, we have two additional parameters: save = 2, adjacent_column = "rainfall"

Thank you as always @MeSophie for being patient and asking good questions - it's definitely clearing things up for me too that I didn't know I didn't know! :)

[MeSophie]

@lilyclements , @rdstern thank you for your explanation.

[MeSophie]

@rdstern , @lilyclements The new change is ready to be checked. Please have a look.

[rdstern]

@MeSophie from first glance it seems to be working on the missing values. Great: a) But you seem to have "hard-wired" the column calledwateras the result. It ignores any change of name in the dialog. b) I next checked with an evaporation variable and the problem (-5) on the first day, seems now to be fived. c) Then the reducing evaporation and here is a problem:

I hope you can fix it, because it is just that the reducing goes over 100. i.e. it isn't respecting the maximum value. I can't check how it is coping with missing, until you fix this.

d) Finally, a small "tweak" that I'm sure you can do. The output looks "messy" with all those decimal places. Could you include round(water, 1) in the calculation. It could be added generally as the non-reducing will usually be to 0.1mm anyway. So the result is rounded each day.

Great - I think it is almost there!

[MeSophie]

@lilyclements Please could you have a look at c) of @rdstern comment ? Thank you.

[MeSophie]

@lilyclements Please could you have a look at c) of @rdstern comment ? Thank you.

@lilyclements Please anything new about this?

[rdstern]

@lilyclements I wonder if @jkmusyoka should take this over, if you don't have time now? He is now in UK and I assume you will be meeting on e-PICSA soon? Meantime, @MeSophie your R is preet good. Have you tried to examine the problem yourself?

[lilyclements]

@rdstern I am struggling to replicate your error - I have tried it with the same data, but I am getting it capped at the 100. Can you remember what parameters you used for this, and/or do you have a copy of the R code output?

I am getting:

Putting the R code here, in case it helps (probably future-me will appreciate this)

# Using Dodoma Data:

# Replace Value In Data
data_book$replace_value_in_data(data_name="dodoma", col_name="rain", rows="51", new_value=67.8)

# Code generated by the dialog, Transform
rain_min <- instat_calculation$new(type="calculation", function_exp="ifelse(test=is.na(x=rain), yes=0, no=rain)", result_name="rain_min", calculated_from=list("dodoma"="rain"))

wb_min <- instat_calculation$new(type="calculation", function_exp="purrr::accumulate(.f= ~ pmin(pmax(..1 + ..2, 0), 100), .x=tail(x=rain_min - 5, n=-1), .init=0)", result_name="wb_min", sub_calculations=list(rain_min))

rain_max <- instat_calculation$new(type="calculation", function_exp="ifelse(yes=100, test=is.na(x=rain), no=rain)", result_name="rain_max", calculated_from=list("dodoma"="rain"))

wb_max <- instat_calculation$new(type="calculation", function_exp="purrr::accumulate(.f= ~ pmin(pmax(..1 + ..2, 0), 100), .x=tail(x=rain_max - 5, n=-1), .init=0)", result_name="wb_max", sub_calculations=list(rain_max))

wb <- instat_calculation$new(type="calculation", function_exp="round(ifelse(test=(wb_min != wb_max) | is.na(x=rain), yes=NA, no=wb_min), 1)", result_name="water", adjacent_column="rain", save=2, sub_calculations=list(wb_min, wb_max))

data_book$run_instat_calculation(calc=wb, display=FALSE)

rm(list=c("rain_min", "wb_min", "rain_max", "wb_max", "wb"))

[MeSophie]

@lilyclements you can obtain the result buy using dodoma data Climatic Transform and reducing box.

rain_min <- instat_calculation$new(type="calculation", function_exp="ifelse(test=is.na(x=rain), yes=0, no=rain)", result_name="rain_min", calculated_from=list("dodoma"="rain"))
wb_min <- instat_calculation$new(type="calculation", function_exp="purrr::accumulate(.f= ~ pmin(pmax(..1 + ..2 - WB_evaporation(..1, 0.5, 100, 5, ..2), 0)), .x=tail(x=rain_min, n=-1), .init=0)", result_name="wb_min", sub_calculations=list(rain_min))
rain_max <- instat_calculation$new(type="calculation", function_exp="ifelse(yes=100, test=is.na(x=rain), no=rain)", result_name="rain_max", calculated_from=list("dodoma"="rain"))
wb_max <- instat_calculation$new(type="calculation", function_exp="purrr::accumulate(.f= ~ pmin(pmax(..1 + ..2 - WB_evaporation(..1, 0.5, 100, 5, ..2), 0)), .x=tail(x=rain_max, n=-1), .init=0)", result_name="wb_max", sub_calculations=list(rain_max))
wb <- instat_calculation$new(type="calculation", function_exp="round(ifelse(test=(wb_min != wb_max) | is.na(x=rain), yes=NA, no=wb_min), 1)", result_name="water", adjacent_column="rain", save=2, sub_calculations=list(wb_min, wb_max))
data_book$run_instat_calculation(calc=wb, display=FALSE)
rm(list=c("rain_min", "wb_min", "rain_max", "wb_max", "wb"))

[lilyclements]

@MeSophie thank you for sharing the code.

From what I can see from inspecting your code against mine, the issue is in the wb_min (and similarly for wb_max)

Yours currently runs:

purrr::accumulate(.f= ~ pmin(pmax(..1 + ..2 - WB_evaporation(..1, 0.5, 100, 5, ..2), 0)), .x=tail(x=rain_min, n=-1), .init=0)

But can you try the following here and in wb_max to see if it still has the issue?:

purrr::accumulate(.f= ~ pmin(pmax(..1 + ..2 - WB_evaporation(..1, 0.5, 100, 5, ..2), 0), 100), .x=tail(x=rain_min, n=-1), .init=0)

There’s a lot of brackets, I know. But essentially what you’re saying is to get the minimum out of “pmax(..1 + ………)” and nothing else. So it doesn’t have anything to cap it at.

This change says to get the minimum out of “pmax(..1 + ………)” and “100”. This means that the water balance value caps at 100 if “pmax(……) > 100”

[MeSophie]

@MeSophie thank you for sharing the code.

From what I can see from inspecting your code against mine, the issue is in the wb_min (and similarly for wb_max)

Yours currently runs:

purrr::accumulate(.f= ~ pmin(pmax(..1 + ..2 - WB_evaporation(..1, 0.5, 100, 5, ..2), 0)), .x=tail(x=rain_min, n=-1), .init=0)

But can you try the following here and in wb_max to see if it still has the issue?:

purrr::accumulate(.f= ~ pmin(pmax(..1 + ..2 - WB_evaporation(..1, 0.5, 100, 5, ..2), 0), 100), .x=tail(x=rain_min, n=-1), .init=0)

There’s a lot of brackets, I know. But essentially what you’re saying is to get the minimum out of “pmax(..1 + ………)” and nothing else. So it doesn’t have anything to cap it at.

This change says to get the minimum out of “pmax(..1 + ………)” and “100”. This means that the water balance value caps at 100 if “pmax(……) > 100”

@lilyclements Thank you the result looks good now. Please @lilyclements this 100 does this 100 come from the capacity control or it is invariable? Mean that I will also correct the End of Season dialog code?

[lilyclements]

@MeSophie good question - yes, you're right it comes from the capacity control. And good point about the End of Seasons dialog. You're absolutely right - this should be the case there too!

[MeSophie]

@MeSophie good question - yes, you're right it comes from the capacity control. And good point about the End of Seasons dialog. You're absolutely right - this should be the case there too!

@lilyclements Thank you.

[MeSophie]

@rdstern @jkmusyoka Please could you test the dialog again ? Thank you

[MeSophie]

@rdstern I have made a change to the code and it's now in the form: rain_min <- instat_calculation$new(type="calculation", function_exp="ifelse(test=is.na(x=rain), yes=capacity/20, no=rain)", result_name="rain_min", calculated_from=list("dodoma"="rain")) Before pushing the change there's something I'm curious about, so I preferred to post here so that @lilyclements can give her point of view too, as she understands the code better than I do. I have noticed that the initial value has no impact on the result. So it has no impact if you enter a missing value in rain (which I've checked is the same whether you have 0 or 5). I have even gone up to 100 and the results are still the same. Is this normal?

rain_min <- instat_calculation$new(type="calculation", function_exp="ifelse(test=is.na(x=rain), yes=5, no=rain)", result_name="rain_min", calculated_from=list("dodoma"="rain"))
wb_min <- instat_calculation$new(type="calculation", function_exp="purrr::accumulate(.f= ~ pmin(pmax(..1 + ..2, 0), 100), .x=tail(x=rain_min - 5, n=-1), .init=0)", result_name="wb_min", sub_calculations=list(rain_min))
rain_max <- instat_calculation$new(type="calculation", function_exp="ifelse(yes=100, test=is.na(x=rain), no=rain)", result_name="rain_max", calculated_from=list("dodoma"="rain"))
wb_max <- instat_calculation$new(type="calculation", function_exp="purrr::accumulate(.f= ~ pmin(pmax(..1 + ..2, 0), 100), .x=tail(x=rain_max - 5, n=-1), .init=0)", result_name="wb_max", sub_calculations=list(rain_max))
wb <- instat_calculation$new(type="calculation", function_exp="round(ifelse(test=(wb_min != wb_max) | is.na(x=rain), yes=NA, no=wb_min), 1)", result_name="water1", adjacent_column="rain", save=2, sub_calculations=list(wb_min, wb_max))
data_book$run_instat_calculation(calc=wb, display=FALSE)
rm(list=c("rain_min", "wb_min", "rain_max", "wb_max", "wb"))

[lilyclements]

Two questions

1. @rdstern what are you suggesting we set to be capacity / 20?

So, let's start small instead - I want to keep thinks simple so let's start anyway (maybe reducing or not) with capacity/20 - so 5mm if capacity = 100?

@MeSophie your R code change is good, but the R code now says that if we have a missing value in rain, set it to be capacity / 20? (then to get the wb_min, wb_max, and if they're the same keep it, etc). I just want to double check that is what we want to do from @rdstern's comment above.

2. About the initial value having no impact:

I have noticed that the initial value has no impact on the result. So it has no impact if you enter a missing value in rain (which I've checked is the same whether you have 0 or 5). I have even gone up to 100 and the results are still the same. Is this normal?

Good point @MeSophie, I will wait to see @rdstern's comment on this.

[rdstern]

@lilyclements I noticed, when the water balance was recovering from a misssing value it seemed to be working really well in the non-reducing case. So what to do after missing values in the reducing case. The method now seemed not to find the values equal, when both (starting full and starting empty) at the low end. It only seems to work when both are full. I assume that's because the formula at the empty end is complicated, because of the reducing nature? I was wondering how to solve that?

I'm not sure my starting at max/20 will help. Instead, could we try checking each day whether the abs(difference) (between starting full and starting empty) is less than a small amount - say 0.5mm rather than when they are equal?

If we get something that works, we will then also need to use the code in the end of season dialog?

Once this is finshed I would like us to write to CIMH to suggest we write this up jointly - not in the current contract, but nice as a paper. This could be you and Sophie here and Adrian and Lisa there?

[lilyclements]

@rdstern I just want to check I understand we're discussing the right bit:

Currently, when a rainfall value is missing:

1. Set rainfall to be the minimum value (5), and calculate the water balance
2. Set rainfall to be the maximum value (the max. capacity), and calculate the water balance.
3. If they're the same, take the water balance value. Otherwise, set as NA.

You're suggesting instead that when a rainfall value is missing:

1. Set rainfall to be the minimum value (5), and calculate the water balance
2. Set rainfall to be the maximum value (the max. capacity), and calculate the water balance.
3. If the absolute difference is less than 0.5mm, take the water balance value. Otherwise, set as NA.

Am I understanding this correctly?

Thanks!

[rdstern]

@lilyclements that sounds fine, though if it started at the lower limit previously, then maybe we keep it as before. I was under the impression it started at zero before as the lower point, and though this resulted in it never being equal, when both were at the lower limit. But it is also such a detail I'd mainly like to be able to close this. Though I am keen we contact Lisa and Adrian, to write up this stuff in a paper!