#########
#
#  1. Read in the CMIP6 climate projection layers, and do exploratory plots
#  2. Plot Fig. 4
#
library(ncdf4) # package for netcdf manipulation
library(raster) # package for raster manipulation
library(rgdal) # package for geospatial analysis
library("tidyverse")

#Read in the netCDF file contents

tc_sum <- nc_open('CMIP/CMIP6 - Sea Surface Temperature (SST) Change deg C - Medium Term (2041-2060) SSP3 7 1981-2010 - June to August (mean of 24 models).nc')

lon <- ncvar_get(tc_sum, "lon")
lat <- ncvar_get(tc_sum, "lat", verbose = F)

val1 <- as.data.frame(ncvar_get(tc_sum))



dsum <- as.data.frame(matrix(NA,nrow = 360*180, ncol = 3))
colnames(dsum) <- c("lon","lat","val")
for(i in 1:180) {
  dsum$lon[(i-1)*360+(1:360)] <- lon
  dsum$lat[(i-1)*360+(1:360)] <- lat[i]
  dsum$val[(i-1)*360+(1:360)] <- val1[,i]
}

# ggplot(data = dsum, aes(x = lon, y = lat, colour = val)) + geom_point()


tc_win <- nc_open('CMIP/CMIP6 - Sea Surface Temperature (SST) Change deg C - Medium Term (2041-2060) SSP3 7 1981-2010 - December to February (mean of 24 models).nc')
lon <- ncvar_get(tc_win, "lon")
lat <- ncvar_get(tc_win, "lat", verbose = F)

val2 <- as.data.frame(ncvar_get(tc_win))



dwin <- as.data.frame(matrix(NA,nrow = 360*180, ncol = 3))
colnames(dwin) <- c("lon","lat","val")
for(i in 1:180) {
  dwin$lon[(i-1)*360+(1:360)] <- lon
  dwin$lat[(i-1)*360+(1:360)] <- lat[i]
  dwin$val[(i-1)*360+(1:360)] <- val2[,i]
}
# windows(8,5)
# ggplot(data = dwin, aes(x = lon, y = lat, colour = val)) + geom_point()

ddiff <- dwin
ddiff$val <- dsum$val - dwin$val
# windows(8,5)
# ggplot(data = ddiff, aes(x = lon, y = lat, colour = val)) + geom_point()

tc_an <- nc_open('CMIP/CMIP6 - Sea Surface Temperature (SST) Change deg C - Medium Term (2041-2060) SSP3 7 1981-2010 - Annual (mean of 24 models).nc')
lon <- ncvar_get(tc_an, "lon")
lat <- ncvar_get(tc_an, "lat", verbose = F)

val2 <- as.data.frame(ncvar_get(tc_an))



dan <- as.data.frame(matrix(NA,nrow = 360*180, ncol = 3))
colnames(dan) <- c("lon","lat","val")
for(i in 1:180) {
  dan$lon[(i-1)*360+(1:360)] <- lon
  dan$lat[(i-1)*360+(1:360)] <- lat[i]
  dan$val[(i-1)*360+(1:360)] <- val2[,i]
}
# windows(8,5)
# ggplot(data = dan, aes(x = lon, y = lat, colour = val)) + geom_point()

tc_win_abs <- nc_open("CMIP/CMIP6 - Sea Surface Temperature (SST) deg C - Medium Term (2041-2060) SSP3 7 1981-2010 - December to February (mean of 24 models).nc")
lon <- ncvar_get(tc_win_abs, "lon")
lat <- ncvar_get(tc_win_abs, "lat", verbose = F)

val2 <- as.data.frame(ncvar_get(tc_win_abs))



win_abs <- as.data.frame(matrix(NA,nrow = 360*180, ncol = 3))
colnames(win_abs) <- c("lon","lat","val")
for(i in 1:180) {
  win_abs$lon[(i-1)*360+(1:360)] <- lon
  win_abs$lat[(i-1)*360+(1:360)] <- lat[i]
  win_abs$val[(i-1)*360+(1:360)] <- val2[,i]
}
# windows(8,5)
# ggplot(data = win_abs, aes(x = lon, y = lat, colour = val)) + geom_point()


tc_sum_abs <- nc_open("CMIP/CMIP6 - Sea Surface Temperature (SST) deg C - Medium Term (2041-2060) SSP3 7 1981-2010 - June to August (mean of 24 models).nc")
lon <- ncvar_get(tc_sum_abs, "lon")
lat <- ncvar_get(tc_sum_abs, "lat", verbose = F)

