#setwd("/Users/new_user/Documents/For_Work/PDSI")
library(ncdf4); library(ggplot2)
library(reshape2); library(factoextra)
#library(NbClust); library(corrplot); library(RColorBrewer)
nc <- nc_open("/home/tanpv/SEACLID_New/RES/PDSI_Res/PDSI_Stn_No_Rain_Center.nc")
var.name <- "T2m"
lon <- round(ncvar_get(nc,varid = "lon"),digits=2)
lat <- round(ncvar_get(nc,varid = "lat"),digits=2)
reg <- round(ncvar_get(nc,varid = "reg"),digits=0)
ivar <- ncvar_get(nc,varid = var.name)

# for k-mean clustering
nsta <- dim(lon)
sta.id <- seq(nsta)
odat <- data.frame(sta.id,lon,lat,reg)
for (iclus in seq(10)) {
     #max.clus <- 4
     km.res <- kmeans(ivar, iclus, nstart = 25)
     odat <- data.frame(odat,km.res$cluster)
}

colnames(odat) <- c("id","lon","lat","reg",seq(10))

write.table(odat,paste("kmean_",var.name,".txt",sep=""),
            col.names = T,row.names = F,append = F,sep = "\t")