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

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

library("class") # for knn classif

## mapping
library(sf)
library(maptools)
library(raster)   ## To convert an "Extent" object to a "SpatialPolygons" object.
library(rgeos)
library(mapsFinland)
library(maps)
library(mapdata)



## .......................................................  Finland clusters
setwd("D:/Helsinki/Moth_project/Samples from whole Finland")
load('PropFINall_Jan23.RDATA')  # not prep

head(Prop00)
length(Prop00)
#Prop06_0.5
#Prop12_0.5
#Prop18_0.5

load('Finland_HABall.RDATA') # prep
head(HabAll)
length(HabAll)
head(HabAll[[1]])

### info clusters Finland and mapData
load('clusterFininfo.RDATA')
head(Ck)
head(mapData3)

Ck$centers

## .......................................................  Moth sampling sites
setwd("D:/Helsinki/Moth_project/Res_jan23")

# Moth data
Moth <- read.csv("Copy_mothMonitoring.csv", header=T)
head(Moth)


Moth.m= Moth %>%
  dplyr::select(SiteID, YearStart, MonthStart, Species,YEUREF, XEUREF, Abundance, Effort) %>%
  group_by(SiteID, YearStart, MonthStart, Species) %>%
  dplyr::summarise(Abundance = sum(Abundance), Effort = sum(Effort),YEUREF=mean(YEUREF), XEUREF=mean(XEUREF),)
head(Moth.m)


Moth.xy = (unique(Moth[,c(5:7)]))
Moth.xy = Moth.xy[order(Moth.xy$SiteID),]

Site = Moth.xy$SiteID
head(Site)

setwd("D:/Helsinki/Moth_project/ResultsAug23/For submission/1.Data/1.EnvironmentalInfo")
load('Moth_hab0.5km.RDATA')
head(MatAllmean)

### get the moth sample under the same cluster types as whole Finland
#install.packages("proxy")
#library("proxy")

# FINcentroid = Ck$centers
# Mothsite.i = MatAllmean[1,]
# colnames(Mothsite.i) = colnames(FINcentroid)
# 
# M_to_FINc = rbind(Mothsite.i, FINcentroid)
# dist(M_to_FINc, method="euclidean")
# dist(M_to_FINc, method="manhattan")
# dist(M_to_FINc, method="canberra")
# dist(M_to_FINc, method="minkowski", p=3) # all are quite similar
# dist(M_to_FINc, method="minkowski", p=2) # is like euclidean dist
# dist(M_to_FINc, method="minkowski", p=5) #somehow similar
# 
# Dist_clusterFin = list()
# for (i in 1:length(Site)) {
#   
#   Mothsite.i = MatAllmean[i,]
#   colnames(Mothsite.i) = colnames(FINcentroid)
#   
#   M_to_FINc.i = rbind(Mothsite.i, FINcentroid)
#   dM_FIN.i = dist(M_to_FINc.i, method="canberra")
#   
#   Dist_clusterFin[[i]] = as.matrix(dM_FIN.i)[-1, 1]
# }
# 
# DistM_clusterFin = do.call(rbind, Dist_clusterFin)
# rownames(DistM_clusterFin) = Site
# 
# table(apply(DistM_clusterFin, 1, which.min))
# 
# 
# ####  with knn
# Moth_site = MatAllmean
# colnames(Moth_site) = colnames(FINcentroid)
# 
# Fin_site = mapData3[, c(1:4, 6, 5)]  # make sure we have the same order
# 
# # classify using knn
# Class_knn = knn(Fin_site, Moth_site, mapData3$col, k=4, prob=TRUE)
# table(Class_knn)
# 
# table(attr(Class_knn,"prob"))
# Moth_classif = as.data.frame(Class_knn)




## plotting the clusters Finland vs Moth sampling
table(Ck$cluster)

#colors = c("cyan3", 'forestgreen',"green3", "brown2")
colors = c("#66CCEE", "#228833",  "#CCBB44", "#EE6677")
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))





#### all

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

setwd("D:/Helsinki/Moth_project/Samples from whole Finland")
Fin_grid = st_read('3km_gridpoints.shp')
class(Fin_grid)

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

## Clip the map
#library(mapdata)
Finmap = maps::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?

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

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

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)]
#Hab00aft = Hab00