val2 <- as.data.frame(ncvar_get(tc_sum_abs))



sum_abs <- as.data.frame(matrix(NA,nrow = 360*180, ncol = 3))
colnames(sum_abs) <- c("lon","lat","val")
for(i in 1:180) {
  sum_abs$lon[(i-1)*360+(1:360)] <- lon
  sum_abs$lat[(i-1)*360+(1:360)] <- lat[i]
  sum_abs$val[(i-1)*360+(1:360)] <- val2[,i]
}
# windows(8,5)
# ggplot(data = sum_abs, aes(x = lon, y = lat, colour = val)) + geom_point()



d_abs_diff <- sum_abs
d_abs_diff$val <- abs(sum_abs$val - win_abs$val)
# windows(8,5)
# ggplot(data = d_abs_diff, aes(x = lon, y = lat, colour = val)) + geom_point() + scale_colour_gradientn(colours = c("blue","green","red"))


min_abs_new <- sum_abs
min_abs_new$val <- sapply(1:nrow(sum_abs), function(x) min(c(sum_abs$val[x],win_abs$val[x])))
max_abs_new <- sum_abs
max_abs_new$val <- sapply(1:nrow(sum_abs), function(x) max(c(sum_abs$val[x],win_abs$val[x])))
# plot(min_abs_new$val,max_abs_new$val)
# abline(a=0,b=1)

win_abs_old <- win_abs
win_abs_old$val <- win_abs$val - dwin$val
# plot(win_abs_old$val,win_abs$val)

sum_abs_old <- sum_abs
sum_abs_old$val <- sum_abs$val - dsum$val
# plot(sum_abs_old$val,sum_abs$val)

min_abs_old <- sum_abs
min_abs_old$val <- sapply(1:nrow(sum_abs), function(x) min(c(sum_abs_old$val[x],win_abs_old$val[x])))
max_abs_old <- sum_abs
max_abs_old$val <- sapply(1:nrow(sum_abs), function(x) max(c(sum_abs_old$val[x],win_abs_old$val[x])))
# plot(min_abs_old$val,max_abs_old$val)
# points(min_abs_new$val,max_abs_new$val, col = rgb(1,0,0,0.5))
# abline(a=0,b=1)

dat <- as.data.frame(cbind(min_abs_old$val,max_abs_old$val))
dat2 <- as.data.frame(cbind(min_abs_new$val,max_abs_new$val))

colnames(dat) <- c("x","y")
colnames(dat2) <- c("x","y")



########  Long-term change

### Absolute T

tc_abs_an2 <- nc_open('CMIP/CMIP6 - Sea Surface Temperature (SST)  deg C - Long Term (2081-2100) SSP3 7 1981-2010 - Annual (mean of 24 models).nc')


lon <- ncvar_get(tc_abs_an2, "lon")
lat <- ncvar_get(tc_abs_an2, "lat", verbose = F)

val1 <- as.data.frame(ncvar_get(tc_abs_an2))

abs_an2 <- as.data.frame(matrix(NA,nrow = 360*180, ncol = 3))
colnames(abs_an2) <- c("lon","lat","val")
for(i in 1:180) {
  abs_an2$lon[(i-1)*360+(1:360)] <- lon
  abs_an2$lat[(i-1)*360+(1:360)] <- lat[i]
  abs_an2$val[(i-1)*360+(1:360)] <- val1[,i]
}

# ggplot(data = abs_an2, aes(x = lon, y = lat, colour = val)) + geom_point()

## Summer Change

tc_sum2 <- nc_open('CMIP/CMIP6 - Sea Surface Temperature (SST) Change deg C - Long Term (2081-2100) SSP3 7 1981-2010 - June to August (mean of 24 models).nc')



lon <- ncvar_get(tc_sum2, "lon")
lat <- ncvar_get(tc_sum2, "lat", verbose = F)

val1 <- as.data.frame(ncvar_get(tc_sum2))

dsum2 <- as.data.frame(matrix(NA,nrow = 360*180, ncol = 3))
colnames(dsum2) <- c("lon","lat","val")
for(i in 1:180) {
  dsum2$lon[(i-1)*360+(1:360)] <- lon
  dsum2$lat[(i-1)*360+(1:360)] <- lat[i]
  dsum2$val[(i-1)*360+(1:360)] <- val1[,i]
}

# ggplot(data = dsum2, aes(x = lon, y = lat, colour = val)) + geom_point()

## Winter Change

