### Figures for Bivalve Manuscript
### Fig. 1, 3a
library(dplyr)

obiv <- readRDS("OBIS_bivalves_epifaunal_cleaned_2021_08_16_realms_and_grid.rds")


#### Everything with Cells 6
gdat <- as.data.frame(matrix(NA,nrow=length(unique(obiv$cells6)), ncol = 1))
gdat <- 
  obiv %>%
  group_by(cells6) %>%
  dplyr::summarise(n = n(), nara = sum(isaragonite), lon = unique(gridlon6), lat = unique(gridlat6),
            sstmean = mean(sstmean), sstrange = mean(sstrange),
            sstmin = mean(sstmin), sstmax = mean(sstmax), realm = unique(realm),

            n_ara_species = length(unique(species[which(isaragonite==1)])),n_species = length(unique(species))
  )

#table(gdat$realm)

evenness <- sapply(gdat$cells6, function(x) {
  spec = c(table(obiv$species[which(obiv$cells6 == x)]))
  specp = spec/sum(spec)
  -sum(specp*log(specp))/log(length(spec))
} )

gdat <- gdat %>% mutate(evenness = NA)
gdat$evenness <- evenness

gdat %>% filter(realm %in% 1:12, n_species >= 2, evenness >= 0.5)

###########
# World map mineralogy

### Plot cells
library(rgdal)
library(raster)
library(plotrix)

# Natural Earth shape files -- global (Robinson) projections
# get shapefiles from http://www.naturalearthdata.com
shape_path <- "map_data/"
coast_shapefile <- paste(shape_path, "ne_10m_coastline.shp", sep="")
ocean_shapefile <- paste(shape_path, "ne_10m_ocean.shp", sep="")
land_shapefile <- paste(shape_path, "ne_10m_land.shp", sep="")

# find out kind of shapefile (lines vs. polygons)
layer <- ogrListLayers(coast_shapefile)
ogrInfo(coast_shapefile, layer=layer)
coast_lines <- readOGR(coast_shapefile, layer=layer)

layer <- ogrListLayers(ocean_shapefile)
ogrInfo(ocean_shapefile, layer=layer)
ocean <- readOGR(ocean_shapefile, layer=layer)

layer <- ogrListLayers(land_shapefile)
ogrInfo(land_shapefile, layer=layer)
land <- readOGR(land_shapefile, layer=layer)


crs(ocean) <- crs(coast_lines)

# try to plot
outline <- bbox(coast_lines)
outline <- data.frame(xmin=outline["x","min"],
                      xmax=outline["x","max"],
                      ymin=outline["y","min"],
                      ymax=outline["y","max"])

x1x2 <- seq(180,-180,-0.05)
y1y2 <- seq(-90,90,0.05)
outline <- data.frame(cbind(c(x1x2,rep(180,length(y1y2)), -x1x2, rep(-180,length(y1y2))),
                            c(rep(-90,length(x1x2)), y1y2, rep(90,length(x1x2)), -y1y2)))
coordinates(outline) <- ~ X1 + X2

outline <- Polygon(outline)
outline <- Polygons(list(outline),1)
outline <- SpatialPolygons(list(outline))
crs(outline) <- crs(land)


### Realm colours:
alph <- 1*255
realmcolours1 <- (c(rgb(201,87,109,alph, maxColorValue = 255),
                   rgb(225,148,66,alph, maxColorValue = 255),
                   
                  rgb(129,94,205,alph, maxColorValue = 255),
                  rgb(174,173,67,alph, maxColorValue = 255),
                  rgb(81,180,141,alph, maxColorValue = 255),
                  
                  rgb(81,127,56,alph, maxColorValue = 255),
                  rgb(75,175,208,alph, maxColorValue = 255),
                  rgb(101,182,75,alph, maxColorValue = 255),
                  
                  rgb(200,101,185,alph, maxColorValue = 255),
                  
                  rgb(93,129,216,alph, maxColorValue = 255),
                  rgb(153,114,56,alph, maxColorValue = 255),
                  rgb(205,81,62,alph, maxColorValue = 255))[12:1])[c(1,7,6,3,11,10,
                                                                     12,2,8,9,5,4)]

