

library(dplyr)

library(raster)
#library(bench)
library(plyr)
library(ggplot2)
library(ggpubr)
#library(rgdal)
library(tidyr)
library(dplyr)
#library(rgeos)
library(sf)

# Moth data
## Old sites to hget the names
setwd("D:/Helsinki/Moth_project/ResultsAug23/RDATA/2.ObsInfo")

Moth <- read.csv("Copy_mothMonitoring.csv", header=T)
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),]
Year = as.factor(sort(unique(Moth$YearStart)))
Site = sort(unique(Moth$SiteID))

Sit = Site



## Prepare the unique site year combination with new dataset
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.m2 = Moth.m[Moth.m$SiteID %in% Sit,]
dim(Moth.m2)


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

Year = as.factor(1998:2020)


# prepare the different buffer size coordinates
Moth.xy2 = Moth.xy[Moth.xy$SiteID %in% Sit,]


## ---- Cluster info

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


## ---- Traits info

## Traits info
Traits = read.csv('Traits_Pasi.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)
levels(Traits.sub$HostPlantForm) <- c('Herbaceous', 'Graminoids', 'Others', 'Others', 'Dwarf shrubs', 
                                      'Decideous', 'Conifer')


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)) # 1457 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)

Traitsp = Traits.sub[Traits.sub$Species %in% Moth.m2$Species,]
table(Traitsp$HostPlantForm)



## Whole period / whole space
load('spprevname.RDATA')  # 10% prevalence
#load('species5.RDATA')   # 5% prevalence

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

# sp5All = unique(unlist(species5))
# sp5All<-gsub(".", " ", sp5All, fixed=TRUE)
# 
# sp5All

# 10% prevalent species
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



# 10% prevalent species
Traits.sub2 = Traits.sub[Traits.sub$Species %in% unique(sp.p$Species),]  ## from preMod23 script
table(Traits.sub2$HostPlantForm)

# Traits.sub2$WL = NA
# Traits.sub2$WL[Traits.sub2$WingspanF<28] = 'W<28'
# Traits.sub2$WL[Traits.sub2$WingspanF>28 & Traits.sub2$WingspanF<34] = '28<W<34'
# Traits.sub2$WL[Traits.sub2$WingspanF>34] = 'W>34'
# 
# levels(Traits.sub2$WL) <- c('W<28', '28<W<34', 'W>34')
# table(Traits.sub2$WL)


Traits.sub2$WL = NA
Traits.sub2$WL[Traits.sub2$WingspanF<10] ='w = 0'
Traits.sub2$WL[Traits.sub2$WingspanF>=17.9 & Traits.sub2$WingspanF<24.1] ='17.9<w<24.1'
Traits.sub2$WL[Traits.sub2$WingspanF>=24.1 & Traits.sub2$WingspanF<30.3] ='24.1<w<30.3'
Traits.sub2$WL[Traits.sub2$WingspanF>=30.3 & Traits.sub2$WingspanF<36.5] ='30.3<w<36.5'
Traits.sub2$WL[Traits.sub2$WingspanF>=36.5 & Traits.sub2$WingspanF<43] ='36.5<w<43'
Traits.sub2$WL[Traits.sub2$WingspanF>=48] ='w=55.5'
table(Traits.sub2$WL)

Traits.sub2$WL = factor(Traits.sub2$WL, 
                                    levels = c('w = 0', '17.9<w<24.1', '24.1<w<30.3', '30.3<w<36.5',
                                               '36.5<w<43', 'w=55.5'))





remotes::install_github("ewenme/ghibli")
library(ghibli)
ghibli_palettes -> ghibli_palettes_list
ghibli_palettes$MononokeMedium

## histogram
par(mfrow=c(1,1))

Traits.sub2.0 = Traits.sub2#[Traits.sub2$WingspanF > 0,]