tc_win2 <- nc_open('CMIP/CMIP6 - Sea Surface Temperature (SST) Change deg C - Long Term (2081-2100) SSP3 7 1981-2010 - December to February (mean of 24 models).nc')
lon <- ncvar_get(tc_win2, "lon")
lat <- ncvar_get(tc_win2, "lat", verbose = F)

val2 <- as.data.frame(ncvar_get(tc_win2))

dwin2 <- as.data.frame(matrix(NA,nrow = 360*180, ncol = 3))
colnames(dwin2) <- c("lon","lat","val")
for(i in 1:180) {
  dwin2$lon[(i-1)*360+(1:360)] <- lon
  dwin2$lat[(i-1)*360+(1:360)] <- lat[i]
  dwin2$val[(i-1)*360+(1:360)] <- val2[,i]
}
# windows(8,5)
# ggplot(data = dwin2, aes(x = lon, y = lat, colour = val)) + geom_point()


## Annual change

tc_an2 <- nc_open('CMIP/CMIP6 - Sea Surface Temperature (SST) Change deg C - Long Term (2081-2100) SSP3 7 1981-2010 - Annual (mean of 24 models).nc')
lon <- ncvar_get(tc_an2, "lon")
lat <- ncvar_get(tc_an2, "lat", verbose = F)

val2 <- as.data.frame(ncvar_get(tc_an2))



dan2 <- as.data.frame(matrix(NA,nrow = 360*180, ncol = 3))
colnames(dan2) <- c("lon","lat","val")
for(i in 1:180) {
  dan2$lon[(i-1)*360+(1:360)] <- lon
  dan2$lat[(i-1)*360+(1:360)] <- lat[i]
  dan2$val[(i-1)*360+(1:360)] <- val2[,i]
}
# windows(8,5)
# ggplot(data = dan2, aes(x = lon, y = lat, colour = val)) + geom_point()


## Winter absolute
tc_win_abs2 <- nc_open("CMIP/CMIP6 - Sea Surface Temperature (SST)  deg C - Long Term (2081-2100) SSP3 7 1981-2010 - December to February (mean of 24 models).nc")
lon <- ncvar_get(tc_win_abs2, "lon")
lat <- ncvar_get(tc_win_abs2, "lat", verbose = F)

val2 <- as.data.frame(ncvar_get(tc_win_abs2))



win_abs2 <- as.data.frame(matrix(NA,nrow = 360*180, ncol = 3))
colnames(win_abs2) <- c("lon","lat","val")
for(i in 1:180) {
  win_abs2$lon[(i-1)*360+(1:360)] <- lon
  win_abs2$lat[(i-1)*360+(1:360)] <- lat[i]
  win_abs2$val[(i-1)*360+(1:360)] <- val2[,i]
}
# windows(8,5)
# ggplot(data = win_abs2, aes(x = lon, y = lat, colour = val)) + geom_point()


## Summer absolute
tc_sum_abs2 <- nc_open("CMIP/CMIP6 - Sea Surface Temperature (SST)  deg C - Long Term (2081-2100) SSP3 7 1981-2010 - June to August (mean of 24 models).nc")
lon <- ncvar_get(tc_sum_abs2, "lon")
lat <- ncvar_get(tc_sum_abs2, "lat", verbose = F)

val2 <- as.data.frame(ncvar_get(tc_sum_abs2))


sum_abs2 <- as.data.frame(matrix(NA,nrow = 360*180, ncol = 3))
colnames(sum_abs2) <- c("lon","lat","val")
for(i in 1:180) {
  sum_abs2$lon[(i-1)*360+(1:360)] <- lon
  sum_abs2$lat[(i-1)*360+(1:360)] <- lat[i]
  sum_abs2$val[(i-1)*360+(1:360)] <- val2[,i]
}
# windows(8,5)
# ggplot(data = sum_abs2, aes(x = lon, y = lat, colour = val)) + geom_point()




min_abs_new2 <- sum_abs2
min_abs_new2$val <- sapply(1:nrow(sum_abs2), function(x) min(c(sum_abs2$val[x],win_abs2$val[x])))
max_abs_new2 <- sum_abs2
max_abs_new2$val <- sapply(1:nrow(sum_abs2), function(x) max(c(sum_abs2$val[x],win_abs2$val[x])))
# plot(min_abs_new2$val,max_abs_new2$val)
# abline(a=0,b=1)

win_abs_old2 <- win_abs2
win_abs_old2$val <- win_abs2$val - dwin2$val
# plot(win_abs_old2$val,win_abs2$val)

