df_env <- import("https://bit.ly/loc_info_wheat")
ENV <- df_env$ENV |> as_character()
LAT <- df_env$LAT
LON <- df_env$LON
ALT <- df_env$ALT
START <- df_env$START
END <- df_env$END

# see more at https://github.com/allogamous/EnvRtype
df_climate <-
  get_weather(env.id = ENV,
              lat = LAT,
              lon = LON,
              start.day = START,
              end.day = END)
# GDD: Growing Degree Day (oC/day)
# FRUE: Effect of temperature on radiation use efficiency (from 0 to 1)
# T2M_RANGE: Daily Temperature Range (oC day)
# SPV: Slope of saturation vapour pressure curve (kPa.Celsius)
# VPD: Vapour pressure deficit (kPa)
# ETP: Potential Evapotranspiration (mm.day)
# PEPT: Deficit by Precipitation (mm.day)
# n: Actual duration of sunshine (hour)
# N: Daylight hours (hour)
# RTA: Extraterrestrial radiation (MJ/m^2/day)
# SRAD: Solar radiation (MJ/m^2/day)
# T2M: Temperature at 2 Meters
# T2M_MAX: Maximum Temperature at 2 Meters
# T2M_MIN: Minimum Temperature at 2 Meters
# PRECTOT: Precipitation
# WS2M: Wind Speed at 2 Meters
# RH2M: Relative Humidity at 2 Meters
# T2MDEW: Dew/Frost Point at 2 Meters
# ALLSKY_SFC_LW_DWN: Downward Thermal Infrared (Longwave) Radiative Flux
# ALLSKY_SFC_SW_DWN: All Sky Insolation Incident on a Horizontal Surface
# ALLSKY_TOA_SW_DWN: Top-of-atmosphere Insolation
# [1] "env"       "ETP"       "GDD"         "PETP"      "RH2M"      "SPV"
# [8] "T2M"       "T2M_MAX"   "T2M_MIN"   "T2M_RANGE" "T2MDEW"    "VPD"

# Compute other parameters
env_data <-
  df_climate %>%
  as.data.frame() %>%
  param_temperature(Tbase1 = 5, # choose the base temperature here
                    Tbase2 = 28, # choose the base temperature here
                    merge = TRUE) %>%
  param_atmospheric(merge = TRUE) %>%
  param_radiation(merge = TRUE)

env_data <- import("https://bit.ly/df_climate_tidy")
str(env_data)
## 'data.frame':    991 obs. of  23 variables:
##  $ env          : chr  "E1" "E1" "E1" "E1" ...
##  $ LON          : num  -46.1 -46.1 -46.1 -46.1 -46.1 ...
##  $ LAT          : num  -19.4 -19.4 -19.4 -19.4 -19.4 ...
##  $ YEAR         : int  2018 2018 2018 2018 2018 2018 2018 2018 2018 2018 ...
##  $ MM           : int  5 5 5 5 5 5 5 5 5 5 ...
##  $ DD           : int  2 3 4 5 6 7 8 9 10 11 ...
##  $ DOY          : int  122 123 124 125 126 127 128 129 130 131 ...
##  $ YYYYMMDD     : chr  "05/02/2018" "05/03/2018" "05/04/2018" "05/05/2018" ...
##  $ daysFromStart: int  1 2 3 4 5 6 7 8 9 10 ...
##  $ tmean        : num  19.5 20 20.8 20.9 21.2 ...
##  $ tmax         : num  27.9 28.5 29.1 28.8 28.3 ...
##  $ tmin         : num  12.4 11.9 13.5 14.4 15 ...
##  $ prec         : num  0 0 0 0.03 0.15 0.09 0 0.02 0.02 0 ...
##  $ ws           : num  1.81 1.13 1.5 1.74 1.5 1.96 2.19 2.66 2.7 1.71 ...
##  $ rh           : num  64.3 67.7 65.8 71.6 70.3 ...
##  $ tdew         : num  11.5 12.9 13.2 14.8 14.9 ...
##  $ gdd          : num  15.1 14.9 15.8 16.2 16.5 ...
##  $ frue         : num  0.721 0.723 0.777 0.79 0.793 ...
##  $ trange       : num  15.5 16.6 15.6 14.5 13.3 ...
##  $ vpd          : num  1.24 1.15 1.27 1.11 1.08 ...
##  $ spv          : num  0.146 0.146 0.156 0.158 0.158 ...
##  $ etp          : num  9.34 8.71 9.12 9.09 8.05 ...
##  $ rta          : num  28.7 28.6 28.4 28.3 28.1 ...
id_var <- names(env_data)[10:23]

mod <-
  waasb(df_traits,
        env = ENV,
        gen = LINHAGEM,
        rep = BLOCO,
        resp = GY,
        wresp = 65) # maior peso para performance
## Evaluating trait GY |============================================| 100% 00:00:02 
## ---------------------------------------------------------------------------
## P-values for Likelihood Ratio Test of the analyzed traits
## ---------------------------------------------------------------------------
##     model       GY
##  COMPLETE       NA
##       GEN 3.42e-03
##   GEN:ENV 1.60e-12
## ---------------------------------------------------------------------------
## All variables with significant (p < 0.05) genotype-vs-environment interaction

waasb_env <-
  mod$GY$model %>%
  select_cols(type, Code, Y, WAASB) %>%
  subset(type == "ENV") %>%
  remove_cols(type) %>%
  rename(env = Code)

blupge <- gmd(mod, "blupge") 

blupge |> 
  make_mat(GEN, ENV, GY) |> 
  kableExtra::kable()

indexes <- blup_indexes(mod)
kableExtra::kable(indexes$GY)

mod_gge <- gge(blupge, ENV, GEN, GY, svp = 1)
p1 <- 
  plot(mod_gge,
       size.text.gen = 2.5,
       size.text.env = 2.5) +
  my_theme

p2 <- 
  plot(mod_gge,
       type = 2,
       size.text.gen = 2.5,
       size.text.env = 2.5) +
  my_theme

p3 <- 
  plot(mod_gge,
       type = 3,
       size.text.gen = 2.5,
       size.text.env = 2.5,
       size.text.win = 3.5) +
  my_theme

p4 <-
  plot(mod_gge,
       type = 4,
       size.text.gen = 2.5,
       size.text.env = 2.5,
       size.text.win = 3.5) +
  my_theme

p5 <-
  plot(mod_gge,
       type = 6,
       size.text.gen = 2.5,
       size.text.env = 2.5,
       size.text.win = 3.5) +
  my_theme

p6 <-
  plot(mod_gge,
       type = 8,
       size.text.gen = 2.5,
       size.text.env = 2.5,
       size.text.win = 3.5) +
  my_theme

arrange_ggplot(p1, p2, p3, p4, p5, p6,
               ncol = 3,
               tag_levels = "a",
               guides = "collect")

ggsave("figs/fig5_gge.png",
       width = 12,
       height = 9)

plot_waasby(mod, size.tex.lab = 6) +
  my_theme +
  theme(legend.position = c(0.8, 0.1))

ggsave("figs/fig6_waasby.png",
       width = 5,
       height = 6)