prepare_data_exp.R

# =============================================================================
# IPCC Atlas Data Acquisition & Processing
# =============================================================================
#~/Documents/GitHub_mac/IPCCAtlas2R/scripts/prepare_data.R
# DATA SOURCE:
#   IPCC WGI AR6 Interactive Atlas – Gridded Monthly Dataset (CMIP6)
#   DOI: 10.24381/cds.5292a2b0
#   Hosted on: Copernicus Climate Data Store (CDS)
#   Digital.CSIC handle: https://hdl.handle.net/10261/332744
#
# SETUP (one-time):
#   1. Create a free CDS account: https://cds.climate.copernicus.eu
#   2. Accept the dataset licence:
#      https://cds.climate.copernicus.eu/datasets/projections-climate-atlas
#   3. Copy your UID and API key from https://cds.climate.copernicus.eu/profile
#   4. Run:
# ecmwfr::wf_set_key(
#   key = "PASTE_YOUR_PERSONAL_ACCESS_TOKEN_HERE", # copy paste API Key from here: https://cds.climate.copernicus.eu/profile
#   user = "ecmwfr"
# )
#      (saves credentials to ~/.cdsapirc – only needed once)
# if that doesn't work do the steps listed here: https://cds.climate.copernicus.eu/how-to-api

# PACKAGES:
# install.packages(c("ecmwfr","ncdf4","terra","raster","dplyr","tidyr",
#                    "readr","lubridate","ggplot2","httr2"))
#        install.packages("ecmwfr")
#        ecmwfr::wf_set_key(user = "<UID>", key = "<API_KEY>", service = "cds")
# library(ecmwfr)
#
# ecmwfr::wf_set_key(
#   key = "PASTE_YOUR_PERSONAL_ACCESS_TOKEN_HERE", # copy paste API Key from here: https://cds.climate.copernicus.eu/profile
#   user = "ecmwfr"
# )
#
# ecmwfr::wf_get_key(user = "ecmwfr")

# =============================================================================

suppressPackageStartupMessages({
  library(ecmwfr)
  library(ncdf4)
  library(terra)
  library(raster)
  library(dplyr)
  library(tidyr)
  library(readr)
  library(lubridate)
})

# ── Configuration ─────────────────────────────────────────────────────────────

DATA_DIR <- "YOURCHOSENFOLDERNAME"
dir.create(DATA_DIR, showWarnings = FALSE, recursive = TRUE)

RAST_FN <- "YOURFILENAME"

#Ecoregion bounding box - ENTER YOUR LAT/LON POINTS
BOUNDING_BOX <- list(xmin = , xmax = , ymin = , ymax = )

# CDS dataset identifier
CDS_DATASET <- "projections-climate-atlas"

SCENARIOS <- c("historical", "ssp1_2_6", "ssp2_4_5", "ssp3_7_0", "ssp5_8_5")
SCEN_KEY  <- c(          # map CDS names → short keys used in app
  historical = "historical",
  ssp1_2_6   = "ssp126",
  ssp2_4_5   = "ssp245",
  ssp3_7_0   = "ssp370",
  ssp5_8_5   = "ssp585"
)
SCENARIO_LABELS <- c(
  historical = "Historical",
  ssp126     = "SSP1-2.6",
  ssp245     = "SSP2-4.5",
  ssp370     = "SSP3-7.0",
  ssp585     = "SSP5-8.5"
)

BASELINE_PERIOD <- c(1985, 2014)
FUTURE_NEAR     <- c(2021, 2040)
FUTURE_MID      <- c(2041, 2060)
FUTURE_FAR      <- c(2081, 2100)

# ── CDS credentials check ─────────────────────────────────────────────────────

# check_cds_credentials <- function() {
#   cds_rc <- path.expand("~/.cdsapirc")
#   if (!file.exists(cds_rc)) {
#     stop(
#       "\nCDS API credentials not found.\n",
#       "Please follow the SETUP instructions at the top of this script:\n",
#       "  1. Create account at https://cds.climate.copernicus.eu\n",
#       "  2. Accept dataset licence\n",
#       "  3. Run: ecmwfr::wf_set_key(user='<UID>', key='<KEY>', service='cds')\n"
#     )
#   }
#   # Read UID from ~/.cdsapirc
#   rc <- readLines(cds_rc)
#   uid_line <- grep("^key:", rc, value = TRUE)
#   if (!length(uid_line)) {
#     stop("~/.cdsapirc found but malformed. Re-run wf_set_key().")
#   }
#   invisible(TRUE)
# }
check_cds_credentials <- function() {
  key <- ecmwfr::wf_get_key()
  if (is.null(key) || !nzchar(key)) {
    stop(
      "\nCDS API key not found.\n",
      "Run:\n",
      "  ecmwfr::wf_set_key(key = '<PERSONAL_ACCESS_TOKEN>')\n",
      "Use the token from: https://cds.climate.copernicus.eu/how-to-api\n"
    )
  }
  invisible(TRUE)
}

