Packages

library("raster")  
library("RStoolbox")
library("ggplot2")
library("gridExtra")
library("ggpubr")
library("rgdal")
library("rasterVis")

Loading the rasterstack

vars <- stack(list.files(path="../input", pattern='.tif$', full.names=T),RAT=F)
names(vars) <- c("Bathymetry", "Current speed", "OCAR", "OC density", "OPD", "Orb Velocity")
vars
## class      : RasterStack 
## dimensions : 2686, 2559, 6873474, 6  (nrow, ncol, ncell, nlayers)
## resolution : 500, 500  (x, y)
## extent     : 3230952, 4510452, 3023466, 4366466  (xmin, xmax, ymin, ymax)
## crs        : +proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs 
## names      :    Bathymetry, Current.speed,          OCAR,    OC.density,           OPD,  Orb.Velocity 
## min values : -857.74938965,    0.00000000,    0.01553161,    1.15797031,    0.93182999,    0.00000000 
## max values :     50.689999,      2.025045,     66.181312,     13.595368,      2.928485,      2.764714
plot(vars)

Normalise raster layers prior to PCA

Performs normalisation (scaling) and centering, i.e. divide by standard deviation

vars.norm <- normImage(vars, norm = TRUE)
vars.norm
## class      : RasterBrick 
## dimensions : 2686, 2559, 6873474, 6  (nrow, ncol, ncell, nlayers)
## resolution : 500, 500  (x, y)
## extent     : 3230952, 4510452, 3023466, 4366466  (xmin, xmax, ymin, ymax)
## crs        : +proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs 
## source     : memory
## names      : Bathymetry, Current.speed,       OCAR, OC.density,        OPD, Orb.Velocity 
## min values : -8.8266827,    -1.2770513, -0.3782746, -1.3081598, -1.2606624,   -1.3171126 
## max values :   1.581071,     11.282040,   9.006722,   3.706218,   1.590953,     8.828313
hist(vars.norm)

Run PCA

https://www.rdocumentation.org/packages/RStoolbox/versions/0.2.3/topics/rasterPCA

set.seed(25)
rpc <- rasterPCA(vars.norm, nSamples = NULL, nComp = nlayers(vars.norm), spca = TRUE, maskCheck = TRUE)
rpc
## $call
## rasterPCA(img = vars.norm, nSamples = NULL, nComp = nlayers(vars.norm), 
##     spca = TRUE, maskCheck = TRUE)
## 
## $model
## Call:
## princomp(cor = spca, covmat = covMat[[1]])
## 
## Standard deviations:
##    Comp.1    Comp.2    Comp.3    Comp.4    Comp.5    Comp.6 
## 1.9554992 0.8714080 0.8033553 0.7099885 0.4073997 0.3181713 
## 
##  6  variables and  6873474 observations.
## 
## $map
## class      : RasterBrick 
## dimensions : 2686, 2559, 6873474, 6  (nrow, ncol, ncell, nlayers)
## resolution : 500, 500  (x, y)
## extent     : 3230952, 4510452, 3023466, 4366466  (xmin, xmax, ymin, ymax)
## crs        : +proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs 
## source     : memory
## names      :       PC1,       PC2,       PC3,       PC4,       PC5,       PC6 
## min values : -8.756754, -9.738404, -5.617817, -2.021233, -4.226620, -2.395745 
## max values :  5.620377,  3.758767,  1.541429,  4.243877,  4.101794,  2.037000 
## 
## 
## attr(,"class")
## [1] "rasterPCA" "RStoolbox"
summary(rpc$model)
## Importance of components:
##                           Comp.1    Comp.2    Comp.3     Comp.4     Comp.5
## Standard deviation     1.9554992 0.8714080 0.8033553 0.70998845 0.40739969
## Proportion of Variance 0.6373295 0.1265587 0.1075633 0.08401393 0.02766242
## Cumulative Proportion  0.6373295 0.7638882 0.8714515 0.95546542 0.98312784
##                            Comp.6
## Standard deviation     0.31817128
## Proportion of Variance 0.01687216
## Cumulative Proportion  1.00000000
rpc$model$loadings
## 
## Loadings:
##               Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Comp.6
## Bathymetry     0.466  0.150  0.250  0.124  0.732  0.383
## Current.speed  0.314 -0.834         0.425 -0.143       
## OCAR          -0.408 -0.159 -0.639  0.225  0.565 -0.177
## OC.density    -0.489 -0.103               -0.110  0.859
## OPD            0.373 -0.233 -0.474 -0.745         0.162
## Orb.Velocity   0.374  0.439 -0.552  0.444 -0.336  0.231
## 
##                Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Comp.6
## SS loadings     1.000  1.000  1.000  1.000  1.000  1.000
## Proportion Var  0.167  0.167  0.167  0.167  0.167  0.167
## Cumulative Var  0.167  0.333  0.500  0.667  0.833  1.000
hist(rpc$map)

plot(rpc$map)

ggRGB(rpc$map,1,2,3, stretch="lin", q=0)

Specify number of PCs to be used

n.pc <- 4

kmeans classification

https://www.rdocumentation.org/packages/RStoolbox/versions/0.2.3/topics/unsuperclass

Specify the maximum number of clusters to be analysed

max_cl <- 20

How many clusters?

df_total <- data.frame()
for (i in 1:max_cl) {
  set.seed(25)
  unC <- unsuperClass(rpc$map[[1:n.pc]], nSamples = 4000, nClasses = i, nStarts = 100, 
                      norm = FALSE, clusterMap = TRUE, algorithm = "Hartigan-Wong")
  
  stat <- 100*(unC[["model"]][["betweenss"]]/unC[["model"]][["totss"]])
  df <- data.frame(i,stat)
  df_total <- rbind(df_total,df)
}

Plot graph

ggplot(df_total,aes(i,stat)) + geom_bar(stat = "identity") +
  xlab("Number of clusters") + ylab("between sum of squares/total sum of squares ") +
  scale_x_continuous(breaks = c(1:max_cl))

Specify the number of clusters to be used

based on the graph above

n.cl <- 3

Run k-means clustering

set.seed(25)
unC <- unsuperClass(rpc$map[[1:n.pc]], nSamples = 10000, nClasses = n.cl, nStarts = 100, 
                    norm = FALSE, clusterMap = TRUE, algorithm = "Hartigan-Wong")