hist(Traits.sub2.0$WingspanF,  breaks = c(0,17.9,24.1,30.3,36.5,43,48,60), col = c("#838A90FF",  '#9FA7BEFF', '#B3B8B1FF',
                                                             '#E7A79BFF', '#F2C695FF','#F5EDC9FF'),
     main ='', xlab = 'Wing size', ylab = 'number of species')


hist(Traits.sub2.0$WingspanF,  breaks = c(0,1, seq(from=18, to=43, by=6.2), 55,60), col = c("#838A90FF",  'white','#9FA7BEFF', '#B3B8B1FF',
                                                                                   '#E7A79BFF', '#F2C695FF', 'white','#F5EDC9FF'),
     main ='', xlab = 'Wing size', ylab = 'number of species')



## custom

## a vector with the center of each bin
bin_ctr <- c(0, 21, 27.2, 33.4, 39.6, 55.5)
## a vector with the width of each bin
bin_wdt <- c(1, 6.2, 6.2, 6.2, 6.2, 1)

props = table(Traits.sub2$WL)

## put everything in a dataframe
bin_prop <- data.frame(bin_ctr, bin_wdt, props)

library(ggplot2)

ggplot(bin_prop, aes(x=bin_ctr, y=props, width=bin_wdt)) + 
  geom_bar(stat="identity", position="identity", 
           fill = c("#838A90FF",  '#9FA7BEFF', '#B3B8B1FF',
                    '#E7A79BFF', '#F2C695FF','#F5EDC9FF'), color = "black") +
  theme_minimal() + 
  theme(axis.title.x = element_text(size=18),
        axis.title.y = element_text(size=18),
        axis.text.x = element_text(size=16),
        axis.text.y = element_text(size=16)) + 
  ylab('Number of species') + xlab('Wing size')
















Moth.m3 = Moth.m2[Moth.m2$Species %in% sp10All,]  


#####  Presence
par(mfrow = c(1,2))

DF1 = merge(Moth.m3, Traits.sub2, by='Species')
DF2 = merge(DF1, DF, by='SiteID')


table(DF2$HostPlantForm, DF2$cluster)
table(DF2$WL, DF2$cluster)

DF.m= DF2 %>%
  dplyr::select(SiteID, Year, Species,Abundance, Effort, HostPlantForm, WL, cluster) %>%
  group_by(SiteID, Year, Species, HostPlantForm, WL, cluster) %>%
  dplyr::summarise(Abundance = sum(Abundance), Effort = sum(Effort))
head(DF.m)



par(mfrow = c(2,2))
# 1 hostplant form
### getting number of observation
tab.df2 = table(DF2$HostPlantForm, DF2$cluster)

x = barplot(tab.df2, col= c('salmon4', 'turquoise', 'purple', 'pink3', 'orange', 'forestgreen'), 
        cex.axis = 1.6, cex.names = 1.8, cex.lab = 1.8, names=c('', '', '', ''),
        legend.text = c('Herbaceous n=22', 'Graminoids n=8', 'Others n=7', 'Dwarf shrubs n=16', 
                       'Decideous n=22', 'Conifer n=3'), ylab = 'Number of observations', 
        args.legend = list(x = "topright", title= 'Species hostplant form',
                           inset = c(0, 0.1), cex=1.4))
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)
 
# DF.m= DF2 %>%
#   dplyr::select(SiteID, YearStart, Species,Abundance, Effort, HostPlantForm, cluster) %>%
#   group_by(SiteID, YearStart, Species, HostPlantForm, cluster) %>%
#   dplyr::summarise(Abundance = sum(Abundance), Effort = sum(Effort))
# head(DF.m)
# 
# Tab = table(DF.m$HostPlantForm, DF.m$cluster)
# 
# 
# barplot(Tab, col= c('salmon4', 'turquoise', 'purple', 'pink3', 'orange', 'forestgreen'), 
#         legend.text = TRUE, main = 'year -Site', 
#         args.legend = list(x = "topright",
#                            inset = c(0, -0.1)))



### with Abundance

