generalizacion fechas

[VeruGHub]

@ajpelu ¿podrías pasarme un ejemplo de script para generalizar las fechas y utilizar los "lag"? Gracias :)

[ajpelu]

Hola @VeruGHub
A ver si lo encuentro. Mientras te cuento:

• En este antiguo repo https://github.com/ajpelu/qpyr_dendro (hace mas de tres años) hice algo parecido para computar resiliencia de diversos índices de satélite y también de BAI (dendro). Las funciones eviResilience y baiResilience son las que usé. Estas son anuales y tendrían que adaptarse (tiene mucho código antiguo)
• Hace poco estuve trabajando en algo parecido en base a esas funciones para computo de resiliencia de índices espectrales (NDVi, EVI, etc) a escalas temporales mas detalladas (meses, quincenas). Ese es el código que no encuentro.

Déjame que lo busque y te digo.

[ajpelu]

Hola @VeruGHub,
sobre lo que te dije que iba a buscar, te paso aquí una modificación de la función que te comenté. Quizá no sea las mas operativa, pero a mi me sirvió.

library(lubridate) 
compute_resilience <- function(data,
                               events,
                               event_type,
                               lags_pre,
                               lags_post,
                               var_interes, avg = c("mean", "median")) {
  
resilience <- c()
  
  for (i in 1:length(events)) {
    if (event_type[i] == "monthly") {
      # genera event as date
      date_event <- as.Date(paste0(events[i], "-01"))

      # subsets
      df_event <- subset(data, date == date_event)
      df_pre <- subset(data, date >= (date_event %m-% months(lags_pre[i])) & date < date_event)
      df_post <- subset(data, date <= date_event %m+% months(lags_post[i]) & date > date_event)

      # compute averages
      pre <- apply(df_pre[, c(var_interes)], MARGIN = 2, FUN = avg)
      post <- apply(df_post[, c(var_interes)], MARGIN = 2, FUN = avg)
      ev <- get(var_interes, df_event)
    } else if (event_type[i] == "yearly") {
      date_event <- as.Date(paste0(events[i], "-01-01"))

      # subsets
      df_event <- subset(data, grepl(events[i], date))
      df_pre <- subset(data, as.numeric(format(date, "%Y")) %in%
        seq(from = (as.numeric(events[i]) - lags_pre[i]), to = as.numeric(events[i]) - 1))

      df_post <- subset(data, as.numeric(format(date, "%Y")) %in%
        seq(
          from = (as.numeric(events[i]) + 1),
          to = (as.numeric(events[i]) + lags_post[i])
        ))

      # compute averages
      pre <- apply(df_pre[, c(var_interes)], MARGIN = 2, FUN = avg)
      post <- apply(df_post[, c(var_interes)], MARGIN = 2, FUN = avg)
      ev <- apply(df_event[, c(var_interes)], MARGIN = 2, FUN = avg)
    } else {
      stop("Invalid event type")
    }

    # compute indices
    rt <- ev / pre
    rc <- post / ev
    rs <- post / pre
    rrs <- (post - ev) / pre


    out <- data.frame(cbind(
      event = events[i],
      var = var_interes,
      type = event_type[i],
      lags_pre = lags_pre[i],
      lags_post = lags_post[i],
      rt = rt,
      rs = rs,
      rc = rc,
      rrs = rrs
    ))
    rownames(out) <- NULL

    resilience <- rbind(resilience, out)
  }

  return(resilience)
}

Esta función necesita un data.frame de al menos dos columnas: date y variable de interés (EVI, temp, etc), y los argumentos son:

• events: un vector de los eventos en los que se quiere computar la resiliencia, e.g. : events <- c("2018-01", "2017-01", "2014", "2012")
• events_type: un vector del tipo de eventos (monthly, yearly). e.g. event_type <- c("monthly","monthly", "yearly", "yearly")
• lags_pre y lags_post: un vector con los el número de lags previos o posteriores al evento. Por ejemplo si quiero calcular la resiliencia de un evento mensual tres meses antes y dos meses despues del evento.
• var_interes el nombre de la variable de interés sobre la cual se quiere computar la resiliencia. Esto lo incluía por si en el data.frame tenía varias columnas al que le quisiera computar las métricas de resiliencia.
• avg sería la función usada para promediar los datos de los meses/años previos al evento y posteriores. Actualmente se implementa la media y la mediana, pero se podría incluir otras. De hecho con la dendro hice algo de biweight robust means, pero no está implementado en esta función.

Finalmente comentarte que está función no está optimizada para ser una función estable y robusta, solamente me ha servido en un par de casos, pero quizá se pueda optimizar y coger ideas. De hecho, aunque permite el cómputo de diferentes métricas de resiliencia a diferentes escalas temporales (e.g meses, años), no se si es muy buena idea que estén mezcladas las diferentes escalas en una sola función. También comentarte que tiene la dependencia de lubridate, pero creo que se puede quitar.

Bueno espero que te sirva,

Te pongo un ejemplo de uso:

