library(tidyr)
library(dplyr)
library(RColorBrewer)

# Environmental data
## CLIM

setwd("D:/Helsinki/Moth_project/ResultsAug23/For submission/1.Data/1.EnvironmentalInfo")

load('ClimSite080623.RDATA')
head(ClimSite)
ClimSite = ClimSite[,-c(2:3)]
colnames(ClimSite)[1:2] = c('SiteID', 'YearStart')

ClimSite2 = ClimSite[, c(1,2,15, 5,9,13,16)] ## with gdd5

Clim.y = ClimSite2 %>%
  group_by(SiteID, YearStart) %>%
  dplyr::summarise(precip = mean(precip), GDD5w = mean(GDD5w), CDD5w = mean(CDD5w), 
                   GDD5s = mean(GDD5s), snow = mean(snow),)
head(Clim.y)
dim(Clim.y) ## 

Clim.y = Clim.y[Clim.y$YearStart > 1997,]
Clim.y = Clim.y[Clim.y$YearStart < 2021,]


# setwd("D:/Helsinki/Moth_project/ResultsAug23/For submission/1.Data/2.ObsInfo")
## load unique sites (109)
load('MothSites109.RDATA')

MothSite = data.frame(SiteID = rep(Moth.m2$SiteID, 29),
                       YEUREF= rep(Moth.m2$YEUREF, 29),
                       XEUREF = rep(Moth.m2$XEUREF, 29))#,
                       #Year = rep(1998:2021, each=109))
head(MothSite)

Site = sort(unique(MothSite$SiteID))
Sit = Site


## Prepare the unique site year combination
Moth <- read.csv("Moth_Monitoring_toRecbase_19072023.csv", header=T)
head(Moth)
dim(Moth)

Moth = Moth[Moth$Month > 3,]
Moth = Moth[Moth$Month < 11,]
Mothsub = Moth[Moth$Month == 10 & Moth$Day > 16,]
M = subset(Moth, !(rownames(Moth) %in% rownames(Mothsub)))

Moth = M
dim(Moth)

Moth = Moth[Moth$Year < 2021,]
Moth = Moth[Moth$Year > 1997,]
dim(Moth)

Moth$Abundance[is.null(Moth$Abundance)] = 0
Moth$Abundance = as.numeric(Moth$Abundance)
Moth$Abundance[is.na(Moth$Abundance)] = 0

Moth.m= Moth %>%
  dplyr::select(SiteID, Year, Species, Abundance, Effort) %>%
  group_by(SiteID, Year, Species) %>%
  dplyr::summarise(Abundance = sum(Abundance), Effort = sum(Effort))
head(Moth.m)

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

Moth.m = Moth.m[Moth.m$SiteID %in% Sit,]
dim(Moth.m)


MothSY = unique(Moth.m[, 1:2])
dim(MothSY)
MothSY$SY = paste(MothSY$SiteID, '', MothSY$Year)

Year = as.factor(1998:2020)

# Subset clim data

Clim.y$SY = paste(Clim.y$SiteID, '', Clim.y$YearStart)
Clim.y2 = Clim.y[Clim.y$SY %in% MothSY$SY,]

#.......................................................................
## HAB
load('HabSites1108.RDATA')

#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

## Merging information
### Adding clim data

dim(Clim.y2)

### Adding CLC class_0.5es info
### Adding fragmentation measure

Mclc.all = merge(Metricsclc, Mclc, by=c('SiteID', 'YearStart'), all.x=T)
summary(Mclc.all)
dim(Mclc.all)

Moth.env = merge(Clim.y2, Mclc.all, by=c('SiteID', 'YearStart'), all.x=T)
summary(Moth.env)
dim(Moth.env)

Moth.envall = Moth.env
summary(Moth.envall)
dim(Moth.envall)


# Boreal regions

load('BorealRegions.RDATA')
NB$BR = rep('NB', length(NB$SiteID))
CB$BR = rep('CB', length(CB$SiteID))
SB$BR = rep('SB', length(SB$SiteID))
HB$BR = rep('HB', length(HB$SiteID))

BR = rbind(NB, CB, SB, HB)
BR2 = BR[BR$SiteID %in% Sit,]
Env.br = merge(Moth.envall, BR2[,c(1,5)], by=c('SiteID'), all.y=T)

head(Env.br)
tail(Env.br)
dim(Env.br)

##...............................................................................................
###  Prepare occurences data

Moth.m= Moth %>%
  dplyr::select(SiteID, Year, Species, Abundance, Effort) %>%
  group_by(SiteID, Year, Species) %>%
  dplyr::summarise(Abundance = sum(Abundance), Effort = sum(Effort))
head(Moth.m)
dim(Moth.m)


Moth.m2 = Moth.m[Moth.m$SiteID %in% Sit,]

## Traits info
Traits = read.csv('Traits_Moths.csv', header=T)
head(Traits)

