# PLOT, MAP AND ANALYSE PRESENT-DAY RICHNESS AND NPP DATA (Figures S6 and S7)

# Contents
# line 11: Functions and packages, data preparation
# line 140: Generalized linear models for Fig. S7
# line 165: Check and account for spatial autocorrelation
# line 337: Plots and maps

##############################################################################

#Functions and packages

#function to plot maps of climate and richness
plotmap <- function(y, x, add.map = TRUE, head = FALSE, xlimits = c(-180,180),
                    ylimits = c(-55,90), zmin=0, zmax=max(dat[,x], na.rm=T),
                    mapint=F, mapcol="black", head_line=1, legendmar=4,
                    colscale = "NORM") 
  #y is a gridded dataset, x is the variable to be plotted 
{
  dat <- as.data.frame(y)
  if(colscale == "NORM")
  {
    if(zmax >= 20) 
    {plotcol <- rev(rainbow(zmax, end=0.7))} else 
    {plotcol <- rev(rainbow(100, end=0.7))}
  } else
  {
    if(zmax >= 20) 
    {plotcol <- rainbow(zmax, end=0.7)} else 
    {plotcol <- rainbow(100, end=0.7)}
  }    
  #plotcol <- rev(rainbow(max(abs(dat[,x]), na.rm=T), end=0.7))
  image.plot(y[,x], zlim=c(zmin, zmax), legend.width=0.5, 
             legend.shrink=0.6, 
             legend.mar=legendmar, col = plotcol, axes=F, 
             horizontal=T,
             xlim=xlimits, ylim=ylimits)
  if(head == FALSE) title(names(dat)[which(names(dat) == x)], 
                          line=head_line) else title(head, line=head_line)
  if(add.map == "TRUE") map(add=T, int=mapint, xlim=xlimits, ylim=ylimits,
                            col=mapcol)
}

#function to calculate NPP from mat (= mean annual temperature) and 
# map (annual precipitation) estimates

npp <- function(mat, map) {
  nppT <- 3000/(1+exp(1.315-0.119*mat))
  nppP <- 3000*(1-exp(-0.000664*map))
  NPP <- nppT
  NPP[which(nppP < nppT)] <- nppP[which(nppP < nppT)]
  NPP[which(is.na(nppT) | is.na(nppP))] <- NA
  return(NPP)
}

#EXPLANATION HABERL DATA (http://www.uni-klu.ac.at/socec/inhalt/1191.htm)
#npp_act_haberl (=tnap_all_gcm): NPP of the actual vegetation in [gC/mē/yr]
#hanpp_haberl (=thanpppallgcm): human appropriation of NPP in [gC/mē/yr]
#hanpp_perc_haberl (=thanpppall_pc): total HANPP as a percentage of NPP_0
#delta_npp_perc_haberl (=tdnpplcpallpc): land-use induced reductions in NPP
# as percentage of NPP_0
#npp_t_haberl (=tntp_all_gcm): NPP remaining in ecosystems after harvest 
# in [gC/mē/yr]
#npp_0_haberl (=tn0_all_gcm): NPP of the potential vegetation in [gC/mē/yr]


#function to calculate species richness per category
length.unique  <- function(x) length(unique(x))

#load packages

library(fields)
library(sp)
library(foreign)
library(lme4)
library(mgcv)
library(rgdal)
library(spdep)
library(raster)


#Load and check data

#set working directory
setwd("D:/Gammadiv ms/FilesForUpload")

#global gridded dataset with species richness, geographical and environmental 
# variables
clim_rich_global <- read.delim("Data_PresentDay_Global.txt", as.is = T)

#gridded dataset for the study regions
clim_rich_globOURreg <- read.delim("Data_PresentDay_StudyRegions.txt", 
                                   as.is = T)

#unique(clim_rich_globreg$Region)

#names(clim_rich_global)

#check spatial setup of different datasets
plot(clim_rich_global$XCoord, clim_rich_global$YCoord)
plot(clim_rich_globOURreg$XCoord, clim_rich_globOURreg$YCoord)


#prepare global grids as grids

