library(Hmsc)
library(ggplot2)
library(abind)
library(dplyr)
library(ggpubr)
library(patchwork)

setwd("E:/Helsinki/ResultsAug23/For submission/1.Data/1.EnvironmentalInfo")
load('PropFINall_Jan23.RDATA')

head(Prop00)
length(Prop00)

library(purrr)
library(tidyr)
library(dplyr)

Prop00prep = Prop00

Mat00.0.5 = plyr::ldply(Prop00prep, rbind)
#rownames(Mat00.0.5) = 1:43453
Mat00.0.5[is.na(Mat00.0.5)] = 0

Mat00.0.5$'100' = Mat00.0.5$'100' + Mat00.0.5$'133'
Mat00.0.5$'200' = Mat00.0.5$'200' + Mat00.0.5$'222'

Mat00.0.5B = Mat00.0.5[, 1:6]

# only use the small Finland one
seq1 = seq(from=1, to=43453-1, by=3)
Hab00 = Mat00.0.5B#[seq1,]
head(Hab00)
dim(Hab00)

Hab00 = Hab00[,order(colnames(Hab00), decreasing = FALSE)]
sum(apply(Hab00, 2, mean, na.rm=T)*100)



Prop06prep = Prop06_0.5

Mat06.0.5 = plyr::ldply(Prop06prep, rbind)
Mat06.0.5[is.na(Mat06.0.5)] = 0
#rownames(Mat06.0.5) = 1:43453
Mat06.0.5$'100' = Mat06.0.5$'100' + Mat06.0.5$'133'
Mat06.0.5$'200' = Mat06.0.5$'200' + Mat06.0.5$'222'

Mat06.0.5B = Mat06.0.5[, 1:6]

# only use the small Finland one
seq1 = seq(from=1, to=43453-1, by=3)
Hab06 = Mat06.0.5B#[seq1,]
head(Hab06)
dim(Hab06)

Hab06 = Hab06[,order(colnames(Hab06), decreasing = FALSE)]
sum(apply(Hab06, 2, mean, na.rm=T)*100)



Prop12prep = Prop12_0.5

Mat12.0.5 = plyr::ldply(Prop12prep, rbind)
Mat12.0.5[is.na(Mat12.0.5)] = 0

#rownames(Mat12.0.5) = 1:43453
Mat12.0.5$'100' = Mat12.0.5$'100' + Mat12.0.5$'133'
Mat12.0.5$'200' = Mat12.0.5$'200' + Mat12.0.5$'222'

Mat12.0.5B = Mat12.0.5[, 1:6]

# only use the small Finland one
seq1 = seq(from=1, to=43453-1, by=3)
Hab12 = Mat12.0.5B#[seq1,]
head(Hab12)
dim(Hab12)

Hab12 = Hab12[,order(colnames(Hab12), decreasing = FALSE)]
sum(apply(Hab12, 2, mean, na.rm=T)*100)


Prop18prep = Prop18_0.5

Mat18.0.5 = plyr::ldply(Prop18prep, rbind)
Mat18.0.5[is.na(Mat18.0.5)] = 0

#rownames(Mat18.0.5) = 1:43453
Mat18.0.5$'100' = Mat18.0.5$'100' + Mat18.0.5$'133'
Mat18.0.5$'200' = Mat18.0.5$'200' + Mat18.0.5$'222'

Mat18.0.5B = Mat18.0.5[, 1:6]

# only use the small Finland one
seq1 = seq(from=1, to=43453-1, by=3)
Hab18 = Mat18.0.5B#[seq1,]
head(Hab18)
dim(Hab18)

Hab18 = Hab18[,order(colnames(Hab18), decreasing = FALSE)]
sum(apply(Hab18, 2, mean, na.rm=T)*100)


###########################################################################

#### all

HabAll = list(Hab00, Hab06, Hab12, Hab18)

HabAllmean = Reduce('+', HabAll)/length(HabAll)
sum(apply(HabAllmean, 2, mean, na.rm=T)*100)
apply(HabAllmean, 2, mean, na.rm=T)*100


