#############################################################################
# This script contains the calibration of weahter effects as well as 
# the extaction of 10 years of average weekly absolute temperature and
# relative humidity for selected locations.
#############################################################################

# Start with a clean session ----
rm(list=ls())

# Import libraries and functions ----
library(tidyverse)
library(doFuture)
library(nasapower)
library(maps)
source("f-compute_beta_var.r")
source("f-optimise_for_location.r")

# First import dataset with cities' coordinates ----
data(world.cities)

# Calibration with one single effect for both variables ----
## Calibration of weather on specific locations using the results from Yaari et al., 2013, with identified variability of 15% -----
Tel_Aviv <- data.frame(loc = "Tel_Aviv-Yafo, Israel",
  pop = 471000,
  lat = world.cities |> dplyr::filter(name == "Tel Aviv-Yafo" & country.etc == "Israel") |> dplyr::pull(lat),
  long = world.cities |> dplyr::filter(name == "Tel Aviv-Yafo" & country.etc == "Israel") |> dplyr::pull(lat),
  capital = 1,
  CZ = "Csa",
  hemisphere = "N")

Tel_Aviv_climate = nasapower::get_power(
  community = "ag",
  lonlat = c(Tel_Aviv$long, Tel_Aviv$lat),
  pars = c("RH2M", "T2M"),
  dates = c("1998-07-01", "2009-06-30"), # Using the same time period of the original paper
  temporal_api = "daily"
) |> 
  dplyr::mutate(week_no = week(as.Date(YYYYMMDD)),
                loc = Tel_Aviv$loc) |> 
  dplyr::select(-YEAR, -MM, -DD, -DOY) |> 
  dplyr::mutate(yr = year(as.Date(YYYYMMDD)), T2M = T2M + 273.15) |> 
  dplyr::group_by(loc, yr, week_no) |> 
  dplyr::summarise(dplyr::across(c(RH2M, T2M), \(x) mean(x, na.rm = TRUE))) |> 
  dplyr::ungroup() |> 
  dplyr::group_by(loc) |> 
  dplyr::mutate(week_no = row_number() - 1,
                RH_norm = RH2M - mean(RH2M),
                Te_norm = T2M - mean(T2M),
                Te_norm_Lag = lag(Te_norm, order_by = week_no),
                RH_norm_Lag = lag(RH_norm, order_by = week_no)) |> 
  dplyr::rename(Te = T2M, RH = RH2M) |> ungroup()

optimised_params_Picked.15 = optimise_for_location(location_name = Tel_Aviv$loc, df = Tel_Aviv_climate, objective_var = 15)

## Calibration for the upper bound of variability to 37.5% according to Weber et al., 2001. Taking as reference climate Turku ----
Turku <- data.frame(loc = "Turku, Finland",
  pop = 194391,
  lat = world.cities |> dplyr::filter(name == "Turku" & country.etc == "Finland") |> dplyr::pull(lat),
  long = world.cities |> dplyr::filter(name == "Turku" & country.etc == "Finland") |> dplyr::pull(long),
  capital = 0,
  CZ = "Dfb",
  hemisphere = "N")
  
Turku_climate <- nasapower::get_power(
  community = "ag",
  lonlat = c(Turku$long, Turku$lat),
  pars = c("RH2M", "T2M"),
  dates = c("1981-01-01", "1990-03-01"), # Using the same time period of the original paper
  temporal_api = "daily"
) |> 
  dplyr::mutate(week_no = week(as.Date(YYYYMMDD)),
                loc = Turku$loc) |> 
  dplyr::select(-YEAR, -MM, -DD, -DOY) |> 
  dplyr::mutate(yr = year(as.Date(YYYYMMDD)), T2M = T2M + 273.15) |> 
  dplyr::group_by(loc, yr, week_no) |> 
  dplyr::summarise(dplyr::across(c(RH2M, T2M), \(x) mean(x, na.rm = TRUE))) |> 
  dplyr::ungroup() |> 
  dplyr::group_by(loc) |> 
  dplyr::mutate(week_no = row_number() - 1,
                RH_norm = RH2M - mean(RH2M),
                Te_norm = T2M - mean(T2M),
                Te_norm_Lag = lag(Te_norm, order_by = week_no),
                RH_norm_Lag = lag(RH_norm, order_by = week_no)) |> 
  dplyr::rename(Te = T2M, RH = RH2M) |> ungroup()

optimised_params_Picked.375 <- optimise_for_location(location_name = Turku$loc, df = Turku_climate, objective_var = 37.5)

# Calibration of two effect, one per the concomitant and lagged effect of variables ----
## First for variability of 15 ----
optim_Lagged.15 <- optimise_for_location_TWO_Deltas(location_name = Tel_Aviv$loc, df = Tel_Aviv_climate, objective_var = 15)

