Cacher is the code snippet organizer for pro developers

library(maps)
library(mapproj)
library(mapplots)
library(maptools)
library(mapdata)
library(RNetCDF)
library(fields)
library(RColorBrewer)
library(colorRamps)
library(SDMTools)
library(raster)
library(ncdf)
library(tidyr)
library(dplyr)
library(ggplot2)
library(reshape2)
library(magrittr)
library(igraph)
library(stringr)


# Ajout surface climatological dynamic height (from CARS2009) 
fname <- ("~/Downloads/hgt2000_cars2009a.nc")
dat <- open.nc(fname)
fid <- read.nc(dat)
mean <- fid$mean
mean[mean==-32766] <- NA
lat <- fid$lat
lon <- fid$lon
close.nc(dat)

fid$mean[fid$mean==-32766] <- NA
zmin <- min(fid$mean[!is.na(fid$mean)])
zmax <- 3


k <- which(lon==180)
px <- c(k:length(lon),1:(k-1))
lon2 <- seq(-180,180,0.5)
mean2 <- mean[px,]
#zmin=0.

#Ajout des coordonnees, polar circle ou pas
coord <- read.table("~/Downloads/Station_coordinates_polarcircle.csv",header=T,stringsAsFactors=FALSE,row.names=1)

#Ajout des tables station/temperature
a <- read.table("~/Downloads/stations_highway_NA.data",row.names=1,sep='\t')


#cls <- a %>%
#  merge(coord, by="row.names") %>%
#  mutate(scl = (V2 - min(V2, na.rm=T))/(max(V2, na.rm=T) - min(V2, na.rm=T)))


a <- rownames(a) %>% 
  as.data.frame() %>% 
  set_colnames("X") %>% 
  mutate(num = str_extract(X, "^([0-9])+"), layer =  str_extract(X, "([A-Z])+"), V2 = a$V2) %>%
  unite(station, c(num, layer), remove = TRUE) %>% 
  set_rownames(rownames(a)) %>%
  dplyr::select(-X)

coord <- rownames(coord) %>% 
  as.data.frame() %>% 
  set_colnames("X") %>% 
  mutate(num = str_extract(X, "^([0-9])+"), layer =  str_extract(X, "([A-Z])+"), Mean_latitude = coord$Mean_latitude, Mean_longitude = coord$Mean_longitude) %>%
  unite(station, c(num, layer), remove = TRUE) %>% 
  set_rownames(rownames(coord)) %>%
  dplyr::select(-X)

dimnames(mean2) <- list(long = lon2, lat = lat)
ret <- reshape::melt(mean2)

# read files
shared_otu <- tbl_df(read.table("~/Desktop/forDaniel/data/north_atlantic_pairwise_OTUs.2016_04_24.txt", header = T, sep = "\t")) %>%
  dplyr::select(station_1, station_2, common_OTUs) %>%
  set_colnames(c("from", "to", "weight"))

shared_otu <- shared_otu %>% 
  graph_from_data_frame() %>%
  simplify(remove.loops = TRUE)

shared_otu_vids <- V(shared_otu)$name


metadata <- as.data.frame(shared_otu_vids) %>%
  set_colnames("station") %>%
  left_join(coord) %>%
  dplyr::select(station, Mean_latitude, Mean_longitude) %>%
  set_colnames(c("station", "latitude", "longitude"))


