library(ncdf4)
library(ggplot2)
library(reshape2)
library(factoextra)
library(cluster)
library(NbClust)
library(corrplot)
library(RColorBrewer)

# Input/output Dir:
#Indir <- "/home/tanpv/Progs/Dr_SEA/EOF_New/"
#Outnc <- "/home/tanpv/Progs/Dr_SEA/Cluster/NC/"
indir <- "/home/tanpv/Progs/Dr_SEA/Cluster/"
Outnc <- "/home/tanpv/Progs/Dr_SEA/Cluster/NC1/"
Outpdf <- "/home/tanpv/Progs/Dr_SEA/PDF/"
#indir <- ""
#Outnc <- ""
#Outpdf <- ""

# read input data
#in.file <- paste(Indir,"pca_APH.nc",sep="")
in.file <- paste("Dr_Chars_4_Clus.nc",sep="")
innc <- nc_open(in.file)
#ivar.name <- "spe12"
ivar.name <- "para"
idat <- ncvar_get(innc,ivar.name)
lon <- ncvar_get(innc,"lon")
nlon <- dim(lon)
lat <- ncvar_get(innc,"lat")
nlat <- dim(lat)
#mod <- ncvar_get(innc,"mode")
mod <- ncvar_get(innc,"par")
nmod <- dim(mod)
max.clus <- 13 # so cum lon nhat
cluster <- array(NA,dim=c(nlon,nlat))

npoint <- nlon*nlat
itmp.lon <- array(NA,dim=c(npoint,1)) 
itmp.lat <- array(NA,dim=c(npoint,1))
r2dat <- array(NA,dim=c(npoint,nmod))



# convert 3d-2d data
itmp <- 0
for (ilon in seq(nlon)) {
     for (ilat in seq(nlat)) {
          check.dat <- idat[ilon,ilat,]
          na.flag <- nrow(as.matrix(which(is.na(check.dat))))
          if (na.flag==0) {
              itmp <- itmp + 1
              r2dat[itmp,] <- idat[ilon,ilat,]
              itmp.lon[itmp] <- ilon
              itmp.lat[itmp] <- ilat
          }
     }
}

ntmp <- itmp
r2dat.new <- r2dat[1:ntmp,]
# for k-means clustering



for (iclus in seq(max.clus)) {
     km.res <- kmeans(r2dat.new, iclus, nstart = 25)
     sta.data <- seq(nrow(r2dat.new))
     dd <- cbind(sta.data,cluster = km.res$cluster)


     # put result data

     for (itmp in seq(ntmp)) {
          cluster[itmp.lon[itmp],itmp.lat[itmp]] <- dd[itmp,2]
     }
     
     # 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(Outnc,"kmeans_",iclus,".nc",sep="")
     ncid_out <- nc_create(file=ofile,list(var_out))
     ncvar_put(ncid_out,var_out,cluster)
     nc_close(ncid_out)
}

idata <- r2dat.new
p <- fviz_nbclust(idata, 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(Outpdf,"new_optimal_k_Elbow_method_1.pdf",sep=""),height = 4, width = 4)
print(p)
dev.off()


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

#p <- fviz_nbclust(idata, kmeans, nstart = 25,  method = "gap_stat", nboot = 25)+
#  labs(subtitle = "Gap statistic method")
#pdf(paste(Outpdf,"optimal_k_gap_stat_method_1.pdf",sep=""),height = 4, width = 4)
#print(p)
#dev.off()