sbegueria / SPEI

An R package for computing the Standardized Precipitation-Evapotranspiration Index (SPEI), Penman-Monteith and other reference evapotranspiration, SPI, etc.
78 stars 34 forks source link

[enhancement] Add events characteristics to SPEI output #41

Open aldotapia opened 2 years ago

aldotapia commented 2 years ago

Based on the package's example:

library(SPEI)

data("wichita")

wichita$PET <- thornthwaite(wichita$TMED, 37.6475)
wichita$BAL <- wichita$PRCP-wichita$PET
wichita <- ts(wichita[,-c(1,2)], end=c(2011,10), frequency=12)

spei12 <- spei(wichita[,'BAL'], 12)

Since events duration, peak and severity can be obtain from spei output, it would be great to attach a table with the summary of these.

Some thoughts using dplyr:

data.frame(year = c(time(spei12$fitted)) %/% 1,
           month = c(cycle(spei12$fitted)),
           values = c(spei12$fitted)) %>%
  filter(complete.cases(.)) %>%
  mutate(group = with(rle(sign(values)),rep(seq_along(values), lengths))) %>% 
  group_by(group) %>% 
  summarise(duration = n(),
            peak = ifelse(max(values)<0,min(values),max(values)),
            severity = sum(values),
            start = paste0(first(year),'-',sprintf("%02d",first(month))),
            end = paste0(last(year),'-',sprintf("%02d",last(month))))

# A tibble: 43 × 6
   group duration    peak severity start   end    
   <int>    <int>   <dbl>    <dbl> <chr>   <chr>  
 1     1       18 -1.92   -20.2    1980-12 1982-05
 2     2        4  0.596    1.76   1982-06 1982-09
 3     3        6 -0.730   -2.64   1982-10 1983-03
 4     4        1  0.283    0.283  1983-04 1983-04
 5     5        2 -0.187   -0.267  1983-05 1983-06
 6     6       11  0.924    6.41   1983-07 1984-05
 7     7       15 -1.23   -10.3    1984-06 1985-08
 8     8        6  0.764    2.36   1985-09 1986-02
 9     9        2 -0.0715  -0.0887 1986-03 1986-04
10    10        2  0.0430   0.0739 1986-05 1986-06
# … with 33 more rows

What do you think?

sbegueria commented 2 years ago

Hi, we have used a similar function in some of our works. It'd be nice to include the routine in a function, with the ability to define a custom threshold for drought other than 0.

aldotapia commented 2 years ago

Excellent, I'm gonna refactoring it with base R, since dpylr isn't a SPEI's dependence