#remove NA values
clim_rich_global_grid  <- clim_rich_global[-which(is.na(
  clim_rich_global$rich_sp_glob)),]
clim_rich_globOURreg_grid  <- clim_rich_globOURreg[-which(is.na(
  clim_rich_globOURreg$rich_sp_glob)),]

#transform gridded datasets to spatial objects
coordinates(clim_rich_global_grid)  <- ~XCoord + YCoord
gridded(clim_rich_global_grid)  <- TRUE
coordinates(clim_rich_globOURreg_grid)  <- ~XCoord + YCoord
gridded(clim_rich_globOURreg_grid)  <- TRUE


#define the limits of the NPP gradient estimated from the 
#Neogene paleobotanical data

#minimum value
NPPmin_HD <- 582.3 
#maximum value
NPPmax_HD <- 2021

limitsHD_rem <- which(clim_rich_global$nppMIA_mean_remain > NPPmin_HD 
                      & clim_rich_global$nppMIA_mean_remain < NPPmax_HD)

limitsHD_ourreg_rem <- which(
  clim_rich_globOURreg$nppMIA_mean_remain > NPPmin_HD 
  & clim_rich_globOURreg$nppMIA_mean_remain < NPPmax_HD)

#grids with genus richness > 0
limits_mam <- which(clim_rich_global$rich_gen_glob > 0)
limitsHD_rem_mam <- intersect(limitsHD_rem, limits_mam)

##############################################################################

# GLMs for Fig. S7

#global grid (full gradient)
glm_glob <- glm(rich_gen_glob ~ nppMIA_mean_remain, family=poisson, 
                data=clim_rich_global)
summary(glm_glob)

#study regions (full gradient)
glm_ourreg <- glm(rich_gen_glob ~ nppMIA_mean_remain, family=poisson, 
                  data=clim_rich_globOURreg)
summary(glm_ourreg)

#global grid (fossil gradient)
glm_glob_fossilgrad <- glm(rich_gen_glob ~ nppMIA_mean_remain, family=poisson, 
                           data=clim_rich_global[limitsHD_rem,])
summary(glm_glob_fossilgrad)

#study regions (fossil gradient)
glm_ourreg_fossilgrad <- glm(rich_gen_glob ~ nppMIA_mean_remain, 
                             family=poisson, 
                             data=clim_rich_globOURreg[limitsHD_ourreg_rem,])
summary(glm_ourreg_fossilgrad)

##############################################################################

# CHECKING AND ACCOUNTING FOR SPATIAL AUTOCORRELATION

#global model, full gradient

#create vector with row IDs where NPP or richness values are not NA
length(which(is.na(clim_rich_global$nppMIA_mean_remain)))
length(which(is.na(clim_rich_global$rich_gen_glob)))

notNA_glob <- intersect(which(!is.na(clim_rich_global$nppMIA_mean_remain)),
                 which(!is.na(clim_rich_global$rich_gen_glob)))

#create new dataframe without NA values in NPP or richness
clim_rich_glob_withoutNA <- clim_rich_global[notNA_glob,]
#names(clim_rich_glob_withoutNA)

#global glm of richness in relation to NPP
glm_glob <- glm(rich_gen_glob ~ nppMIA_mean_remain, family=poisson, 
                data=clim_rich_glob_withoutNA)

#create vector with residuals of global glm
res_glm_glob <- glm_glob$residuals

#create matrix of coordinates (for spatial autocorrelation analyses)
xy_glob <- cbind(clim_rich_glob_withoutNA$XCoord, 
                 clim_rich_glob_withoutNA$YCoord)

#create neighbour file (class: nb)
nb_glob <- tri2nb(xy_glob)

#check neighbour values
summary(unlist(nbdists(nb_glob, xy_glob)))

#create list with spatial weights
listw_glob <- nb2listw(nb_glob)
length(listw_glob)

#test for residual spatial autocorrelation in using Moran's I
moran.test(x = res_glm_glob, listw = listw_glob)

