library(dplyr)
library(tidyr)


##----------------------------------------------------------------------------------------
#
#			Preparation of VP cond - Trait information
#					Finnish grid sites
##----------------------------------------------------------------------------------------

# Predictions from model - grid info
setwd("E:/Helsinki/ResultsAug23/For submission/1.Data/3.ModelRun")
load('Predsmall2408.RDATA')

hM = modelyear[[1]]

######++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# for multiple species
# looking at traits
spNames = hM$spNames

## Traits info
setwd("E:/Helsinki/ResultsAug23/For submission/1.Data/2.ObsInfo")
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)

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

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

Moth.TR = Traits.sub[ Traits.sub$Species %in% spNames,]
table(Moth.TR$HostPlantForm)
table(Moth.TR$Diet)

## order from species hmsc order
Moth.TR = Moth.TR %>%
  slice(match(spNames, Species))


### Finland .........................................................
## VP 

setwd("E:/Helsinki/ResultsAug23/For submission/1.Data/3.ModelRun")
load('VPbwFIN_PA01408.RDATA')
load('VPcondAFIN100cdd_PA_1408.RDATA')
VP.a.PA = VP.a
length(VP.a.PA)
dim(VP.a.PA[[1]])

load('VPbwFIN_AB01608.RDATA')
load('VPcondAFIN100cdd_AB_1408.RDATA')
VP.a.AB = VP.a
dim(VP.a.AB[[1]])


## why is different from HMsc calc?
apply(VP.a.PA[[1]], c(1), mean)
apply(VP.a.AB[[1]], c(1), mean)

apply(VP.a.PA[[1]], c(1), sd)
apply(VP.a.AB[[1]], c(1), sd)

###
VP.a.PA2 <- aperm(VP.a.PA[[1]], c(2,3,1))
dim(VP.a.PA2)

VP_PA.prep = list()
for (i in 1:2) {
  VP.a.PA2_i = VP.a.PA2[,,i]
  dim(VP.a.PA2_i)
  vec_pa_i = c(VP.a.PA2_i)
  VP_PA.prep[[i]] = data.frame(VP = vec_pa_i, Species = rep(rownames(VP.a.PA2_i), 100))
}


VP_PA.prep[[1]]$Var = 'Clim'
VP_PA.prep[[2]]$Var = 'Hab'

VP_PA.p = do.call(rbind, VP_PA.prep)
VP_PA.p = data.frame(VP_PA.prep[[1]], VPhab = VP_PA.prep[[2]]$VP)
                     
##

VP.a.AB2 <- aperm(VP.a.AB[[1]], c(2,3,1))
dim(VP.a.AB2)

VP_AB.prep = list()
for (i in 1:2) {
  VP.a.AB2_i = VP.a.AB2[,,i]
  dim(VP.a.AB2_i)
  vec_AB_i = c(VP.a.AB2_i)
  VP_AB.prep[[i]] = data.frame(VP = vec_AB_i, Species = rep(rownames(VP.a.AB2_i), 100))
}


VP_AB.prep[[1]]$Var = 'Clim'
VP_AB.prep[[2]]$Var = 'Hab'

VP_AB.p = do.call(rbind, VP_AB.prep)
VP_AB.p = data.frame(VP_AB.prep[[1]], VPhab = VP_AB.prep[[2]]$VP)
#########


VP_PA.all = merge(VP_PA.p, Moth.TR, by='Species')
VP_AB.all = merge(VP_AB.p, Moth.TR, by='Species')

levels(VP_PA.all$HostPlantForm) = c('Herbaceous', 'Graminoids', 'Others', 'Others',
                                       'Shurbs', 'Decideous', 'Coniferous')
levels(VP_PA.all$HostPlantForm) 

levels(VP_AB.all$HostPlantForm) = c('Herbaceous', 'Graminoids', 'Others', 'Others',
                                    'Shurbs', 'Decideous', 'Coniferous')
levels(VP_AB.all$HostPlantForm) 

save(VP_PA.all, VP_AB.all, file='VPall_traits2FIN.RDATA')



# Sites ................................................................................
## VP

load('VPycdd-PAclust_1108.RDATA')
VP.a.PA = VP.a
length(VP.a.PA)
dim(VP.a.PA[[1]])

load('VPycdd-ABclust_1108.RDATA')
VP.a.AB = VP.a
dim(VP.a.AB[[1]])


## why is different from HMsc calc?
apply(VP.a.PA[[1]], c(1), mean)
apply(VP.a.AB[[1]], c(1), mean)

apply(VP.a.PA[[1]], c(1), sd)
apply(VP.a.AB[[1]], c(1), sd)

###
VP.a.PA2 <- aperm(VP.a.PA[[1]], c(2,3,1))
dim(VP.a.PA2)

VP_PA.prep = list()
for (i in 1:2) {
  VP.a.PA2_i = VP.a.PA2[,,i]
  dim(VP.a.PA2_i)
  vec_pa_i = c(VP.a.PA2_i)
  VP_PA.prep[[i]] = data.frame(VP = vec_pa_i, Species = rep(rownames(VP.a.PA2_i), 1000))
}


VP_PA.prep[[1]]$Var = 'Clim'
VP_PA.prep[[2]]$Var = 'Hab'

VP_PA.p = do.call(rbind, VP_PA.prep)
VP_PA.p = data.frame(VP_PA.prep[[1]], VPhab = VP_PA.prep[[2]]$VP)

##

VP.a.AB2 <- aperm(VP.a.AB[[1]], c(2,3,1))
dim(VP.a.AB2)

VP_AB.prep = list()
for (i in 1:2) {
  VP.a.AB2_i = VP.a.AB2[,,i]
  dim(VP.a.AB2_i)
  vec_AB_i = c(VP.a.AB2_i)
  VP_AB.prep[[i]] = data.frame(VP = vec_AB_i, Species = rep(rownames(VP.a.AB2_i), 1000))
}


VP_AB.prep[[1]]$Var = 'Clim'
VP_AB.prep[[2]]$Var = 'Hab'

VP_AB.p = do.call(rbind, VP_AB.prep)
VP_AB.p = data.frame(VP_AB.prep[[1]], VPhab = VP_AB.prep[[2]]$VP)
#########


VP_PA.all = merge(VP_PA.p, Moth.TR, by='Species')
VP_AB.all = merge(VP_AB.p, Moth.TR, by='Species')

levels(VP_PA.all$HostPlantForm) = c('Herbaceous', 'Graminoids', 'Others', 'Others',
                                    'Shurbs', 'Decideous', 'Coniferous')
levels(VP_PA.all$HostPlantForm) 

levels(VP_AB.all$HostPlantForm) = c('Herbaceous', 'Graminoids', 'Others', 'Others',
                                    'Shurbs', 'Decideous', 'Coniferous')
levels(VP_AB.all$HostPlantForm) 

save(VP_PA.all, VP_AB.all, file='VPall_traits_Sitesprep.RDATA')
