preanalysis.r

## Paper "Climate shift", R. Fillon, M. Linsenmeier, G. Wagner #####################
####################################################################################
####################################################################################
### This file uses large temperature outputs to create temperature distributions ###
####################################################################################
####################################################################################
#This dataset has three major steps
#1 create file a for geographical weight (KOPPEN & DOSE-KOPPEN levels)
#2 create SSP distribution of 1° temperatures for 2050 and save mean temperature 2050
#3 create synthetic distribution of 1° temperatures from historical (rounded 1°)
#4 create historical distribution of 1° temperatures from historical 

libraries=c("ncdf4","ggpubr","lwgeom","gam","raster","caret","sys","sf","terra","stats","exactextractr","spaMM","metR","rgdal","dplyr","tidyr","base","geifer","data.table","hrbrthemes","pROC","DescTools","rgeos","maptools")
lapply(libraries, require, character.only = TRUE) #load libraries
shape_file="KOPPEN5" # 5 climate zones Koppen
list_SSP=c("ssp126","ssp370","ssp585") #3 SSP
floor_value=0.01 #meaning, we work at 0.001° 

subsample_run=0
subsample_creation=0

#list of climate models
list_huge2=c("UK","MRI","IPSL","MPI","GFDL")

#set working directory
file_path_reproducibility <- file.path("")

setwd("")
mydirection <- getwd()

####################################################################################
####################################################################################
## STEP 1 - LOAD GENERAL DATASETS AND SHAPEFILES  ##################################
####################################################################################
####################################################################################

if (shape_file=="KOPPEN5"){
boundaries <- rast(file.path(mydirection,"/data_additional/koppen_geiger_0p5.tif")) 
name_regions_complete=c(1:5)
name_regions=c(1:5)
for (k in c(2,3)){boundaries <- subst(boundaries, k, 1)}
for (k in c(5,6,7)){boundaries <- subst(boundaries, k, 4)}
for (k in c(9,10,11,12,13,14,15,16)){boundaries <- subst(boundaries, k, 8)}
for (k in c(18,19,20,21,22,23,24,25,26,27,28)){boundaries <- subst(boundaries, k, 17)}
boundaries <- subst(boundaries, 30, 29)
boundaries = as.polygons(boundaries)
boundaries=st_as_sf(boundaries)}

#Keep only DOSE v2 regions
DOSE_shapefiles <- st_read(file.path(mydirection,"/data_additional/DOSE/geometries_all_simplified.shp"))

resolution <- c(0.5, 0.5)
n_rows <- 360 / resolution[2]  
n_cols <- 360 / resolution[1]  
nc_weight <- raster(nrows = n_rows, ncols = n_cols, xmn = -180, xmx = 180, ymn = -90, ymx = 90)
res(nc_weight) <- resolution
crs(nc_weight) <- "+proj=longlat +datum=WGS84 +no_defs"
z <- exact_extract(nc_weight, DOSE_shapefiles, coverage_area = T, include_xy=TRUE)

sf_use_s2(FALSE)

#We want share of each polygon DOSE of use here in each raster 0.5 cell
c<-data.frame()
for (j in c(1:length(z))){
  z[[j]]$DOSE_region=j
  c=rbind(c,z[[j]])}
c=c[,c("x","y","coverage_area","DOSE_region")]

##Treating each DOSE as a unit, take the distribution of degree-days (gridded) at the DOSE level
c <- c %>%
  group_by(DOSE_region) %>%
  mutate(weight_ID=coverage_area/sum(coverage_area,na.rm=T))

if (!all(st_is_valid(DOSE_shapefiles))) {DOSE_shapefiles <- st_make_valid(DOSE_shapefiles)}