#calculate spatial covariate
clim_rich_glob_withoutNA$spatcov_glob <- autocov_dist(
  res_glm_glob, xy_glob, nbs=max(unlist(nbdists(nb_glob, xy_glob))),
  longlat=TRUE, zero.policy=TRUE, style="W")


#global glm of richness in dependency of NPP, including spatial covariate to
# account for spatial autocorrelation
glm_glob_SAC <- glm(rich_gen_glob ~ spatcov_glob + nppMIA_mean_remain, 
                family=poisson, data=clim_rich_glob_withoutNA)
#vector of residuals
res_glm_glob_SAC <- glm_glob_SAC$residuals

#test for residual spatial autocorrelation using Moran's I
moran.test(x = res_glm_glob_SAC, listw = listw_glob)



# global model, only for fossil gradient
# (exact same setup as for full gradient)

notNA_glob_fossilgrad <- intersect(notNA_glob, limitsHD_rem)

#create new dataframe without NA values in NPP or richness
clim_rich_glob_fossilgrad_withoutNA <- 
  clim_rich_global[notNA_glob_fossilgrad,]

#global glm of richness in dependency of NPP
glm_glob_fossilgrad <- glm(rich_gen_glob ~ nppMIA_mean_remain, family=poisson, 
                           data=clim_rich_glob_fossilgrad_withoutNA)

#create vector with residuals
res_glm_glob_fossilgrad <- glm_glob_fossilgrad$residuals

#create matrix of coordinates (for spatial autocorr analyses)
xy_glob_fossilgrad <- cbind(clim_rich_glob_fossilgrad_withoutNA$XCoord, 
                 clim_rich_glob_fossilgrad_withoutNA$YCoord)

#create neighbour file (class: nb)
nb_glob_fossilgrad<- tri2nb(xy_glob_fossilgrad)

#check neighbour values
summary(unlist(nbdists(nb_glob_fossilgrad, xy_glob_fossilgrad)))

#create list with spatial weights
listw_glob_fossilgrad <- nb2listw(nb_glob_fossilgrad)

#test for spatial autocorrelation using Moran's I
moran.test(x = res_glm_glob_fossilgrad, listw = listw_glob_fossilgrad, 
           zero.policy = TRUE)

#calculate spatial covariate
clim_rich_glob_fossilgrad_withoutNA$spatcov_glob_fossilgrad <- autocov_dist(
  res_glm_glob_fossilgrad, xy_glob_fossilgrad, 
  nbs=max(unlist(nbdists(nb_glob_fossilgrad, xy_glob_fossilgrad))),
  longlat=TRUE, zero.policy=TRUE, style="W")

#global glm of richness in dependency of NPP, including spatial covariate to
# account for spatial autocorrelation
glm_glob_fossilgrad_SAC <- glm(rich_gen_glob ~ spatcov_glob_fossilgrad + 
                                 nppMIA_mean_remain, family=poisson, 
                               data=clim_rich_glob_fossilgrad_withoutNA)
summary(glm_glob_fossilgrad_SAC)

#vector of residuals
res_glm_glob_fossilgrad_SAC <- glm_glob_fossilgrad_SAC$residuals

#test for spatial autocorrelation using Moran's I
moran.test(x = res_glm_glob_fossilgrad_SAC, listw = listw_glob_fossilgrad)



# model just for study regions

#create vector with row IDs where NPP or richness values are not NA
notNA_globOURreg <- intersect(
  which(!is.na(clim_rich_globOURreg$nppMIA_mean_remain)),
  which(!is.na(clim_rich_globOURreg$rich_gen_glob)))

#create new dataframe without NA values in NPP or richness
clim_rich_globOURreg_withoutNA <- clim_rich_globOURreg[notNA_globOURreg,]
names(clim_rich_globOURreg_withoutNA)

#global glm of richness in dependency of NPP
glm_ourreg <- glm(rich_gen_glob ~ nppMIA_mean_remain, family=poisson, 
                      data=clim_rich_globOURreg_withoutNA)
summary(glm_ourreg)

#create vector with residuals
res_glm_ourreg <- glm_ourreg$residuals