alph <- 0.35*255

realmcolours2 <- (c(rgb(201,87,109,alph, maxColorValue = 255),
                   rgb(225,148,66,alph, maxColorValue = 255),
                   
                   rgb(129,94,205,alph, maxColorValue = 255),
                   rgb(174,173,67,alph, maxColorValue = 255),
                   rgb(81,180,141,alph, maxColorValue = 255),
                   
                   rgb(81,127,56,alph, maxColorValue = 255),
                   rgb(75,175,208,alph, maxColorValue = 255),
                   rgb(101,182,75,alph, maxColorValue = 255),
                   
                   rgb(200,101,185,alph, maxColorValue = 255),
                   
                   rgb(93,129,216,alph, maxColorValue = 255),
                   rgb(153,114,56,alph, maxColorValue = 255),
                   rgb(205,81,62,alph, maxColorValue = 255))[12:1])


### Cell colour preparation
library(icosa)
library(maptools)
library(sp)
grid6 <- hexagrid(c(6,6), sp = F) # 

dat <- gdat %>% filter(evenness>=0.5,n_species>1)
col_cell6 <- color.scale(dat$nara/dat$n,cs1=c(0.2,1),cs2=c(0.2,0.2),cs3=c(1,0.2),alpha=1)
grid6sp <- newsp(grid6)
spgrid6poly <- icosa::SpPolygons(grid6sp, res=50)
spgrid6polytrans <- sp::spTransform(spgrid6poly, CRS("+proj=moll +lon_0=0 +x_0=0 +y_0=0") )

### Run if not already loaded from the preparation file:
regions <- readOGR("Meow/meow_ecos.shp")
realms <- unionSpatialPolygons(regions, regions$REALM)
###

realms6polytrans <- sp::spTransform(realms, CRS("+proj=moll +lon_0=0 +x_0=0 +y_0=0") )


polycol <- rep(rgb(0,0,0,0),length(spgrid6polytrans@polygons))

cells6all <- dat$cells6
for(i in 1:length(dat$lon)) polycol[which(rownames(grid6@faces) == cells6all[i])] <- col_cell6[i]

#############################################################################
####
#### Fig. 3a

windows(8,4)
# png("Figs_06_2021/realms_map.png", units = "in", width = 12, height = 6, res = 600)
# pdf("realms_map.pdf", width = 8, height = 4)

par(mar = c(0,0,0,0))
plot(spTransform(land, CRS=("+proj=moll +lon_0=0 +x_0=0 +y_0=0")), col = "grey75", border = NA, add = F, xlim = c(-18000000,18000000), ylim = c(-9000000,9000000))


plot(realms6polytrans, add = T, border = NA,col = realmcolours2, lwd = 2)
plot(spTransform(outline, CRS=("+proj=moll +lon_0=0 +x_0=0 +y_0=0")), col = NA, border = rgb(0,0,0,0.5), add = T, lwd = 0.5)

text(-1.8*10^6,5700000, expression(bold("3")), col = realmcolours1[3], cex = 2) #
text(-0*10^6,9000000, expression(bold("1")), col = realmcolours1[1], cex = 2) #
text(13.1*10^6,2000000, expression(bold("8")), col = realmcolours1[8], cex = 2) #
text(13.2*10^6,5000000, expression(bold("2")), col = realmcolours1[2], cex = 2) # 
text(-11.5*10^6,5000000, expression(bold("2")), col = realmcolours1[2], cex = 2) #
text(-14.8*10^6,700000, expression(bold("4")), col = realmcolours1[4], cex = 2) #4
text(18.1*10^6,1200000, expression(bold("4")), col = realmcolours1[4], cex = 2) # 4
text(1.6*10^6,-6500000, expression(bold("12")), col = realmcolours1[12], cex = 2) #
text(10.75*10^6,-5100000, expression(bold("11")), col = realmcolours1[11], cex = 2)
text(-7*10^6,-5300000, expression(bold("9")), col = realmcolours1[9], cex = 2) #
text(0.2*10^6,-3500000, expression(bold("10")), col = realmcolours1[10], cex = 2)
text(-3.9*10^6,1300000, expression(bold("6")), col = realmcolours1[6], cex = 2)
text(-10.1*10^6,350000, expression(bold("5")), col = realmcolours1[5], cex = 2)
text(6.2*10^6,430000, expression(bold("7")), col = realmcolours1[7], cex = 2)