d <- data.frame()
for (koppen in c(1:nrow(boundaries))){
   for (DOSE in c(1:nrow(DOSE_shapefiles))){
     d3=st_intersection(boundaries[koppen,]$geometry, DOSE_shapefiles[DOSE,]$geometry)
      if (!all(st_is_valid(d3))) {
          d3 <- st_make_valid(d3)}
      d2 = st_area(d3)/st_area(DOSE_shapefiles[DOSE,]$geometry)
      d2=as.data.frame(as.numeric(d2[1]))
      d2$DOSE_region=DOSE
      d2$koppen=koppen
      d=rbind(d,d2)}}
colnames(d)=c("share_DOSE_in_koppen","DOSE_region","region")

##################################################################################################
d=d[!is.na(d$share_DOSE_in_koppen),]

#make sure we take as a share (i.e. the sum of individual shares of Koppen region in DOSE region is)
d <- d %>%  
 group_by(DOSE_region) %>%
 mutate(share_DOSE_in_koppen=share_DOSE_in_koppen/sum(share_DOSE_in_koppen,na.rm=TRUE))
a=merge(d,c,by=c("DOSE_region"))

a$weight_final=a$weight_ID * a$share_DOSE_in_koppen
a$weight_DOSElevel=a$weight_final#

a <- a %>% 
  group_by(region) %>%
  mutate(weight_sum=sum(weight_final,na.rm=TRUE))

a$weight_final=a$weight_final/a$weight_sum
a$x=round(a$x,2)
a$y=round(a$y,2)
    
file_path <- file.path(mydirection,"/climate_weighting/a.csv")
write.csv(a, file = file_path, row.names = FALSE)

#################################################
##### Step 2 Define some functions
#################################################

#Change a to prepare for merging and reduce size
a=a[,c("DOSE_region","region","x","y","weight_ID","weight_DOSElevel","weight_final")]
a_dt <- as.data.table(a)[, .(x, y, weight_DOSElevel, DOSE_region, region)]

#create raster mask to drop unused data
r <- rast(ncols = 360/0.5, nrows = 180/0.5, xmin = -180, xmax = 180, ymin = -90, ymax = 90)
values(r) <- NA
cells <- cellFromXY(r, a[, c("x", "y")])
r[cells] <- 1
rm(cells)

#This function quickly computes mean over SSP for each raster [without bin]
#Then function also gives temp distribution for each day [with bin] for SSP
processing_quick <-function(dataset) {
  dataset <- as.data.table(as.data.frame(dataset, xy = TRUE))
  setnames(dataset, c("x", "y", "value"))
  dataset[, value := value - 273.15]
  dataset <- merge(dataset, a_dt, by = c("x", "y"))
  return(dataset)}

#This function takes the SSP data to 1 bin (day level)
#For different levels (Koppen, Koppen Dose)
#For projections or synthetic
function_round_tobin <- function(dataset){
  dataset <- as.data.table(dataset)
  #dataset[, value := floor(value / floor_value) * floor_value]
  dataset[, value := round(value / floor_value) * floor_value]
  #koppen level
  datasetKOPPEN <- dataset[, .(frequency = sum(weight_DOSElevel, na.rm = TRUE)), by = .(region, value, type)]
  datasetKOPPEN <- datasetKOPPEN[, .(region, type, value, frequency)]

  #dose koppen level
  datasetDOSEKOPPEN <- dataset[, .(frequency = sum(weight_DOSElevel, na.rm = TRUE)), by = .(DOSE_region, region, value, type)]
  datasetDOSEKOPPEN <- datasetDOSEKOPPEN[, .(DOSE_region, region, value, type, frequency)]

rm(dataset)
return(list(datasetKOPPEN,datasetDOSEKOPPEN))}