# ── Download one scenario from CDS ────────────────────────────────────────────

download_cds_sst <- function(scenario) {
  #
  # download_cds_sst <- function(scenario, uid) {
  # Map scenario to CDS period parameter
  # NB: valid period depends on the variable. For sst the historical block is
  # 1850-2014
  period <- if (scenario == "historical") "1850-2014" else "2015-2100"
  short  <- SCEN_KEY[scenario]
  dest   <- file.path(DATA_DIR, sprintf("tos_cmip6_%s.nc", short))

  if (file.exists(dest)) {
    message("  Cached: ", basename(dest))
    return(dest)
  }

  zip_dest <- file.path(DATA_DIR, sprintf("tos_cmip6_%s.zip", short))
  message("  Requesting from CDS: ", scenario, " (", period, ") ...")

  req <- list(
    dataset_short_name = CDS_DATASET,
    origin             = "cmip6",
    experiment         = scenario,
    domain             = "global",
    variable           = "monthly_mean_of_sea_surface_temperature",
    #variable          = 'monthly_mean_of_acidity_of_seawater',
    period             = period,
    target             = basename(zip_dest)
  )

  tryCatch({
    wf_request(
      request  = req,
      path     = DATA_DIR,
      transfer = TRUE,
      verbose  = TRUE
    )
  }, error = function(e) {
    stop(
      "CDS request failed for scenario '", scenario, "'.\n",
      "  CDS message: ", conditionMessage(e), "\n",
      "  The job reached the server and was rejected there. The exact reason is\n",
      "  shown in your request queue: https://cds.climate.copernicus.eu/requests?tab=all\n",
      "  Common causes: an invalid parameter value (variable / period / domain /\n",
      "  experiment), or the dataset licence not yet accepted at:\n",
      "  https://cds.climate.copernicus.eu/datasets/projections-climate-atlas?tab=download",
      call. = FALSE
    )
  })

  # Unzip – CDS returns a zip with one .nc inside
  if (file.exists(zip_dest)) {
    nc_files <- unzip(zip_dest, list = TRUE)$Name
    nc_name  <- nc_files[grepl("\\.nc$", nc_files)][1]
    unzip(zip_dest, files = nc_name, exdir = DATA_DIR)
    file.rename(file.path(DATA_DIR, nc_name), dest)
    file.remove(zip_dest)
  }

  if (!file.exists(dest)) stop("NC file not found after extraction: ", dest)
  dest
}

# ── Clip NC to Bounding Box and extract time series ────────────────────────────

