Blog ODES - Ejemplo de cálculo índice de sequía SPI-12

Acumulación de precipitación para 12 meses

Para el ejemplo se utilizan datos de precipitación acumulada mensual desde enero de 1981 a julio 2023 para la cuenca del río Bueno.

El set de datos contiene los valores de precipitacion obtenidos de la aplicación ODES Unidades. Tiene dos columnas, date para las fechas y pre para los valores de precipitación mensual.

Código

d |> 
  gt() |> 
  opt_interactive()

Código

# library(zoo)
# library(dygraphs)
# d_zoo <- zoo(d$pre,d$date)
# dygraph(d_zoo,main = 'Precipitación mensual (mm)')

highchart2() |>
  hc_add_series(data =d, type = 'column',name = 'Precipitación mensual (mm)',hcaes(x= date, y = round(pre))) |> 
  hc_xAxis(dateTimeLabelFormats = list(day = '%m  %Y'), type = "datetime")

El Índice Estanadrizado de Precipitación (SPI) considera la precipitación acumulada de \(n\) meses, para el ejemplo se considerará la escala de tiempo de \(n=12\), para doce meses. Esto quiere decir que para el cálculo del SPI-12 se acumulará la precipitación de los últimos 12 meses, partiendo desde el último mes de la serie, en este caso julio 2023 y yendo mes a mes hacia el primer mes de la serie que corresponde a enero de 1981. Por lo que la precipitación acumulada del mes de julio 2023 será.

\[PA12_{julio_{2023}} = \sum_{agosto_{2022}}^{julio_{2023}} P\]

dónde, \(PA12_{julio_{2023}}\) corresponde a la precipitación acumulada de los últimos doce meses para julio 2023; y \(P\) corresponde a la precipitación acumulada mensual. Esto se repite mes a mes hasta llegar a diciembre de 1981 (\(PA12_{dic_{1981}}\)).

\[PA12_{dic_{1981}} = \sum_{enero_{1981}}^{dic_{1981}} P\] Por lo que el primer mes de la serie será diciembre 2012, ya que los meses anteriores no alzanzan a acumular doce meses de precipitación.

La siguiente animación ejemplifica el proceso de acumulación de 12 meses desde julio 2023 hasta diciembre 1981.

SPI-12 para Julio

Del proceso anterior se obtiene la siguiente tabla con los valores acumulados de doce meses entre diciembre 1981 a julio 2023.

Código

library(zoo)
x.t <- zoo(rev(d$pre))
xt_sum <- rollapply(x.t,12,sum)

df_12 <- tibble(dates = dates[1:500],pre_12 = as.numeric(xt_sum))

df_12 |> 
    gt() |> 
  opt_interactive()

La serie de tiempo de los valores de precipitación acumulada de 12 mese es la siguiente:

Código

df_12 |> 
  ggplot(aes(dates,pre_12)) +
  geom_col(fill='blue',alpha = .6) + 
  scale_x_date(date_breaks = "2 years",date_labels = "%Y",expand =c(0,0)) +
  ylab('Precipitación acumulada 12 meses (mm)') +
  theme_bw() +
  theme(axis.title.x = element_blank())

Ahora para el ejemplo continuaré realizando el proceso para los meses de julio, sin embargo, para el resto de los once meses el procedimiento es el mismo.

Entonces, filtrando los julios tenemos:

Código

df_12_jul <- df_12 |> filter(month(dates) == 7)

df_12_jul |> 
  gt() |> 
  opt_interactive()

Código

df_12_jul |> 
  ggplot(aes(dates,pre_12)) + 
  geom_col(fill = 'blue') +
  scale_x_date(date_breaks = "2 years",date_labels = "%Y",expand = c(0,0)) +
  scale_y_continuous(expand = c(0,0)) +
  labs(y = 'Precipitación acumulada \n doce meses (julios) (mm)') +
  theme_bw() +
  theme(axis.title.x = element_blank())

Los datos de precipitación acumulada mensual de doce meses para todos los julios tiene la siguiente distribución de probabilidades.

Código

df_12_jul |> 
  ggplot(aes(pre)) +
  geom_density(aes(pre_12),fill = 'blue',alpha=.6) +
  theme_bw()

Esta distribución se normaliza (en caso de ser necesario), es decir se transforma a una distribución normal. Existen diferentes técnicas para normalizar una distribución, hay métodos paramétricos (ej: L-momentos) y no paramétricos. En ODES para el cálculo de todos los indicadores de sequía se utilizan métodos no paramétricos. Por simplicidad, en este caso de ejemplo se realizó el procedimiento de forma paramétrica utilizando {bestNormalize}.