DF.m2= DF2 %>%
  dplyr::select(Abundance, HostPlantForm, cluster) %>%
  group_by(HostPlantForm, cluster) %>%
  dplyr::summarise(Abundance = sum(Abundance))
head(DF.m2)

DFf = as.data.frame(DF.m2)

RSH = reshape(DFf, idvar = "HostPlantForm", timevar = "cluster", direction = "wide")

RSH2 = RSH[-7, -1]
rownames(RSH2) = RSH$HostPlantForm[-7]
colnames(RSH2) = c('Forest', 'Heterogeneous', 'Water', 'Wetland - vegetation')



barplot((as.matrix(RSH2)), col= c('salmon4', 'turquoise', 'purple', 'pink3', 'orange', 'forestgreen'), 
        cex.axis = 1.6, cex.names = 1.8, cex.lab = 1.8, names=c('', '', '', ''),
        legend.text = TRUE, ylab = 'Abundance', 
        args.legend = list(x = "topright",
                           inset = c(0, -0.4)))
axis(1, at = x, cex.axis=1.8, c('Forest', 'Heterogeneous', 'Water', 'Wetland -\n Vegetation'), 
     xpd=TRUE, tick = FALSE, padj=0.9)
title("b)", cex.main = 2,   font.main= 4, adj=0)



# 2 wing size
### getting number of observation
levels(DF2$WL) <- c('w = 0', '17.9<w<24.1', '24.1<w<30.3', '30.3<w<36.5',
                    '36.5<w<43', 'w=55.5')
tab.dfW = table(DF2$WL, DF2$cluster)


barplot(tab.dfW, col= c("#838A90FF",  '#9FA7BEFF', '#B3B8B1FF',
                        '#E7A79BFF', '#F2C695FF', '#F5EDC9FF'), 
        legend.text = c('w = 0, n=4', '17.9<w<24.1, n=16', '24.1<w<30.3, n=16', '30.3<w<36.5, n=30',
                        '36.5<w<43, n=11', 'w=55.5, n=1'), 
        cex.axis = 1.6, cex.names = 1.8, cex.lab = 1.8, 
        names=c('', '', '', ''), ylab = 'Number of observations', 
        args.legend = list(x = "topright", title= 'Species wing size',
                           inset = c(0, 0.1), cex=1.4))
axis(1, at = x, cex.axis=1.8, c('Forest', 'Heterogeneous', 'Water', 'Wetland -\n Vegetation'), 
     xpd=TRUE, tick = FALSE, padj=0.9)
title("c)", cex.main = 2,   font.main= 4, adj=0)


### with Abundance

DF.m2= DF2 %>%
  dplyr::select(Abundance, WL, cluster) %>%
  group_by(WL, cluster) %>%
  dplyr::summarise(Abundance = sum(Abundance))
head(DF.m2)

DFf = as.data.frame(DF.m2)

RSH = reshape(DFf, idvar = "WL", timevar = "cluster", direction = "wide")

RSH2 = RSH[-7, -1]
rownames(RSH2) = RSH$HostPlantForm[-7]
colnames(RSH2) = c('Forest', 'Heterogeneous', 'Water', 'Wetland - vegetation')

barplot((as.matrix(RSH2)), col= c("#838A90FF", '#BA968AFF', '#9FA7BEFF', '#B3B8B1FF',
                                  '#E7A79BFF', '#F2C695FF', '#F5EDC9FF'), names=c('', '', '', ''),
        legend.text = TRUE, ylab = 'Abundance', cex.axis = 1.6, cex.names = 1.8, cex.lab = 1.8,
        args.legend = list(x = "topright", 
                           inset = c(0, -0.4)))
axis(1, at = x, cex.axis=1.8, c('Forest', 'Heterogeneous', 'Water', 'Wetland -\n Vegetation'), 
     xpd=TRUE, tick = FALSE, padj=0.9)
title("d)", cex.main = 2,   font.main= 4, adj=0)