# dev.off()

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

#############################################################################
####
#### Fig. 1a

windows(12,6)
# png("Figs_06_2021/grid6_mineral_map.png", units = "in", width = 12, height = 6, res = 600)
plot(spTransform(land, CRS=("+proj=moll +lon_0=0 +x_0=0 +y_0=0")), col = "grey75", border = NA, add = F)


plot(realms6polytrans, add = T, border = realmcolours2,col = NA, lwd = 2)

plot(spgrid6polytrans, add = T, col = polycol,border = NA)
plot(spTransform(outline, CRS=("+proj=moll +lon_0=0 +x_0=0 +y_0=0")), col = NA, border = rgb(0,0,0,0.5), add = T, lwd = 0.5)
# dev.off()

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



plot_grid <- function(datout, gridres,xbreaks, ybreaks, colors, xlable, ylable, ncolour, legendat, legend, legendlegend) {
  
  layout(matrix(1:2,nrow=1),widths=c(3,1))
  par(las = 1, mgp=c(2.5,0.75,0), mar = c(4,4,1,0))
  
  plot(0,0, type = "n", xlim = xbreaks[c(1,length(xbreaks))], ylim = ybreaks[c(1,length(ybreaks))], xaxs = "i", yaxs = "i", mar = c(4,4,1,0),
       xlab = xlable, ylab = ylable)
  
  for(i in 1:gridres[1]) {
    for (j in 1:gridres[2]) {
      
      if (datout[i,j] == 0)  polygon(xbreaks[c(i,i,i+1,i+1)], ybreaks[c(j,j+1,j+1,j)], border = NA, col = NA)
      if (datout[i,j] > 0)  polygon(xbreaks[c(i,i,i+1,i+1)], ybreaks[c(j,j+1,j+1,j)], border = NA, col = colors[i+(j-1)*gridres[1]])
      
    }
  }
  
  abline(h = ybreaks[c(1,length(ybreaks))], v = xbreaks[c(1,length(xbreaks))])
  
  par(las = 1, mgp=c(2.5,0.75,0), mar = c(4,0.75,1,0))
  
  plot.new()
  
  sapply(1:(ncolour),function(x) polygon(c(-0,-0,0.33,0.33),c((x-1)/(ncolour),(x)/(ncolour),(x)/(ncolour),(x-1)/(ncolour)),
                                         border = NA, col = colorRampPalette(colpoints, alpha=T)(ncolour)[x]))
  sapply(1:length(legendat),function(x) points(c(0.33,0.42),rep(legendat[x],2), type = "l"))
  sapply(1:length(legendat),function(x) text(0.45,legendat[x], adj = c(0,0.5), labels = legend[x]))
  text(0.1667,1.05,legendlegend, xpd = T, adj = c(0.5,0.5))
}