function_round_tobin_synthetic <- function(dataset, modelname, ssp){
    dataset <- as.data.table(dataset)
  #koppen
  datasetKOPPEN <- dataset[, .(frequency = sum(weight_DOSElevel, na.rm = TRUE)), by = .(region, value, type)]
  datasetKOPPEN <- datasetKOPPEN[, .(region, type, value, frequency)]

  #dose koppen
  datasetDOSEKOPPEN <- dataset[, .(frequency = sum(weight_DOSElevel, na.rm = TRUE)), by = .(DOSE_region, region, value, type)]
  datasetDOSEKOPPEN <- datasetDOSEKOPPEN[, .(DOSE_region, region, value, type, frequency)]

  # Synthetic adjustments
  datasetDOSEKOPPEN1 <- copy(datasetDOSEKOPPEN)
  datasetKOPPEN1 <- copy(datasetKOPPEN)
 
  datasetDOSEKOPPEN1[, value := value + get(paste0(model, "_mean_", ssp))[[1]] - get(paste0("picontrol_", model, "_mean"))[[1]]]
  datasetKOPPEN1[, value := value + get(paste0(model, "_mean_", ssp))[[1]] - get(paste0("picontrol_", model, "_mean"))[[1]]]
  datasetKOPPEN1[, type := paste0("sm_", ssp)]
  datasetDOSEKOPPEN1[, type := paste0("sm_", ssp)]

  datasetKOPPEN[, type := "picontrol"]
  datasetDOSEKOPPEN[, type := "picontrol"]

  datasetKOPPEN <- rbind(datasetKOPPEN, datasetKOPPEN1)
  datasetDOSEKOPPEN <- rbind(datasetDOSEKOPPEN, datasetDOSEKOPPEN1)

  # Binning final
  datasetKOPPEN[, value := round(value / floor_value) * floor_value]
  datasetDOSEKOPPEN[, value := round(value / floor_value) * floor_value]


  # Recalcul des fréquences après shift
  datasetKOPPEN <- datasetKOPPEN[, .(frequency = sum(frequency, na.rm = TRUE)), by = .(region, value, type)]
  datasetKOPPEN <- datasetKOPPEN[, .(region, type, value, frequency)]

  datasetDOSEKOPPEN <- datasetDOSEKOPPEN[, .(frequency = sum(frequency, na.rm = TRUE)), by = .(DOSE_region, region, value, type)]
  datasetDOSEKOPPEN <- datasetDOSEKOPPEN[, .(DOSE_region, region, value, type, frequency)]

  # Ajout du modèle et du SSP
  datasetKOPPEN[, `:=`(model = modelname, SSP = ssp)]
  datasetDOSEKOPPEN[, `:=`(model = modelname, SSP = ssp)]
  rm(dataset, datasetDOSEKOPPEN1, datasetDOSEKOPPEN2, datasetKOPPEN1, datasetKOPPEN2)
return(list(datasetKOPPEN,datasetDOSEKOPPEN))}

# This loop creates complete distribution [1°]
# Loop over days in 2045-2055
# This loop also saves complete mean [no bin]
# For projections
process_loop <- function(dataset, modelname, typename, sspname) {
  dataset <- as.data.frame(dataset, xy = TRUE)
  setDT(dataset)
  value_cols <- setdiff(names(dataset), c("x", "y"))
  if (length(value_cols) == 0) stop("Pas de couches raster dans le dataset")
  dataset <- merge(dataset, a_dt[, .(x, y, weight_DOSElevel, region, DOSE_region)], by = c("x", "y"))
  dataset <- melt(dataset,
                  id.vars = c("x", "y", "weight_DOSElevel", "region", "DOSE_region"),
                  measure.vars = value_cols,
                  variable.name = "time",
                  value.name = "value")
  dataset[, c("x", "y", "time") := NULL]
  dataset[, value := value - 273.15]
  weight_total <- dataset[, sum(value * weight_DOSElevel, na.rm = TRUE)]
  total <- dataset[, sum(weight_DOSElevel, na.rm = TRUE)]
  dataset[, type := typename]
  rounded_list <- function_round_tobin(dataset)
  dataKOPPEN <- as.data.table(rounded_list[[1]])
  dataDOSEKOPPEN <- as.data.table(rounded_list[[2]])
  dataKOPPEN[, SSP := sspname]
  dataDOSEKOPPEN[, SSP := sspname]
  return(list(dataKOPPEN, dataDOSEKOPPEN, weight_total, total))
}