#dev.off()
## Finland
par(mar=c(2,1,3,2))
zones=matrix(c(9,1,2, 9,3,4, 10,5,6, 10,7,8), ncol=3, nrow=4, 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=4, cex.axis=1.8)
axis(1, at=1:6, labels = c('urban', 'agri', 'SNH', 'forest', 'water', 'wetland'), cex.axis=1.8)
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=4, cex.axis=1.8)
axis(1, at=1:6, labels = c('urban', 'agri', 'SNH', 'forest', 'water', 'wetland'), cex.axis=1.8)
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=4, cex.axis=1.8)
axis(1, at=1:6, labels = c('urban', 'agri', 'SNH', 'forest', 'water', 'wetland'), cex.axis=1.8)
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=4, cex.axis=1.8)
axis(1, at=1:6, labels = c('urban', 'agri', 'SNH', 'forest', 'water', 'wetland'), cex.axis=1.8)
box(col=colors[4], lwd=2)

## Sites
setwd("D:/Helsinki/Moth_project/Res_jan23")
load('MothclusterFIN.RDATA')
head(Moth_knn)
table(Moth_knn$Class_knn)/sum( table(Moth_knn$Class_knn))

Moth_classif = Moth_knn$Class_knn
Moth_site = Moth_knn[1:6]
C1m = Moth_site[Moth_classif == 'Water' ,]
vioplot(C1m, col=colB, plotCentre  = 'line',
        main='Water (14.7%)', ylim=c(0,1),
        xaxt = 'n', cex=4, cex.axis=1.8)
axis(1, at=1:6, labels = c('urban', 'agri', 'SNH', 'forest', 'water', 'wetland'), cex.axis=1.8)
box(col=colors[1], lwd=2)

C2m = Moth_site[Moth_classif == 'Forest',]
vioplot(C2m, col=colB, plotCentre  = 'line',
        main='Forest (42.2%)', ylim=c(0,1),
        xaxt = 'n',  cex=4, cex.axis=1.8)
axis(1, at=1:6, labels = c('urban', 'agri', 'SNH', 'forest', 'water', 'wetland'), cex.axis=1.8)
box(col=colors[2], lwd=2)

C3m = Moth_site[Moth_classif == 'Wetland - vegetation' ,]
vioplot(C3m, col=colB, plotCentre  = 'line',
        main='Wetland - vegetation (2.8%)', ylim=c(0,1),
        xaxt = 'n',  cex=4, cex.axis=1.8)
axis(1, at=1:6, labels = c('urban', 'agri', 'SNH', 'forest', 'water', 'wetland'), cex.axis=1.8)
box(col=colors[3], lwd=2)

C4m = Moth_site[Moth_classif == 'Heterogeneous',]
vioplot(C4m, col=colB, plotCentre  = 'line',
        main='Heterogeneous (40.4%)', ylim=c(0,1),
        xaxt = 'n',  cex=4, cex.axis=1.8)
axis(1, at=1:6, labels = c('urban', 'agri', 'SNH', 'forest', 'water', 'wetland'), cex.axis=1.8)
box(col=colors[4], lwd=2)

##


par(mar=c(2,1,3,2))
zones=matrix(c(1,2, 3,4, 5,6, 7,8), ncol=2, nrow=4, byrow=TRUE)
layout(zones, widths=c(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=4, cex.axis=1.8, cex.main = 2.2)
axis(1, at=1:6, labels = c('urban', 'agri', 'SNH', 'forest', 'water', 'wetland'), cex.axis=1.8)
box(col=colors[1], lwd=2)


#table(Class_knn)/sum( table(Class_knn))

C1m = Moth_site[Moth_classif == 'Water' ,]
vioplot(C1m, col=colB, plotCentre  = 'line',
        main='Water (14.7%)', ylim=c(0,1),
        xaxt = 'n', cex=4, cex.axis=1.8, cex.main = 2.2)
axis(1, at=1:6, labels = c('urban', 'agri', 'SNH', 'forest', 'water', 'wetland'), cex.axis=1.8)
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=4, cex.axis=1.8, cex.main = 2.2)
axis(1, at=1:6, labels = c('urban', 'agri', 'SNH', 'forest', 'water', 'wetland'), cex.axis=1.8)
box(col=colors[2], lwd=2)


C2m = Moth_site[Moth_classif == 'Forest',]
vioplot(C2m, col=colB, plotCentre  = 'line',
        main='Forest (42.2%)', ylim=c(0,1),
        xaxt = 'n',  cex=4, cex.axis=1.8, cex.main = 2.2)
axis(1, at=1:6, labels = c('urban', 'agri', 'SNH', 'forest', 'water', 'wetland'), cex.axis=1.8)
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=4, cex.axis=1.8, cex.main = 2.2)
axis(1, at=1:6, labels = c('urban', 'agri', 'SNH', 'forest', 'water', 'wetland'), cex.axis=1.8)
box(col=colors[3], lwd=2)