#load('Coord_small_fin.RDATA')
library(sf)
Fin_grid = st_read('3km_gridpoints.shp')
class(Fin_grid)

Grid_sf = as(Fin_grid, "Spatial")
plot(Grid_sf)

## Clip the map
library(maptools)
library(raster)   ## To convert an "Extent" object to a "SpatialPolygons" object.
library(rgeos)

Finmap = map('worldHires','Finland', plot=FALSE, fill = TRUE)

Fin.sp = map2SpatialPolygons(map = Finmap, IDs = Finmap$names, 
                             proj4string = CRS("+init=epsg:4326"), 
                             checkHoles=FALSE)
newMFin.sp <- spTransform(Fin.sp, CRS('+proj=utm +zone=35 +ellps=GRS80 +units=m +no_defs'))

out <- gIntersection(newMFin.sp, Grid_sf, byid=TRUE)
plot(out)
extcrop = extent(62000, 6606750, 7700000, 7773800)
out2 = crop(out, extcrop)
plot(out2)

CoordFinall = as.data.frame(Fin_grid)[, 1:2]

# MAKE A DATAFRAME OF THE DATA TO BE PLOTTED
mapData = data.frame(CoordFinall, HabAllmean)
colnames(mapData)[1:2] = c('X', 'Y')

# get rid of extra popints in the south in the ocean?
library(mapsFinland)
library(maps)
library(mapdata)

par(mfrow=c(1,1), oma= c(2,2,2,2))
Finmap = map('worldHires','Finland', plot=FALSE, fill = TRUE)
class(Finmap)

Finmap$names[2] = "Finland1:?" 
Finmap$names[5] = "Finland2:?" 
Finmap$names[6] = "Finland3:?" 

library(maptools)
Fin.sp = map2SpatialPolygons(map = Finmap, IDs = Finmap$names, 
                             proj4string = CRS("+init=epsg:4326"), 
                             checkHoles=FALSE)
newMFin.sp <- spTransform(Fin.sp, CRS("+proj=utm +zone=35 +ellps=GRS80 +units=m +no_defs"))
class(newMFin.sp)

# prepare coordinates, data, and proj4string
coords <- mapData[ , c("X", "Y")]   # coordinates
data   <- mapData[ , -c(1:2)]          # data
crs    <- CRS("+proj=utm +zone=35 +ellps=GRS80 +units=m +no_defs") # proj4string of coords

# make the SpatialPointsDataFrame object
spdf <- SpatialPointsDataFrame(coords      = coords,
                               data        = data, 
                               proj4string = crs)


spdf_cor = spdf[newMFin.sp,]
mapData2 = as.data.frame(spdf_cor)

Hab00aft = mapData2[, -c(7,8)]
Hab00aft = Hab00aft[,order(colnames(Hab00aft), decreasing = FALSE)]

Hab00aft = Hab00aft[, c(1,2,3,4,6,5)]


## test with best method
Ck = kmeans(Hab00aft, 4, iter.max = 1000, nstart = 50)
#Ck$cluster
table(Ck$cluster)

library(vioplot)
colors = c("darkcyan", 'forestgreen',"green3", "brown2")
colB = c(alpha('gray50', 0.5), alpha('brown', 0.5), alpha('green', 0.5), alpha('forestgreen', 0.5),
         alpha('cornflowerblue', 0.5), alpha('lightblue', 0.5))

#dev.off()
zones=matrix(c(5,1,2, 5,3,4), ncol=3, byrow=TRUE)
layout(zones, widths=c(3,5,5))
#par(mfrow=c(2,2), mar=c(2,2,2,2))
C1 = Hab00aft[Ck$cluster==1 ,]
vioplot(C1, col=colB, plotCentre  = 'line',
        main='Water (9.8%)', ylim=c(0,1),
        xaxt = 'n', cex=3, cex.axis=1.4)
axis(1, at=1:6, labels = c('urban', 'agri', 'SNH', 'forest', 'water B', 'wetland'), cex.axis=1.4)
box(col=colors[1], lwd=2)