process_nc <- function(nc_path, scenario) {
  message("  Processing: ", basename(nc_path))
  nc <- nc_open(nc_path)
  on.exit(nc_close(nc))

  lon  <- ncvar_get(nc, "lon")   # expect 0–360
  lat  <- ncvar_get(nc, "lat")
  time_raw   <- ncvar_get(nc, "time")
  time_units <- ncatt_get(nc, "time", "units")$value

  # This CDS Atlas file stores the RAW ensemble, not pre-computed stats:
  #   sst : [lon, lat, time, member]   (member = individual CMIP6 models)
  # We compute the ensemble mean/percentiles across the member dimension here.

  if (!"sst" %in% names(nc$var)) {
    stop("Expected variable 'sst' not found in ", basename(nc_path),
         ". Variables present: ", paste(names(nc$var), collapse = ", "))
  }


  # Spatial indices for the ecoregion
  lon_idx <- which(lon >= BOUNDING_BOX$xmin & lon <= BOUNDING_BOX$xmax)
  lat_idx <- which(lat >= BOUNDING_BOX$ymin & lat <= BOUNDING_BOX$ymax)


  if (!length(lon_idx) || !length(lat_idx)) {
    stop("No grid cells found in Ecoregion domain for: ", nc_path)
  }


  # Read & clip the box only (avoids loading the global 4-D array).
  # Result dims: [lon', lat', time, member]
  start <- c(lon_idx[1], lat_idx[1], 1, 1)
  count <- c(length(lon_idx), length(lat_idx), -1, -1)
  sst   <- ncvar_get(nc, "sst", start = start, count = count)
  if (length(dim(sst)) == 3) dim(sst) <- c(dim(sst), 1L)   # guard: single member

  n_lon <- dim(sst)[1]; n_lat <- dim(sst)[2]
  n_time <- dim(sst)[3]; n_member <- dim(sst)[4]

  fill <- ncatt_get(nc, "sst", "_FillValue")$value
  if (!is.null(fill) && !is.na(fill)) sst[sst == fill] <- NA
  # Convert K → °C if needed
  if (mean(sst, na.rm = TRUE) > 200) sst <- sst - 273.15

  # Parse times (assume "<unit> since YYYY-MM-DD ..."; CDS Atlas uses days)
  origin <- regmatches(time_units, regexpr("\\d{4}-\\d{1,2}-\\d{1,2}", time_units))
  origin <- if (length(origin)) as.Date(origin) else as.Date("1850-01-01")
  dates  <- tryCatch(
    as.Date(time_raw, origin = origin),
    error = function(e) origin + as.integer(time_raw)
  )

  lat_crop <- lat[lat_idx]
  lon_crop <- lon[lon_idx]


  # # ── Area-weighted regional mean per (time, member) ──────────────────────────
  # # cosine-latitude weights over the clipped box: [lon', lat'] -> vector
  wt_mat <- outer(rep(1, n_lon), cos(lat_crop * pi / 180))
  wt_vec <- as.vector(wt_mat)
  # collapse space -> rows, (time x member) -> cols, then weighted column means
  sst_mat <- matrix(sst, nrow = n_lon * n_lat, ncol = n_time * n_member)
  num     <- colSums(sst_mat * wt_vec, na.rm = TRUE)
  den     <- colSums((!is.na(sst_mat)) * wt_vec)
  reg     <- matrix(num / den, nrow = n_time, ncol = n_member)   # [time, member]
  reg[!is.finite(reg)] <- NA_real_

  # ── Ensemble statistics across members, per time step ───────────────────────
  ens_q   <- function(p) apply(reg, 1, quantile, probs = p, na.rm = TRUE)
  ts_mean <- rowMeans(reg, na.rm = TRUE)
  ts_p10  <- ens_q(0.10)
  ts_p25  <- ens_q(0.25)
  ts_p75  <- ens_q(0.75)
  ts_p90  <- ens_q(0.90)
  n_models <- sum(apply(reg, 2, function(col) any(!is.na(col))))

  # Annual aggregation
  yrs <- year(dates)
  yr_uniq <- sort(unique(yrs))

  ann_agg <- function(vec) tapply(vec, yrs, mean, na.rm = TRUE)[as.character(yr_uniq)]

  ts_df <- tibble(
    year           = yr_uniq,
    scenario       = SCEN_KEY[scenario],
    mean           = as.numeric(ann_agg(ts_mean)),
    p10            = as.numeric(ann_agg(ts_p10)),
    p25            = as.numeric(ann_agg(ts_p25)),
    p75            = as.numeric(ann_agg(ts_p75)),
    p90            = as.numeric(ann_agg(ts_p90)),
    n_models       = n_models
  )

  # Save clipped spatial raster (annual ensemble mean) for map tab.
  # Ensemble mean over members -> [lon', lat', time], then annual mean per year
  # so every (lon', lat') cell keeps one sst value for each year & scenario.
  ens_spatial <- apply(sst, c(1, 2, 3), mean, na.rm = TRUE)   # [lon', lat', time]
  annual_spatial <- vapply(
    yr_uniq,
    function(y) apply(ens_spatial[, , yrs == y, drop = FALSE], c(1, 2),
                      mean, na.rm = TRUE),
    matrix(0, nrow = n_lon, ncol = n_lat)
  )                                                            # [lon', lat', year]

  # Build a raster with -180:180 lon convention for leaflet.
  # Long form: one row per (lon, lat, year); expand.grid varies lon fastest,
  # then lat, then year -> matches as.vector() order of annual_spatial.
  lon_display <- ifelse(lon_crop > 180, lon_crop - 360, lon_crop)
  rast_df <- expand.grid(lon = lon_display, lat = lat_crop, year = yr_uniq)
  rast_df$scenario <- SCEN_KEY[scenario]
  rast_df$sst <- as.vector(annual_spatial)
  rast_df <- rast_df[!is.na(rast_df$sst), ]

  rast_path <- file.path(DATA_DIR, sprintf("raster_%s.rds", SCEN_KEY[scenario]))
  saveRDS(rast_df, rast_path)
  message("    Saved spatial raster: ", basename(rast_path))

  ts_df
}

# ── Warming summary table ──────────────────────────────────────────────────────

