Adding tricot rainfall indices into R-Instat

[lilyclements]

The rainfall indices used in Tricot is "Total Rainfall", "SDII", and/or "Rx5day" "Max. 5 Days Precipitation" so we want to be able to get these indices:

• Rx5day: maximum 5-day precipitation (mm)
• Rtotal: total precipitation (mm) in wet days, rain >= 1 (mm)
• SDII: simple daily intensity index, total precipitation divided by the number of wet days (mm/days)

This is run for each ID (essentially a "station")
This is run from a different start day depending on the "station" (defined by the planting_date variable in the covar data set)
All run for a length of 45 days

@rdstern are you able to outline to me if there is a way at the moment for us to calculate:

1. The total precipitation on days with at least 1mm of rainfall from the planting_date variable in the covar data set to 45 days later (different planting_date for each "station")
   ( I can see how to do this with the "Climatic Summaries" dialog for the Rainfall Total, but not how to run this for rainfall above a threshold amount - unless I convert all rainfall amounts less than 1 to 0, which is easy enough!)
2. The number of rainy days from the planting_date variable in the covar data set to 45 days later (different planting_date for each "station")
3. The maximum 5 day precipitation from the planting_date variable in the covar data set to 45 days later (different planting_date for each "station")

Note that I am able to do 2. ordinarily, but, we want to look at it from the planting_date variable in the covar data set to 45 days later. I think you could probably quite quickly tell me if this is possible or not at the moment!

I've attached the CHIRPS data and corresponding data which defines the start year here.

They do have a few other rainfall indices which they look at. We can look at getting these into R-Instat too if we want!

[rdstern]

@lilyclements thanks for sending the message. In principle we can do everything, but I think there will be quite a lot of "devil in the detail", so we may have a bit of work to do to make it all easy.

Here are the rainfall indices we can get easily - from climdex package:

And in our climatic summaries we can get the summaries for 45 days from varying starting days:

Perfect - except those two figures are from different dialogs.

So, if we have to use the current R-Instat, then we could either use the calculator first, to get the variables, or hopefully use the Climatic > Prepare > Transform dialog to get the variables we need.

I would strongly prefer the transform dialog, because it (I think) knows about multiple stations, and the calculator does not.

I suspect we may need to amend the Prepare > Transform dialog a bit. We may also wish to consider adding a Climatic tab with these summaries to the Climatic > Prepare > Summarise sub-dialog, to make the calculations an easier one-step operation?

In the dataset did you get the chirps data - the file just had the day numbers. It would be nice to try to see what we can do now?

[lilyclements]

Ooo look at those indices! Nice! Can we add it into the Climatic > Prepare > Transform, perhaps under a general "Summaries" heading? Or "Indices"? Not sure! But, this is great! Thank you for sharing.

Apologies for that on the data. I'm not sure how I lost the rainfall variable on the CHIRPS data. I've reattached it now to this message

tricot_climatic_data.zip

[rdstern]

@lilyclements Thanks I'll have a look. There was one further point with R-Instat and statistical aspects in the data. We want some indices for the 45 days following the date of planting.

Now what are the dates of planting? I noticed for some farmers they were in the November-January "season". But when I look further down there are some planting dates more in (maybe) June/July. I need to look more carefully.
a) R-Instat point - we shift the season to cope easily with the November-January planting - which is like Zambia. We leave it starting in January for (say) Ghana, where the season is UK Summer months. There are places where there are two planting seasons - maybe like these data - but this is the first time that we want to deal with the two seasons together. There is no fundamental reason why our "events" need to stop at the end of each year/season, (e,g, data starting in January, and we look for the rainfall totals in 45 days, when planting is any time from November to January.) But I'm not sure if our current code copes yet.
b) Statistical point: I'm not comfortable just looking at a few climatic covariates when there are really different seasons in different farms in the study. That's where there is no overlap in the sets of 45 days for different sets of farms. I would also want to include co-factors in the study. So, is there evidence of a seasonal effect in the data and possibly an interaction between the seasonal and treatment effects? If so, then we move to the climatic covariates to see how much of the seasonal effect can be explained by the climatic indices?

Let me give you an example from India. The monsoon season is (say) June to September and that's when lots is grown. That's on all soil types. Then, they mainly clay soils hold much more water than the more sandy soils. So there is a second season from (say) October to January on the clay soils. So, any seasonal effects could either be soil or climate related. If data spans both seasons, then nice to see first is there is any seasonal effect. If so, then how much, and then why. We often get back to levels of variation - here there is variation in the climatic indices, both between and within seasons. This is a big study, so nice not to confuse them. In addition, often some questions are answered more simply if you consider this as two studies. So, if you are growing these crops in the main season, then which is best. Rather than always lumping all together!