C3m = Moth_site[Moth_classif == 'Wetland - vegetation' ,]
vioplot(C3m, col=colB, plotCentre  = 'line',
        main='Wetland - vegetation (2.8%)', ylim=c(0,1),
        xaxt = 'n',  cex=4, cex.axis=1.8, cex.main = 2.2)
axis(1, at=1:6, labels = c('urban', 'agri', 'SNH', 'forest', 'water', 'wetland'), cex.axis=1.8)
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=4, cex.axis=1.8, cex.main = 2.2)
axis(1, at=1:6, labels = c('urban', 'agri', 'SNH', 'forest', 'water', 'wetland'), cex.axis=1.8)
box(col=colors[4], lwd=2)


C4m = Moth_site[Moth_classif == 'Heterogeneous',]
vioplot(C4m, col=colB, plotCentre  = 'line',
        main='Heterogeneous (40.4%)', ylim=c(0,1),
        xaxt = 'n',  cex=4, cex.axis=1.8, cex.main = 2.2)
axis(1, at=1:6, labels = c('urban', 'agri', 'SNH', 'forest', 'water', 'wetland'), cex.axis=1.8)
box(col=colors[4], lwd=2)

mtext('Whole Finland', side=3, outer=TRUE, line = 0, adj = 0, cex=1.7)
mtext('Moth sampling sites', side=3, outer=TRUE, line = 0, adj = 0.98, cex=1.7)
###



## plotting Moth sampling sites clusters
Moth.spdf <- SpatialPointsDataFrame(Moth.xy[,c(3,2)], Moth.xy)
proj4string(Moth.spdf) <-CRS("+proj=utm +zone=35 +ellps=GRS80 +units=m +no_defs")

newMoth<-spTransform(Moth.spdf, CRS("+init=epsg:4326"))
#plot(newMoth, main="Map of Moths Locations")
# make the SpatialPointsDataFrame object
spdf.moth <- SpatialPointsDataFrame(coords      = Moth_knn[, c(11:12)],
                               data        = Moth_knn, 
                               proj4string = CRS("+init=epsg:4326"))


# plotting

Moth_knn = data.frame(Moth_site, Moth_classif, newMoth)
head(Moth_knn)

colors = data.frame(color=c("#66CCEE", "#228833",  "#CCBB44", "#EE6677"), Ck = c('Water', 'Forest', 'Wetland - vegetation', 'Heterogeneous'))
colors


Fin.spdf <- SpatialPointsDataFrame(mapData3[,c(7,8)], mapData3)
proj4string(Fin.spdf) <-CRS("+proj=utm +zone=35 +ellps=GRS80 +units=m +no_defs")
newFin <- spTransform(Fin.spdf, CRS("+init=epsg:4326"))

#general map
par(mfrow=c(1,2), mar=c(6,1,3,6), oma=c(2,2,2,2), xpd = T)
#zones=matrix(c(1,3, 2,3), ncol=2, nrow=2, byrow=FALSE)
#layout(zones, widths=c(5,5), heights = c(5,2))

maps::map(Finmap, fill=TRUE, col = 'gray80')
points(newFin, col=colors$color[match(mapData3$col, colors$Ck)], pch=16, cex=0.65)
legend(x=24, y=59.5, legend = c('Water', 'Forest', 'Wetland - vegetation', 'Heterogeneous'), 
       col = c("#66CCEE", "#228833",  "#CCBB44", "#EE6677"), pch = 16, bty = 'n', cex=1.4, xpd = NA, ncol=2)

map = maps::map(Finmap, fill=TRUE, col = 'gray80')

color.choice = colors$color[match(as.data.frame(spdf.moth)$Moth_classif, colors$Ck)]
points(spdf.moth, col=color.choice, pch=16, cex=1.8)


# Other plotting
par(mfrow=c(2,1), mar=c(0,0,0,0), oma=c(4,0,0,0), xpd = T)

maps::map(Finmap, fill=TRUE, col = 'gray80')
points(newMoth, col=colors$color[match(Moth_knn$Class_knn, colors$Ck)], pch=16, cex=1.2)

maps::map(Finmap, fill=TRUE, col = 'gray80')
points(newFin, col=colors$color[match(mapData3$col, colors$Ck)], pch=16, cex=0.65)

legend(x=16, y=59.5, legend = c('Water', 'Forest', 'Wetland - vegetation', 'Heterogeneous'), 
       col = c("#66CCEE", "#228833",  "#CCBB44", "#EE6677"), pch = 16, bty = 'n', cex=1.4, 
       xpd = NA, ncol=2)



save(Moth_knn, file='MothclusterFIN.RDATA')