#create matrix of coordinates (for spatial autocorr analyses)
xy_globOURreg <- cbind(clim_rich_globOURreg_withoutNA$XCoord, 
                       clim_rich_globOURreg_withoutNA$YCoord)

#create neighbour file (class nb)
nb_globOURreg<- tri2nb(xy_globOURreg)

#plot(nb_glob, coords=xy_globOURreg)

#check neighbour values
summary(unlist(nbdists(nb_globOURreg, xy_globOURreg)))

#create list with spatial weights
listw_globOURreg <- nb2listw(nb_globOURreg)

#test for spatial autocorrelation using Moran's I
moran.test(x = res_glm_ourreg, listw = listw_globOURreg)

#calculate spatial covariate
clim_rich_globOURreg_withoutNA$spatcov_globOURreg <- autocov_dist(
  res_glm_ourreg, xy_globOURreg, 
  nbs=max(unlist(nbdists(nb_globOURreg, xy_globOURreg))),
  longlat=TRUE, zero.policy=TRUE, style="W")

#global glm of richness in dependency of NPP, including spatial covariate to
# account for spatial autocorrelation
glm_ourreg_SAC <- glm(rich_gen_glob ~ spatcov_globOURreg + 
                            nppMIA_mean_remain, family=poisson, 
                          data=clim_rich_globOURreg_withoutNA)
summary(glm_ourreg_SAC)

#vector of residuals
res_glm_ourreg_SAC <- glm_ourreg_SAC$residuals

#test for spatial autocorrelation using Moran's I
moran.test(x = res_glm_ourreg_SAC, listw = listw_globOURreg)


##############################################################################

# PLOTS AND MAPS

#create folder for figures and maps
dir.create("figures_maps/")

# create new x data for model predictions of GLMs

#global dataset, full gradient
newdata_x <- seq(0, max(clim_rich_global$nppMIA_mean_remain, na.rm=T), by=1)

#global dataset, fossil gradient
newdata_x_fossilgrad <- seq(
  min(clim_rich_global$nppMIA_mean_remain[limitsHD_rem], na.rm=T),
  max(clim_rich_global$nppMIA_mean_remain[limitsHD_rem], na.rm=T), by=1)

#dataset for study regions
newdata_x_ourreg <- seq(min(clim_rich_globOURreg$nppMIA_mean_remain, na.rm=T),
                        max(clim_rich_globOURreg$nppMIA_mean_remain, na.rm=T), 
                        by=1)

newdata_x_ourreg_fossilgrad <- seq(
  min(clim_rich_globOURreg$nppMIA_mean_remain[limitsHD_ourreg_rem], na.rm=T),
  max(clim_rich_globOURreg$nppMIA_mean_remain[limitsHD_ourreg_rem], na.rm=T),
  by=1)

#define limits of x any y axes
x.lim <- c(0,max(clim_rich_global$nppMIA_mean_remain, na.rm=T))
y.lim <- c(0,max(clim_rich_global$rich_gen_glob, na.rm=T))

#define colours
MyGrey1 <- rgb(t(col2rgb("grey90")), alpha=100, maxColorValue=255)
MyGrey2 <- rgb(t(col2rgb("grey60")), alpha=80, maxColorValue=255)
MyBlue1 <- rgb(t(col2rgb("dodgerblue")), alpha=80, maxColorValue=255)
MyBlue2 <- rgb(t(col2rgb("dodgerblue4")), alpha=80, maxColorValue=255)


# Fig. S6 A

par(mfcol=c(1,1))

pdf("figures_maps/FigS6NEW_NPP_maps_A.pdf", width=7, height=5)

par(mfcol=c(1,1), mar=c(  3,1,2,1), oma=c(2,1,1,1))

#A: Miami model, REMAINING
clim_rich_global_grid <- clim_rich_global
coordinates(clim_rich_global_grid)  <- ~XCoord + YCoord
gridded(clim_rich_global_grid)  <- TRUE

plotmap(clim_rich_global_grid, "nppMIA_mean", head="", 
        ylimits=c(-55,90), legendmar=8)