unC
## unsuperClass results
## 
## *************** Map ******************
## $map
## class      : RasterLayer 
## dimensions : 2686, 2559, 6873474  (nrow, ncol, ncell)
## resolution : 500, 500  (x, y)
## extent     : 3230952, 4510452, 3023466, 4366466  (xmin, xmax, ymin, ymax)
## crs        : +proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs 
## source     : memory
## names      : layer 
## values     : 1, 3  (min, max)
unC$model
## K-means clustering with 3 clusters of sizes 4401, 4543, 1056
## 
## Cluster means:
##          PC1        PC2        PC3          PC4
## 1 -0.4722387  0.1900850  0.6077706 -0.031614993
## 2  1.5152840 -0.1132621 -0.3727417 -0.008600807
## 3 -4.4595705 -0.3890747 -0.9077637  0.155511737
## 
## Clustering vector:
##     [1] 1 1 1 2 1 1 1 3 1 3 1 1 2 1 1 1 2 2 2 2 1 1 2 2 3 3 3 1 2 1 2 1 1 3 1 3
##    [37] 2 1 1 2 1 3 1 1 1 1 2 2 1 2 1 2 2 2 2 1 1 2 1 2 1 3 1 2 1 1 2 2 1 3 2 1
##    [73] 1 2 3 2 2 1 1 3 1 3 1 1 2 1 1 2 1 2 2 1 1 1 3 2 1 2 1 1 2 1 2 2 3 1 1 2
##   [109] 2 1 2 2 1 1 1 2 2 2 1 3 1 1 2 3 2 2 1 2 2 2 1 2 1 2 2 1 2 2 1 2 2 2 3 1
##   [145] 2 2 1 1 2 2 1 2 2 1 1 2 1 1 1 2 2 1 1 2 1 2 1 3 1 2 1 1 2 2 2 2 1 2 1 1
##   [181] 2 3 2 1 1 2 2 2 1 2 2 1 1 2 1 2 2 1 1 1 1 2 3 1 1 2 1 3 2 3 3 2 1 2 1 2
##   [217] 3 2 3 2 1 1 1 2 2 2 2 1 2 1 3 1 2 1 2 1 1 2 1 3 1 1 1 2 1 3 2 3 2 3 1 2
##   [253] 2 1 1 1 2 2 2 1 1 2 1 1 3 2 3 2 1 2 1 1 2 2 1 2 3 2 2 2 1 2 2 2 2 2 1 2
##   [289] 1 1 3 3 1 1 1 1 2 2 3 2 1 2 2 2 1 2 2 2 1 2 2 1 1 2 2 2 3 1 3 1 1 2 1 2
##   [325] 1 1 1 2 1 3 3 1 1 2 1 1 1 2 2 1 2 3 2 2 3 1 2 2 1 2 1 1 1 2 2 1 2 1 1 1
##   [361] 2 1 2 2 3 2 2 1 2 2 2 2 2 2 3 1 2 1 1 1 2 2 2 2 1 2 1 2 2 1 1 1 2 1 1 1
##   [397] 3 1 3 1 1 2 1 2 1 2 1 1 2 1 1 1 2 2 2 1 1 1 3 2 1 1 1 2 1 1 2 2 2 1 2 2
##   [433] 1 1 2 1 2 2 2 2 2 1 1 1 3 1 3 1 1 1 1 1 1 1 2 2 2 2 2 1 1 3 1 1 2 2 1 1
##   [469] 3 1 1 2 2 1 2 1 1 2 2 1 2 2 2 1 1 2 2 2 2 2 1 2 3 1 2 2 2 3 3 1 3 2 1 2
##   [505] 1 2 3 2 1 2 1 1 1 1 1 1 2 2 2 2 1 2 1 2 2 1 3 1 3 2 1 1 1 2 2 2 3 3 2 3
##   [541] 1 1 2 1 1 1 1 2 1 3 1 1 1 1 2 2 2 1 3 2 1 2 1 1 1 2 1 1 1 2 3 2 2 2 2 1
##   [577] 2 1 2 1 2 2 1 2 1 1 2 1 2 2 2 2 2 1 1 2 2 3 2 1 3 1 2 1 2 2 1 1 3 1 3 2
##   [613] 1 1 1 1 1 2 1 2 1 3 1 1 2 2 1 2 2 3 2 1 1 2 2 1 1 2 1 3 2 2 2 3 2 1 1 1
##   [649] 1 1 2 2 2 2 1 2 1 1 1 1 1 1 1 2 2 2 1 2 2 2 2 2 2 3 1 1 1 1 3 1 2 1 1 1
##   [685] 1 2 2 2 2 2 1 2 2 2 1 1 1 1 1 1 2 1 1 2 2 3 3 1 1 1 1 2 1 2 2 2 2 1 1 1
##   [721] 1 2 2 3 1 2 2 1 2 2 3 1 2 1 1 1 3 2 2 2 2 1 2 3 2 2 2 1 2 1 1 1 2 2 2 1
##   [757] 1 1 1 1 1 1 1 1 1 2 1 2 3 1 1 2 2 2 2 3 2 1 2 2 1 1 1 2 3 1 1 2 1 2 3 1
##   [793] 2 3 1 1 2 2 2 1 2 1 1 2 1 2 2 1 1 2 1 2 3 1 2 1 1 1 2 1 1 1 2 1 1 1 2 2
##   [829] 3 1 1 1 2 2 1 3 1 3 1 1 1 2 2 2 1 1 3 2 2 1 1 1 2 2 1 1 2 2 2 1 1 1 1 1
##   [865] 1 1 2 2 1 1 1 2 1 2 1 2 2 2 1 3 2 1 1 1 2 2 1 1 2 1 2 2 2 1 3 2 2 1 2 1
##   [901] 2 2 3 1 1 2 3 1 1 2 2 2 2 1 1 1 3 1 2 1 1 3 3 1 2 3 2 2 1 1 1 1 1 1 1 1
##   [937] 2 2 1 3 1 2 2 3 2 2 2 2 2 1 1 2 1 1 2 2 3 1 1 1 3 2 1 2 2 2 2 2 1 1 2 2
##   [973] 2 2 2 1 2 2 2 1 1 1 1 2 1 3 1 1 1 1 1 1 3 2 2 1 2 2 1 1 2 1 2 1 3 2 1 3
##  [1009] 2 2 3 2 2 2 1 1 1 1 2 2 1 2 3 2 1 3 1 3 1 1 1 2 1 2 1 2 1 3 2 2 2 1 2 1
##  [1045] 1 2 1 2 1 2 3 2 1 2 1 2 1 1 2 2 2 2 1 2 2 3 1 1 2 3 1 2 2 2 3 1 2 1 2 2
##  [1081] 2 2 3 1 2 2 2 1 1 1 3 3 2 1 1 3 2 2 3 2 1 1 1 2 1 2 1 3 1 2 1 1 2 1 2 2
##  [1117] 1 1 1 1 1 1 1 1 2 3 2 1 2 1 3 2 1 2 2 2 2 2 1 2 1 1 1 3 1 2 2 3 1 2 2 1
##  [1153] 1 1 2 2 1 1 2 2 3 1 2 2 2 1 1 1 2 1 1 2 1 1 2 2 2 2 2 1 2 1 1 2 2 1 3 1
##  [1189] 2 2 3 3 1 2 2 2 2 1 2 1 2 2 2 2 2 3 2 2 1 1 2 2 2 2 1 2 1 2 1 2 1 2 3 3
##  [1225] 2 1 1 2 1 3 1 1 2 1 2 1 3 2 2 2 1 1 2 2 2 2 1 2 1 2 2 2 1 2 2 2 2 1 1 2
##  [1261] 2 2 2 2 2 2 1 1 2 1 1 2 2 1 2 1 3 2 3 1 1 1 2 2 2 1 1 3 1 2 1 2 2 1 2 2
##  [1297] 2 3 2 2 1 2 1 2 1 1 2 2 1 2 3 1 2 1 2 2 2 1 3 2 2 2 2 3 1 2 2 1 2 3 1 1
##  [1333] 3 2 1 2 1 1 2 2 2 1 1 1 1 2 2 2 2 2 1 1 1 2 2 1 1 3 2 2 2 3 2 1 2 3 3 1
##  [1369] 1 1 2 2 2 1 1 1 2 2 2 2 2 3 1 2 2 1 2 1 1 2 2 2 2 2 1 2 2 1 2 2 2 3 3 1
##  [1405] 1 1 2 2 1 2 1 2 1 2 2 1 2 2 1 2 2 1 1 2 2 1 2 1 1 2 1 1 1 1 2 2 1 1 1 1
##  [1441] 2 2 1 1 2 1 2 2 3 2 1 1 1 2 1 2 1 2 3 2 2 2 1 1 1 1 1 3 1 2 3 2 1 2 2 1
##  [1477] 1 2 3 2 3 1 2 2 2 1 2 2 2 1 2 1 2 2 1 2 2 2 2 2 1 3 2 1 2 2 2 1 2 1 1 3
##  [1513] 3 2 2 2 2 1 1 2 1 2 2 1 1 1 2 1 3 2 1 1 1 2 2 2 1 1 1 3 2 1 1 2 2 2 1 2
##  [1549] 1 1 2 1 2 1 2 2 2 2 2 2 1 1 1 1 2 3 1 2 2 3 2 3 2 2 2 2 1 2 2 1 1 3 1 3
##  [1585] 1 1 2 2 2 2 2 1 1 1 1 2 1 1 2 1 3 1 3 1 2 1 1 2 2 1 2 1 1 1 2 2 1 3 2 2
##  [1621] 2 1 1 2 1 1 2 1 1 2 2 1 2 2 2 1 3 1 3 1 2 2 1 2 2 2 2 2 3 2 1 2 3 2 1 2
##  [1657] 2 2 1 1 1 1 2 1 2 1 2 2 2 1 2 2 3 2 1 1 2 2 3 2 2 3 2 2 1 1 1 1 1 1 2 1
##  [1693] 2 1 2 1 2 2 1 1 1 1 2 3 1 2 3 1 1 1 2 3 3 2 2 1 3 2 2 1 2 2 2 2 2 2 1 2
##  [1729] 1 1 2 1 1 2 2 1 2 1 2 1 2 1 1 3 2 1 1 1 1 2 2 1 1 2 2 1 3 1 2 1 1 1 2 2
##  [1765] 1 2 1 2 2 1 1 1 1 2 2 2 1 2 1 3 2 2 1 1 1 1 1 2 2 2 2 1 2 2 2 1 2 2 1 3
##  [1801] 2 1 2 1 2 2 1 2 1 2 2 1 2 1 1 1 2 1 2 1 1 2 2 2 1 2 2 2 1 1 2 1 1 2 1 2
##  [1837] 1 1 3 1 2 2 1 3 1 1 2 1 1 1 2 1 2 1 1 1 2 2 1 1 1 1 2 1 1 1 2 2 1 1 1 3
##  [1873] 1 2 1 2 1 2 2 2 2 3 2 1 1 2 1 3 2 3 2 2 1 3 1 1 2 1 1 1 1 1 1 3 3 1 2 1
##  [1909] 1 2 2 2 2 1 1 1 1 2 1 2 2 2 1 1 2 2 2 3 2 2 2 1 2 2 2 2 3 1 1 2 1 2 3 1
##  [1945] 2 2 2 3 2 2 1 1 2 3 3 1 1 2 1 2 1 1 1 1 2 2 2 1 2 1 2 2 2 2 3 2 2 3 3 2
##  [1981] 1 2 1 2 2 2 2 2 2 2 3 1 3 1 1 2 1 1 1 1 2 1 2 2 2 1 1 1 3 1 2 2 2 3 1 1
##  [2017] 2 2 1 1 1 1 2 1 2 3 2 2 2 2 1 3 1 1 3 2 1 1 3 2 1 3 2 2 2 1 1 2 2 1 3 2
##  [2053] 1 3 1 1 1 1 2 1 1 3 1 1 2 1 3 2 2 2 2 2 2 1 2 2 3 1 2 2 2 1 2 2 2 2 2 2
##  [2089] 1 2 3 2 2 1 3 2 2 1 1 1 1 1 1 3 1 1 2 2 1 1 2 2 2 2 2 1 2 2 2 2 2 1 2 2
##  [2125] 1 1 2 1 1 1 2 2 1 2 3 1 1 2 3 1 1 3 3 3 1 2 1 2 3 2 1 1 1 1 2 3 1 2 3 2
##  [2161] 2 2 2 2 1 2 1 3 2 3 2 1 3 2 1 1 1 1 2 2 1 3 3 1 1 1 1 3 1 1 1 2 1 1 3 2
##  [2197] 1 2 1 2 1 1 1 3 1 1 2 1 1 3 1 3 2 1 3 1 1 2 2 1 1 2 2 2 2 2 2 1 1 2 1 2
##  [2233] 2 1 2 1 2 1 1 1 2 2 1 2 1 2 2 1 2 1 1 3 3 2 2 3 1 2 2 3 2 2 3 2 1 1 1 2
##  [2269] 1 2 2 1 2 1 1 2 1 3 1 2 1 1 2 3 2 1 1 3 1 1 1 3 3 3 1 1 1 2 3 1 1 1 2 2
##  [2305] 1 1 1 2 2 2 1 3 1 2 2 2 2 2 1 1 1 2 2 1 1 3 1 2 2 1 2 2 3 1 1 2 3 1 1 2
##  [2341] 2 1 1 1 2 1 1 1 1 2 2 1 2 2 1 2 1 3 3 2 1 2 1 2 1 3 2 1 1 3 3 2 1 1 3 1
##  [2377] 3 1 2 1 2 1 1 1 1 2 1 2 2 3 2 1 1 2 1 2 1 1 2 1 2 1 1 2 3 1 1 2 1 2 1 2
##  [2413] 2 2 1 3 2 2 1 2 1 2 2 2 2 2 3 2 1 3 1 1 2 1 2 2 1 2 1 1 1 2 2 1 2 1 1 2
##  [2449] 1 2 1 2 1 2 1 1 2 2 2 1 2 2 2 1 1 1 2 2 2 2 2 2 1 2 2 2 2 2 2 2 2 1 1 2
##  [2485] 2 3 2 1 1 2 2 2 1 2 2 1 1 1 1 2 1 2 3 1 2 2 2 2 2 2 1 1 2 1 1 2 1 1 2 1
##  [2521] 2 1 1 1 1 2 3 1 2 3 2 2 2 1 3 1 1 2 1 2 1 1 2 2 2 2 2 2 1 2 1 2 2 1 3 2
##  [2557] 3 2 2 3 1 2 2 1 2 3 2 1 2 2 1 1 2 1 2 2 2 2 1 2 1 3 1 1 1 1 2 1 2 3 1 1
##  [2593] 3 2 2 2 2 1 2 1 3 3 2 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 1 1 2 3 1 2 3
##  [2629] 2 2 3 2 1 2 2 1 1 2 2 2 2 3 2 1 2 1 2 1 2 1 2 2 2 3 3 1 1 2 1 1 1 2 1 1
##  [2665] 1 1 2 2 2 1 1 3 1 2 3 2 2 2 2 2 2 1 2 3 2 3 2 3 1 3 2 1 2 1 3 1 2 1 2 2
##  [2701] 2 1 2 2 2 2 2 1 2 1 2 2 1 2 3 1 3 3 2 1 1 2 1 1 1 2 1 3 1 1 1 3 1 3 2 1
##  [2737] 1 1 1 1 2 3 3 1 1 2 2 1 1 1 1 1 1 2 2 2 2 1 1 1 1 2 2 1 1 2 2 3 2 1 2 3
##  [2773] 2 1 1 2 2 1 2 2 1 1 1 1 2 2 1 1 2 1 2 2 2 1 3 1 2 2 2 2 2 2 2 1 1 2 2 3
##  [2809] 2 1 2 2 1 2 1 2 2 2 2 3 1 2 2 1 2 3 1 2 2 3 2 1 2 2 1 1 2 1 1 3 2 1 2 2
##  [2845] 1 1 1 1 1 2 2 1 2 2 3 1 2 2 2 1 1 1 2 2 1 2 1 2 1 1 2 1 2 2 1 2 1 1 2 1
##  [2881] 1 1 2 2 1 2 1 2 1 2 1 2 1 1 3 1 2 2 1 2 1 2 2 1 1 2 1 3 3 1 3 3 1 1 3 2
##  [2917] 1 3 2 3 2 2 1 2 1 1 2 1 2 1 1 1 2 2 1 1 2 2 2 2 1 1 1 1 1 2 1 1 2 1 1 1
##  [2953] 1 3 1 2 1 1 3 1 3 2 1 2 1 1 2 2 1 1 2 1 2 1 2 3 2 3 3 2 2 2 2 2 2 2 2 2
##  [2989] 1 1 1 1 1 3 1 2 2 1 2 1 2 1 2 3 1 2 2 2 2 1 1 2 3 1 1 1 2 2 2 2 2 1 2 2
##  [3025] 1 2 1 2 2 2 1 2 1 1 1 1 1 2 2 1 2 1 2 1 1 2 1 3 2 1 1 1 1 1 3 2 3 1 3 2
##  [3061] 1 2 1 3 1 3 2 1 1 2 3 1 2 1 2 2 2 1 1 1 1 2 2 2 1 3 1 1 2 1 3 1 2 1 2 2
##  [3097] 3 2 3 3 1 1 1 1 3 2 1 1 1 1 1 2 1 2 2 2 3 1 2 2 2 3 1 2 1 2 2 1 1 2 2 2
##  [3133] 1 3 1 1 3 2 1 2 2 1 1 1 2 2 3 2 2 1 1 1 2 1 3 1 2 1 3 1 2 2 2 2 2 1 1 1
##  [3169] 2 3 1 3 1 2 1 3 1 1 1 2 3 2 2 1 1 2 1 1 3 1 1 2 2 3 1 2 1 1 1 2 1 3 1 1
##  [3205] 1 1 2 2 2 1 2 1 1 2 2 2 1 1 2 2 2 1 2 2 1 1 1 2 1 3 1 1 2 2 2 1 1 2 2 3
##  [3241] 2 1 1 1 2 1 3 3 2 3 2 1 1 2 1 2 1 1 2 2 2 1 1 2 2 3 1 1 2 2 2 1 1 2 1 2
##  [3277] 2 1 1 1 1 1 1 3 1 2 2 2 2 1 1 1 1 2 2 1 2 1 3 1 2 2 2 2 2 1 2 3 2 1 1 3
##  [3313] 3 2 1 2 1 1 1 1 3 2 1 2 2 1 1 3 2 1 1 1 1 2 2 3 1 1 2 1 1 3 1 3 2 2 2 2
##  [3349] 2 1 2 2 2 1 1 1 2 2 1 3 3 1 2 2 2 2 2 1 2 2 1 2 3 1 1 1 1 2 2 2 2 2 2 2
##  [3385] 2 2 2 2 1 2 2 1 2 2 2 2 1 1 2 2 2 2 1 2 2 2 1 3 2 1 1 1 1 1 1 3 2 2 3 1
##  [3421] 1 1 3 1 3 3 2 2 2 2 2 1 1 1 1 2 2 2 1 1 2 3 1 3 1 2 2 1 2 1 2 1 2 2 1 2
##  [3457] 2 1 1 1 1 1 2 3 3 2 1 3 2 2 2 1 2 1 1 1 3 2 2 1 2 2 1 2 1 2 2 1 2 1 2 1
##  [3493] 1 2 2 1 1 1 1 2 2 1 1 2 2 1 1 1 1 3 2 2 2 1 2 1 2 2 1 3 1 2 1 2 3 1 1 2
##  [3529] 1 1 2 2 2 1 2 1 1 2 1 1 2 2 1 1 2 2 1 1 3 1 1 1 2 1 3 1 1 1 2 2 2 1 1 2
##  [3565] 2 1 1 2 2 2 1 1 3 1 1 2 1 2 2 2 1 1 2 2 2 1 2 1 2 2 2 3 1 2 1 2 1 1 1 1
##  [3601] 2 3 2 1 2 2 1 2 2 1 2 2 1 2 2 3 2 1 1 1 2 2 1 2 1 2 1 1 2 1 1 2 2 3 1 1
##  [3637] 2 3 2 2 2 1 2 3 2 1 1 2 1 2 2 2 2 1 1 2 1 1 3 2 2 3 2 2 1 1 1 2 1 2 2 2
##  [3673] 1 1 1 1 1 1 2 2 1 1 2 2 1 1 2 2 1 2 2 1 1 1 3 2 2 2 2 2 3 1 2 2 2 2 1 2
##  [3709] 1 2 1 3 1 2 2 1 2 1 2 2 2 2 2 3 1 1 1 1 3 1 1 1 3 2 2 2 1 1 1 1 2 1 1 1
##  [3745] 2 1 1 1 2 1 2 1 1 2 2 2 2 1 1 1 2 2 2 1 2 3 2 1 1 3 2 3 1 1 2 2 2 1 1 1
##  [3781] 1 1 1 1 2 3 1 2 1 3 2 3 1 2 1 2 3 2 2 2 1 2 2 3 2 1 1 1 1 1 2 3 2 1 3 2
##  [3817] 3 1 2 1 2 1 2 2 2 1 2 1 2 1 1 2 3 2 1 1 1 1 2 2 1 2 1 1 1 2 1 1 1 2 2 2
##  [3853] 1 3 3 1 2 2 1 1 2 3 2 2 1 2 2 3 1 1 2 2 1 1 1 1 1 2 2 2 1 3 1 2 3 1 2 2
##  [3889] 1 1 2 2 1 1 2 3 3 3 2 1 2 2 1 1 1 2 3 2 2 3 1 3 1 2 1 1 1 2 2 2 1 2 2 3
##  [3925] 2 1 2 2 2 2 1 1 1 2 1 2 1 2 1 2 2 2 1 2 1 2 1 1 3 1 1 2 1 1 1 1 1 1 1 1
##  [3961] 3 2 2 2 2 1 2 3 2 2 2 1 1 1 2 3 2 2 2 2 2 1 1 1 2 1 2 2 2 1 1 2 1 1 2 1
##  [3997] 2 2 2 2 3 2 2 1 1 1 2 1 1 2 3 1 2 2 1 1 2 1 2 2 2 2 1 2 2 1 1 3 2 1 2 1
##  [4033] 2 1 2 2 3 1 2 2 1 2 1 1 2 1 2 1 1 1 3 2 1 2 2 1 2 2 1 1 1 1 1 2 2 1 1 2
##  [4069] 1 1 2 2 1 2 2 2 1 3 2 1 3 2 2 2 1 2 2 3 1 2 2 2 1 2 2 1 1 2 1 2 3 2 2 1
##  [4105] 2 2 2 1 2 2 1 2 2 1 1 1 1 2 1 2 2 2 3 1 1 1 1 1 2 1 1 2 2 1 1 1 1 1 2 1
##  [4141] 2 1 3 1 2 1 1 3 2 1 1 1 1 2 2 2 2 2 2 3 1 3 3 1 2 2 3 1 2 1 1 1 2 1 2 1
##  [4177] 2 1 1 2 2 3 2 2 2 3 1 1 2 2 2 1 1 2 2 1 2 2 3 1 1 2 2 1 2 2 1 1 2 2 2 3
##  [4213] 1 2 2 3 2 1 2 2 1 2 1 1 1 2 2 1 1 1 1 2 1 1 2 1 2 1 1 1 2 2 1 3 2 2 2 1
##  [4249] 2 1 3 2 1 1 1 3 2 2 2 1 2 1 1 2 2 1 1 1 1 1 1 1 2 3 1 2 1 3 2 1 2 2 1 1
##  [4285] 2 1 2 2 2 1 1 1 2 2 2 1 2 2 1 2 1 2 3 2 1 1 2 1 2 1 1 1 2 2 3 2 2 1 2 2
##  [4321] 1 2 2 3 2 2 1 1 2 2 1 1 2 1 2 1 1 1 2 2 1 2 1 2 1 2 1 1 2 3 2 1 2 2 2 1
##  [4357] 2 3 2 1 2 1 3 1 3 2 2 1 1 2 1 1 1 2 1 2 2 1 1 1 1 1 1 2 2 2 1 1 2 1 2 2
##  [4393] 2 1 2 2 3 2 2 1 2 1 3 1 1 1 2 2 2 3 1 2 2 1 2 1 1 2 2 2 1 2 2 1 3 2 2 3
##  [4429] 2 1 2 1 1 1 2 2 2 2 1 1 2 2 2 2 1 3 2 2 1 2 1 2 3 1 2 1 1 3 3 2 2 2 1 3
##  [4465] 2 2 1 2 1 2 2 1 2 1 2 1 1 1 3 2 1 1 2 1 1 1 1 2 2 1 2 2 1 1 3 2 2 1 1 2
##  [4501] 2 2 2 2 1 2 1 2 2 1 1 1 2 2 1 3 1 3 2 2 1 1 2 2 1 1 1 1 2 2 1 1 1 1 3 1
##  [4537] 1 2 1 2 2 1 2 1 2 3 3 2 2 1 1 1 3 2 2 2 3 1 1 2 1 2 2 1 2 2 1 2 1 1 2 1
##  [4573] 1 2 2 2 3 1 2 1 2 3 2 1 1 2 2 2 1 1 2 2 1 1 3 1 2 2 1 1 2 2 2 2 1 1 1 1
##  [4609] 1 2 1 2 2 1 3 1 2 1 2 1 2 1 3 2 1 1 2 1 3 1 1 2 3 2 3 2 2 2 1 2 1 2 1 2
##  [4645] 1 1 2 2 2 1 1 2 2 2 1 3 1 1 2 2 2 1 1 1 1 1 3 3 3 1 2 1 2 1 1 1 1 3 2 2
##  [4681] 2 2 1 1 1 2 1 3 1 1 1 1 2 3 1 1 1 2 3 1 3 2 2 1 1 2 3 1 1 2 1 1 2 1 2 2
##  [4717] 2 2 1 1 2 2 3 1 1 1 2 2 2 1 3 2 2 2 3 1 1 1 1 1 1 2 3 2 2 1 1 1 1 2 1 1
##  [4753] 1 2 2 2 2 3 3 2 3 1 2 2 1 2 2 3 1 2 2 1 2 2 1 1 2 2 3 2 2 2 2 2 2 2 2 1
##  [4789] 3 1 1 1 1 2 1 3 2 1 1 1 3 2 2 1 3 1 1 2 2 1 1 2 1 2 3 2 2 3 2 3 1 2 3 1
##  [4825] 1 1 2 1 2 3 2 3 3 1 2 2 1 1 3 2 2 3 2 1 2 3 2 1 2 2 2 3 2 3 2 2 1 1 3 1
##  [4861] 2 1 2 1 1 1 2 2 3 3 2 1 1 1 3 1 1 2 1 1 2 2 2 1 1 2 2 2 2 2 2 1 2 2 2 1
##  [4897] 1 1 2 1 2 1 2 2 1 1 2 2 1 2 1 2 2 2 1 1 2 1 3 3 1 2 2 1 1 1 2 3 1 1 1 1
##  [4933] 1 1 2 2 1 2 2 2 1 1 2 1 2 3 3 1 2 2 1 2 1 3 1 1 1 2 3 2 1 1 1 2 2 2 2 2
##  [4969] 1 1 1 2 2 1 2 2 2 2 1 1 2 1 3 3 2 2 3 1 2 1 1 2 1 2 1 1 1 2 2 2 1 1 2 2
##  [5005] 1 2 1 3 2 2 2 1 1 1 2 2 2 2 1 2 2 2 2 2 2 1 1 2 1 2 1 1 2 1 1 2 2 3 2 2
##  [5041] 2 2 2 2 2 1 2 1 1 1 1 2 1 1 1 1 3 1 1 1 2 2 1 1 1 1 2 1 1 1 2 1 2 2 2 1
##  [5077] 3 3 1 2 2 1 1 3 1 1 1 2 2 1 2 1 2 1 1 2 2 2 3 2 2 1 2 2 1 2 1 1 2 2 2 2
##  [5113] 1 2 3 3 1 3 1 2 2 2 2 1 2 2 1 2 2 2 1 1 1 1 1 2 2 3 1 2 2 2 2 1 1 3 1 1
##  [5149] 1 1 1 3 2 1 1 1 1 1 1 1 1 2 1 1 1 2 1 2 2 1 3 3 2 2 1 2 2 1 2 2 1 1 1 2
##  [5185] 1 3 2 1 1 1 1 1 1 1 2 2 2 2 1 2 1 1 3 1 1 2 1 2 1 3 2 2 2 3 2 1 2 1 2 1
##  [5221] 3 3 2 3 1 1 3 1 1 2 2 3 3 1 1 3 1 1 3 2 1 1 2 2 3 2 2 2 2 2 3 2 1 2 2 3
##  [5257] 1 2 1 1 1 2 2 2 2 2 3 1 2 2 2 1 1 2 3 3 2 2 2 1 2 2 2 2 2 2 2 1 2 1 3 1
##  [5293] 2 2 3 1 2 1 2 2 1 2 2 1 2 2 1 3 1 1 1 2 1 2 1 1 2 1 2 2 2 1 1 1 2 3 2 1
##  [5329] 1 2 2 1 2 1 1 2 1 1 1 2 3 2 1 3 1 2 2 2 2 2 1 2 2 1 2 3 1 1 1 2 1 2 1 1
##  [5365] 2 3 3 1 1 1 2 2 1 1 1 1 1 1 1 3 2 3 2 1 3 2 3 1 1 1 2 1 3 1 1 1 2 1 2 2
##  [5401] 2 3 2 1 1 2 2 3 2 2 2 3 1 1 2 2 2 1 1 1 2 2 2 3 1 1 1 1 2 2 3 2 1 2 1 2
##  [5437] 1 1 1 1 1 2 2 2 2 1 1 2 1 1 1 1 1 1 2 2 3 1 1 2 2 2 2 2 2 1 1 1 1 3 2 1
##  [5473] 2 1 1 2 1 2 1 1 2 1 2 2 1 3 2 2 2 1 2 2 1 2 2 1 2 2 1 1 2 1 2 2 1 2 1 1
##  [5509] 2 1 1 1 1 1 1 2 1 1 2 1 1 3 2 2 3 3 3 1 2 1 1 2 3 1 1 2 1 1 2 2 2 3 2 1
##  [5545] 1 2 1 1 1 2 3 1 2 2 2 3 1 1 1 1 2 2 2 1 2 2 2 2 2 3 1 1 1 3 2 1 1 2 1 2
##  [5581] 3 1 1 2 2 2 1 2 2 1 2 2 2 1 1 1 1 1 1 1 1 1 3 1 1 2 1 1 3 1 2 2 2 2 1 1
##  [5617] 2 1 2 1 1 2 1 1 2 2 2 1 2 1 1 1 2 3 3 2 2 2 1 2 2 2 2 1 2 1 2 1 2 1 2 1
##  [5653] 3 2 2 1 1 1 1 2 1 1 2 1 1 2 2 2 1 2 1 2 2 1 2 1 1 3 3 1 1 1 2 1 1 2 1 1
##  [5689] 2 1 3 2 1 2 3 1 1 2 3 1 2 1 2 1 2 3 1 1 1 2 2 3 2 1 2 2 2 2 2 2 1 2 3 1
##  [5725] 3 2 2 1 2 1 2 2 3 1 1 2 2 2 2 3 1 1 1 2 2 1 1 1 2 2 2 2 3 2 2 1 1 2 2 2
##  [5761] 2 3 1 2 1 3 1 3 1 1 3 1 2 2 1 1 1 1 1 3 2 2 1 2 1 1 1 2 1 2 2 2 2 1 2 2
##  [5797] 2 2 2 1 2 1 2 2 1 2 2 1 2 2 1 1 2 2 1 1 1 1 3 1 1 2 2 1 1 1 3 1 2 1 1 2
##  [5833] 3 1 2 2 1 2 1 1 2 3 2 2 1 2 2 1 3 2 1 1 2 2 1 1 1 2 2 2 3 1 2 2 1 2 2 1
##  [5869] 1 2 1 3 1 2 1 1 2 1 3 1 1 2 2 2 2 2 2 1 3 3 1 1 1 2 1 1 3 3 2 2 1 2 1 1
##  [5905] 2 1 1 2 1 2 2 2 2 1 1 2 1 1 3 1 2 2 2 2 3 2 1 1 2 2 1 1 2 2 1 2 2 1 3 2
##  [5941] 3 3 2 1 3 2 1 1 1 2 1 2 3 2 1 2 2 2 3 1 1 2 3 2 2 1 2 2 2 2 2 2 2 2 2 2
##  [5977] 3 2 2 3 2 1 1 2 2 1 2 2 2 1 1 2 1 3 3 1 1 2 2 2 1 2 2 2 2 3 1 2 1 1 2 1
##  [6013] 2 2 2 1 1 2 1 2 2 2 2 2 2 1 1 2 3 1 1 1 2 2 2 1 1 2 2 2 2 1 3 2 2 2 3 1
##  [6049] 1 3 2 2 1 2 1 1 1 2 2 2 2 2 2 1 2 1 1 2 1 3 1 1 1 1 1 1 2 2 3 2 1 2 1 1
##  [6085] 2 1 2 2 1 2 1 1 2 2 2 1 1 2 1 3 1 2 1 2 1 2 1 1 2 1 2 1 1 2 2 2 1 1 1 3
##  [6121] 2 1 2 1 1 2 3 1 2 2 1 2 3 2 2 1 1 1 1 2 3 1 2 3 2 1 1 2 2 2 1 1 2 2 2 2
##  [6157] 3 1 1 3 2 2 1 1 3 2 2 2 1 2 2 2 2 1 2 1 1 1 1 3 2 2 1 2 1 1 3 1 1 2 1 2
##  [6193] 2 1 2 2 1 1 2 1 1 1 2 2 1 3 2 2 2 2 2 3 2 1 2 2 1 3 2 2 2 1 3 2 3 2 3 2
##  [6229] 1 1 1 1 1 1 2 1 2 3 2 1 2 2 1 3 2 1 3 1 1 1 1 1 2 3 2 1 1 2 2 3 1 2 1 3
##  [6265] 2 2 1 1 2 1 2 1 1 1 2 2 1 3 2 2 3 1 3 2 1 2 1 2 1 1 1 1 1 1 1 1 2 2 2 2
##  [6301] 2 2 2 1 2 3 2 2 1 2 1 2 2 2 2 1 3 2 1 2 2 1 1 1 1 2 1 2 1 2 2 2 2 1 1 1
##  [6337] 1 2 2 3 2 1 1 1 1 2 1 1 1 2 1 2 1 1 2 1 3 1 1 2 1 2 2 2 1 1 1 2 1 3 1 3
##  [6373] 1 2 2 1 2 2 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 1 2 1 2 1 1 2 2 2 1 1 2 3 1
##  [6409] 2 2 2 1 1 2 1 1 1 2 2 1 2 1 2 1 1 1 1 2 1 2 1 2 2 1 3 3 1 1 1 2 1 2 1 1
##  [6445] 1 2 2 2 3 1 3 1 2 3 2 3 1 1 2 1 3 2 2 2 1 2 2 2 2 1 1 1 2 1 1 1 1 2 1 1
##  [6481] 3 1 3 1 1 3 1 1 1 1 1 1 1 2 1 1 2 2 2 2 2 1 1 2 1 2 2 1 3 2 2 1 2 2 1 1
##  [6517] 2 2 3 2 1 2 1 1 1 1 1 2 1 1 2 2 2 1 1 2 1 1 1 1 1 1 2 2 2 2 2 1 2 2 2 1
##  [6553] 2 3 1 1 1 1 1 3 1 2 3 2 1 1 2 2 2 1 1 2 1 1 2 1 1 2 3 1 2 1 2 1 2 1 1 1
##  [6589] 2 2 1 2 1 1 2 1 1 3 2 1 1 2 2 1 1 2 2 3 2 3 2 2 1 2 2 1 1 3 2 3 3 1 1 1
##  [6625] 2 1 1 2 2 2 2 2 2 2 3 1 2 2 1 1 1 1 1 2 2 2 1 2 1 2 1 2 3 3 2 2 2 2 3 3
##  [6661] 2 3 2 1 1 2 2 2 2 2 1 1 2 1 1 1 1 1 2 2 1 2 3 2 3 2 1 1 1 1 1 1 3 2 2 1
##  [6697] 1 1 2 2 1 1 2 1 1 2 2 2 1 3 2 1 1 1 1 3 1 2 1 1 1 1 1 2 1 1 2 2 1 2 1 2
##  [6733] 2 2 2 3 2 3 2 2 2 2 2 2 1 2 1 2 2 2 2 2 1 2 1 1 1 2 2 2 2 1 1 1 1 2 3 2
##  [6769] 1 2 2 1 2 1 1 2 2 1 1 1 1 2 2 1 2 1 2 1 2 2 2 2 1 3 2 3 1 2 2 1 2 1 1 1
##  [6805] 1 1 2 3 2 1 1 1 1 1 2 2 1 3 3 2 2 1 2 1 1 1 1 1 2 3 1 1 2 2 2 3 3 1 2 1
##  [6841] 1 1 3 2 1 1 1 2 2 1 2 1 1 2 1 1 1 3 2 1 2 1 1 1 1 2 3 1 1 3 1 2 1 1 1 2
##  [6877] 2 2 2 2 2 2 2 1 2 2 1 1 3 1 1 1 2 3 1 1 2 1 1 1 3 3 3 1 1 1 3 1 1 2 1 1
##  [6913] 1 2 1 2 2 2 1 2 2 2 1 2 3 1 2 1 2 2 1 2 2 1 3 2 3 1 2 2 1 2 2 1 1 2 1 2
##  [6949] 2 2 1 2 2 1 1 2 2 2 1 2 3 1 3 2 1 2 1 2 3 2 1 2 2 1 1 1 1 3 1 2 1 1 2 1
##  [6985] 1 1 1 2 1 2 1 1 1 1 1 2 2 2 1 2 1 1 1 2 2 2 1 2 1 1 3 2 1 1 1 2 2 2 2 2
##  [7021] 1 1 2 1 1 1 2 2 1 3 1 1 2 1 1 1 1 1 2 3 1 1 2 1 1 1 1 1 2 2 1 2 1 2 1 2
##  [7057] 1 2 2 1 1 2 3 2 1 1 2 3 2 2 2 2 3 1 1 1 1 1 2 1 1 1 2 1 1 2 2 2 2 1 2 2
##  [7093] 2 2 1 1 2 1 2 2 1 2 1 1 2 1 1 1 2 1 2 2 2 2 1 2 1 2 1 2 2 1 1 1 2 1 3 1
##  [7129] 1 2 1 2 2 2 1 2 3 2 1 1 2 2 1 1 2 1 1 3 1 1 2 2 2 2 2 2 1 2 1 2 1 1 1 1
##  [7165] 2 3 1 1 2 1 2 1 1 1 1 1 1 1 1 2 3 1 2 2 2 1 2 1 2 2 2 1 3 2 2 1 2 1 2 1
##  [7201] 1 3 1 1 2 1 1 1 1 1 1 3 2 1 2 1 2 2 2 2 1 1 1 2 2 1 2 2 2 2 2 3 2 2 3 2
##  [7237] 2 2 2 1 2 1 3 1 3 2 2 2 1 2 2 3 2 3 2 2 2 1 1 1 2 3 2 3 3 1 1 1 2 1 2 1
##  [7273] 2 2 1 2 2 3 2 2 3 2 1 3 2 2 2 2 1 2 1 1 2 2 1 2 2 1 1 1 3 1 1 1 2 2 1 1
##  [7309] 2 1 2 1 2 2 2 1 1 2 1 2 1 2 2 3 2 2 2 3 2 1 2 1 2 1 1 1 2 1 1 1 1 1 1 1
##  [7345] 2 3 2 1 2 1 3 2 1 1 2 1 2 1 1 2 1 1 2 1 2 3 2 2 2 1 1 2 3 1 1 2 1 2 2 2
##  [7381] 1 1 2 1 1 2 2 2 2 1 1 1 1 1 1 1 2 2 1 2 1 1 2 1 1 3 2 1 1 2 1 2 1 1 2 1
##  [7417] 1 2 1 3 2 1 1 3 1 2 2 1 1 1 2 1 2 1 2 1 1 3 1 2 2 2 2 3 2 1 1 1 2 2 1 2
##  [7453] 2 1 1 3 1 2 2 1 2 3 2 1 2 1 2 1 1 2 3 2 1 2 1 2 1 1 1 2 1 1 1 2 2 1 1 1
##  [7489] 1 1 2 2 2 1 2 1 2 2 1 3 2 1 1 1 1 2 1 3 1 2 1 2 1 1 1 2 2 1 1 2 1 1 3 1
##  [7525] 1 2 2 3 3 1 1 1 2 2 3 2 2 1 2 3 2 1 2 2 2 2 1 1 1 1 2 1 1 1 2 2 2 3 1 2
##  [7561] 2 1 1 2 1 1 1 2 2 1 2 2 3 3 1 3 2 2 1 2 2 1 2 1 1 2 1 2 2 2 2 3 2 2 2 2
##  [7597] 2 1 1 2 2 2 2 2 1 2 2 1 1 1 1 3 2 1 2 1 1 1 1 3 2 2 2 2 2 1 1 2 1 2 1 2
##  [7633] 1 2 1 1 3 2 2 2 3 2 2 3 2 1 1 2 1 2 2 2 2 3 2 1 1 2 2 1 2 3 3 2 1 1 1 2
##  [7669] 2 1 1 2 1 2 3 2 1 1 1 1 1 1 1 2 1 2 1 2 1 3 2 3 2 2 2 3 1 3 3 2 2 3 1 2
##  [7705] 1 1 2 3 2 2 2 1 1 2 3 1 2 1 2 2 2 1 1 2 2 2 1 1 2 2 2 3 1 1 3 1 3 1 2 2
##  [7741] 1 2 2 2 1 1 3 2 1 1 2 2 2 1 1 2 2 1 1 1 2 1 2 3 1 3 2 2 2 3 2 1 1 1 1 3
##  [7777] 1 2 2 1 1 1 1 1 2 2 1 1 1 1 1 2 1 2 2 1 3 1 1 3 2 2 2 1 2 3 2 2 1 1 2 1
##  [7813] 2 2 1 2 1 1 3 1 1 1 2 1 2 2 3 1 2 1 3 2 1 2 2 2 1 2 2 2 3 2 2 2 2 3 1 1
##  [7849] 2 2 1 1 2 1 1 1 2 1 1 1 2 3 1 2 2 2 1 2 1 2 2 2 2 1 1 2 1 1 2 1 2 3 1 3
##  [7885] 2 1 2 2 1 2 1 2 1 2 3 1 3 2 1 3 1 2 1 3 1 2 2 1 1 2 1 1 2 1 1 2 2 2 2 1
##  [7921] 1 2 2 3 1 2 2 1 1 1 2 2 2 1 3 2 2 1 1 3 2 2 2 1 1 3 2 1 3 2 2 1 1 2 2 3
##  [7957] 2 2 2 2 2 1 1 1 2 1 1 2 1 2 1 1 2 2 2 2 1 1 2 1 2 2 2 1 2 2 2 2 1 1 1 1
##  [7993] 1 1 1 1 1 3 1 2 3 1 2 1 1 1 1 1 1 1 3 1 1 1 1 1 1 1 1 1 1 2 1 2 1 1 1 1
##  [8029] 1 2 2 3 1 1 2 2 1 3 1 2 2 2 2 2 3 2 2 1 2 2 1 3 2 1 1 2 2 2 3 2 1 2 2 1
##  [8065] 2 3 3 1 2 3 1 3 3 2 1 2 1 1 3 2 2 1 2 1 2 2 2 1 2 2 1 2 2 2 1 1 3 1 3 3
##  [8101] 2 1 2 1 2 1 1 1 2 1 2 2 2 2 1 1 2 2 3 2 1 1 3 2 1 1 2 2 2 1 2 1 3 2 1 3
##  [8137] 2 2 1 2 1 2 2 1 1 2 2 1 1 1 1 2 2 2 1 2 1 1 2 1 2 3 2 3 1 2 1 1 2 1 1 2
##  [8173] 1 2 1 1 2 2 2 2 1 2 1 1 1 2 2 2 2 2 2 3 2 2 1 2 1 2 2 2 2 3 3 1 2 1 2 2
##  [8209] 1 1 2 2 2 2 2 1 2 1 3 1 3 1 1 1 3 1 2 1 3 1 1 1 2 2 1 1 2 3 1 2 2 1 2 2
##  [8245] 1 1 1 3 3 1 1 1 2 2 2 3 1 1 1 3 1 2 3 2 2 2 2 2 2 2 1 2 2 1 3 3 2 1 2 1
##  [8281] 2 2 1 2 2 1 3 1 2 1 1 1 3 2 3 1 1 1 2 3 2 1 2 3 1 2 1 2 1 2 2 1 2 2 2 1
##  [8317] 2 3 3 1 2 1 2 1 2 2 2 1 1 2 2 3 1 2 1 2 2 1 2 1 2 3 2 1 2 1 1 3 1 1 1 1
##  [8353] 2 2 3 1 2 2 1 1 3 2 2 2 1 2 2 2 1 1 1 1 1 2 1 2 3 1 3 2 1 3 1 2 2 2 1 2
##  [8389] 3 3 2 1 2 2 2 2 1 1 2 1 1 1 1 2 1 2 1 1 1 2 2 2 1 2 1 2 1 1 2 2 1 1 1 2
##  [8425] 2 1 1 1 1 2 2 3 2 2 2 1 2 2 2 2 1 1 2 1 1 1 1 2 2 3 2 2 2 2 1 1 2 1 2 2
##  [8461] 2 2 2 2 3 2 1 1 3 1 2 1 2 2 1 2 1 1 3 1 2 2 1 1 1 1 3 1 2 2 2 1 3 1 2 2
##  [8497] 2 1 3 2 1 1 2 2 1 2 2 2 1 3 1 2 2 1 1 2 2 2 2 2 1 1 2 2 1 1 3 1 1 2 3 1
##  [8533] 1 2 2 1 1 2 3 1 2 1 2 1 1 1 1 2 3 1 2 1 2 2 1 2 1 1 1 2 1 1 1 2 2 2 1 1
##  [8569] 2 1 3 1 1 1 1 3 2 1 2 2 1 2 3 2 2 2 1 1 3 1 1 1 2 3 2 2 1 3 1 2 2 3 2 1
##  [8605] 2 1 1 3 2 1 1 1 2 2 1 1 2 3 2 1 1 2 2 3 1 1 1 2 1 1 2 2 3 1 2 2 2 1 1 1
##  [8641] 1 2 1 2 2 1 1 2 2 1 1 1 2 1 2 1 2 2 1 1 2 3 1 2 2 3 2 1 2 3 1 1 2 1 2 1
##  [8677] 2 1 2 2 1 1 2 2 2 2 2 2 2 1 2 3 2 1 2 3 1 2 2 2 1 1 1 2 1 1 1 1 1 2 1 2
##  [8713] 1 1 3 1 2 2 2 1 2 2 1 3 1 3 1 2 1 1 1 1 2 2 1 3 1 2 2 2 3 2 2 1 2 1 1 2
##  [8749] 1 3 1 1 1 1 2 2 1 1 1 3 1 2 2 1 1 1 1 2 1 2 1 3 2 3 2 2 1 2 1 1 1 3 3 3
##  [8785] 2 1 1 1 2 2 2 2 1 1 2 2 2 1 3 1 2 2 2 2 2 2 2 2 2 1 1 1 3 1 2 2 1 2 2 1
##  [8821] 2 2 1 2 1 3 2 1 1 2 2 1 1 1 2 2 2 1 2 1 2 1 2 2 1 2 2 1 2 2 2 2 1 1 2 1
##  [8857] 3 1 3 2 2 2 1 2 2 1 1 1 2 2 2 2 1 2 2 1 2 1 1 2 2 2 3 1 1 3 2 2 2 1 1 1
##  [8893] 1 1 1 1 3 2 1 2 3 1 1 1 1 1 2 1 2 1 1 1 2 2 2 1 3 2 1 1 1 2 2 2 2 2 1 1
##  [8929] 2 3 1 2 1 2 2 1 1 2 2 2 2 1 1 1 1 1 2 3 1 2 2 2 1 1 3 2 1 1 1 2 2 2 2 2
##  [8965] 2 2 2 2 2 1 2 2 2 2 1 3 2 1 1 1 1 2 2 1 1 2 2 2 2 2 1 3 2 1 1 1 1 3 2 1
##  [9001] 3 2 1 2 1 2 2 1 2 2 1 1 1 1 1 2 2 1 2 2 1 1 1 2 2 2 1 2 1 1 2 1 3 1 1 1
##  [9037] 1 2 1 1 1 1 1 2 2 2 2 1 1 1 1 1 2 2 2 2 2 2 2 1 1 1 3 1 2 2 1 1 1 2 1 1
##  [9073] 2 1 1 2 2 2 1 2 1 2 1 1 1 3 1 1 1 1 2 2 1 1 2 1 1 1 3 2 1 3 2 2 1 3 1 1
##  [9109] 1 2 3 2 1 1 2 1 3 1 3 1 2 3 1 2 2 1 2 2 1 2 3 2 1 1 2 2 1 3 1 2 1 1 3 2
##  [9145] 2 1 1 2 2 2 2 2 1 2 2 2 1 1 2 2 2 1 1 1 2 1 2 2 2 2 2 1 1 1 2 1 1 2 2 1
##  [9181] 3 2 3 1 1 2 2 2 1 1 2 3 1 2 1 1 3 1 2 2 1 3 1 2 1 2 2 2 2 1 3 2 1 2 2 2
##  [9217] 1 1 1 2 2 2 3 2 1 1 2 3 2 2 1 1 1 2 3 2 1 2 2 2 2 2 2 3 2 2 1 1 2 1 2 2
##  [9253] 1 1 2 2 3 1 3 2 3 2 1 1 2 2 2 1 2 2 1 1 1 2 2 1 2 2 2 1 2 1 3 3 2 1 2 3
##  [9289] 1 2 1 2 1 2 1 2 2 3 1 1 1 1 2 2 1 2 2 2 2 1 2 1 3 2 1 1 2 1 1 2 1 1 2 2
##  [9325] 1 2 2 1 1 2 1 1 1 2 2 2 2 1 2 2 1 1 3 3 1 1 2 2 2 2 3 2 1 1 1 1 1 2 1 1
##  [9361] 1 2 2 1 1 2 1 2 1 1 3 2 2 1 2 2 3 2 2 2 2 2 2 2 2 1 1 1 2 1 1 1 1 2 1 1
##  [9397] 1 2 2 2 1 2 2 2 1 2 2 2 1 1 3 1 1 3 2 2 2 2 1 2 2 2 2 3 2 2 2 1 1 1 1 2
##  [9433] 1 2 1 2 2 2 1 1 1 1 1 1 1 1 2 2 2 3 2 2 1 3 2 1 2 2 1 2 3 1 2 2 2 1 1 1
##  [9469] 1 1 1 2 1 3 1 1 2 1 2 2 1 1 2 3 3 2 2 2 3 2 1 2 1 2 1 2 1 2 1 1 2 2 3 2
##  [9505] 2 2 2 1 1 2 1 2 2 1 2 3 1 2 2 2 2 1 3 2 1 2 2 1 1 2 1 2 2 2 2 1 1 1 3 2
##  [9541] 1 2 2 1 2 2 1 2 1 1 2 3 2 2 1 2 2 1 1 2 3 1 2 2 2 3 1 2 1 2 2 2 1 2 2 2
##  [9577] 1 2 2 2 2 2 1 1 3 1 2 2 1 1 2 1 2 2 2 1 2 2 1 2 2 1 2 2 2 2 2 2 2 1 2 2
##  [9613] 2 2 2 1 1 2 2 1 1 1 2 2 3 1 3 2 2 1 1 2 2 1 1 2 2 2 2 2 1 2 2 2 2 2 2 3
##  [9649] 1 3 3 2 1 2 1 3 2 2 2 2 2 2 1 1 2 2 2 1 2 1 2 1 1 2 2 2 1 1 2 3 2 2 3 2
##  [9685] 1 2 2 2 2 1 1 2 2 2 2 1 2 2 1 1 2 1 1 2 1 2 2 3 1 2 1 1 2 2 1 2 1 2 2 2
##  [9721] 2 1 2 1 2 2 1 2 2 1 2 1 1 2 2 2 1 1 2 1 2 2 2 2 1 1 2 1 2 2 1 1 2 1 2 2
##  [9757] 2 2 3 2 2 2 2 2 2 1 2 1 2 3 2 1 1 1 2 2 2 2 1 2 3 1 2 1 2 2 1 2 3 2 1 1
##  [9793] 2 1 2 2 1 2 1 2 2 2 2 2 1 1 1 2 3 2 1 2 1 1 1 2 1 2 1 2 2 2 1 3 3 2 2 1
##  [9829] 1 1 3 1 2 1 2 2 1 1 3 2 2 2 2 1 1 1 2 1 1 1 2 1 1 2 3 2 1 2 2 2 2 2 1 2
##  [9865] 2 2 2 3 2 2 1 1 1 2 1 2 1 1 2 2 2 2 1 2 2 1 1 1 2 1 1 1 2 2 1 2 2 2 1 1
##  [9901] 1 2 2 1 2 1 3 1 1 2 1 1 2 1 2 1 1 1 2 3 2 2 2 1 2 1 3 2 1 1 1 1 1 1 2 2
##  [9937] 1 1 2 2 2 2 2 2 1 1 3 1 1 1 1 2 1 3 2 2 1 3 1 2 1 3 2 2 2 2 1 1 3 2 1 1
##  [9973] 1 2 1 3 1 1 2 1 2 2 2 1 1 2 1 1 1 1 3 2 1 2 1 1 3 2 1 1
## 
## Within cluster sum of squares by cluster:
## [1]  6103.296 12369.985  2895.329
##  (between_SS / total_SS =  62.7 %)
## 
## Available components:
## 
## [1] "cluster"      "centers"      "totss"        "withinss"     "tot.withinss"
## [6] "betweenss"    "size"         "iter"         "ifault"
col.pal <- c('#a6cee3','#b2df8a','#1f78b4')
levels(unC$map) <- data.frame(ID=1:3, code=1:3)
levelplot(unC$map, col.regions = col.pal)