grid_density <- function(x,y,gridres = c(10,10), lowbreak = 1) {
  xlimit = minmax(x)
  ylimit = minmax(y)
  
  xspace = xlimit[2]-xlimit[1]
  yspace = ylimit[2] - ylimit[1]
  xbreaks = seq(xlimit[1]-0.05*xspace/gridres[1], xlimit[2]+0.05*xspace/gridres[1], length.out = gridres[1]+1)
  ybreaks = seq(ylimit[1]-0.05*yspace/gridres[2], ylimit[2]+0.05*yspace/gridres[2], length.out = gridres[2]+1)
  
  gridbins = gridres
  
  datout = matrix(NA, nrow = gridbins[1], ncol = gridbins[2])
  
  for(i in 1:gridres[1]) {
    for (j in 1:gridres[2]) {
      
      datout[i,j] <- length(which(x >= xbreaks[i] & x<xbreaks[i+1] & y >= ybreaks[j] & y < ybreaks[j+1]))
      
    }
  }
  
  
  out <- list(datout,
              xbreaks,
              ybreaks
              
  )
  
  return(out)
}

minmax <- function(x) c(min(x),max(x))
#############################################################################
####
#### Fig. 1b and c

# crossplot T - Tr - colour
library(plyr)
library(dplyr)

# png("Mineralogy_heatmap.png", units = "in", width = 3.5, height = 3.5, res = 600)

gridres <- c(8,8)


out <- grid_density(dat$sstmean,dat$sstrange, gridres = gridres)
colpoints <- c(rgb(0,0,0,0.09), rgb(0,0,0,1))  

ii = cut(c(log(out[[1]])), breaks = seq(0, max(log(out[[1]]), na.rm = T), len = 400), 
         include.lowest = TRUE)


#colpoints <- c(rgb(0,0,1,0.72), rgb(1,0,0,0.72))  

colors <- colorRampPalette(colpoints, alpha=T)(400)[ii]


percara = matrix(NA, nrow = gridres[1], ncol = gridres[2])

for(i in 1:gridres[1]) {
  for (j in 1:gridres[2]) {
    
    percara[i,j] <- mean((dat$nara/dat$n)[which(dat$sstmean >= out[[2]][i] & dat$sstmean <out[[2]][i+1] & dat$sstrange >= out[[3]][j] & dat$sstrange < out[[3]][j+1])]
                           )
    
  }
}

colorsrgb <- matrix(0,nrow = 4, ncol = gridres[1]*gridres[2])
colorsrgb[,which(!(is.na(colors)))] <- sapply(which(!(is.na(colors))), function(x) col2rgb(colors[x], alpha = T))

colorsrgb[1,] <- c(percara)
colorsrgb[3,] <- c(1-percara)


colorsname <- rep(NA,gridres[1]*gridres[2])

colorsname[which(!(is.na(colors)))] <- rgb(colorsrgb[1,which(!(is.na(colors)))],colorsrgb[2,which(!(is.na(colors)))],
                                           colorsrgb[3,which(!(is.na(colors)))], colorsrgb[4,which(!(is.na(colors)))]/255)

#plot_grid(out[[1]],gridres,out[[2]],out[[3]], colorsname, expression("Temperature ("*degree*"C)"), expression("Temperature range ("*degree*"C)"), gridres[1]*gridres[2])



windows(5,4)
png("Grid6_t_trange_density.png", units = "in", width = 4.4, height = 3.5, res = 600)
colpoints <- c(rgb(0,0,0,0.09), rgb(0,0,0,1))  
plot_grid(out[[1]],gridres,out[[2]],out[[3]], colors, expression("mean temperature ("*degree*"C)"), expression("seasonality ("*degree*"C)"), gridres[1]*gridres[2],
          legend = c(1,3,10,30,84,"1,000","10,000"), legendat = log(c(1,3,10,30,84,1000,10000))/log(max(out[[1]],na.rm=T)),legendlegend = "N")
dev.off()

windows(5,4)
svg("Grid6_t_trange_para_density.svg",  width = 4.4, height = 3.5)
colpoints <- c(rgb(0,0,1,0.72), rgb(1,0,0,0.72))  
plot_grid(out[[1]],gridres,out[[2]],out[[3]], colorsname, expression("mean temperature ("*degree*"C)"), expression("seasonality ("*degree*"C)"), gridres[1]*gridres[2],
          legend = seq(0,1,0.2), legendat = seq(0,1,0.2),legendlegend = "P arag.")
dev.off()


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