mtext("A", 3, cex=1.5, adj=0, line=-2, font=2)
mtext("Potential NPP [Miami, g dry matter/mē/year]", 1, line=5.5, cex=1.2)

dev.off()


# Fig. S6 B

pdf("figures_maps/FigS6NEW_NPP_maps_B.pdf", width=7, height=5)
par(mar=c(3,1,2,1), oma=c(2,1,1,1))
#C: Potential NPP, Haberl
clim_rich_global_grid  <- clim_rich_global[-which(is.na(
  clim_rich_global$npp_0_haberl)),]
coordinates(clim_rich_global_grid)  <- ~XCoord + YCoord
gridded(clim_rich_global_grid)  <- TRUE

plotmap(clim_rich_global_grid, "npp_0_haberl", 
        head="", legendmar=8)
mtext("B", 3, cex=1.5, adj=0, line=-2, font=2)
mtext("Potential NPP [DGVM, g C/mē/year]", 1, line=5.5, cex=1.2)

dev.off()


# Fig. S6 C

pdf("figures_maps/FigS6NEW_NPP_maps_C.pdf", width=7, height=5)
par(mar=c(3,1,2,1), oma=c(2,1,1,1))
#B: HANPP (correction factor)
clim_rich_global_grid  <- clim_rich_global[-which(is.na(
  clim_rich_global$MYhanpp_factor)),]
coordinates(clim_rich_global_grid)  <- ~XCoord + YCoord
gridded(clim_rich_global_grid)  <- TRUE

plotmap(clim_rich_global_grid, "MYhanpp_factor", 
        head="", colscale="R", legendmar=8)
mtext("C", 3, cex=1.5, adj=0, line=-2, font=2)
mtext("HANPP correction factor [proportion]", 1, line=5.5, cex=1.2)

dev.off()


# Fig. S6 D

pdf("figures_maps/FigS6NEW_NPP_maps_D.pdf", width=7, height=5)
par(mar=c(3,1,2,1), oma=c(2,1,1,1))
#D: Remaining (A * B)
clim_rich_global_grid  <- clim_rich_global[-which(is.na(
  clim_rich_global$nppMIA_mean_remain)),]
coordinates(clim_rich_global_grid)  <- ~XCoord + YCoord
gridded(clim_rich_global_grid)  <- TRUE

plotmap(clim_rich_global_grid, "nppMIA_mean_remain", 
        head="", legendmar=8)
mtext("D", 3, cex=1.5, adj=0, line=-2, font=2)
mtext("Remaining NPP [g dry matter/mē/year]", 1, line=5.5, cex=1.2)

dev.off()


# Fig. S6 E

#scatter plot: NPP Miami ~ NPP 0 Haberl

#linear model: NPP Miami ~ NPP 0 Haberl
lm_NPP <- lm(nppMIA_mean ~ npp_0_haberl, data=clim_rich_global_grid)
#summary(lm_NPP)

#new x variables for plotting the regression line
x_Haberl <- seq(min(clim_rich_global_grid$npp_0_haberl, na.rm=T),
                max(clim_rich_global_grid$npp_0_haberl, na.rm=T), by=1)

pdf("figures_maps/FigS6NEW_NPP_HaberlMiami_globally_E.pdf", width=6, height=6)

par(mfrow=c(1,1), mar=c(4,4,1,1), oma=c(1,1,1,1))
plot(clim_rich_global_grid$nppMIA_mean ~ 
       clim_rich_global_grid$npp_0_haberl, 
     xlab="Potential NPP [DGVM, g C/mē/year]", 
     ylab="Potential NPP [Miami, g dry matter/mē/year]",
     pch=16, col=MyGrey2, ylim=c(-100,2800))

lines(x_Haberl, predict(lm_NPP, data.frame(npp_0_haberl = x_Haberl),
                        type="response"), lwd=3, col="grey35")

dev.off()


# Fig. S7 A (maps of present-day genus richness)

#global
pdf("figures_maps/GenRichness_globally.pdf", width=8, 
    height=4.5)