x <- data.frame(
  stringsAsFactors = FALSE,
              site = c("APH","APH","APH","APH",
                       "APH","APH","APH","APH","APH","APH","APH","APH",
                       "APH","APH","APH","APH","APH","APH","APH","APH","APH",
                       "APH","APH","APH","APH","APH","APH","APH","APH",
                       "APH","APH","APH","APH","APH","APH","APH","APH",
                       "APH","APH","APH","APH","APH","APH","APH","APH",
                       "APH","APH","APH","APH","APH","APH","APH","APH",
                       "APH","APH","APH","APH","APH","APH","APH","APH",
                       "APH","APH","APH","APH","APH","APH","APH","APH",
                       "APH","APH","APH","APH","APH","APH","APH","APH","APH",
                       "APH","APH","APH","APH","APH","APH","APH","APH",
                       "APH","APH","APH","APH","APH","APH","APH","APH",
                       "APH","APH","APH","APH","APH","APH","APH","APH",
                       "APH","APH","APH","APH","APH","APH","APH","APH",
                       "APH","APH","APH","APH","APH","APH","APH","APH",
                       "APH","APH","APH","APH","APH","APH","APH","APH",
                       "APH","APH","APH","APH","APH","APH","APH","APH","APH",
                       "APH","APH","APH","APH","APH","APH","APH","APH",
                       "APH","APH","APH","APH","APH","APH","APH","APH",
                       "APH","APH","APH","APH","APH","APH","APH","APH",
                       "APH","APH","APH","APH","APH","APH","APH","APH",
                       "APH","APH","APH","APH","APH","APH","APH","APH",
                       "APH","APH","APH","APH"),
              tmin = c(3.09999990463257,
                       2.40000009536743,4,7.90000009536743,11.6000003814697,
                       15.6999998092651,15.3999996185303,12,8.19999980926514,
                       4.59999990463257,2,2.70000004768372,3.40000009536743,
                       3.79999995231628,5.40000009536743,7.19999980926514,
                       12.3000001907349,15.8000001907349,16.5,12.3999996185303,8,2,1,
                       -0.100000001490116,1.5,4.69999980926514,3.70000004768372,
                       9.39999961853027,13.6000003814697,17.2000007629395,16,
                       11.6999998092651,11.5,7.09999990463257,3,
                       1.29999995231628,2.09999990463257,3.09999990463257,
                       6.09999990463257,7.30000019073486,11.3999996185303,17.1000003814697,
                       17.2000007629395,12.6999998092651,7.80000019073486,
                       4.5,3.20000004768372,2.29999995231628,1.79999995231628,
                       2.5,7.59999990463257,9.10000038146973,
                       12.8999996185303,15.1000003814697,16.5,13.1000003814697,
                       9.60000038146973,4.90000009536743,2.59999990463257,
                       1.79999995231628,-1.39999997615814,3.5,3.79999995231628,
                       9.80000019073486,14.5,16,18.1000003814697,12.5,8.60000038146973,
                       4.69999980926514,3.20000004768372,2.20000004768372,
                       0.400000005960464,3.20000004768372,5.30000019073486,
                       6.90000009536743,11.1999998092651,15.5,16.7999992370605,
                       12.6999998092651,10.8000001907349,3.70000004768372,
                       2.20000004768372,3.09999990463257,1.29999995231628,
                       2.90000009536743,7.30000019073486,9.10000038146973,
                       11.8000001907349,14.8000001907349,15.8000001907349,
                       12.6999998092651,11,5.5,1.60000002384186,1.79999995231628,0,
                       4.40000009536743,6.59999990463257,11,12.6000003814697,
                       18.6000003814697,16.3999996185303,11.8000001907349,
                       9.30000019073486,5.90000009536743,5.40000009536743,3,
                       2.20000004768372,1.89999997615814,5,8,13.3000001907349,
                       16.7999992370605,16.2999992370605,13.3999996185303,
                       10.3000001907349,3.90000009536743,3.09999990463257,
                       0.400000005960464,2.79999995231628,3.5,6.59999990463257,
                       9.69999980926514,14.8999996185303,16.1000003814697,16.5,
                       12.3999996185303,10.8999996185303,5.40000009536743,
                       1.29999995231628,1.20000004768372,-0.300000011920929,1.5,
                       4.5,6.40000009536743,11.6999998092651,15,
                       16.2000007629395,13.3000001907349,7.90000009536743,3.59999990463257,
                       4.30000019073486,1.39999997615814,2.20000004768372,
                       3.70000004768372,3.70000004768372,9.10000038146973,12,
                       15.6000003814697,16.2000007629395,12.8999996185303,
                       9.39999961853027,3.79999995231628,4.30000019073486,
                       1.60000002384186,4.69999980926514,4.19999980926514,
                       5.09999990463257,9.69999980926514,12.1000003814697,
                       16.7000007629395,16.3999996185303,12.6999998092651,7.59999990463257,
                       6.19999980926514,2.29999995231628,1.70000004768372,
                       3.5,3.40000009536743,5.19999980926514,8.89999961853027,
                       12,16.2000007629395,17.2000007629395,12.8999996185303,
                       9.89999961853027,3.29999995231628,4.5),
              date = c("2007-02-01","2007-03-01",
                       "2007-04-01","2007-05-01","2007-06-01","2007-07-01",
                       "2007-08-01","2007-09-01","2007-10-01","2007-11-01",
                       "2007-12-01","2008-01-01","2008-02-01","2008-03-01",
                       "2008-04-01","2008-05-01","2008-06-01","2008-07-01",
                       "2008-08-01","2008-09-01","2008-10-01","2008-11-01",
                       "2008-12-01","2009-01-01","2009-02-01","2009-03-01",
                       "2009-04-01","2009-05-01","2009-06-01","2009-07-01",
                       "2009-08-01","2009-09-01","2009-10-01","2009-11-01",
                       "2009-12-01","2010-01-01","2010-02-01","2010-03-01","2010-04-01",
                       "2010-05-01","2010-06-01","2010-07-01","2010-08-01",
                       "2010-09-01","2010-10-01","2010-11-01","2010-12-01",
                       "2011-01-01","2011-02-01","2011-03-01","2011-04-01",
                       "2011-05-01","2011-06-01","2011-07-01","2011-08-01",
                       "2011-09-01","2011-10-01","2011-11-01","2011-12-01",
                       "2012-01-01","2012-02-01","2012-03-01","2012-04-01",
                       "2012-05-01","2012-06-01","2012-07-01","2012-08-01",
                       "2012-09-01","2012-10-01","2012-11-01","2012-12-01",
                       "2013-01-01","2013-02-01","2013-03-01","2013-04-01",
                       "2013-05-01","2013-06-01","2013-07-01","2013-08-01",
                       "2013-09-01","2013-10-01","2013-11-01","2013-12-01",
                       "2014-01-01","2014-02-01","2014-03-01","2014-04-01",
                       "2014-05-01","2014-06-01","2014-07-01","2014-08-01",
                       "2014-09-01","2014-10-01","2014-11-01","2014-12-01","2015-01-01",
                       "2015-02-01","2015-03-01","2015-04-01","2015-05-01",
                       "2015-06-01","2015-07-01","2015-08-01","2015-09-01",
                       "2015-10-01","2015-11-01","2015-12-01","2016-01-01",
                       "2016-02-01","2016-03-01","2016-04-01","2016-05-01",
                       "2016-06-01","2016-07-01","2016-08-01","2016-09-01",
                       "2016-10-01","2016-11-01","2016-12-01","2017-01-01",
                       "2017-02-01","2017-03-01","2017-04-01","2017-05-01",
                       "2017-06-01","2017-07-01","2017-08-01","2017-09-01",
                       "2017-10-01","2017-11-01","2017-12-01","2018-01-01",
                       "2018-02-01","2018-03-01","2018-04-01","2018-05-01",
                       "2018-06-01","2018-07-01","2018-08-01","2018-09-01",
                       "2018-10-01","2018-11-01","2018-12-01","2019-01-01",
                       "2019-02-01","2019-03-01","2019-04-01","2019-05-01",
                       "2019-06-01","2019-07-01","2019-08-01","2019-09-01","2019-10-01",
                       "2019-11-01","2019-12-01","2020-01-01","2020-02-01",
                       "2020-03-01","2020-04-01","2020-05-01","2020-06-01",
                       "2020-07-01","2020-08-01","2020-09-01","2020-10-01",
                       "2020-11-01","2020-12-01","2021-01-01","2021-02-01",
                       "2021-03-01","2021-04-01","2021-05-01","2021-06-01",
                       "2021-07-01","2021-08-01","2021-09-01","2021-10-01",
                       "2021-11-01","2021-12-01")
)