Traits.sub= Traits %>%
  dplyr::select(Genus.and.species, female.size..average.wingspan..mm.,
                HostPlantForm, Diet, male.size..average.wingspan..mm.)

Traits.sub$HostPlantForm = as.factor(Traits.sub$HostPlantForm)

head(Traits.sub)
colnames(Traits.sub)[1] = 'Species'
colnames(Traits.sub)[2] = 'WingspanF'
colnames(Traits.sub)[5] = 'WingspanM'

Moth.TR = merge(Moth.m2, Traits.sub, by=c('Species'), all=T)
head(Moth.TR)

Moth.TR2 = Moth.TR[!is.na(Moth.TR$SiteID),] # some extra species that we don't have info on
summary(Moth.TR2)
length(unique(Moth.TR2$Species)) # 1441 species

Moth.TR = Moth.TR2[!is.na(Moth.TR2$WingspanF),] # some extra species that we don't have info on
summary(Moth.TR)
length(unique(Moth.TR$Species)) # 503 species

dim(Moth.m2)
dim(Moth.TR)


# 10% prevalent species
## Get occurences in Shape for the model
load('spprevname.RDATA') ## Focal species 10% prevalence

sp10All<-gsub(".", " ", sp10All, fixed=TRUE)
sp10All


sp.p = Moth.TR[Moth.TR$Species %in% sp10All,]  
head(sp.p)
dim(sp.p)

sp40.brY=sp.p %>%
  dplyr::select(SiteID, Year, Species, Abundance)%>%
  group_by(SiteID, Year, Species) %>%
  dplyr::summarise(Abundance= sum(Abundance))
dim(sp40.brY) # same dim ok but less species because redundant ones or not present during that new period?

sp40.brY_spread <- sp40.brY %>%
  spread(key = Species, value = Abundance)
dim(sp40.brY_spread)
colnames(sp40.brY_spread)  # Still with SiteID and YearStart



### Get Effort in shape for data
spEff=sp.p %>%
  dplyr::select(SiteID, Year, Effort)%>%
  group_by(SiteID, Year) %>%
  dplyr::summarise(Effort= sum(Effort))
dim(spEff)
colnames(spEff) ## getting Effort
head(spEff)

sp40Y_spread = sp40.brY_spread


## if all species
sp40.brYA=Moth.TR %>%
  dplyr::select(SiteID, Year, Species, Abundance)%>%
  group_by(SiteID, Year, Species) %>%
  dplyr::summarise(Abundance= sum(Abundance))
dim(sp40.brYA) # same dim ok but less species because redundant ones or not present during that new period?

sp40.brY_spreadA <- sp40.brYA %>%
  spread(key = Species, value = Abundance)
dim(sp40.brY_spreadA) # same because same site/year
colnames(sp40.brY_spreadA)  # Still with SiteID and YearStart

# total abundance
sp.p %>%
  dplyr::select(Abundance)%>%
  #group_by(SiteID, YearStart, Species) %>%
  dplyr::summarise(Abundance= sum(Abundance))


# Overall abundance
Moth.TR %>%
  dplyr::select(Abundance)%>%
  #group_by(SiteID, YearStart, Species) %>%
  dplyr::summarise(Abundance= sum(Abundance))

##
# prepare trait info
spT = sp.p %>%
  dplyr::select(Species, WingspanF, HostPlantForm, WingspanM, Diet)%>% 
  group_by(Species) 

spT= unique(spT)
summary(spT)
dim(spT)

par(mfrow=c(2,2))
hist(spT$WingspanF)
hist(spT$WingspanF, breaks=20)
hist(spT$WingspanM)
hist(spT$WingspanM, breaks=20)

wingDiff = spT$WingspanF[spT$WingspanF>0]-spT$WingspanM[spT$WingspanF>0]
summary(wingDiff)


##### --------------------------------------------------------------------
# to plot only
Moth.br = merge(sp.p, BR2[,c(1,5)], by=c('SiteID'), all.y=T)

m.env = Moth.br %>%
  dplyr::select(SiteID, Year, Effort, Species, Abundance,
                WingspanF, HostPlantForm, Diet, BR)%>%
  group_by(SiteID, Year, BR) %>%
  dplyr::summarise(Effort= sum(Effort), Abundance = sum(Abundance))


dim(m.env)
summary(m.env)


load('MothclusterFIN.RDATA')
head(Moth_knn)
DF = data.frame(SiteID = Moth_knn$SiteID, cluster = Moth_knn$Class_knn)
head(DF)


m.env2 = merge(m.env, DF, by=c('SiteID'), all.x=T)
head(m.env2)