C2 = Hab00aft[Ck$cluster==2 ,]
vioplot(C2, col=colB,plotCentre  = 'line',
        main='Forest (63.4%)', ylim=c(0,1),
        xaxt = 'n',  cex=3, cex.axis=1.4)
axis(1, at=1:6, labels = c('urban', 'agri', 'SNH', 'forest', 'water B', 'wetland'), cex.axis=1.4)
box(col=colors[2], lwd=2)

C3 = Hab00aft[Ck$cluster==3 ,]
vioplot(C3, col=colB,plotCentre  = 'line',
        main='Wetland - vegetation (14.6%)', ylim=c(0,1),
        xaxt = 'n',  cex=3, cex.axis=1.4)
axis(1, at=1:6, labels = c('urban', 'agri', 'SNH', 'forest', 'water B', 'wetland'), cex.axis=1.4)
box(col=colors[3], lwd=2)

C4 = Hab00aft[Ck$cluster==4 ,]
vioplot(C4, col=colB,plotCentre  = 'line',
        main='Heterogeneous (12.1%)', ylim=c(0,1),
        xaxt = 'n',  cex=3, cex.axis=1.4)
axis(1, at=1:6, labels = c('urban', 'agri', 'SNH', 'forest', 'water B', 'wetland'), cex.axis=1.4)
box(col=colors[4], lwd=2)


mapData3 = data.frame(mapData2, Ck$cluster)
mapData3$col[mapData3$Ck.cluster == 1] = 'Water'
mapData3$col[mapData3$Ck.cluster == 2] = 'Forest'
mapData3$col[mapData3$Ck.cluster == 3] = 'Wetland - vegetation'
mapData3$col[mapData3$Ck.cluster == 4] = 'Heterogeneous'

mapData3$Ck.cluster = as.factor(mapData3$Ck.cluster)
ggplot(data = mapData3, aes(x=X, y=Y, col=col))+geom_point(size=1.8) + 
  ggtitle("Clusters") +
  xlab("East coordinate (km)") + ylab("North coordinate (km)") +
  scale_colour_manual(values = c("Water" = "cyan3","Wetland - vegetation" = "green3",
                                 "Heterogeneous" = "brown2", "Forest" = "forestgreen")) + theme_minimal() + 
  theme(plot.title = element_text(size=0), axis.text=element_text(size=0),
        axis.title = element_text(size = 0), legend.text = element_text(size = 16),
        legend.title = element_text(size = 0)) + 
  theme(legend.position = c(0.9, 0.1)) +
        guides(colour = guide_legend(override.aes = list(size=10)))



save(Ck, C1, C2, C3, C4, mapData3, file='clusterFininfo.RDATA')

###--------------------------------------------------------------------------------------------------



## compare clustering algorithm
library(clValid)
library(mclust)
library(cluster)
# metric "euclidean", "correlation", and "manhattan"
# validity "internal", "stability", and "biological"
# clMethods "hierarchical", "kmeans", "diana", "fanny", "som", "model", "sota", "pam", "clara", and "agnes"
# method "ward", "single", "complete" and "average"

# test = clValid(Hab00, 4, clMethods = c('kmeans'), 
#                validation = c("stability", 'internal'), 
#                metric = 'manhattan', method = "ward")
# 
# summary(test)
# res <- getRanksWeights(test)  # does not work for manahattan
# print(res$ranks[,1:4], quote=FALSE)
# 
# op <- par(no.readonly=TRUE)
# par(mfrow=c(2,2),mar=c(4,4,3,1))
# plot(test, legend=FALSE)
# plot(nClusters(test),measures(test,"Dunn")[,,1],type="n",axes=F,
#      xlab="",ylab="")
# legend("center", clusterMethods(test), col=1:8, lty=1:98, pch=paste(1:8))
# par(op)