#We loop over SSP
for (SSP in list_SSP){

#Load SSP datasets
GFDL_4050=rast(file.path(mydirection, paste0("/SSP/gfdl-esm4_r1i1p1f1_w5e5_",SSP,"_tas_global_daily_2041_2050.nc")))
GFDL_5060=rast(file.path(mydirection, paste0("SSP/gfdl-esm4_r1i1p1f1_w5e5_",SSP,"_tas_global_daily_2051_2060.nc")))
IPSL_4050=rast(file.path(mydirection, paste0("SSP/ipsl-cm6a-lr_r1i1p1f1_w5e5_",SSP,"_tas_global_daily_2041_2050.nc")))
IPSL_5060=rast(file.path(mydirection, paste0("SSP/ipsl-cm6a-lr_r1i1p1f1_w5e5_",SSP,"_tas_global_daily_2051_2060.nc")))
UK_4050=rast(file.path(mydirection, paste0("SSP/ukesm1-0-ll_r1i1p1f2_w5e5_",SSP,"_tas_global_daily_2041_2050.nc")))
UK_5060=rast(file.path(mydirection, paste0("SSP/ukesm1-0-ll_r1i1p1f2_w5e5_",SSP,"_tas_global_daily_2051_2060.nc")))
MPI_4050=rast(file.path(mydirection, paste0("SSP/mpi-esm1-2-hr_r1i1p1f1_w5e5_",SSP,"_tas_global_daily_2041_2050.nc")))
MPI_5060=rast(file.path(mydirection, paste0("SSP/mpi-esm1-2-hr_r1i1p1f1_w5e5_",SSP,"_tas_global_daily_2051_2060.nc")))
MRI_4050=rast(file.path(mydirection, paste0("SSP/mri-esm2-0_r1i1p1f1_w5e5_",SSP,"_tas_global_daily_2041_2050.nc")))
MRI_5060=rast(file.path(mydirection, paste0("SSP/mri-esm2-0_r1i1p1f1_w5e5_",SSP,"_tas_global_daily_2051_2060.nc")))

MPI_4050=MPI_4050[[(dim(MPI_4050)[3]-6*365):(dim(MPI_4050)[3])]] #1 bissextile
MPI_5060=MPI_5060[[1:4*365+1]] #1 bissextile
MPI=c(MPI_4050,MPI_5060)
MPI=MPI*r 

UK_4050=UK_4050[[(dim(UK_4050)[3]-6*365):(dim(UK_4050)[3])]] #1 bissextile
UK_5060=UK_5060[[1:4*365+1]] #1 bissextile
UK=c(UK_4050,UK_5060)
UK=UK*r

MRI_4050=MRI_4050[[(dim(MRI_4050)[3]-6*365):(dim(MRI_4050)[3])]] #1 bissextile
MRI_5060=MRI_5060[[1:4*365+1]] #1 bissextile
MRI=c(MRI_4050,MRI_5060)
MRI=MRI*r 

GFDL_4050=GFDL_4050[[(dim(GFDL_4050)[3]-6*365):(dim(GFDL_4050)[3])]] #1 bissextile
GFDL_5060=GFDL_5060[[1:4*365+1]] #1 bissextile
GFDL=c(GFDL_4050,GFDL_5060)
GFDL=GFDL*r

IPSL_4050=IPSL_4050[[(dim(IPSL_4050)[3]-6*365):(dim(IPSL_4050)[3])]] #1 bissextile
IPSL_5060=IPSL_5060[[1:4*365+1]] #1 bissextile
IPSL=c(IPSL_4050,IPSL_5060)
IPSL=IPSL*r

#first wekk of january 2050
if (subsample_creation == 1) {
  select_january_2050 <- function(dataset) {
    layer_names <- time(dataset)
    layer_dates <- as.Date(layer_names)
    start_date <- as.Date("2050-01-01")
    end_date <- as.Date("2050-01-07")
    first_week_indices <- which(layer_dates >= start_date & layer_dates <= end_date)
    dataset <- dataset[[first_week_indices]]
    return(dataset)}
  GFDL_2050 <- select_january_2050(GFDL_4050)*r
  IPSL_2050 <- select_january_2050(IPSL_4050)*r
  UK_2050 <- select_january_2050(UK_4050)*r
  MPI_2050 <- select_january_2050(MPI_4050)*r
  MRI_2050 <- select_january_2050(MRI_4050)*r

writeRaster(GFDL_2050, file.path(file_path_reproducibility,paste0("GFDL_",SSP,"_subsample.tif")),  overwrite=TRUE)
writeRaster(IPSL_2050, file.path(file_path_reproducibility,paste0("IPSL_",SSP,"_subsample.tif")),  overwrite=TRUE)
writeRaster(UK_2050, file.path(file_path_reproducibility,paste0("UK_",SSP,"_subsample.tif")),  overwrite=TRUE)
writeRaster(MPI_2050, file.path(file_path_reproducibility,paste0("MPI_",SSP,"_subsample.tif")),  overwrite=TRUE)
writeRaster(MRI_2050, file.path(file_path_reproducibility,paste0("MRI_",SSP,"_subsample.tif")),  overwrite=TRUE)}

if (subsample_run == 1) {
  IPSL=rast(file.path(file_path_reproducibility,paste0("IPSL_",SSP,"_subsample.tif")))
  GFDL=rast(file.path(file_path_reproducibility,paste0("GFDL_",SSP,"_subsample.tif")))
  UK=rast(file.path(file_path_reproducibility,paste0("UK_",SSP,"_subsample.tif")))
  MRI=rast(file.path(file_path_reproducibility,paste0("MPI_",SSP,"_subsample.tif")))
  MPI=rast(file.path(file_path_reproducibility,paste0("MRI_",SSP,"_subsample.tif")))
}

#This loop computes mean for each projection
#This loop also computes distribution for projection 
for (model in list_huge2){
  dataset=get(model)
  dimension=dim(dataset)[3]
      midpoint <- ceiling(nlyr(dataset)/ 2)
      dataseta=process_loop(dataset[[1:midpoint]],model,"projections",SSP)
      datasetb=process_loop(dataset[[(midpoint + 1):nlyr(dataset)]],model,"projections",SSP)
      dataset1 <- rbindlist(list(as.data.table(dataseta[[1]]),as.data.table(datasetb[[1]])))
      dataset2 <- rbindlist(list(as.data.table(dataseta[[2]]),as.data.table(datasetb[[2]])))
      dataset3=(dataseta[3][[1]]+datasetb[3][[1]])/(dataseta[4][[1]]+datasetb[4][[1]])
      print("dataset3")
      print(dataset3)

  saveRDS(dataset1, file = file.path(mydirection,paste0("/climate_outputs_pretreated/",model,"_KOPPEN_",shape_file,"_",SSP,".rds")))
  saveRDS(dataset2, file = file.path(mydirection,paste0("/climate_outputs_pretreated/",model,"_DOSEKOPPEN_",shape_file,"_",SSP,".rds")))
  assign(paste0(model,"_KOPPEN"),dataset[1])
  assign(paste0(model,"_DOSEKOPPEN"),dataset[2])
  assign(paste0(model,"_mean_",SSP),dataset3)
rm(dataset,dataseta,datasetb,dataset1,dataset2)}}