Export Geotiff

writeRaster(unC$map,paste("../output/kmeans_", n.cl, "cl.tif", sep = ""), format="GTiff", overwrite = TRUE)

Boxplots

r <- stack(unC$map, vars)
smp <- as.data.frame(sampleRandom(x = r, size = 10000))
smp[,1] <- as.factor(smp[,1])
names(smp)[names(smp)=="layer"] <- "Cluster"

for (i in 2:ncol(smp)) {
  
  print(ggplot(smp, aes(x = Cluster, y = smp[,i],fill = Cluster)) +
 geom_boxplot() +
 scale_y_continuous(name = names(smp[i])) +
 theme(axis.text.x = element_blank(),
                axis.title.x = element_blank(),
                axis.ticks.x = element_blank()))
  
  }

Figure for publication

AoI <- extent(3450000,4450000,3100000,4300000)
land <- readOGR(dsn = "../input", layer = "land")
## OGR data source with driver: ESRI Shapefile 
## Source: "N:\Prosjekter\311700_MAREANO\311720_Havbunnskart\c_Metodeprosjekter\xxx_Karbonlagring\clustering\input", layer: "land"
## with 8186 features
## It has 1 fields
clmap <- crop(unC$map, AoI)
land <- crop(land, AoI)

p <- ggR(clmap, geom_raster = TRUE, forceCat = TRUE) + 
  geom_polygon(data = land, aes(x = long, y = lat, group = group)) + 
  scale_fill_manual(values = col.pal) +
  theme(axis.title = element_blank()) +
  theme(legend.position = "none" ) + 
  ggpubr::rotate_y_text()
