Additional tests TotalTimeExceedingThreshold S2S4E cases
Hi @amanriqu and @llledo, @lpalma and @ngonzal2 may also interested in this issue.
The function TotalTimeExceedingThreshold
has already been included in the master branch. I did basic tests and Chung did tests for his seasonal case.
Since this function computes the number of days, weeks or months exceeding a threshold, I think it is similar to the computation that you did in the S2S4E for extremes. Therefore, if I am right, I will need that you test the function for your cases.
In the sub-seasonal case, I understood that you do the weekly average. Let say that you have two data array with dimensions: member, sweek, lweek (lead time mean for each week), lat and lon. The first array will be a reference, maybe observations (dimension member of length 1) for which you want to know the absolute temperature that corresponds to the 90th percentile. It could also be the hindcast.
ref <- array(rnorm(5*4*10*2*3), c(member = 1, sweek = 10, lweek = 4, lat = 2, lon = 3))
p90 <- Threshold(ref, threshold = 0.9, sdate_dim = 'sweek')
dim(p90)
#lweek lat lon
# 4 2 3
For each gridpoint and lead time (weeks), Threshold() returns the corresponding value (e.g.: temperature) of the 90th percentile. It is possible to get the percentile using all start dates and lead times as samples if you wish:
p90 <- Threshold(ref, threshold = 0.9, sdate_dim = c('sweek', 'lweek'))
dim(p90)
dim(p90)
lat lon
2 3
Once you have the threshold translated to the same units of the data, you can use the function TotalTimeExceedingThreshold() to actually compute the indicator on the second data that in the case of S2S4E was a forecast, so, a dimension of sweek with length 1.
#using threshold depending on lweek:
fcst <- array(rnorm(5*4*2*3), c(member = 5, sweek = 1, lweek = 4, lat = 2, lon = 3))
index_weeksexceeding90p <- TotalTimeExceedingThreshold(fcst, p90, time_dim = 'sweek')
dim(index_weeksexceeding90p)
#member lweek lat lon
# 5 4 2 3
range(index_weeksexceeding90p)
#[1] 0 1
Given that you may want to know the result for each week separately, I have defined time_dim
as start date week and not lweek but it is flexible.
The result indicates for each gridpoint, member and week if the 90th percentile has been exceeded with 1 and if not with 0.
I know that in S2S4E do you provide the probability. Then, there is a missing function as we commented that could be called something like Alert(), or whatever, to explore from the result the number of members that detect the extreme to be exceeded and compare it with a probability threshold: if the 25 % of the members detect the extreme then shows a triangle.
I think that a similar approach applies to the seasonal case.
Then, if you agree that the S2S4E is performing the analysis of extremes in this way, could you add tests for your cases to the Threshold() and TotalTimeExceedingThreshold() functions?
There is also the possibility to create the function Alert if you wish.
Please let me know if I misunderstood something, maybe the functions are already able to fulfill your needs but I haven't reflected properly. On the contrary, let me know if you have doubts about creating the tests.
Thanks in advance,
Núria