ls()
rm(dataset3, GFDL,GFDL_4050,GFDL_5060, GFDL_KOPPEN, GFDL_DOSEKOPPEN, IPSL, IPSL_4050,IPSL_5060, IPSL_DOSEKOPPEN, IPSL_KOPPEN, MPI, MPI_4050, MPI_5060, MPI_KOPPEN, MPI_DOSEKOPPEN)
rm(MRI, MRI_4050, MRI_5060, MRI_KOPPEN, MRI_DOSEKOPPEN, UK, UK_4050, UK_5060, UK_KOPPEN, UK_DOSEKOPPEN)

                   
#Step 2 

#This function takes projections raster (much bigger)
#And create synthetic (model & general) with raster
#bit complicated because must add DOSE/Koppen/DOSE-Koppen anomaly to large raster without first aggregation
#Either at Koppen, DOSE, DOSE-Koppen levels
#Taking the mean based on coverage_area
#Directly translate to bin 1° to save memory

#data DOSE, Koppen, Dose Koppen (synthetic model, synthetic general, picontrol)
process_loop_picontrol <- function(dataset,modelname,sspname){
  dataKOPPEN <- data.table()
  dataDOSEKOPPEN <- data.table()

  dataset=dataset-273.15
  dataset=as.data.frame(dataset, xy=TRUE)
  setDT(dataset)
  dataset <- melt(dataset, id.vars = c("x", "y"))

  dataset <- merge(dataset, a_dt, by = c("x", "y"), allow.cartesian = TRUE)
  dataset[, c("x", "y") := NULL]

  data_synthKOPPEN <- data.table()
  data_synthDOSEKOPPEN <- data.table()
  
  dataset[, type := "pi"]

  for (SSP in list_SSP) {
    data_synth2 <- function_round_tobin_synthetic(dataset, modelname, SSP)
    data_synthKOPPEN <- rbindlist(list(data_synthKOPPEN, as.data.table(data_synth2[[1]])))
    data_synthDOSEKOPPEN <- rbindlist(list(data_synthDOSEKOPPEN, as.data.table(data_synth2[[2]])))
  }

  dataKOPPEN <- rbindlist(list(dataKOPPEN, data_synthKOPPEN))
  dataDOSEKOPPEN <- rbindlist(list(dataDOSEKOPPEN, data_synthDOSEKOPPEN))
  
  dataKOPPEN[, model := modelname]
  dataDOSEKOPPEN[, model := modelname]
  
  rm(data_synthDOSEKOPPEN, data_synth2, data_synthKOPPEN, dataset)
  return(list(dataKOPPEN,dataDOSEKOPPEN))}