sum_abs_old2 <- sum_abs2
sum_abs_old2$val <- sum_abs2$val - dsum2$val
# plot(sum_abs_old2$val,sum_abs2$val)

min_abs_old2 <- sum_abs2
min_abs_old2$val <- sapply(1:nrow(sum_abs2), function(x) min(c(sum_abs_old2$val[x],win_abs_old2$val[x])))
max_abs_old2 <- sum_abs2
max_abs_old2$val <- sapply(1:nrow(sum_abs2), function(x) max(c(sum_abs_old2$val[x],win_abs_old2$val[x])))
# plot(min_abs_old2$val,max_abs_old2$val)
# points(min_abs_new2$val,max_abs_new2$val, col = rgb(1,0,0,0.5))

dan_baseline <- abs_an2 - dan2

dat <- as.data.frame(cbind(max_abs_old$val,(max_abs_old$val-min_abs_old$val)))
dat2 <- as.data.frame(cbind(max_abs_new$val,(max_abs_new$val-min_abs_new$val)))
dat3 <- as.data.frame(cbind(max_abs_new2$val,(max_abs_new2$val-min_abs_new2$val)))

colnames(dat) <- c("x","y")
colnames(dat2) <- c("x","y")
colnames(dat3) <- c("x","y")

loess_fit <- loess(y ~ x , dat, span = 0.5)
loess_fit2 <- loess(y ~ x , dat2, span = 0.5)
loess_fit3 <- loess(y ~ x , dat3, span = 0.5)


#########################################################################################
###
###   Fig. 4a
###
###

### mean temperature


an_baseline <- abs_an2 - dan2

dat <- as.data.frame(cbind(max_abs_old$val,(an_baseline$val)))
#dat2 <- as.data.frame(cbind(max_abs_new$val,(max_abs_new$val-min_abs_new$val)))
dat3 <- as.data.frame(cbind(max_abs_new2$val,(abs_an2$val)))

colnames(dat) <- c("x","y")
#colnames(dat2) <- c("x","y")
colnames(dat3) <- c("x","y")

loess_fit <- loess(y ~ x , dat, span = 0.5)
#loess_fit2 <- loess(y ~ x , dat2, span = 0.5)
loess_fit3 <- loess(y ~ x , dat3, span = 0.5)

windows(3.4,3.3)

#png("Figs_06_2021/CMIP_warm_mean_annual.png", units = "in", width = 3.4, height = 3.3, res = 600)
#pdf("CMIP_warm_mean_annual.pdf",  width = 3.4, height = 3.3)

par(las = 1, mar = c(4.1,4.1,1.1,1.1), mgp = c(2.5,0.75,0))

plot(dat$x,dat$y, bg = rgb(0.5,0.5,0.5,0.2), pch = 21, col = NA, xlim = c(15,35), ylim = c(15,35), xaxs = "i", yaxs = "i", 
     xlab = expression("temperature warm season ("*degree~"C)"), ylab = expression("mean annual temperature ("*degree~"C)"), cex = 0.7, cex.lab = 0.9, cex.axis = 0.9)

#points(dat3$x,dat3$y, bg = rgb(1,0,0,0.01), pch = 21, col = NA,  cex = 0.7, cex.lab = 0.9, cex.axis = 0.9)


points(seq(min(dat$x, na.rm = T), max(dat$x, na.rm = T),0.1) ,
       stats::predict(loess_fit, newdata = seq(min(dat$x, na.rm = T), max(dat$x, na.rm = T),0.1)), col = rgb(0,0,0,1), cex = 0.5, type = "l", lwd = 2)

points(seq(min(dat2$x, na.rm = T), max(dat2$x, na.rm = T),0.1) ,
       stats::predict(loess_fit2, newdata = seq(min(dat2$x, na.rm = T), max(dat2$x, na.rm = T),0.1)), col = "orange", cex = 0.5, type = "l", lwd = 2)

points(seq(min(dat3$x, na.rm = T), max(dat3$x, na.rm = T),0.1) ,
       stats::predict(loess_fit3, newdata = seq(min(dat3$x, na.rm = T), max(dat3$x, na.rm = T),0.1)), col = "red", cex = 0.5, type = "l", lwd = 2)


legend("topleft", legend = c("Baseline (1981 - 2010)", "SSP3.7 (2081 - 2100)"), col = c("black", "red"), lwd = 2, cex = 0.75, bty = "n")

dev.off()