par(mfrow=c(1,1), mar=c(1,1,1,1), oma=c(0,0,3,0))
plotmap(clim_rich_global_grid, "rich_gen_glob", head="", ylimits=c(-55,90),
        zmax=63)
dev.off()

#study regions
pdf("figures_maps/GenRichness_ourregions_globscale.pdf", width=8, 
    height=4.5)
par(mfrow=c(1,1), mar=c(1,1,1,1), oma=c(0,0,3,0))
plotmap(clim_rich_globOURreg_grid, "rich_gen_glob", head="", ylimits=c(15,75),
        xlimits=c(-125,53), zmax=63, mapint=FALSE)
dev.off()


# Fig. S7 B

pdf("figures_maps/plot_glob_ourreg_greyblue_TRANSPARENT.pdf", 
    width=6, height=6)

pointsize <- 0.8

par(mfcol=c(1,1), oma=c(1,1,1,1), mar=c(4,4,2,2), las=1, xpd=FALSE)

plot(clim_rich_global$rich_gen_glob ~ clim_rich_global$nppMIA_mean_remain,
     xlab="Net primary production (g dry matter/mē/year)", 
     ylab="Present-day genus richness", 
     xlim=x.lim, ylim=y.lim,
     main="", pch=16, cex=pointsize, col=MyGrey1)

abline(v=NPPmax_HD, lty=2, lwd=2, col="grey20", ylim=y.lim)
abline(v=NPPmin_HD, lty=2, lwd=2, col="grey20", ylim=y.lim)

par(xpd=TRUE)

points(clim_rich_global$rich_gen_glob[limitsHD_rem] ~ 
         clim_rich_global$nppMIA_mean_remain[limitsHD_rem], pch=16,
       col=MyGrey2, cex=pointsize)

points(clim_rich_globOURreg$rich_gen_glob ~ 
         clim_rich_globOURreg$nppMIA_mean_remain, pch=16,
       col=MyBlue1, cex=0.6)

points(clim_rich_globOURreg$rich_gen_glob[limitsHD_ourreg_rem] ~ 
         clim_rich_globOURreg$nppMIA_mean_remain[limitsHD_ourreg_rem], pch=16,
       col=MyBlue2, cex=0.6)

#add glm lines

lines(newdata_x, 
      predict(glm_glob, 
              data.frame(nppMIA_mean_remain = newdata_x), 
              type="response"), lwd=3, col="grey40")

lines(newdata_x_fossilgrad, 
      predict(glm_glob_fossilgrad, 
              data.frame(nppMIA_mean_remain = newdata_x_fossilgrad), 
              type="response"), lwd=3, col="black")


lines(newdata_x_ourreg, 
      predict(glm_ourreg, 
              data.frame(nppMIA_mean_remain = newdata_x_ourreg), 
              type="response"), lwd=3, col="red")

dev.off()


# Fig. S7 C (map illustrating location of coloured points in Fig. S7B)

#jpeg("figures_maps/map_glob_fossilgrad_ourreg_greyblue.jpeg", res=300, 
#     width=1500, height=1000)
pdf("figures_maps/map_glob_fossilgrad_ourreg_greyblue.pdf", 
    width=10, height=7.5)

map(int=F, ylim=c(-55, 90))

#points(clim_rich_global$XCoord, clim_rich_global$YCoord, pch=16, cex=0.3, 
#     col="grey90")

points(clim_rich_global$XCoord[limits_mam], 
       clim_rich_global$YCoord[limits_mam], 
       pch=16, cex=0.4, 
       col="grey90")

points(clim_rich_global$XCoord[limitsHD_rem_mam], 
       clim_rich_global$YCoord[limitsHD_rem_mam], pch=16, cex=0.4, 
       col="grey60")

points(clim_rich_globOURreg$XCoord, 
       clim_rich_globOURreg$YCoord, pch=16, cex=0.4, 
       col="lightblue")

points(clim_rich_globOURreg$XCoord[limitsHD_ourreg_rem], 
       clim_rich_globOURreg$YCoord[limitsHD_ourreg_rem], pch=16, cex=0.4, 
       col="dodgerblue4")

map(add=T, int=F)

dev.off()
