# =========================================================
# R version 4.3.2 (2023-10-31)  #
# Operating System: arm64    #
# Code for statistical analysis   #
# Supplement to: Gabriel Fénisse*, Manuel Chevalier, Odile Peyron, David Vincent Bekaert, Pierre-Henri Blard: Advancing Last Glacial Maximum paleoclimate reconstructions in Europe using pollen data: a multi-method (mega)biomization approach, EGUsphere [preprint], https://doi.org/10.5194/egusphere-2026-1590, 2026. #
# Contact: *Gabriel Fénisse (gabriel.fenisse@univ-lorraine.fr)[Centre de Recherches Pétrographiques et Géochimiques, Vandoeuvre-les-Nancy, 2025] #
# =========================================================

# Load required libraries
library(ggplot2)
library(plotly)
library(dplyr)
library(sf)
library(rnaturalearth)
library(stringr)

# -------------------------------------------------------
# EMPD2 calibration dataset (Davis et al., 2000)
# and associated site names
# -------------------------------------------------------

# Classification scheme
PA=read.csv2("Taxpftpi_Class1.csv",sep=';',dec=",")
rownames(PA)=PA[,1]
PA=PA[,-1]
PA1=PA
BPFT=read.csv2("biopftpi_Class1.csv",sep=';',dec=".")
Name_Biome=BPFT[,1]
row.names(BPFT)=Name_Biome
BPFT=BPFT[,-1]
BPFT1=BPFT

# -------------------------------------------------------
# Modern dataset preprocessing and taxonomic harmonization
# -------------------------------------------------------

Table=read.csv2("harmonization_table.txt",sep='\t',dec=',')
IntEMPD=read.csv2("IntEMPD2_ModernDatasets.txt",sep='\t')

# Harmonize taxon names across datasets
for (a in Table$acc_varname)
{if (a %in% IntEMPD$acc_varname)
{IntEMPD$acc_varname[IntEMPD$acc_varname==a]=Table$TaxonName[Table$acc_varname==a]}}
# Check taxonomic consistency
setdiff(rownames(PA),unique(IntEMPD$acc_varname))

# -------------------------------------------------------
# Normalisation of calibration dataset (relative abundances)
# -------------------------------------------------------

IntEMPD3=IntEMPD
for(i in 1:(length(IntEMPD[,3])))
{IntEMPD[i,3]=(100*IntEMPD[i,3])/(sum(IntEMPD3[,3][(IntEMPD[,1]==IntEMPD3[i,1])]))}
write.csv2(IntEMPD,"EMPD2_Normalized_v4.csv")

# =======================================================
# MODERN DATA PREPARATION
# =======================================================
# Load processed modern dataset
EMPD=read.csv2("EMPD2_Normalized_v4.csv",sep=';',dec=".")
Loc=read.csv2("PositionSiteCal.csv",sep=';',dec=".")
EMPD=subset(EMPD,EMPD$SampleName %in% Loc$SampleName)
SiteNameCal=Loc$SampleName
NameCal=Loc$SampleName
namefile="OutModern_Biomization2_v5"
# Taxonomic consistency check
setdiff(rownames(PA),unique(EMPD[,2]))

#Lat constraints-> Europe
#Loc<-Loc[Loc[, 'Latitude'] > 0, ]
#Loc<-Loc[Loc[, 'Latitude'] < 90, ]
#Loc<- Loc[Loc[, 'Longitude'] > 0, ]
#Loc<- Loc[Loc[, 'Longitude'] < 90, ]
#AssPol=AssPol[,-c(grep("Delete",Name$Names))]
#Name=Name[-c(grep("Delete",Name$Names)),]

# =======================================================
# FOSSIL DATA PREPARATION
# =======================================================

# Harmonized fossil dataset
# Age model (example LDB)
AgeCal=read.csv2("ModeleAge_LDB.txt",sep='\t',dec='.')
Ages=AgesCal$X2/1E3
#Ages=AssPol[,1]/1E3#Eifel_2

