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

# read input data
in.file <- paste("pca_APH.nc",sep="")
innc <- nc_open(in.file)
ivar.name <- "spe12"
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")
nmod <- dim(mod)
max.clus <- 11 # 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)) {
          if (is.na(idat[ilon,ilat,1])==FALSE) {
              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("kmeans_",iclus,".nc",sep="")
     ncid_out <- nc_create(file=ofile,list(var_out))
     ncvar_put(ncid_out,var_out,cluster)
     nc_close(ncid_out)
}