## Second for variability of 37.5% ----
optim_Lagged.375 <- optimise_for_location_TWO_Deltas(location_name = Turku$loc, df = Turku_climate, objective_var = 37.5)

# Downloading weather data of 4 picked locations ----
picked_Locations <- rbind(
  data.frame(
    loc = "Dubai, United Arab Emirates",
    pop = 1182439,
    lat = world.cities |> dplyr::filter(name == "Dubai" & country.etc == "United Arab Emirates") |> dplyr::pull(lat),
    long = world.cities |> dplyr::filter(name == "Dubai" & country.etc == "United Arab Emirates") |> dplyr::pull(long),
    capital = 0,
    CZ = "BWh",
    hemisphere = "N",
    stringsAsFactors = FALSE
  ),
  data.frame(
    loc = "Rio_de_Janeiro, Brazil",
    pop = 6055582,
    lat = world.cities |> dplyr::filter(name == "Rio de Janeiro" & country.etc == "Brazil") |> dplyr::pull(lat),
    long = world.cities |> dplyr::filter(name == "Rio de Janeiro" & country.etc == "Brazil") |> dplyr::pull(long),
    capital = 0,
    CZ = "Aw",
    hemisphere = "S",
    stringsAsFactors = FALSE
  ),
  data.frame(
    loc = "Rome, Italy",
    pop = 2561181,
    lat = world.cities |> dplyr::filter(name == "Rome" & country.etc=="Italy") |> dplyr::pull(lat),
    long = world.cities |> dplyr::filter(name == "Rome" & country.etc=="Italy") |> dplyr::pull(long),
    capital = 1,
    CZ = "Csa",
    hemisphere = "N",
    stringsAsFactors = FALSE
  ),
  data.frame(
    loc = "Toronto, Canada",
    pop = 3026000,
    lat = world.cities |> dplyr::filter(name == "Toronto" & country.etc=="Canada") |> dplyr::pull(lat),
    long = world.cities |> dplyr::filter(name == "Toronto" & country.etc=="Canada") |> dplyr::pull(long),
    capital = 0,
    CZ = "Dfb",
    hemisphere = "N",
    stringsAsFactors = FALSE
  )
)

# If necessary, create and save weather data time series for selected locations
out_file <- "weather_data/picked_Locations_weather.rds"

if (file.exists(out_file) && file.access(out_file, mode = 4) == 0) {

  message("Reading weather data")
  picked_Locations_weather <- readRDS(out_file)

} else {

  message("File not found — pulling weather data and creating file")

  ## Parallel query for weather data
  n.cores <- nrow(picked_Locations)
  future::plan(multisession, workers = n.cores)

  tmp <- foreach(
    idx = seq_len(nrow(picked_Locations)),
    .combine = "comb_rbind",
    .options.future = list(seed = TRUE)
  ) %dofuture% {

    roi <- picked_Locations[idx, ] |> dplyr::select(loc, lat, long)

    final_clim <- nasapower::get_power(
      community = "ag",
      lonlat = c(roi$long, roi$lat),
      pars = c("RH2M", "T2M"),
      dates = c("2011-01-01", "2020-12-31"),
      temporal_api = "daily"
    ) |> 
     dplyr::mutate(
  date = as.Date(YYYYMMDD),
  loc = roi$loc
) |> 
dplyr::select(-YEAR, -MM, -DD, -DOY) |> 
dplyr::mutate(
  yr = year(date),
  week_in_year = week(date),
  T2M = T2M + 273.15
) |> 
dplyr::filter(week_in_year <= 52) |>
dplyr::group_by(loc, yr, week_in_year) |> 
dplyr::summarise(
  dplyr::across(c(RH2M, T2M), \(x) mean(x, na.rm = TRUE)),
  .groups = "drop"
) |> 
dplyr::arrange(loc, yr, week_in_year) |>
dplyr::group_by(loc) |> 
dplyr::mutate(
  week_no = row_number() - 1,
  RH_norm = RH2M - mean(RH2M),
  Te_norm = T2M - mean(T2M),
  Te_norm_lag = lag(Te_norm),
  RH_norm_lag = lag(RH_norm)
) |> 
dplyr::rename(Te = T2M, RH = RH2M) |> 
dplyr::ungroup()

    cat("Iteration:", idx, "for", roi$loc, "\n")

    list(final_clim)
  }

  future::plan(sequential)

  picked_Locations_weather <- tmp[[1]]

  ## Ensure directory exists
  dir.create(dirname(out_file), recursive = TRUE, showWarnings = FALSE)

  saveRDS(picked_Locations_weather, out_file)
}