# Age model (example LDB)
AssPol=read.csv2("CoolingPollen.csv",sep=';',dec=",")
Name=read.csv2("Name_CoolingPollen_Harmonized.txt",sep=';',dec=".")
Ages=AssPol$SiteName
AssPol=AssPol[,-c(1)]
AssPol=AssPol*100
colnames(AssPol)=Name$Names

# Taxonomic consistency check
setdiff(names(AssPol),rownames(PA))
namefile="OutCooling_Biomisation_Classe2"

# =======================================================
# FUNCTION: BIOME INDEX (Prentice PFT scheme)
# =======================================================

EuropeBiomeIndexModern<-function(EMPD,NameCal,namefile) {
library(tidyverse)
# Initialise output table
Output=data.frame(matrix(NA,ncol=3,nrow=length(NameCal)))
dimnames(Output)=list(c(NameCal),c("Biome","Score","Std"))
# Progress bar
pb <- txtProgressBar(min = 0, 
                       max = length(NameCal),
                       style = 3, 
                       width = 50,
                       char = "=")   
n=1
for (a in NameCal)
{PA=PA1
PA=PA[which(row.names(PA) %in% EMPD[,2][EMPD[,1]==a]),] 
# Adjust PFT contributions based on observed abundances
for (i in EMPD[,2][EMPD[,1]==a])
{PA[(row.names(PA)==i)&(PA>0)&(PA>EMPD[,3][(EMPD[,1]==a)&(EMPD[,2]==i)])]=0
PA[(row.names(PA)==i)&(PA>0)&(PA<sum(EMPD[,3][(EMPD[,1]==a)&(EMPD[,2]==i)]))]=sqrt(sum(EMPD[,3][(EMPD[,1]==a)&(EMPD[,2]==i)])-PA[(row.names(PA)==i)&(PA>0)&(PA<sum(EMPD[,3][(EMPD[,1]==a)&(EMPD[,2]==i)]))])}
BPFT=BPFT1
  
# PFT scoring
ScorePFT=data.frame(matrix(NA,nrow=1,ncol=ncol(PA)))
colnames(ScorePFT)=colnames(BPFT)
BPFT=BPFT1
for(i in colnames(PA))
{ScorePFT[colnames(ScorePFT)==i]=sum(PA[i])}

# Aggregation into biome scores
for (j in rownames(BPFT))
{B=BPFT[j,]
for(i in 1:ncol(B))
{if ((BPFT[j,i]!='')&(BPFT[j,i]!="+")&(BPFT[j,i]!="-"))
{if ((BPFT[j,i]==names(ScorePFT[i])))
{BPFT[j,i]=as.numeric(ScorePFT[,i])}}}}

for(j in rownames(BPFT))
{B=BPFT[j,]
for (i in 1:ncol(B))
{if ((BPFT[j,i]!='')&(BPFT[j,i]!="+")&(BPFT[j,i]!="-"))
{if ((BPFT[j,i]!=names(ScorePFT[i]))&(length(BPFT[j,BPFT[j,i]])!=0))
{BPFT[j,BPFT[j,i]]=as.numeric(BPFT[j,BPFT[j,i]])+as.numeric(ScorePFT[,colnames(B[i])])
BPFT[j,colnames(B[i])]=''}}}}
  
BPFT2=BPFT
BPFT2[BPFT2=='+']=''
BPFT2[BPFT2=='-']=''
BPFT3=BPFT2
for(i in colnames(PA))
{BPFT2[BPFT2==i]=sum(ScorePFT[i])}
  
BPFT[BPFT=='']=NA
BPFT3[BPFT3=='']=NA
  
SumBiome=data.frame(matrix(NA,ncol=1,nrow=nrow(BPFT)))
for(i in 1:(nrow(SumBiome)))
{SumBiome[i,]<-sum(as.numeric(BPFT2[i,]),na.rm=TRUE)}
  
for(i in colnames(PA))
{BPFT2[BPFT2==i]=sum(ScorePFT[i])}
BPFT[BPFT=='']=NA
BPFT3[BPFT3=='']=NA

BPFT4=BPFT
for(i in 1:ncol(BPFT))
{for(j in 1:nrow(BPFT))
{y=j-1
z=j+1
if (is.na(BPFT4[j,i])==FALSE)
{if((BPFT4[j,i]=="+")&(sum(as.numeric(BPFT2[j,]),na.rm=TRUE)==0))
{BPFT[j,i]=as.numeric(ScorePFT[names(BPFT4[i])])}
if((BPFT4[j,i]=="+")&(sum(as.numeric(BPFT2[j,]),na.rm=TRUE)!=0))
{BPFT[z,i]=as.numeric(ScorePFT[names(BPFT4[i])])
BPFT[j,i]=0}
if((BPFT4[j,i]=="-")&(sum(as.numeric(BPFT2[j,]),na.rm=TRUE)==0))
{BPFT[j,i]=as.numeric(ScorePFT[names(BPFT4[i])])}
if((BPFT4[j,i]=="-")&(sum(as.numeric(BPFT2[j,]),na.rm=TRUE)!=0))
{BPFT[y,i]=as.numeric(ScorePFT[names(BPFT4[i])])
BPFT[j,i]=0}}
}}
BPFT[(BPFT=='+')|(BPFT=='-')]=NA

# Compute biome scores
BiomeIndexResult=data.frame(matrix(NA,ncol=1,nrow=length(Name_Biome)))
row.names(BiomeIndexResult)=Name_Biome
for(i in 1:(nrow(BiomeIndexResult)))
{BiomeIndexResult[i,]<-sum(as.numeric(BPFT[i,]),na.rm=TRUE)}
NameScore=Name_Biome[which(BiomeIndexResult==BiomeIndexResult[BiomeIndexResult==max(BiomeIndexResult,na.rm=TRUE)])]
Max=max(BiomeIndexResult,na.rm=TRUE)
Sd=sd(BiomeIndexResult$matrix.NA..ncol...1..nrow...length.Name_Biome..,na.rm=TRUE)
if (length(max(BiomeIndexResult,na.rm=TRUE)!=1))
{NameScore=NameScore[1]
Max=Max[1]
Sd=Sd[1]}

# Best biome selection
Output[n,]=c(NameScore,Max,Sd)
n=n+1
Sys.sleep(0.1)
setTxtProgressBar(pb, n)
if(n==length(NameCal)/2){print("Eso eso que")}}
nameFile=str_c("./",namefile,".csv") 
return(write.csv2(Output,nameFile))}

