5. Self-Organizing Maps (SOM)
Apply self-organizing maps (SOM) to the winter 3-day maximum precipitation over Southwest US. Try 3x3 or 4x4 node configuration. For each node composite the rainfall and the SSTs and show the maps. This will provide insights into the relationship between the SSTs and their signature on the rainfall extremes.

#### CLEAR MEMORY
rm(list=ls())
#### Prevent Warnings
options(warn=-1)
#### Initialize Graphics
graphics.off()
.pardefault <- par()

#### Source Libraries and set working directory
mainDir="/Users/annastar/OneDrive - UCB-O365/CVEN 6833 - Advanced Data Analysis Techniques/Homework/HW 02/"
setwd(mainDir)
suppressPackageStartupMessages(source("hw2_library.R"))
source("hw2_library.R")
# Unload packages that create a conflict with map("world2", add = T)
detach(package:mclust)
detach(package:sm)

nrows=72
ncols=36
ntime = 118    #Nov-Mar 1901 - Nov-Mar 2018

nyrs = ntime    #1902 - 2016

nglobe = nrows*ncols
N = nrows*ncols

### Lat - Long grid..
locs=matrix(scan("http://civil.colorado.edu/~balajir/CVEN6833/R-sessions/session2/files-4HW2/sst-lat-long.txt"), ncol=2, byrow=T)
ygrid=seq(-87.5,87.5,by=5)
ny=length(ygrid)

xgrid=seq(27.5,382.5,by=5)
#xgrid[xgrid > 180]=xgrid[xgrid > 180]-360  #longitude on 0-360 grid if needed
xgrid[xgrid > 180]=xgrid[xgrid > 180]
nx=length(xgrid)

xygrid=matrix(0,nrow=nx*ny,ncol=2)

i=0
for(iy in 1:ny){
  for(ix in 1:nx){
    i=i+1
    xygrid[i,1]=ygrid[iy]
    xygrid[i,2]=xgrid[ix]
  }
  
}

Read Kaplan SST data

data=readBin(paste(mainDir, "data/prob4/data.r4", sep = ""),what="numeric", n=( nrows * ncols * ntime), size=4,endian="swap")

data <- array(data = data, dim=c( nrows, ncols, ntime ) )

data1=data[,,1]

# the lat -long data grid..

index=1:(nx*ny)

index1=index[data1 < 20 & data1 != "NaN"]   # only non-missing data.
xygrid1=xygrid[index1,]

x1=xygrid1[,2]
#x1[x1 < 0]= x1[x1 < 0] + 360
#xygrid1[,2]=x1

nsites=length(index1)
data2=data1[index1]

### SSTdata matrix - rows are seasonal (i.e. one value per year)
## and columns are locations
sstdata=matrix(NA,nrow=nyrs, ncol=nsites)

for(i in 1:nyrs){
  data1=data[,,i]
  index1=index[data1 < 20 & data1 != "NaN"]
  data2=data1[index1]
  sstdata[i,]=data2
}

sstannavg = sstdata
indexgrid = index1
rm("data")  #remove the object data to clear up space

lon_sst = sort(unique(xygrid[,2]))
lat_sst = sort(unique(xygrid[,1]))

Read winter 3-day maximum precipitation over Southwest US

data = read.table(paste(mainDir, "data/Winter_temporal/Max_Winter_Seas_Prec.txt", sep = ""), header = T)
# Omit extreme values
non_values <- which(data[,-1] > 1000, arr.ind = T)
data[non_values] = NaN
data = na.omit(data)
loc = read.table(paste(mainDir, "data/Precipitaton_data.txt", sep = ""), header = T)
lon = loc[,1]
lat = loc[,2]

# # Convert precipitation data into spatial object
# coordinates(data) <- -lon+lat
# gridded(data) = TRUE
# precipstack <- stack(data)

Perform 3x3 SOM

n = 3
n_node = n*n
precip = data[,-1]
#index = complete.cases(precip)
#vals = precip[index,]
begin <- Sys.time()
SOM <- som(scale(precip), grid = somgrid(n,n,"rectangular"))
end <- Sys.time()
end-begin
## Time difference of 0.05239391 secs

Plot number of observations in each node

colors <- function(n, alpha = 1) {
  rev(heat.colors(n,alpha))
}
plot(SOM, type = "counts", palette.name = colors, heatkey = TRUE)

Display composited precipitation data into nodes

prec_nodes = as.data.frame(array(NaN, dim = c(dim(precip)[2],n_node)))  # Array: Location x Node
for (i in 1:n_node) {
  index = which(SOM$unit.classif == i)  # Determine which years belong to which node
  tmp_df = precip[index,]
  if (length(index) > 1) {
    prec_nodes[,i] = colMeans(tmp_df, na.rm = T)
  }
  else {
    prec_nodes[,i] = tmp_df
  }
}

zr = range(prec_nodes)
# Plot
nx = dim(precip)[2]
par(oma = c(0,2,3,6))
par(mfrow = c(n,n))
par(mar = c(3, 3, 2, 1))
for (i in 1:n_node) {
  quilt.plot(lon, lat, prec_nodes[,i], main = paste("Node", i, sep = " "), zlim = zr, add.legend = F, nx = nx, ny = nx)
  US(add = T, col = "gray30", lwd = 1)
}
mtext("Winter 3-Day Maximum Precipitation SOM, 3x3", outer = T, side = 3, cex = 1.2, line = 1)
par(oma = c(1,0,0,1))
par(mar = c(0,0,0,0))
image.plot(zlim = zr, legend.only = T, legend.shrink = 1, legend.width = 3)

Display composited SST data for each node

# unload package kohonen that creates a conflict with map("world2",add=T)
detach(package:kohonen)
# Display composited SST days for each node case: 3x3
SST_nodes=as.data.frame(array(NaN, dim = c(dim(sstannavg)[2],n_node)))
for(i in 1:n_node){
  index=which(SOM$unit.classif==i)
  tmp_df=sstannavg[index,]
  if(length(index)>1){
    SST_nodes[,i]=colMeans(tmp_df,na.rm=T)
  }else{
    SST_nodes[,i]=tmp_df
  }
}

zr=range(SST_nodes)
Nx=50
par( oma=c(0,0,2.5,1)) # save some room for the legend
par(mfrow = c(n, n))
#set.panel(4,4)
par(mar = c(2, 3, 2, 1))
for(i in 1:n_node){
  zfull = rep(NaN,nglobe)   #also equal 72*36
  zfull[indexgrid]=SST_nodes[,i]
  zmat = matrix(zfull,nrow=nrows,ncol=ncols)
  image.plot(lon_sst,lat_sst,zmat,ylim=range(-40,70),lim=zr,main = paste("Node ",i,sep=""))
  contour(lon_sst,lat_sst,(zmat),ylim=range(-40,70),add=TRUE,nlev=6,lwd=1)
  map("world2",add=T)
}
mtext("Composite SST SOM Maps 3x3", outer = TRUE, side = 3, cex = 1.2, line = 1)

Discussion