compute_warming_table <- function(ts_data) {
  hist_mean <- ts_data |>
    filter(scenario == "historical",
           year >= BASELINE_PERIOD[1], year <= BASELINE_PERIOD[2]) |>
    summarise(ref = mean(mean, na.rm = TRUE)) |>
    pull(ref)

  periods <- list(
    "Near-term (2021-2040)" = FUTURE_NEAR,
    "Mid-term (2041-2060)"  = FUTURE_MID,
    "Long-term (2081-2100)" = FUTURE_FAR
  )

  bind_rows(lapply(names(periods), function(pname) {
    p <- periods[[pname]]
    ts_data |>
      filter(scenario != "historical", year >= p[1], year <= p[2]) |>
      group_by(scenario) |>
      summarise(
        delta_mean = mean(mean, na.rm = TRUE) - hist_mean,
        delta_p10  = mean(p10,  na.rm = TRUE) - hist_mean,
        delta_p90  = mean(p90,  na.rm = TRUE) - hist_mean,
        .groups = "drop"
      ) |>
      mutate(period = pname)
  })) |>
    mutate(scenario_label = factor(SCENARIO_LABELS[scenario],
                                   levels = SCENARIO_LABELS[-1]))
}

# ── Main pipeline ─────────────────────────────────────────────────────────────

message("=== IPCC Atlas sst Processing: Ecoregion ===\n")
message("Domain: ", BOUNDING_BOX$xmin, "–", BOUNDING_BOX$xmax,
        "°E  /  ", BOUNDING_BOX$ymin, "–", BOUNDING_BOX$ymax, "°N\n")


# Check credentials before starting
check_cds_credentials()

# Read UID from ~/.cdsapirc
# cds_uid <- trimws(
#   sub("key:\\s*", "",
#       grep("^key:", readLines(path.expand("~/.cdsapirc")), value = TRUE)[1])
# )
# Note: ecmwfr uses "key" field as the user ID for CDS API routing
# The actual API key is stored internally by wf_set_key()

all_ts <- list()

for (scen in SCENARIOS) {
  short <- SCEN_KEY[scen]
  message("\n[", which(SCENARIOS == scen), "/", length(SCENARIOS), "] Scenario: ", scen)

  nc_path <- tryCatch(
    #download_cds_sst(scen, cds_uid),
    download_cds_sst(scen),
    error = function(e) {
      message("  ERROR: ", e$message)
      NULL
    }
  )
  if (is.null(nc_path)) next

  ts_df <- tryCatch(
    process_nc(nc_path, scen),
    error = function(e) {
      message("  ERROR processing NC: ", e$message)
      NULL
    }
  )
  if (!is.null(ts_df)) all_ts[[short]] <- ts_df
}

if (length(all_ts) == 0) {
  stop(
    "\nNo data was successfully downloaded.\n",
    "Check CDS credentials, licence acceptance, and internet connectivity.\n"
  )
}

ts_data <- bind_rows(all_ts) |>
  mutate(scenario_label = factor(SCENARIO_LABELS[scenario],
                                 levels = SCENARIO_LABELS))

warming <- compute_warming_table(ts_data)

# Baseline spatial raster (from historical)
hist_raster_path <- file.path(DATA_DIR, "raster_historical.rds")
spatial_baseline <- if (file.exists(hist_raster_path)) {
  readRDS(hist_raster_path)
} else NULL

# Future spatial rasters: combine all SSP scenarios (2015-2100), excl. historical
# Long form: lon, lat, year, scenario, sst
future_scen <- setdiff(names(all_ts), "historical")
spatial_future <- bind_rows(lapply(future_scen, function(short) {
  rp <- file.path(DATA_DIR, sprintf("raster_%s.rds", short))
  if (file.exists(rp)) readRDS(rp) else NULL
}))
if (nrow(spatial_future) == 0) spatial_future <- NULL

sst_data <- list(
  timeseries       = ts_data,
  warming_summary  = warming,
  spatial_baseline = spatial_baseline,
  spatial_future   = spatial_future,
  source           = "IPCC Atlas CMIP6 (Copernicus CDS – doi:10.24381/cds.5292a2b0)"
)


saveRDS(sst_data, file.path(DATA_DIR, RAST_FN))
write_csv(ts_data, file.path(DATA_DIR, "timeseries.csv"))
write_csv(warming, file.path(DATA_DIR, "warming_summary.csv"))

message("\n✓ Done. Scenarios loaded: ", paste(names(all_ts), collapse = ", "))
message("  bering_sst.rds  – main data object")
message("  timeseries.csv  – annual ensemble stats")
message("  warming_summary.csv – delta sst by period & scenario")
message("  raster_*.rds    – clipped spatial rasters for map tab")
message("  sst_data$spatial_future – combined SSP spatial sst (2015-2100)")
message("\nRun the app: shiny::runApp('app.R')")