[rdstern]

@lily, thanks for the data. Here are the tasks, copied from above:

"The rainfall indices used in Tricot is "Total Rainfall", "SDII", and/or "Rx5day" "Max. 5 Days Precipitation" so we want to be able to get these indices:

Rx5day: maximum 5-day precipitation (mm)
Rtotal: total precipitation (mm) in wet days, rain >= 1 (mm)
SDII: simple daily intensity index, total precipitation divided by the number of wet days (mm/days)
This is run for each ID (essentially a "station")
This is run from a different start day depending on the "station" (defined by the planting_date variable in the covar data set)
All run for a length of 45 days"

Let me take them in turn, starting from the chirps rainfall as provided:
a) Checked the 3 variables - Used Climatic > Examine > Summarise: About 1.5 years, and looks ok
b) Climatic > Dates > Fill date Gaps - nothing to fill - good
c) Climatic > Dates > Use Date - add year, month, doy, day of month.
d) Climatic > Define Climatic Data

Seems fine so far.

e) Rx5day. I think that is easy as now. First calculate the moving 5-day sums:

f) Rtotal and Rcount.

Rcount is easy already, with >=1mm using the Count option below:

I suggest we include 1 in the drop-down list. (Still let's keep 0.85 as the default.)
Also add 0.5 which is used in Zambia, and change 4.85, 9.85. and 24.85 to 4.95 etc.

(Note we don;t need to calculate it as an extra variable - but often simple and convenient.)

Rtotal with a threshold we don't yet have on this dialog.

It is easy to do with the calculator - so I'll do it that way, to continue checking. But we should add an option into this dialog, to make it easier. We could add an extra button at the top, but I wonder instead of adding an extra checkbox to the Moving button, see below:

In the figure above, when we do the Moving summary over 1 day, then all the summaries (Mean, Sum, etc) are the same and so are the following controls. So we could disable them and add a new checkbox Move zero to xxx. The xxx is the same drop-down as in the count, but with the default of 1 (but see below) and no control after - it is always effectively >=. And the default name changes from Moving to adding a 1 to the variable name.

I have now done that with the calculator. I used simply:

rain1 <- ifelse(rain>=1, rain, 0) which works fine, but with 2 questions:

1. What does that formula do with missing values? We want to keep them as missing - does it do that?
2. Do we really want 1 in the formula: Here are the frequencies of the actual data:

Interesting that there are 46 rows of data with 0.96 in the rain column and nothing else less than 1mm. The frequencies are a little odd anyway, I now suggest we have 0.95 as our default. (And we look at the chirps values in detail more generally.) If we do use 0.95 should we do a more complicated transformation and make those values into 1 for simplicity? Or doesn't it matter?

So we now have the daily variables needed for the indices:

g) I add an index, so we now need to get:

Rx5day, Rtotal, Rcount and SDII.

I've added Rcount, because I think it is useful anyway. If I have understood correctly then we get the first 3 from the daily data, and then SDII = Rtotal/Rcount?

So now we use our climatic summaries dialog:

It is almost there.

1. We use the Station button at the top, because we want summaries by station.
2. Then we choose the planting date and add 44 days.
3. Sum and Count first , then we'll get the maximum for the Rx5day!

Let's see if it works, then I'll show the problem we need to solve: (What a nice example this is!)

There I've done it. @lilyclements I hope you are so impressed you don't look too closely, because then you see that the answers are all wrong! There is a bit of work needed for this method to work, but I think it can.

Three further points at this stage.

1. These resaults highlight an inherent limitation of the chirps data, namely there are a lot of exact repeats. That's because multiple stations are close to the same grid-point, and hence exactly the same data for the whole period. Note that feature will be hidden from the correct results, because of the different planting dates.
2. The whole process of getting these indices will be fairly easy - for this group - but nowhere near a single button. That's R-Instat. Nothing is difficult, but we want to combine maximum flexibility with pretty easy. We can make exactly these steps easier, by writing a script. Then you have just one button (Run All) and hey presto. Now you have very easy, but zero flexibility. (For example all the columns names must remain exactly the same.) So you are fine, until something goes wrong!
3. I bet this exercise in their current training is one where they need a merge. That's to get these indices added to the covar data frame.

Now I have actually done it, and I think here are the correct results with the existing R-Instat:

So we have two methods and I suggest we might allow both eventually.

1. We can also choose whether to put the summaries (indices) into a new data frame and then use Prepare > Data Frame > Add (Merge) Variables - which works quite well. Or we can add a key (ID) to the covar data frame and then link it and the chirps long data. I chose to do that. One limitation is that the linking only works if the linking variables have the same names in the two data frames. So I had to change planting_date to date in covar. Then the summaries go straight into covar.
2. The obvious method - which doesn't work as we saw earlier - is to use the Day Range to select the 45 days for each station. But there are 2 current problems for this Station option. The first is that it uses doy as the variable and here - with data for more than 1 year, this is not sufficient. We either need to add the year, or preferably use the date instead. Why not. The second problem is that the current calculation for the station option seems to ignore that filter anyway. I suggest making these edits should be done?
3. So instead I filtered the data first, and that worked fine. I first added the planting data (Add(Merge) dialog, to the Chirps data frame. Then I couldn't see how to use the filter dialog directly, (can anyone?), so I used the calculator as follows:

into <- ifelse((date <planting.date) | (date >= planting.date+45) ,FALSE ,TRUE )

Then I just filtered on the dates that into is TRUE.

This is turning into quite a nice example. If, instead, I make diff <-date - planting.date, which is easy with the calculator, then it is also easy to use the filter on the diff variable. I'm liking these options, in relation to different users of R-Instat.

Now our summaries - at the station level - work fine on these filtered data and give the results!

I'm liking the idea that, in our teaching, once we can do both methods, and we can link or not, then this might be a good example to show different ways of this data wrangling step - because it is often the harder step than the subsequent analyses. So, in the future we want different staff to choose the method that they find most suitable for them. That will give them more time to concentrate on becoming a good data detective and a powerful data analyst. If you are struggling with the computing steps, then you often try blindly with a single method, and then may get stuck if the data are at all awkward!

Here are the data again with the rainfall indices included:

tricot_climatic_data_indices.zip

[lilyclements]

@rdstern I will dig into this and reply properly to it later today (and compare these results to their function too) - but you're right, I am impressed - this R-Instat software has had quite a bit of thought behind it, I guess!

[rdstern]

@lilyclements as you dig, I found your initial message to be a good starting point. No hurry but you also said:

"They do have a few other rainfall indices which they look at. We can look at getting these into R-Instat too if we want!"

I found your initial issue to be a very useful starter. Happy for you to decide whether to discuss those in a new issue, and/or do a similar starting issue on their temperaure indices?

[lilyclements]

@rdstern I've finally read through and processed this! This is great. And a very nice example as you say - perhaps would be a good example to showcase for Tricot training videos on using the dialogs, especially since this follows the vignette so there is some alignment there.

I agree it would be good to show these different options on handling the day range in training, and to have the different options.

To answer some of your questions:

1. How do they cope with NAs in their vignette example? They just treat NAs as if they're 0's.
2. We can ask Kaue about the reason for using 1 as the default for rainy days.
3. We can discuss the temperature indices in a new issue perhaps. This is not covered in the vignette so I will need to check and familiarise myself with what indices they tend to use in the training. I also need to get an example with temperature data (chirts?). Although, come to think of it, what do we want in our first "draft" of this to show the tricots team? I was thinking we want to be able to show they can run their two vignettes in R-Instat, so temperature indices could wait a few months. You know how the rainfall and temperature indices intersect much better than me -- does it make more sense to dig into that now as well?
4. Same point to 3. but with the remaining rainfall indices. It could be quite good to see if we can replicate them in R-Instat, and where we can't, see if this is something to help enhance the dialogs further! This might be relatively straightforward since we already have the data for them. The definitions (which I believe you are familiar with anyway!) are the following, which I assume follow a similar structure of using the Transform then Climatic Summaries dialogs:

MLDS | maximum length of consecutive dry day, rain < 1 mm (days)
MLWS | maximum length of consecutive wet days, rain >= 1 mm (days)
R10mm | number of heavy precipitation days 10 >= rain < 20 mm (days)
R20mm | number of very heavy precipitation days rain >= 20 (days)
Rx1day | maximum 1-day precipitation (mm)
Rx5day | maximum 5-day precipitation (mm)
R95p | total precipitation when rain > 95th percentile (mm)
R99p | total precipitation when rain > 99th percentile (mm)
Rtotal | total precipitation (mm) in wet days, rain >= 1 (mm)
SDII | simple daily intensity index, total precipitation divided by the number of wet days (mm/days)

5. To help with 4., in case you are interested, I've attached the

vignette1_rainfall_indices.xlsx

Or if you wanted to try it yourself, the code is

library(chirps)
library(gosset)
library(climatrends)

data(nicabean)
covar <- nicabean$covar

# 1. Getting the Data frame (File) - Note this takes a while to run!
dates = c(min(covar[, "planting_date"]),
         max(covar[, "planting_date"]) + 70)

chirps = get_chirps(covar[, c("longitude","latitude")],
                   dates = as.character(dates),
                   as.matrix = TRUE,
                   server = "ClimateSERV")

# rename the date variable
newnames = dimnames(chirps)[[2]]
newnames = gsub("chirps-v2.0.", "", newnames)
newnames = gsub("[.]", "-", newnames)
dimnames(chirps)[[2]] = newnames

rain <- rainfall(chirps, day.one = covar$planting_date, span = 45)

Next steps to dialog changes (from what I can see - but please correct me if I'm wrong!)
a) On dlgTransform, under "Count", amend the drop-down list to include 0.5 and 1. And change 4.85, 9.85. and 24.85 to 4.95 etc.
b) On dlgTransform, under "Moving", add a new checkbox with an associated dropdown. This is default 0, but have in there options and we do the moving sum where the "Element" variable is >= the value in this dropdown.
c) On dlgTransform, under "Moving", change the new default name from "moving" to "<element_variable> _1"
d) On the Climatic Summary dialog, the current calculation for the station option seems to ignore the filter with dates (presumably a bug!)