#Large loop to handle many files
#Unfortunately, we cant work without bin + with aggregate data in df format
#While df needed (can't stay to raster with DOSE, Koppen, DOSEKoppen format)
#Thus, trade off between loops (computation time) and memory 
process_loop_picontrol_decades <- function(dataset) {
  
  # Liste des fichiers NetCDF
  file_list <- list.files(
    file.path(mydirection, paste0("/picontrol/picontrol_", dataset, "/")), 
    full.names = TRUE)
  
  # Si subsample_run == 1, on charge un fichier sp�cial
  if (subsample_run == 1) {
    nc <- rast(file.path(file_path_reproducibility, paste0(dataset, "_hist_subsample.tif")))
    file_list <- NULL }

  # Fonction interne : traite un fichier netCDF en chunks
  process_file_chunks <- function(nc_file, trim_start = TRUE) {
    nc_data <- rast(nc_file)
    
    if (trim_start) {
      start_idx <- 5 * 365 + 2
      end_idx <- dim(nc_data)[3]
      nc_data <- nc_data[[start_idx:end_idx]]}

    total_layers <- nlyr(nc_data)

    chunk_size <- 10
    merge_every <- 5

    KOPPEN_list <- list()
    DOSEKOPPEN_list <- list()

    chunk_idx <- 1
    for (i in seq(1, total_layers, by = chunk_size)) {
      chunk_end <- min(i + chunk_size - 1, total_layers)
      chunk <- nc_data[[i:chunk_end]] * r

if (all(is.na(values(chunk)))) {
  cat(paste("Chunk", chunk_idx, "is completely NA\n"))
  next}

      result <- process_loop_picontrol(chunk, model, "picontrol")
if (is.null(result[[1]]) || nrow(result[[1]]) == 0) {
  cat(paste("Chunk", chunk_idx, "-> empty result[[1]]\n"))}

      setcolorder(result[[1]], c("region", "type", "value", "model", "SSP", "frequency"))
      setcolorder(result[[2]], c("DOSE_region", "region", "type", "value", "model", "SSP", "frequency"))

      result[[1]] <- result[[1]][, .(frequency = sum(frequency, na.rm = TRUE)), by = .(region, type, value, model, SSP)]
      result[[2]] <- result[[2]][, .(frequency = sum(frequency, na.rm = TRUE)), by = .(DOSE_region, region, type, value, model, SSP)]

      KOPPEN_list[[chunk_idx]] <- as.data.table(result[[1]])
      DOSEKOPPEN_list[[chunk_idx]] <- as.data.table(result[[2]])
      rm(chunk, result)
      gc()

cat(paste0("== MERGE INTERMEDIAIRE a CHUNK ", chunk_idx - 1, " ==\n"))
for (i in seq_along(KOPPEN_list)) {
  cat(paste0("Chunk ", i, " cols: ", paste(names(KOPPEN_list[[i]]), collapse = ", "), "\n"))}

      if (chunk_idx %% merge_every == 0 || i == total_layers) {
        tempK <- rbindlist(KOPPEN_list)
        tempD <- rbindlist(DOSEKOPPEN_list)
    
if (!exists("ncKOPPEN")) {
  ncKOPPEN <- tempK
} else {
ncKOPPEN <- rbindlist(list(ncKOPPEN, tempK))[
  , .(frequency = sum(frequency, na.rm = TRUE)),
  by = .(region, type, value, model, SSP)
]
}

if (!exists("ncDOSEKOPPEN")) {
  ncDOSEKOPPEN<- tempD
} else {
ncDOSEKOPPEN <- rbindlist(list(ncDOSEKOPPEN, tempD))[
  , .(frequency = sum(frequency, na.rm = TRUE)),
  by = .(DOSE_region, region, type, value, model, SSP)
]
}
        rm(tempK, tempD)
        gc()
       KOPPEN_list <- list()
       DOSEKOPPEN_list <- list()
    	chunk_idx <- 1  # reset index after merge
      } else {chunk_idx <- chunk_idx + 1}
    }

    return(list(ncKOPPEN, ncDOSEKOPPEN))
  }
  if (subsample_run == 1) {
    result <- process_file_chunks(nc_file = nc, trim_start = FALSE)
  } else {
    result <- process_file_chunks(nc_file = file_list[[1]], trim_start = TRUE)
  }

  KOPPEN_combined <- result[[1]]
  DOSEKOPPEN_combined <- result[[2]]

  if (!is.null(file_list) && subsample_run != 1) {
    for (k in 2:length(file_list)) {
      result <- process_file_chunks(file_list[[k]], trim_start = FALSE)

      KOPPEN_combined <- rbindlist(list(KOPPEN_combined, result[[1]]))[
        , .(frequency = sum(frequency, na.rm = TRUE)), by = .(region, type, value, model, SSP)
      ]
      DOSEKOPPEN_combined <- rbindlist(list(DOSEKOPPEN_combined, result[[2]]))[
        , .(frequency = sum(frequency, na.rm = TRUE)), by = .(DOSE_region, region, type, value, model, SSP)
      ]
      gc()
    }
  }
  return(list(KOPPEN_combined, DOSEKOPPEN_combined))
}