lags_pre <- c(3, 2, 3, 2)
lags_post <- c(2, 1, 3, 2)
events <- c("2018-01", "2017-01", "2014", "2012")
event_type <- c("monthly","monthly", "yearly", "yearly")
event_length <- c(1, 1, 12, 12)

ex <- compute_resilience(x, events, event_type, lags_pre, lags_post, var_interes  = "tmin", avg = "mean")

Bueno espero que te sirva, cualquier cosa hablamos la semana próxima
AntonioJ

[VeruGHub]

Hola Antonio,
He estado estos días con tu código y tengo un par de preguntas para cuando tengas un rato:

• ¿Para qué es necesario tener eventos de diferente resolución temporal? Si quieres un evento de año y otro de meses, ¿se podría poner como c("1990-01:1990-12", "2001-04") para mantener siempre la misma resolución O tiene alguna connotación que no estoy pillando?
• ¿Por qué querrías tener lags diferentes para cada evento de perturbación? ¿No hace que los índices sean incomparables entre eventos?
• Discutimos ayer lo de incluir diferentes maneras de agregar los datos para calcular el performance pre y post pero si los índices están definidos de una manera, decicimos que sería mejor ajustarnos a la forma original (p.e. usar solo la media). ¿Que opinas?
  Siento que no pudieras estar ayer! Abrazos