fototermal.R

# https://itsalocke.com/blog/understanding-rolling-calculations-in-r/
# https://inta.gob.ar/sites/default/files/inta_impacto_de_los_factores_ambientales_en_la_definicion_de_los_rendimientos_de_los_cultivos_resultados_de_la_campana_2013_de_trigo.pdf
# http://rafaela.inta.gov.ar/info/miscelaneas/101/trigo2004_n1.pdf 

pacman::p_load(tidyverse, lubridate, zoo)

stations <- tibble::tribble(
              ~station,   ~lat,   ~lon,
       "Balcarce INTA", -37.75,  -58.3,
  "Mar del Plata AERO", -37.93, -57.58,
           "Azul AERO", -36.83, -59.88,
  "Benito Juarez AERO", -37.72, -59.78,
             "Laprida", -37.57, -60.77,
         "Barrow INTA", -38.32, -60.25,
         "TANDIL AERO", -37.23, -59.25
  ) %>% arrange(station)

dates <- seq.Date(lubridate::dmy("1-8-2019"), lubridate::dmy("30-12-2019"), by = "1 day")
temp <- metR::GetSMNData(dates, type = "daily", bar = TRUE)

rad <- nasapower::get_power(
  community = "AG",
  lonlat =  as.vector(unlist(stations[7,-1])),
  pars = c("ALLSKY_SFC_SW_DWN"),
  dates = c(min(dates), max(dates)),
  temporal_average = "DAILY"
) 

tandil <- rad %>% 
  rename(date = YYYYMMDD, rad = ALLSKY_SFC_SW_DWN) %>% 
  left_join(temp %>% 
              filter(station == "TANDIL AERO") %>% 
              mutate(date = as.Date(date))) %>% 
  # mutate(rfa = rad*0.5, 
  #       q = rfa/((tmax + tmin)/2)) %>% 
  select(date, tmin, tmax, rad) 

tandil %>% tail

tandil %>% 
  ggplot(aes(
    # date, 
    seq(1, length(rad)),
    rad))+ 
  geom_line()+ geom_rug(sides="b") -> p1
p1

filtro <- tibble(y = pracma::hampel(tandil$rad, 5, 3)$y)

p1 + geom_point(data = filtro, 
                aes(x = seq(1, length(y)), y = y), 
                col = "darkred")

tandil <- tandil %>% 
  mutate(tmean = (tmax + tmin)/2, 
         rad_y = filtro$y, 
         q = rad_y*0.5/tmean-4.5) 
           
tandil %>% 
  ggplot(aes(date, q))+ 
  # geom_line()+
  geom_quantile(quantiles = c(0.5),
                formula = NULL,
                method = "rqss") 

# geom_line(data = . %>%
#               mutate(q = zoo::rollmean(q, 2, align = "right", fill = NA)),
#             col = "darkred")

tandil %>% 
  filter(month(date) %in%  10:12) %>%
  group_by(date=if_else(day(date) >= 30,
                        floor_date(date, "20 days"),
                        floor_date(date, "10 days"))) %>%
  summarize(q10 = mean(q),
            days = n()) ->tan_q

tandil <- 
  tandil %>% 
  # mutate(q_pc = rollapply(q, width = list(-20:10), mean, 
  #                             align = "center", 
  #                             fill = NA, 
  #                             na.rm = T)) 
  mutate(tmean_pc = rollapply(tmean, width = list(-20:10), mean, 
                              align = "center", 
                              fill = NA, 
                              na.rm = T)-4.5, 
         rad_pc = rollapply(rad_y*0.5, 
                            width = list(-20:10), mean, 
                            align = "center", 
                            fill = NA, 
                            na.rm = T), 
         q = rad_pc/tmean_pc)


# saveRDS(nasa, here::here("meteo", paste0(id_trial,"_nasa.rds")))

tandil %>% 
  rename("q10" = q) %>% 
  ggplot(aes(date, q10))+ 
  geom_line()+
  # geom_quantile(quantiles = c(0.5),
  #               formula = NULL,
  #               method = "rqss") +
  geom_point(data = tan_q, aes(x = date, y = q10))