load('Coord_small_fin.RDATA')
# MAKE A DATAFRAME OF THE DATA TO BE PLOTTED
mapData=data.frame(Coord_small_fin, Hab00)
mapData=data.frame(Coord_small_fin, Hab06)
mapData=data.frame(Coord_small_fin, Hab12)
mapData=data.frame(Coord_small_fin, Hab18)

# get rid of extra popints in the south in the ocean?
library(mapsFinland)
library(maps)
library(mapdata)

par(mfrow=c(1,1), oma= c(2,2,2,2))
Finmap = map('worldHires','Finland', plot=FALSE, fill = TRUE)
class(Finmap)

Finmap$names[2] = "Finland1:?" 
Finmap$names[5] = "Finland2:?" 
Finmap$names[6] = "Finland3:?" 

library(maptools)
Fin.sp = map2SpatialPolygons(map = Finmap, IDs = Finmap$names, 
                             proj4string = CRS("+init=epsg:4326"), 
                             checkHoles=FALSE)
plot(Fin.sp)
newMFin.sp <- spTransform(Fin.sp, CRS("+proj=utm +zone=35 +ellps=GRS80 +units=m +no_defs"))
class(newMFin.sp)

# prepare coordinates, data, and proj4string
coords <- mapData[ , c("X", "Y")]   # coordinates
data   <- mapData[ , -c(1:2)]          # data
crs    <- CRS("+proj=utm +zone=35 +ellps=GRS80 +units=m +no_defs") # proj4string of coords

# make the SpatialPointsDataFrame object
spdf <- SpatialPointsDataFrame(coords      = coords,
                               data        = data, 
                               proj4string = crs)


spdf_cor = spdf[newMFin.sp,]
mapData2 = as.data.frame(spdf_cor)

Hab00aft = mapData2[, -c(7,8)]
## test with best method
Ck = kmeans(Hab00aft, 4, iter.max = 1000, nstart = 50)
#Ck$cluster
table(Ck$cluster)

Hab00 = Hab00aft[,order(colnames(Hab00aft), decreasing = FALSE)]

par(mfrow=c(2,2), mar=c(2,2,2,2))
C1 = Hab00[Ck$cluster==1 ,]
boxplot(C1, col=c('gray50', 'brown','green','forestgreen','cornflowerblue', 'lightblue'),
        main='Water', ylim=c(0,1))
box(col="lightblue3")

C2 = Hab00[Ck$cluster==2 ,]
boxplot(C2, col=c('gray50', 'brown','green','forestgreen','cornflowerblue', 'lightblue'),
        main='low vegetation', ylim=c(0,1))
box(col="green3")

C3 = Hab00[Ck$cluster==3 ,]
boxplot(C3, col=c('gray50', 'brown','green','forestgreen','cornflowerblue', 'lightblue'),
        main='heterogen', ylim=c(0,1))
box(col="orange3")

C4 = Hab00[Ck$cluster==4 ,]
boxplot(C4, col=c('gray50', 'brown','green','forestgreen','cornflowerblue', 'lightblue'),
        main='Forest', ylim=c(0,1))
box(col="forestgreen")

mapData3 = data.frame(mapData2, Ck$cluster)
mapData3$col[mapData3$Ck.cluster == 1] = 'Water'
mapData3$col[mapData3$Ck.cluster == 2] = 'LowVG'
mapData3$col[mapData3$Ck.cluster == 3] = 'Heterogen'
mapData3$col[mapData3$Ck.cluster == 4] = 'Forest'

mapData3$Ck.cluster = as.factor(mapData3$Ck.cluster)
ggplot(data = mapData3, aes(x=X, y=Y, col=col))+geom_point(size=2.9) + 
  ggtitle("Clusters") +
  xlab("East coordinate (km)") + ylab("North coordinate (km)") +
  scale_colour_manual(values = c("Water" = "lightblue3","LowVG" = "green3",
                                 "Heterogen" = "orange3", "Forest" = "forestgreen")) + 
  theme(plot.title = element_text(size=7.5), axis.text=element_text(size=6),
        axis.title = element_text(size = 6), legend.text = element_text(size = 12),
        legend.title = element_text(size = 12))