#########################################################################################
###
###   Fig. 4c
###
##
datb <- as.data.frame(cbind(max_abs_old$val,(min_abs_old$val)))
datb <- as.data.frame(cbind(max_abs_old$val,(min_abs_old$val)))
#dat2 <- as.data.frame(cbind(max_abs_new$val,(max_abs_new$val-min_abs_new$val)))
dat3b <- as.data.frame(cbind(max_abs_new2$val,(min_abs_new2$val)))

colnames(datb) <- c("x","y")
#colnames(dat2) <- c("x","y")
colnames(dat3b) <- c("x","y")

loess_fitb <- loess(y ~ x , datb, span = 0.5)
#loess_fit2 <- loess(y ~ x , dat2, span = 0.5)
loess_fit3b <- loess(y ~ x , dat3b, span = 0.5)

windows(3.4,3.3)
#png("Figs_06_2021/CMIP_warm_cold_mean_diff.png", units = "in", width = 3.4, height = 3.3, res = 600)
#pdf("CMIP_warm_cold_mean_diff.pdf", width = 3.4, height = 3.3)

par(las = 1, mar = c(4.1,4.1,1.1,1.1), mgp = c(2.5,0.75,0))

plot(   seq(min(dat$x, na.rm = T), max(dat$x, na.rm = T),0.1), 
       stats::predict(loess_fit3, newdata = seq(min(dat$x, na.rm = T), max(dat$x, na.rm = T),0.1))-
      stats::predict(loess_fit, newdata = seq(min(dat$x, na.rm = T), max(dat$x, na.rm = T),0.1)),
     type = "l", lwd = 2, xlim = c(-2,33), xaxs = "i", ylim = c(-1.65,0), yaxs = "i", col = "black",
     xlab = expression("constant"~italic("T")~"warm season ("*degree~"C)"), ylab = expression(Delta~italic("T")))

points(seq(min(datb$x, na.rm = T), max(datb$x, na.rm = T),0.1),
       stats::predict(loess_fit3b, newdata = seq(min(datb$x, na.rm = T), max(datb$x, na.rm = T),0.1)) -
         stats::predict(loess_fitb, newdata = seq(min(datb$x, na.rm = T), max(datb$x, na.rm = T),0.1)),
     type = "l", lwd = 2, col = "dodgerblue")

legend("bottomleft", legend = c(expression("mean annual T"),  expression("T cold season")), col = c("black", "dodgerblue"), lwd = 2, cex = 0.75, bty = "n")

dev.off()


#########################################################################################
###
###   Fig. 4b
###
##
windows(3.4,3.3)

#png("Figs_06_2021/CMIP_warm_cold.png", units = "in", width = 3.4, height = 3.3, res = 600)
#pdf("CMIP_warm_cold.pdf", width = 3.4, height = 3.3)

par(las = 1, mar = c(4.1,4.1,1.1,1.1), mgp = c(2.5,0.75,0))

plot(datb$x,datb$y, bg = rgb(0.5,0.5,0.5,0.2), pch = 21, col = NA, xlim = c(15,35), ylim = c(15,34.5), xaxs = "i", yaxs = "i", 
     xlab = expression("temperature warm season ("*degree~"C)"), ylab = expression("temperature cold season ("*degree~"C)"), cex = 0.7, cex.lab = 0.9, cex.axis = 0.9)

#points(dat3$x,dat3$y, bg = rgb(1,0,0,0.01), pch = 21, col = NA,  cex = 0.7, cex.lab = 0.9, cex.axis = 0.9)


points(seq(min(datb$x, na.rm = T), max(datb$x, na.rm = T),0.1) ,
       stats::predict(loess_fitb, newdata = seq(min(datb$x, na.rm = T), max(datb$x, na.rm = T),0.1)), col = rgb(0,0,0,1), cex = 0.5, type = "l", lwd = 2)

#points(seq(min(dat2$x, na.rm = T), max(dat2$x, na.rm = T),0.1) ,
#       stats::predict(loess_fit2, newdata = seq(min(dat2$x, na.rm = T), max(dat2$x, na.rm = T),0.1)), col = "orange", cex = 0.5, type = "l", lwd = 2)

points(seq(min(dat3b$x, na.rm = T), max(dat3b$x, na.rm = T),0.1) ,
       stats::predict(loess_fit3b, newdata = seq(min(dat3b$x, na.rm = T), max(dat3b$x, na.rm = T),0.1)), col = "red", cex = 0.5, type = "l", lwd = 2)
legend("topleft", legend = c("Baseline (1981 - 2010)", "SSP3.7 (2081 - 2100)"), col = c("black", "red"), lwd = 2, cex = 0.75, bty = "n")

dev.off()