# Some accessory functions
get_edges <- function(dist_filt = dist_filt, layout_coordinates = layout_coordinates){
  
  edge_maker <- function(which_row, len = 100, curved = FALSE, layout_coordinates = layout_coordinates, adjacency_list = adjacency_list){
    # Get coordinates for each metagenome
    from_c <- as.numeric(subset(layout_coordinates, station == as.character(adjacency_list[which_row,]$X))[2:3])  # Origin
    to_c <- as.numeric(subset(layout_coordinates, station == as.character(adjacency_list[which_row,]$Y))[2:3])  # Terminus
    
    # Add curve:
    graph_center <- colMeans(layout_coordinates[2:3])  # Center of the overall graph
    bezier_mid <- c(from_c[1], to_c[2])  # A midpoint, for bended edges
    distance1 <- sum((graph_center - bezier_mid)^2)
    if(distance1 < sum((graph_center - c(to_c[1], from_c[2]))^2)){
      bezier_mid <- c(to_c[1], from_c[2])
    }  # To select the best Bezier midpoint
    bezier_mid <- (from_c + to_c + bezier_mid) / 2.5  # Moderate the Bezier midpoint
    if(curved == FALSE){bezier_mid <- (from_c + to_c) / 2}  # Remove the curve
    
    edge <- data.frame(bezier(c(from_c[1], bezier_mid[1], to_c[1]),  # Generate
                              c(from_c[2], bezier_mid[2], to_c[2]),  # X & y
                              evaluation = len))  # Bezier path coordinates
    edge$Sequence <- 1:len  # For size and colour weighting in plot
    edge$Group <- paste(adjacency_list[which_row,]$X,adjacency_list[which_row,]$Y, sep = "#")
    edge$Similarity <- adjacency_list[which_row,]$dist
    return(edge)
  }
  
  dist_filt<-dist_filt %>%
    setnames(c("X", "Y", "dist"))
  
  layout_coordinates<-layout_coordinates %>%
    setnames(c("station", "latitude", "longitude"))
  
  # Generate a (curved) edge path for each pair of connected nodes
  all_edges <- lapply(1:nrow(dist_filt), edge_maker, 
                      len = 100, curved = TRUE, adjacency_list = dist_filt, 
                      layout_coordinates = layout_coordinates)
  all_edges <- rbindlist(all_edges, fill = T)  # a fine-grained path ^, with bend ^
  return(all_edges)
}

plot_on_map_mst <- function(all_edges = all_edges, metadata = metadata, name = name){
  
  xlim=c(-100,45)
  ylim=c(-27,50)
  
  
  worldmap <- map_data("world") %>%
    set_colnames(c("X","Y","PID","POS","region","subregion")) %>%
    clipPolys(xlim=xlim,ylim=ylim, keepExtra=TRUE)
  
  zp1 <- ggplot(all_edges) +
    geom_raster(data = ret %>% filter(value <= 3), 
                aes(x =long, y = lat, fill = value), 
                interpolate = F,na.rm = TRUE) + 
    coord_map(xlim=xlim,ylim=ylim) + 
    geom_polygon(data=worldmap,aes(X,Y,group=PID),
                 fill = "#0a0c2a",color="#0a0c2a") +
    theme_bw()+ #scale_fill_gradient( space = "Lab", low = "#edf8b1", high = "#2c7fb8")  +
    scale_fill_gradientn(colours = (colorRampPalette(brewer.pal(9,"Blues"))(25)), na.value = "white" ) +
  #scale_fill_gradientn(colours = rev(viridis(n = 2, option = "B")), na.value = "white" ) +
    #scale_fill_gradient2() +
    coord_fixed(1.5) + xlim(xlim) + ylim(ylim) 
  # Pretty simple plot code
  zp1 <- zp1 + geom_path(aes(x = y, y = x, group = Group,  # Edges with gradient
                             colour = Similarity, size = -Sequence/200), alpha = 1)  # and taper
  zp1 <- zp1 + scale_size(range = c(0.2, 0.6), guide = "none")  # Customize taper
  zp1 <- zp1 + scale_colour_gradientn(colours = rev(colorRampPalette(brewer.pal(11,"Spectral"))(25)), name="Shared OTUs")
  #scale_colour_gradientn(colours = ())
  zp1 <- zp1 + theme(panel.background = element_rect(fill = "white", colour = "black"), 
                     panel.grid.major = element_blank(),
                     panel.grid.minor = element_blank(),
                     axis.ticks = element_blank(), 
                     axis.text.x = element_text (size = 12, vjust = 0), 
                     axis.text.y = element_text (size = 12, hjust = 1.3)) + 
    coord_cartesian() + 
    labs(y="",x="") 
  zp1 <- zp1
  
  zp1<-zp1 + geom_point(data = metadata ,  
                        aes(x = longitude, y = latitude), size = 1.5, color = "#149067", alpha = 0.8) # Customize gradient v
  zp1 <- zp1 + theme_bw() + coord_equal(ratio=1) + theme(legend.position="right") +
    xlab("") + ylab("")
  ggtitle(bquote(atop(.(name), atop(italic(.(name)), "")))) #
  return(zp1)
}



dist <- induced_subgraph(shared_otu, vids =  metadata$station)
#dist<- delete_vertices(dist, V(dist)[degree(dist) == 0] )
#dist <- delete_edge_attr(dist, name = "weight_orig")
tara_all_edges<-get_edges(dist_filt = dist %>% get.data.frame, layout_coordinates = metadata %>% dplyr::select(station, latitude, longitude))

l<-plot_on_map_mst(tara_all_edges, metadata = metadata, name = "Shared OTUs North Atlantic (0.8-5)") 

print(l)