La siguiente figura en el panel superior, muestra la distribución de valores de precipitación acumulada de 12 meses para todos los junios, normalizada. En el panel de abajo, se presenta la misma distribución pero con los valores estandarizados, es decir con una media 0 y desviación estandar de 1 (z-score).

Ambas corresponde a la misma distribución, pero la figura del panel de arriba tiene los valores absolutos y la de abajo valores relativos de acuerdo a una distribución normal estandarizada.

Código

library(bestNormalize)
set.seed(354)
norm <- bestNormalize(df_12_jul$pre_12,standardize = FALSE)
df_12norm <- tibble(date = df_12_jul$dates,pre_12n = norm$x.t,pre_12ns = as.numeric(scale(norm$x.t)))

vals <- c(mean(df_12norm$pre_12n)-sd(df_12norm$pre_12n),mean(df_12norm$pre_12n),mean(df_12norm$pre_12n)+sd(df_12norm$pre_12n))
plot_d1 <- ggplot() +
  geom_density(data =df_12norm,aes(pre_12n),fill = 'blue',alpha=.6) +
  geom_vline(xintercept = vals[1],col='red',lty = 'dashed',alpha =.4 ) +
  geom_vline(xintercept = vals[3],col='red',lty = 'dashed',alpha =.4  ) +
  geom_vline(xintercept = vals[2],col='red',lty = 'dashed',alpha =.4  ) +
  annotate('text',x = vals[1],y=0.0013,label = paste("mu-~sigma==",round(vals[1],1),"~mm"),parse = TRUE ) +
  annotate('text',x = vals[2],y=0.0013,label = paste("mu==",round(vals[2],1),"~mm"),parse = TRUE ) +
  annotate('text',x = vals[3],y=0.0013,label = paste("mu+~sigma==",round(vals[3],1),"~mm"),parse = TRUE ) +
  xlab('Precipitación acumulada 12 meses (julios) [mm]') +
  theme_bw()

plot_d2 <- ggplot() +
  geom_density(data =df_12norm,aes(pre_12ns),fill = 'blue',alpha=.6) +
  geom_vline(xintercept = -1,col='red',lty = 'dashed',alpha =.4 ) +
  geom_vline(xintercept = 0,col='red',lty = 'dashed',alpha =.4  ) +
  geom_vline(xintercept = 1,col='red',lty = 'dashed',alpha =.4  ) +
  annotate('text',x = -1,y=0.3,label = "mu-~sigma==-1",parse = TRUE ) +
  annotate('text',x = 0,y=0.3,label = "mu==0",parse = TRUE ) +
  annotate('text',x = 1,y=0.3,label = "mu+~sigma==1",parse = TRUE ) +
  xlab('z-score') +
  theme_bw()

library(gridExtra)
grid.arrange(plot_d1,plot_d2,nrow = 2)

Entonces, lo que hacen los índices de sequía y en este caso el SPI-12 para los meses de julio. Es llevar los valores de precipitación acumulado a una distribución normal estandarizada. En este caso el promedio de precipitación acumulada de 12 meses corresponden a \(1600\) yla desviación estandar a \(~310mm\), lo que en la distribución normal estandarizada corresponde a \(0\) y \(1\), respectivamente.

Luego, en las figuras de abajo, se muestran las series de tiempo de los valores normalizados absolutos y relativos. Los índices de sequía utilizan los valores relativos, ya que permiten la comparación con diferentes zonas y momentos del año, ya que no depende de la magnitud de la precipitación, sino que de la distribución de probabilidad.

Finalmente, para calcular el SPI-12 para el resto de los 11 meses, se debe seguir el mismo procedimiento.

Código

plot_ts1 <- df_12norm |> 
  mutate(preC = pre_12n - vals[2]) |>  
  ggplot(aes(date,preC)) +
  geom_col(aes(fill = pre_12n > vals[2])) +
  geom_col(aes(fill = pre_12n < vals[2])) +
  scale_fill_manual(values = c('blue', 'red')) +
  scale_y_continuous(labels = function(x) round(x + vals[2])) +
  ylab('Precipitación acumulada de \n 12 meses normalizada (mm)') +
  theme_bw() + 
  theme(legend.position = 'none',axis.title.x = element_blank())

plot_ts2 <- ggplot(df_12norm,aes(date,pre_12ns)) +
  geom_col(aes(fill = pre_12ns > 0)) +
  geom_col(aes(fill = pre_12ns < 0)) +
  scale_fill_manual(values = c('blue', 'red')) +
  ylab('SPI-12')+
  theme_bw() + 
  theme(legend.position = 'none',axis.title.x = element_blank())

grid.arrange(plot_ts1,plot_ts2,nrow = 2)