#setwd("/Users/new_user/Desktop/for_thay")
library(ggplot2)
library(reshape2)
library(kohonen)
library(RColorBrewer); library(factoextra)
library(NbClust); library(ncdf4); library(tidyr)

#loading the data file, tab-delimited text file
in.file <- "spe12_pca_APH"
innc <- nc_open(paste(in.file,".nc",sep=""))
spe12 <- ncvar_get(innc,"spe12")
lon <- round(ncvar_get(innc,"lon"),digits = 3); rownames(spe12) <- lon; nlon <- dim(lon)
lat <- round(ncvar_get(innc,"lat"),digits = 3); colnames(spe12) <- lat; nlat <- dim(lat)
time <- ncvar_get(innc,"mode"); nmod <- dim(time)

colnames(spe12) <- c(lat)
rownames(spe12) <- c(lon)
new <- melt(spe12,id=colnames(spe12))
#new[is.na(new)] <- -9999
colnames(new) <- c("lon","lat","mod","value")
df <- spread(new,mod,value)
D.input <- na.omit(df)
colnames(D.input) <- c("lon","lat",paste("var",seq(15),sep=""))
nvar <- 15
D <- D.input[,3:(2+nvar)]
D <- data.frame(D.input$lon,D.input$lat,D)
#Z value - scalee values
Z <- scale(D,center=T,scale=T)


p <- fviz_nbclust(Z, kmeans, method = "wss") +
  geom_vline(xintercept = 4, linetype = 2)+
  #geom_vline(xintercept = 6, linetype = 2)+
  labs(subtitle = "Elbow method")+
  theme(                                                                  
        legend.title=element_text(size = 22, face = "bold"))

  
pdf(paste("new_optimal_k_Elbow_method.pdf",sep=""),height = 4, width = 4)
print(p)
dev.off()


p <- fviz_nbclust(Z, kmeans, method = "silhouette") +
  geom_vline(xintercept = 4, linetype = 2)+
  labs(subtitle = "Silhouette method")
pdf(paste("optimal_k_Silhouette_method.pdf",sep=""),height = 4, width = 4)
print(p)
dev.off()


# for k-mean clustering
for (iclus in seq(10)){
remove(out.dat)
max.clus <-iclus
km.res <- kmeans(Z, max.clus, nstart = 25)
ind.groupe <- as.matrix(km.res$cluster)
##print(ind.groupe)
out.dat <- data.frame(D.input$lon,D.input$lat,ind.groupe)
colnames(out.dat) <- c("lon","lat","group")

clus.out <- array(NA,dim=c(nlon,nlat)) 
colnames(clus.out) <- lat
rownames(clus.out) <- lon

for (ilon in seq(nlon)) {
  for (ilat in seq(nlat)){
    tmp.lon <- lon[ilon]
    tmp.lat <- lat[ilat]
    tmp <- subset(out.dat,out.dat$lon==tmp.lon&out.dat$lat==tmp.lat)
    #print(paste(lon[ilon],"_",lat[ilat],sep=""))
    #print(tmp)
    if (nrow(tmp)!=0){
      clus.out[ilon,ilat] <- as.integer(tmp[,3])
      #print(as.integer(tmp[,3]))
    }
  }
}


# write to ncdf file
dim_x <- ncdim_def('lon','degrees_east',lon)
dim_y <- ncdim_def('lat','degrees_north',lat)
var_out <- ncvar_def("cluster",units="",list(dim_x,dim_y),-9999)
ofile=paste("kmeans_",nvar,"vars_",in.file,"_clus",iclus,".nc",sep="")
ncid_out <- nc_create(file=ofile,list(var_out))
ncvar_put(ncid_out,var_out,clus.out)
nc_close(ncid_out)
}