#### Plotting
{
  par(mfrow=c(3,2), mar=c(4,5,2,2), oma=c(2,1,2,2), xpd=NA)
  
  ## Boreal regions -------------------------------------------
  BR_nb = m.env %>%
    dplyr::select(BR, Abundance) %>%
    group_by(BR) %>%
    dplyr::summarize(Abundance = sum(Abundance))
  
  barplot(BR_nb$Abundance[c(3,1,4,2)], col=c("#FDE725FF","#35B779FF", "#31688EFF", "#440154FF"),
          names=c('NB', 'CB', 'SB', 'HB'), ylab = 'Total number of observations', 
          cex.axis = 1.2, cex.names = 1.4, cex.lab = 1.4)
 
  
  table(m.env$BR, m.env$SiteID)
  BR_nb = m.env %>%
    dplyr::select(SiteID, BR) %>%
    group_by(BR) %>%
    dplyr::summarize(siteNB = n_distinct(SiteID))
  
  barplot(BR_nb$siteNB[c(3,1,4,2)], col=c("#FDE725FF","#35B779FF", "#31688EFF", "#440154FF"),
          names=c('NB', 'CB', 'SB', 'HB'), ylab = 'Total number of sites', 
          cex.axis = 1.2, cex.names = 1.4, cex.lab = 1.4)
  mtext('', outer=TRUE, line=-1.5)
  
  
  
  par(mfrow=c(2,2), mar=c(4,5,2,2), oma=c(2,1,2,2), xpd=NA)
  ## clusters sites ----------------------------------------------
  BR_nb = m.env2 %>%
    dplyr::select(cluster, Abundance) %>%
    group_by(cluster) %>%
    dplyr::summarize(Abundance = sum(Abundance))
  
  x <- barplot(BR_nb$Abundance[c(2,1,3,4)], col=c("#228833", "#EE6677", "#66CCEE", "#CCBB44"),
          names=c('', '', '', ''), ylab = 'Total number of observations', 
          cex.axis = 1.6, cex.names = 1.8, cex.lab = 1.8)
  axis(1, at = x, cex.axis=1.8, c('Forest', 'Heterogeneous', 'Water', 'Wetland -\n Vegetation'), 
       xpd=TRUE, tick = FALSE, padj=0.9)
  title("a)", cex.main = 2,   font.main= 4, adj=0)
  
  table(m.env2$BR, m.env2$SiteID)
  BR_nb = m.env2 %>%
    dplyr::select(SiteID, cluster) %>%
    group_by(cluster) %>%
    dplyr::summarize(siteNB = n_distinct(SiteID))
  
  barplot(BR_nb$siteNB[c(2,1,3,4)], col=c("#228833", "#EE6677", "#66CCEE", "#CCBB44"),
          names=c('', '', '', ''), ylab = 'Total number of sites', 
          cex.axis = 1.6, cex.names = 1.8, cex.lab = 1.8)
  axis(1, at = x, cex.axis=1.8, c('Forest', 'Heterogeneous', 'Water', 'Wetland -\n Vegetation'), 
       xpd=TRUE, tick = FALSE, padj=0.9)
  mtext('', outer=TRUE, line=-1.5)
  title("b)", cex.main = 2,   font.main= 4, adj=0)
  
  # years ---------------------------------------------------------------
  library(colorspace)
  q26 <- qualitative_hcl(26, "Dark2")
  
  m.env = m.env[m.env$Year < 2021 ,]
  
  BR_nb = m.env %>%
    dplyr::select(Year, Abundance) %>%
    group_by(Year) %>%
    dplyr::summarize(Abundance = sum(Abundance))
  barplot(BR_nb$Abundance~ BR_nb$Year, col='gray60', las=3, xlab='',
          ylab = '', cex.axis = 1.6, cex.names = 1.8, cex.lab = 1.8)
  title(ylab = "Total number of observations", line = 4.4, cex.lab = 1.8)
  title("c)", cex.main = 2,   font.main= 4, adj=0)
  
  table(m.env$Year, m.env$SiteID)
  BR_nb = m.env %>%
    dplyr::select(SiteID, Year) %>%
    group_by(Year) %>%
    dplyr::summarise(yearNB = n_distinct(SiteID))
  
  barplot(BR_nb$yearNB ~ BR_nb$Year, col='gray60', las=3, xlab='',
          ylab = 'Total number of sites sampled', 
          cex.axis = 1.6, cex.names = 1.8, cex.lab = 1.8)
  title("d)", cex.main = 2,   font.main= 4, adj=0)
  
}





####

Env.br$Effort = spEff$Effort

XData = Env.br[-c(1:2, 8, 16)]
head(XData)

par(mfrow=c(1,1))
library(corrplot)
X.cor = cor(XData)
testRes = cor.mtest(X.cor, conf.level = 0.95)
corrplot(X.cor, p.mat = testRes$p, method = 'square', type = 'lower', insig='blank',
         order = 'AOE', diag = FALSE)$corrPos -> p1
text(p1$x, p1$y, round(p1$corr, 2))



### Save prepData

save(Env.br, spT, sp40Y_spread, BR2, file = 'prep_centCLC_A23.RDATA')

## +++++++++++++++++++++++++++++++++++++++++++++
