Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fixed Reducing Water Issue in Transform Climatic Dialogue #8935

Merged
merged 24 commits into from
Jun 19, 2024

Conversation

MeSophie
Copy link
Contributor

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
Copy link
Collaborator

rdstern commented Apr 17, 2024

@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:

image

@MeSophie MeSophie marked this pull request as draft April 19, 2024 15:18
@MeSophie MeSophie marked this pull request as ready for review April 22, 2024 16:59
@MeSophie
Copy link
Contributor Author

@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
Copy link
Collaborator

rdstern commented Apr 23, 2024

@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:

image

@MeSophie
Copy link
Contributor Author

@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
Copy link
Collaborator

rdstern commented Apr 27, 2024

@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
Copy link
Contributor

@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
Copy link
Contributor

lilyclements commented May 1, 2024

@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
Copy link
Collaborator

rdstern commented May 1, 2024

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

@rdstern
Copy link
Collaborator

rdstern commented May 6, 2024

@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:

image

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
Copy link
Contributor

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
Copy link
Contributor Author

MeSophie commented May 8, 2024

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
Copy link
Contributor Author

MeSophie commented May 8, 2024

@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
Copy link
Contributor

@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
Copy link
Collaborator

rdstern commented May 13, 2024

@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:

image

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
Copy link
Contributor Author

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

@MeSophie
Copy link
Contributor Author

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

@lilyclements Please anything new about this?

@rdstern
Copy link
Collaborator

rdstern commented May 29, 2024

@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
Copy link
Contributor

@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:

image

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
Copy link
Contributor Author

MeSophie commented Jun 4, 2024

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

image

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
Copy link
Contributor

@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
Copy link
Contributor Author

MeSophie commented Jun 5, 2024

@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
Copy link
Contributor

@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
Copy link
Contributor Author

MeSophie commented Jun 5, 2024

@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
Copy link
Contributor Author

MeSophie commented Jun 5, 2024

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

rdstern
rdstern previously approved these changes Jun 5, 2024
Copy link
Collaborator

@rdstern rdstern left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@MeSophie and @lilyclements brilliant - I think this is now working though I have just tested for Dodoma. So very well done both of you.
I am approving, though I still am afraid I have one last suggestion, but @N-thony could you check in parallel, so we can get this merged as soon as possible.

So I will discuss my suggestion with Lily because we will meet tomorrow. We need real payoff from all this work. I suggest it is our most ambitious use of the calculation system, and most people will not notice, because they will not be mirroring the commands. But, for our documentation and papers it would be really good if we added comments to the code! Sophie I presume you now understand enough to do this?

Comment on lines 407 to 409
'Dim strEndofSeason As String = "end_of_season_combined"
'Dim strDoyFilter As String = "doy_filter"
'Dim strConditionsFilter As String = "conditions_filter"
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

remove?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is solved.

Comment on lines 427 to 431
'clsRWaterBalanceFunction = New RFunction
'clsRWaterBalanceFunction1 = New RFunction
'clsRWaterBalanceFunction2 = New RFunction
'clsWBEvaporation = New RFunction
'clsTailFunction = New RFunction
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

remove?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is solved.

Copy link
Collaborator

@rdstern rdstern left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@MeSophie I assume I am approving again. At the same time I record some actions I am taking to test the new dialog.
a) I introduced a missing value in the rainfall to check its action on the water balance:

image

This also confirms my request @MeSophie of the value for you to add comments to the R-code - perhaps in a new pull request.

b) I introduced a column - value 5 with the calculator, called evap - to check it is the same with a column as with a constant value.

c) Oh wow @MeSophie and @lilyclements I am no longer approving! But just request a very small change - which is now clearer to see as you have split up (I think) the R code, so it is clearer. I may be back-tracking on the need for more internal comments in the code - but I would like you to consider a small edit. Here is the result of the one missing value when you have reducing evaporation:

image

So, with reducing evaporation, it only seems to recover 200 days later - when there is a full profile, for the first time.

Now, in addition, here is your code from the script window:

image

That's beautiful, and it uses the R-Instat calculation system which we are now documenting! So this can become an excellent advert for R-Instat. And that's partly because I would like statistics packages to set an excellent example of coping well with missing values.

@MeSophie first, then maybe @lilyclements I am hoping that just a slight tweak in the code for wb_min. Currently I think you start with zero or max. But, with reducing evaporation you never get to zero, so the missing value trigger is only when the bucket is full?
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?

Exciting, because this is timed with including documentation on the R-Instat functions in a linkable way from the help.

And we need to write this up! I'll discuss with @volloholic

@MeSophie
Copy link
Contributor Author

MeSophie commented Jun 7, 2024

@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?

image

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
Copy link
Contributor

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.

  1. 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
Copy link
Collaborator

rdstern commented Jun 15, 2024

@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
Copy link
Contributor

@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
Copy link
Collaborator

rdstern commented Jun 18, 2024

@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!

Copy link
Collaborator

@rdstern rdstern left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@N-thony I am now approving this one, as I suggest it is working sort of ok and we can continue checking once merged. @lilyclements I hope you agree?
So, @N-thony could you check anyway, and then merge unless Lily complains earlier.

@N-thony N-thony merged commit b89a310 into IDEMSInternational:master Jun 19, 2024
2 checks passed
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

Reducing Evaporation in the Transform => Water Balance Dialogue not working as expected
4 participants