Are there any other next steps you can see on the dialog changes?

[lilyclements]

@rdstern I can confirm we're getting the same Rtotal and SDII as them - but our Max 5 Days is different interestingly!

Looking at ID 2121, they got 11.70 and we got 23.57.

Our 23.57 is including rainfall not in the day 1-45 period, but, the rolling 23.57mm does occur in this period (days 1-2).

Since these are indices from Climdex, I looked at how the climdex.rx5day function handles it, and they have a function to decide. They have center.mean.on.last.day, which if TRUE takes goes from day 1 (like us), and if FALSE, goes from day 3.

E.g., see March 1961 from below:

• In Climdex, March 1961 would give 41.6mm (day 2) if center.mean.on.last.day = TRUE
• In Climdex, March 1961 would give 40.3mm (day 3) if center.mean.on.last.day = FALSE
• In Tricot, March 1961 would give 25.6 since it wouldn't look for rainfall that occurred outside of this month-year combination
• For us, March 1961 would currently give 41.6 as it is. If we put in a group_by for year-month when calculating the rolling sum, then it would give the same as tricot.

@rdstern I'm not sure how we want to handle it (or perhaps we already have this option in the dialog too). What do you suggest?

[rdstern]

@lilyclements quick response to some points above. More after we discuss.
a) When getting the 5-day maximum the transform dialog has the option of Right, Centred, Left. The default is currently right, on the assumption that you only know the 5-day total once you have had the 5 days. In this application I suggest there could be a good argument for centering, because it could be good to know about a heavy rainfall in the 2 days before planting?
b) I'd need quite a bit of convincing before including a year/month component into this calculation. It is the 45 days from planting that are important. Why do we need to know the month in addition?
c) They have quite a long list of other rainfall indices. And more could be added, e.g. we suggest 21 days after planting as interesting. Once it is easy to produce indices there can be lots. That's the benefit - and curse of having daily rainfall values to relate to simgle yield orders. Data at different levels is everywhere!
d) Of those you list, I think all except the last 2 are in the current summary or spells dialogs. The exceptions are R95p and R99p. They are in the Climdex dialog, but not in a way that can be related to exactly the 45 days from planting. I discuss what they are in the R-Instat climatic guide.
e) Oh they just maybe quite interesting to produce in the Transform dialog - and now perhaps by just adjusting the Count button. Here is the current dialog:

Currently we just give the 0/1 variable, i.e. the count. That's where we code the R10mm and R20mm, (though we don't yet have 10 to 20mm. I wonder how important that is? It is also a simple calculation if really needed! If we just aim to select one of them, then I'd prefer "ours"!)

We could just add the option to give the value instead of the count!
Then could we just add the 95% and 99% to the pull-down. Maybe others too, or one could type the percentage - so value and percent.
Then we have the option to sum over a number of days - maybe here 45? And we need to add the Left, Centred, Right option here too I guess? Here we would want Left, so we can associate the 45 day value

Interesting! But no. I think instead we just need mode summaries.

To be continued!