#Unfortunately, for synthetic general, we need mean anomaly for all models
#Thus, trade off between loops (computation time) and memory 
process_loop_syntheticgeneral_mean <- function(dataset){
list = as.list(list.files(file.path(mydirection, paste0("/picontrol/picontrol_",dataset,"/"))))
 meandata <- NULL
 chunk_size = 100
 dimension=0
 
if (subsample_creation==1){
      nc=rast(file.path(mydirection, paste0("/picontrol/picontrol_",dataset,"/",list[[1]])))
      layer_names <- time(nc)
  layer_dates <- as.Date(layer_names)  # utile si ce sont des strings ou POSIXct
   start_date <- as.Date("2010-01-01")
  end_date <- as.Date("2010-01-07")
  indices <- which(layer_dates >= start_date & layer_dates <= end_date)
  dataset2 <- nc[[indices]]
    writeRaster(dataset2,
              file.path(file_path_reproducibility, paste0(dataset, "_hist_subsample.tif")),
              overwrite = TRUE)
}

  
  if (subsample_run==0){
  for (k in seq_along(list)){	
    mean=rast(file.path(mydirection, paste0("/picontrol/picontrol_",dataset,"/",list[[k]])))
    if (k == 1) {mean <- mean[[(5 * 365 + 2):(dim(mean)[3])]]}
    dimension=dimension+dim(mean)[3]
    total_layers <- nlyr(mean)
    if (is.null(meandata)) {
        meandata <- mean[[1]] * 0  # Create an empty raster (0 values) with the same extent and resolution
    }

    for (i in seq(1, total_layers, by = chunk_size)) {

      chunk_end <- min(i + chunk_size - 1, total_layers)
      chunk <- mean[[i:chunk_end]]
      chunk <- chunk * r
      chunk <- chunk - 273.15
      chunk_mean <- sum(chunk, na.rm = TRUE)
      meandata <- meandata + chunk_mean
}    
    rm(mean)}
  }

  if (subsample_run==1){
    mean=rast(file.path(file_path_reproducibility,paste0(dataset,"_hist_subsample.tif")))
    dimension=dim(mean)[3]
    mean=mean*r
    mean=mean-273.15
    meandata <- sum(mean, na.rm = TRUE)
  }
  meandata <- meandata / dimension
  meandata <- as.data.table(as.data.frame(meandata, xy = TRUE))
  meandata <- merge(meandata, a_dt, by = c("x", "y"))
  cols_to_melt <- setdiff(names(meandata), c("x", "y", "weight_DOSElevel", "region", "DOSE_region"))
  meandata <- melt(meandata, id.vars = c("x", "y", "weight_DOSElevel", "region", "DOSE_region"), measure.vars = cols_to_melt, value.name = "value")
  meandata[, c("x", "y", "variable") := NULL]
  meandata <- sum(meandata$value * meandata$weight_DOSElevel, na.rm = TRUE) /sum(meandata$weight_DOSElevel, na.rm = TRUE)
  rm(chunk)
  return(meandata)}



for (model in list_huge2){
	mean=process_loop_syntheticgeneral_mean(model)
	print("mean")
	print(mean)
	mean=round(mean/ floor_value) * floor_value 
	print("mean")
	print(mean)
	assign(paste0("picontrol_",model,"_mean"),mean)}


if (subsample_creation==1){stop("subsample_creation = 1 ? Script stopped intentionally")}

#this loop creates distribution 1° synthetic/picontrol
for (model in list_huge2){
    dataset=process_loop_picontrol_decades(model)
    resultKOPPEN=dataset[1]
    resultDOSEKOPPEN=dataset[2]

for (aggregation in c("KOPPEN","DOSEKOPPEN")){
  dataset=get(paste0("result",aggregation))
  if (subsample_run==1){saveRDS(dataset, file = file.path(mydirection, paste0("/climate_outputs_pretreated/",aggregation,model,"_",shape_file,"_picontrol_synthetic_",SSP,"_subsample.rds")))}
  if (subsample_run==0){saveRDS(dataset, file = file.path(mydirection, paste0("/climate_outputs_pretreated/",aggregation,model,"_",shape_file,"_picontrol_synthetic_",SSP,".rds")))}}}