EuropeBiomeIndexFossil<-function(Name,AssPol,Ages,namefile){
# Initialise output table
Output=data.frame(matrix(NA,ncol=nrow(BPFT),nrow=nrow(AssPol)))
colnames(Output)=rownames(BPFT)
rownames(Output)=make.names(Ages, unique=TRUE)

n=1
# Progress bar
pb <- txtProgressBar(min = 0, 
                     max = nrow(AssPol),
                     style = 3, 
                     width = 50,
                     char = "=")   
for (a in 1:nrow(AssPol))
{if(n==1){print("Fasten your seatbelts!")}
PA=PA1
for(i in Name$Names)
# Adjust PFT contributions based on observed abundances
{PA=PA[which(row.names(PA) %in% Name$Names),] 
PA[(row.names(PA)==i)&(PA>0)&(PA>AssPol[a,i])]=0
PA[(row.names(PA)==i)&(PA>0)&(PA<AssPol[a,i])]=sqrt(AssPol[a,i]-PA[(row.names(PA)==i)&(PA>0)&(PA<AssPol[a,i])])}
ScorePFT=data.frame(matrix(NA,nrow=1,ncol=ncol(PA)))
colnames(ScorePFT)=colnames(BPFT)
BPFT=BPFT1
# PFT scoring
for(i in colnames(PA))
{ScorePFT[colnames(ScorePFT)==i]=sum(PA[i])}

for (j in rownames(BPFT))
{B=BPFT[j,]
for(i in 1:ncol(B))
{if ((BPFT[j,i]!='')&(BPFT[j,i]!="+")&(BPFT[j,i]!="-"))
{if ((BPFT[j,i]==names(ScorePFT[i])))
{BPFT[j,i]=as.numeric(ScorePFT[,i])}}}}

for(j in rownames(BPFT))
{B=BPFT[j,]
for (i in 1:ncol(B))
{if ((BPFT[j,i]!='')&(BPFT[j,i]!="+")&(BPFT[j,i]!="-"))
{if ((BPFT[j,i]!=names(ScorePFT[i]))&(length(BPFT[j,BPFT[j,i]])!=0))
{BPFT[j,BPFT[j,i]]=as.numeric(BPFT[j,BPFT[j,i]])+as.numeric(ScorePFT[,colnames(B[i])])
BPFT[j,colnames(B[i])]=''}}}}

BPFT2=BPFT
BPFT2[BPFT2=='+']=''
BPFT2[BPFT2=='-']=''
BPFT3=BPFT2
BPFT2[BPFT2=='']=NA
SumBiome=data.frame(matrix(NA,ncol=1,nrow=nrow(BPFT)))
for(i in 1:(nrow(SumBiome)))
{SumBiome[i,]<-sum(as.numeric(BPFT2[i,]),na.rm=TRUE)}

for(i in colnames(PA))
{BPFT2[BPFT2==i]=sum(ScorePFT[i])}
BPFT[BPFT=='']=NA
BPFT3[BPFT3=='']=NA

BPFT4=BPFT
for(i in 1:ncol(BPFT))
{for(j in 1:nrow(BPFT))
{y=j-1
z=j+1
if (is.na(BPFT4[j,i])==FALSE)
{if((BPFT4[j,i]=="+")&(sum(as.numeric(BPFT2[j,]),na.rm=TRUE)==0))
{BPFT[j,i]=as.numeric(ScorePFT[names(BPFT4[i])])}
  if((BPFT4[j,i]=="+")&(sum(as.numeric(BPFT2[j,]),na.rm=TRUE)!=0))
  {BPFT[z,i]=as.numeric(ScorePFT[names(BPFT4[i])])
  BPFT[j,i]=0}
  if((BPFT4[j,i]=="-")&(sum(as.numeric(BPFT2[j,]),na.rm=TRUE)==0))
  {BPFT[j,i]=as.numeric(ScorePFT[names(BPFT4[i])])}
  if((BPFT4[j,i]=="-")&(sum(as.numeric(BPFT2[j,]),na.rm=TRUE)!=0))
  {BPFT[y,i]=as.numeric(ScorePFT[names(BPFT4[i])])
  BPFT[j,i]=0}}
}}
BPFT[(BPFT=='+')|(BPFT=='-')]=NA
# Biome score
BiomeIndexResult=data.frame(matrix(NA,ncol=1,nrow=length(Name_Biome)))
row.names(BiomeIndexResult)=Name_Biome
for(i in 1:(nrow(BiomeIndexResult)))
{BiomeIndexResult[i,]<-sum(as.numeric(BPFT[i,]),na.rm=TRUE)
Output[n,i]=BiomeIndexResult[i,]}

# Compute biome scores
n=n+1
Sys.sleep(0.1)
setTxtProgressBar(pb, n)
if(n==nrow(AssPol)/2){print("Eso eso que")}
close(pb) }
for(i in 1:nrow(Output)){Output[i,"Best"]=colnames(Output)[which.max(Output[i,])]
}
nameFile=str_c("./",namefile,".csv") 
return(write.csv2(Output,nameFile))
}
# =========================================================
# END
# =========================================================