## Regions defined for each Polygons
p1 <- ggplot(smp, aes(x = Cluster, y = Bathymetry, fill = Cluster)) +
  geom_boxplot() +
 scale_fill_manual(values = col.pal) +
 scale_y_continuous(name = "Bathymetry (m)") +
 theme(axis.title.x = element_blank()) +
  theme(legend.position = "none")
p2 <- ggplot(smp, aes(x = Cluster, y = Current.speed, fill = Cluster)) +
  geom_boxplot()+
  scale_fill_manual(values = col.pal) +
  scale_y_continuous(name = expression(M2~current~speed~(m~s^{-1}))) +
 theme(axis.title.x = element_blank()) +
  theme(legend.position = "none")
p3 <- ggplot(smp, aes(x = Cluster, y = OCAR, fill = Cluster)) +
  geom_boxplot()+
  scale_fill_manual(values = col.pal) +
  scale_y_continuous(name = expression(OCAR~(g~m^{-2}~yr^{-1}))) +
 theme(axis.title.x = element_blank()) +
  theme(legend.position = "none")
p4 <- ggplot(smp, aes(x = Cluster, y = OC.density, fill = Cluster)) +
  geom_boxplot()+
  scale_fill_manual(values = col.pal) +
  scale_y_continuous(name = expression(OC~density~(kg~m^{-3}))) +
 theme(axis.title.x = element_blank()) +
  theme(legend.position = "none")
p5 <- ggplot(smp, aes(x = Cluster, y = OPD, fill = Cluster)) +
  geom_boxplot() +
  scale_fill_manual(values = col.pal) +
  scale_y_continuous(name = "OPD (cm)") +
 theme(axis.title.x = element_blank()) +
  theme(legend.position = "none")
p6 <- ggplot(smp, aes(x = Cluster, y = Orb.Velocity, fill = Cluster)) +
  geom_boxplot()+
  scale_fill_manual(values = col.pal) +
  scale_y_continuous(name = expression(Orbital~velocity~(m~s^{-1}))) +
 theme(axis.title.x = element_blank()) +
  theme(legend.position = "none")

lay <- rbind(c(1,1,1,2,3,4),
             c(1,1,1,5,6,7))

tiff("../figs/OC_processing_regions.tif", width = 30, height = 15, units = "cm", res = 500, compression = "lzw")
grid.arrange(p, p1, p2, p3, p5, p6, p4, nrow=2, layout_matrix=lay)
dev.off()
## png 
##   2

Save sessionInfo to file

sessionInfo = sessionInfo()
save(sessionInfo, file = "sessionInfo.Rdata")
rm("sessionInfo")