#setwd("/Users/new_user/Documents/For_Work/2021/Drought_Thay_Tan")
library(ncdf4)
library(ggplot2)
#library(reshape2)
library(factoextra)
library(cluster)
library(NbClust)
library(corrplot)
library(maps)
library(mapdata);library(dplyr)
library(mapproj)
library(viridis)
library(ggthemes)
#library(RColorBrewer)

# Input/output Dir:
Indir <- "" #"/home/tanpv/Progs/2021/NC/"
Outtxt <- "" #"/home/tanpv/Progs/2021/NC/"
Outpdf <-""  #"/home/tanpv/Progs/2021/PDF/"

# read input data
in.file <- paste(Indir,"HWday.nc",sep="")
innc <- nc_open(in.file)
ivar.name <- "hwday"
idat <- ncvar_get(innc,ivar.name)
lon <- round(ncvar_get(innc,"lon"),digits = 4)
nlon <- dim(lon)
lat <- round(ncvar_get(innc,"lat"),digits = 4)
nlat <- dim(lat)

max.clus <- 7 # so cum lon nhat
cluster <- array(NA,dim=c(nlon,1))

# for k-means clustering

r2dat.new <- idat
for (iclus in seq(max.clus)) {
     km.res <- kmeans(r2dat.new, iclus, nstart = 25)
     sta.data <- seq(nrow(r2dat.new))
     dd <- cbind(sta.data,lon,lat,cluster = km.res$cluster)
 
     # write to ncdf file
     write.table(dd, file=paste(Outtxt,"cluster_",iclus,"groups.txt",sep=""),col.names = T, row.names = F,sep="\t")
 
     # for draw
     map('world',col="white", fill=TRUE, bg="white", lwd=0.05,boundary = TRUE,
         mar=rep(0,4),border=1.2, ylim=c(8,25),xlim = c(100,115))
     
     
     # Get the world polygon and extract UK
     library(maps)
     VN <- map_data("world") %>% filter(region=="Vietnam")
     data <- world.cities %>% filter(country.etc=="Vietnam")
     dd <- as.data.frame(dd)
     p <- ggplot() +
       geom_polygon(data = VN, aes(x=long, y = lat, group = group), fill="grey", alpha=1,color="black") +
       geom_point(data=dd, aes(x=lon, y=lat,shape=as.factor(cluster),col=cluster)) +
       xlab("")+ylab("")+
       #scale_shape_manual(values=c(25, 21, 24))+
       scale_size_continuous(range=c(1.5,3.2)) +coord_map() +
       scale_x_continuous(breaks=seq(102,116,2),limits = c(102,116))+
       scale_y_continuous(breaks=seq(8,24,2))+
       guides(shape = FALSE,size=FALSE)+ ggtitle(paste("cluster",iclus,sep=""))
     pdf(paste("Fig_kmeans_",iclus,"group.pdf",sep=""), width = 4, height = 5)
     print(p)
     dev.off()    
}

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.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.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.pdf",sep=""),height = 4, width = 4)
print(p)
dev.off()