##########################################################################################
#
#		Traits prediction preparation - Finnish grid sites Sites - emprirical
#
##########################################################################################

{
  library(Hmsc)
  library(tidyr)
  library(dplyr)
  library(abind)

}
# Get grid site info
setwd("E:/Helsinki/ResultsAug23/For submission/1.Data/1.EnvironmentalInfo")
#........................................................................................
load('CurrentSmallFIN_1408.RDATA')
tail(mapData2)
tail(mapData2hab)

mapDataALL = cbind(mapData2[, c(1:2, 6, 15, 10:11)], mapData2hab[, c(2:10)])
head(mapDataALL)
dim(mapDataALL)
mapcoord1 = mapDataALL[mapDataALL$Year == 1999, 14:15]
mapDataALL$SiteID = paste(mapDataALL$X, '', mapDataALL$Y)
SiteFIn_un = unique(mapDataALL$SiteID)
MapdataALL.list = split(mapDataALL, f = mapDataALL$Year) 
Fin_coord = MapdataALL.list[[1]][, 14:15]



## Results from predictions
#........................................................................................
setwd("E:/Helsinki/ResultsAug23/For submission/1.Data/4.Predictions/1.PredRich")

load('predPArich_Y1.RDATA')
load('predPArich_Y3.RDATA')
load('predPArich_Y5.RDATA')
load('predPArich_Y7.RDATA')
load('predPArich_Y9.RDATA')
load('predPArich_Y11.RDATA')
load('predPArich_Y13.RDATA')
load('predPArich_Y15.RDATA')
load('predPArich_Y17.RDATA')
load('predPArich_Y19.RDATA')
load('predPArich_Y21.RDATA')
load('predPArich_Y23.RDATA')
load('predPArich_Y2.RDATA')
load('predPArich_Y4.RDATA')
load('predPArich_Y6.RDATA')
load('predPArich_Y8.RDATA')
load('predPArich_Y10.RDATA')
load('predPArich_Y12.RDATA')
load('predPArich_Y14.RDATA')
load('predPArich_Y16.RDATA')
load('predPArich_Y18.RDATA')
load('predPArich_Y20.RDATA')
load('predPArich_Y22.RDATA')

## All

dat.iALLb = c(S_Y1, S_Y2, S_Y3, S_Y4, S_Y5, S_Y6, S_Y7, S_Y8, S_Y9, S_Y10, S_Y11, 
              S_Y12, S_Y13, S_Y14, S_Y15, S_Y16, S_Y17, S_Y18, S_Y19, S_Y20, S_Y21, S_Y22, S_Y23)
dat.iALL = data.frame(S = dat.iALLb, 
                      Year = c(rep(1998, length(S_Y1)), rep(1999, length(S_Y2)), rep(2000, length(S_Y3)), rep(2001, length(S_Y4)),
                               rep(2002, length(S_Y5)),  rep(2003, length(S_Y6)),rep(2004, length(S_Y7)), rep(2005, length(S_Y8)),
                               rep(2006, length(S_Y9)),  rep(2007, length(S_Y10)),rep(2008, length(S_Y11)), rep(2009, length(S_Y12)),
                               rep(2010, length(S_Y13)),  rep(2011, length(S_Y14)),rep(2012, length(S_Y15)), rep(2013, length(S_Y16)),
                               rep(2014, length(S_Y17)),  rep(2015, length(S_Y18)),rep(2016, length(S_Y19)), rep(2017, length(S_Y20)),
                               rep(2018, length(S_Y21)),  rep(2019, length(S_Y22)),rep(2020, length(S_Y23))))



#####
# Trend and residuals
######........................................................................................

##------------------ WL - traits

dat.iW1 = c(predYRW1_Y1, predYRW1_Y2, predYRW1_Y3, predYRW1_Y4, predYRW1_Y5, predYRW1_Y6, predYRW1_Y7, predYRW1_Y8, predYRW1_Y9, predYRW1_Y10, predYRW1_Y11, 
              predYRW1_Y12, predYRW1_Y13, predYRW1_Y14, predYRW1_Y15, predYRW1_Y16, predYRW1_Y17, predYRW1_Y18, predYRW1_Y19, predYRW1_Y20, predYRW1_Y21, predYRW1_Y22, predYRW1_Y23)
dat.iW1 = data.frame(S = dat.iW1, 
                     Year = c(rep(1998, length(S_Y1)), rep(1999, length(S_Y2)), rep(2000, length(S_Y3)), rep(2001, length(S_Y4)),
                               rep(2002, length(S_Y5)),  rep(2003, length(S_Y6)),rep(2004, length(S_Y7)), rep(2005, length(S_Y8)),
                               rep(2006, length(S_Y9)),  rep(2007, length(S_Y10)),rep(2008, length(S_Y11)), rep(2009, length(S_Y12)),
                               rep(2010, length(S_Y13)),  rep(2011, length(S_Y14)),rep(2012, length(S_Y15)), rep(2013, length(S_Y16)),
                               rep(2014, length(S_Y17)),  rep(2015, length(S_Y18)),rep(2016, length(S_Y19)), rep(2017, length(S_Y20)),
                               rep(2018, length(S_Y21)),  rep(2019, length(S_Y22)),rep(2020, length(S_Y23))))



dat.iW2 = c(predYRW2_Y1, predYRW2_Y2, predYRW2_Y3, predYRW2_Y4, predYRW2_Y5, predYRW2_Y6, predYRW2_Y7, predYRW2_Y8, predYRW2_Y9, predYRW2_Y10, predYRW2_Y11, 
              predYRW2_Y12, predYRW2_Y13, predYRW2_Y14, predYRW2_Y15, predYRW2_Y16, predYRW2_Y17, predYRW2_Y18, predYRW2_Y19, predYRW2_Y20, predYRW2_Y21, predYRW2_Y22, predYRW2_Y23)
dat.iW2 = data.frame(S = dat.iW2, 
                     Year = c(rep(1998, length(S_Y1)), rep(1999, length(S_Y2)), rep(2000, length(S_Y3)), rep(2001, length(S_Y4)),
                               rep(2002, length(S_Y5)),  rep(2003, length(S_Y6)),rep(2004, length(S_Y7)), rep(2005, length(S_Y8)),
                               rep(2006, length(S_Y9)),  rep(2007, length(S_Y10)),rep(2008, length(S_Y11)), rep(2009, length(S_Y12)),
                               rep(2010, length(S_Y13)),  rep(2011, length(S_Y14)),rep(2012, length(S_Y15)), rep(2013, length(S_Y16)),
                               rep(2014, length(S_Y17)),  rep(2015, length(S_Y18)),rep(2016, length(S_Y19)), rep(2017, length(S_Y20)),
                               rep(2018, length(S_Y21)),  rep(2019, length(S_Y22)),rep(2020, length(S_Y23))))

dat.iW3 = c(predYRW3_Y1, predYRW3_Y2, predYRW3_Y3, predYRW3_Y4, predYRW3_Y5, predYRW3_Y6, predYRW3_Y7, predYRW3_Y8, predYRW3_Y9, predYRW3_Y10, predYRW3_Y11, 
              predYRW3_Y12, predYRW3_Y13, predYRW3_Y14, predYRW3_Y15, predYRW3_Y16, predYRW3_Y17, predYRW3_Y18, predYRW3_Y19, predYRW3_Y20, predYRW3_Y21, predYRW3_Y22, predYRW3_Y23)
dat.iW3 = data.frame(S = dat.iW3, 
                     Year = c(rep(1998, length(S_Y1)), rep(1999, length(S_Y2)), rep(2000, length(S_Y3)), rep(2001, length(S_Y4)),
                               rep(2002, length(S_Y5)),  rep(2003, length(S_Y6)),rep(2004, length(S_Y7)), rep(2005, length(S_Y8)),
                               rep(2006, length(S_Y9)),  rep(2007, length(S_Y10)),rep(2008, length(S_Y11)), rep(2009, length(S_Y12)),
                               rep(2010, length(S_Y13)),  rep(2011, length(S_Y14)),rep(2012, length(S_Y15)), rep(2013, length(S_Y16)),
                               rep(2014, length(S_Y17)),  rep(2015, length(S_Y18)),rep(2016, length(S_Y19)), rep(2017, length(S_Y20)),
                               rep(2018, length(S_Y21)),  rep(2019, length(S_Y22)),rep(2020, length(S_Y23))))



PArich = data.frame(RichH1=c(dat.iW1$S), RichH2=c(dat.iW2$S), RichH4=c(dat.iW3$S),
                    Richall=c(dat.iALL$S), SiteID = rep(SiteFIn_un, 23),
                    Year = c(rep(1998, length(S_Y1)), rep(1999, length(S_Y2)), rep(2000, length(S_Y3)), rep(2001, length(S_Y4)),
                               rep(2002, length(S_Y5)),  rep(2003, length(S_Y6)),rep(2004, length(S_Y7)), rep(2005, length(S_Y8)),
                               rep(2006, length(S_Y9)),  rep(2007, length(S_Y10)),rep(2008, length(S_Y11)), rep(2009, length(S_Y12)),
                               rep(2010, length(S_Y13)),  rep(2011, length(S_Y14)),rep(2012, length(S_Y15)), rep(2013, length(S_Y16)),
                               rep(2014, length(S_Y17)),  rep(2015, length(S_Y18)),rep(2016, length(S_Y19)), rep(2017, length(S_Y20)),
                               rep(2018, length(S_Y21)),  rep(2019, length(S_Y22)),rep(2020, length(S_Y23))))


summary(PArich)
head(PArich)

#PArich$Year = as.factor(PArich$Year)
PArich$Year = as.integer(PArich$Year)


Trend_site_FIN = Resid_site_FIN = Resid2_site_FIN = coef_trend_FIN = matrix(NA, nrow= 12340, ncol=4)
for (i in 1:length(SiteFIn_un)) {
  MapdataALL_i = PArich[PArich$SiteID == SiteFIn_un[i],]
  
  for (j in colnames(MapdataALL_i)[1:4]) {
    index_j = match(j, colnames(MapdataALL_i)[1:4])
    Fit_site_FIN = lm(reformulate('Year', j), data = MapdataALL_i)
    
    Resid_site_FIN[i, index_j] = var(resid(Fit_site_FIN))
    Resid2_site_FIN[i, index_j] = var(rstandard(Fit_site_FIN), na.rm=TRUE)
    coef_trend_FIN[i, index_j] = Fit_site_FIN$coefficients[2]
  }
}

Resid_site_FIN = as.data.frame(Resid_site_FIN)
colnames(Resid_site_FIN)[1:4] = colnames(MapdataALL_i)[1:4]
Resid_site_FIN$X = Fin_coord$X
Resid_site_FIN$Y = Fin_coord$Y
Resid_site_FIN$spatial = 'FINLAND'

Resid2_site_FIN = as.data.frame(Resid2_site_FIN)
colnames(Resid2_site_FIN)[1:4] = colnames(MapdataALL_i)[1:4]
Resid2_site_FIN$X = Fin_coord$X
Resid2_site_FIN$Y = Fin_coord$Y
Resid2_site_FIN$spatial = 'FINLAND'

coef_trend_FIN = as.data.frame(coef_trend_FIN)
colnames(coef_trend_FIN)[1:4] = colnames(MapdataALL_i)[1:4]
coef_trend_FIN$X = Fin_coord$X
coef_trend_FIN$Y = Fin_coord$Y
coef_trend_FIN$spatial = 'FINLAND'

## 

save(PArich, Resid_site_FIN, Resid2_site_FIN, coef_trend_FIN,
     file='PAtrendWL_FIN_2808.RDATA')




##------------------ HP - trait

dat.iH1 = c(predYRH1_Y1, predYRH1_Y2, predYRH1_Y3, predYRH1_Y4, predYRH1_Y5, predYRH1_Y6, predYRH1_Y7, predYRH1_Y8, predYRH1_Y9, predYRH1_Y10, predYRH1_Y11, 
              predYRH1_Y12, predYRH1_Y13, predYRH1_Y14, predYRH1_Y15, predYRH1_Y16, predYRH1_Y17, predYRH1_Y18, predYRH1_Y19, predYRH1_Y20, predYRH1_Y21, predYRH1_Y22, predYRH1_Y23)
dat.iH1 = data.frame(S = dat.iH1, 
                     Year = c(rep(1998, length(S_Y1)), rep(1999, length(S_Y2)), rep(2000, length(S_Y3)), rep(2001, length(S_Y4)),
                               rep(2002, length(S_Y5)),  rep(2003, length(S_Y6)),rep(2004, length(S_Y7)), rep(2005, length(S_Y8)),
                               rep(2006, length(S_Y9)),  rep(2007, length(S_Y10)),rep(2008, length(S_Y11)), rep(2009, length(S_Y12)),
                               rep(2010, length(S_Y13)),  rep(2011, length(S_Y14)),rep(2012, length(S_Y15)), rep(2013, length(S_Y16)),
                               rep(2014, length(S_Y17)),  rep(2015, length(S_Y18)),rep(2016, length(S_Y19)), rep(2017, length(S_Y20)),
                               rep(2018, length(S_Y21)),  rep(2019, length(S_Y22)),rep(2020, length(S_Y23))))



dat.iH2 = c(predYRH2_Y1, predYRH2_Y2, predYRH2_Y3, predYRH2_Y4, predYRH2_Y5, predYRH2_Y6, predYRH2_Y7, predYRH2_Y8, predYRH2_Y9, predYRH2_Y10, predYRH2_Y11, 
              predYRH2_Y12, predYRH2_Y13, predYRH2_Y14, predYRH2_Y15, predYRH2_Y16, predYRH2_Y17, predYRH2_Y18, predYRH2_Y19, predYRH2_Y20, predYRH2_Y21, predYRH2_Y22, predYRH2_Y23)
dat.iH2 = data.frame(S = dat.iH2, 
                     Year = c(rep(1998, length(S_Y1)), rep(1999, length(S_Y2)), rep(2000, length(S_Y3)), rep(2001, length(S_Y4)),
                               rep(2002, length(S_Y5)),  rep(2003, length(S_Y6)),rep(2004, length(S_Y7)), rep(2005, length(S_Y8)),
                               rep(2006, length(S_Y9)),  rep(2007, length(S_Y10)),rep(2008, length(S_Y11)), rep(2009, length(S_Y12)),
                               rep(2010, length(S_Y13)),  rep(2011, length(S_Y14)),rep(2012, length(S_Y15)), rep(2013, length(S_Y16)),
                               rep(2014, length(S_Y17)),  rep(2015, length(S_Y18)),rep(2016, length(S_Y19)), rep(2017, length(S_Y20)),
                               rep(2018, length(S_Y21)),  rep(2019, length(S_Y22)),rep(2020, length(S_Y23))))


dat.iH4 = c(predYRH4_Y1, predYRH4_Y2, predYRH4_Y3, predYRH4_Y4, predYRH4_Y5, predYRH4_Y6, predYRH4_Y7, predYRH4_Y8, predYRH4_Y9, predYRH4_Y10, predYRH4_Y11, 
              predYRH4_Y12, predYRH4_Y13, predYRH4_Y14, predYRH4_Y15, predYRH4_Y16, predYRH4_Y17, predYRH4_Y18, predYRH4_Y19, predYRH4_Y20, predYRH4_Y21, predYRH4_Y22, predYRH4_Y23)
dat.iH4 = data.frame(S = dat.iH4, 
                     Year = c(rep(1998, length(S_Y1)), rep(1999, length(S_Y2)), rep(2000, length(S_Y3)), rep(2001, length(S_Y4)),
                               rep(2002, length(S_Y5)),  rep(2003, length(S_Y6)),rep(2004, length(S_Y7)), rep(2005, length(S_Y8)),
                               rep(2006, length(S_Y9)),  rep(2007, length(S_Y10)),rep(2008, length(S_Y11)), rep(2009, length(S_Y12)),
                               rep(2010, length(S_Y13)),  rep(2011, length(S_Y14)),rep(2012, length(S_Y15)), rep(2013, length(S_Y16)),
                               rep(2014, length(S_Y17)),  rep(2015, length(S_Y18)),rep(2016, length(S_Y19)), rep(2017, length(S_Y20)),
                               rep(2018, length(S_Y21)),  rep(2019, length(S_Y22)),rep(2020, length(S_Y23))))



dat.iH31 = c(predYRH31_Y1, predYRH31_Y2, predYRH31_Y3, predYRH31_Y4, predYRH31_Y5, predYRH31_Y6, predYRH31_Y7, predYRH31_Y8, predYRH31_Y9, predYRH31_Y10, predYRH31_Y11, 
              predYRH31_Y12, predYRH31_Y13, predYRH31_Y14, predYRH31_Y15, predYRH31_Y16, predYRH31_Y17, predYRH31_Y18, predYRH31_Y19, predYRH31_Y20, predYRH31_Y21, predYRH31_Y22, predYRH31_Y23)
dat.iH31 = data.frame(S = dat.iH31, 
                     Year = c(rep(1998, length(S_Y1)), rep(1999, length(S_Y2)), rep(2000, length(S_Y3)), rep(2001, length(S_Y4)),
                               rep(2002, length(S_Y5)),  rep(2003, length(S_Y6)),rep(2004, length(S_Y7)), rep(2005, length(S_Y8)),
                               rep(2006, length(S_Y9)),  rep(2007, length(S_Y10)),rep(2008, length(S_Y11)), rep(2009, length(S_Y12)),
                               rep(2010, length(S_Y13)),  rep(2011, length(S_Y14)),rep(2012, length(S_Y15)), rep(2013, length(S_Y16)),
                               rep(2014, length(S_Y17)),  rep(2015, length(S_Y18)),rep(2016, length(S_Y19)), rep(2017, length(S_Y20)),
                               rep(2018, length(S_Y21)),  rep(2019, length(S_Y22)),rep(2020, length(S_Y23))))



dat.iH32 = c(predYRH32_Y1, predYRH32_Y2, predYRH32_Y3, predYRH32_Y4, predYRH32_Y5, predYRH32_Y6, predYRH32_Y7, predYRH32_Y8, predYRH32_Y9, predYRH32_Y10, predYRH32_Y11, 
              predYRH32_Y12, predYRH32_Y13, predYRH32_Y14, predYRH32_Y15, predYRH32_Y16, predYRH32_Y17, predYRH32_Y18, predYRH32_Y19, predYRH32_Y20, predYRH32_Y21, predYRH32_Y22, predYRH32_Y23)
dat.iH32 = data.frame(S = dat.iH32, 
                     Year = c(rep(1998, length(S_Y1)), rep(1999, length(S_Y2)), rep(2000, length(S_Y3)), rep(2001, length(S_Y4)),
                               rep(2002, length(S_Y5)),  rep(2003, length(S_Y6)),rep(2004, length(S_Y7)), rep(2005, length(S_Y8)),
                               rep(2006, length(S_Y9)),  rep(2007, length(S_Y10)),rep(2008, length(S_Y11)), rep(2009, length(S_Y12)),
                               rep(2010, length(S_Y13)),  rep(2011, length(S_Y14)),rep(2012, length(S_Y15)), rep(2013, length(S_Y16)),
                               rep(2014, length(S_Y17)),  rep(2015, length(S_Y18)),rep(2016, length(S_Y19)), rep(2017, length(S_Y20)),
                               rep(2018, length(S_Y21)),  rep(2019, length(S_Y22)),rep(2020, length(S_Y23))))



dat.iH33 = c(predYRH33_Y1, predYRH33_Y2, predYRH33_Y3, predYRH33_Y4, predYRH33_Y5, predYRH33_Y6, predYRH33_Y7, predYRH33_Y8, predYRH33_Y9, predYRH33_Y10, predYRH33_Y11, 
              predYRH33_Y12, predYRH33_Y13, predYRH33_Y14, predYRH33_Y15, predYRH33_Y16, predYRH33_Y17, predYRH33_Y18, predYRH33_Y19, predYRH33_Y20, predYRH33_Y21, predYRH33_Y22, predYRH33_Y23)
dat.iH33 = data.frame(S = dat.iH33, 
                     Year = c(rep(1998, length(S_Y1)), rep(1999, length(S_Y2)), rep(2000, length(S_Y3)), rep(2001, length(S_Y4)),
                               rep(2002, length(S_Y5)),  rep(2003, length(S_Y6)),rep(2004, length(S_Y7)), rep(2005, length(S_Y8)),
                               rep(2006, length(S_Y9)),  rep(2007, length(S_Y10)),rep(2008, length(S_Y11)), rep(2009, length(S_Y12)),
                               rep(2010, length(S_Y13)),  rep(2011, length(S_Y14)),rep(2012, length(S_Y15)), rep(2013, length(S_Y16)),
                               rep(2014, length(S_Y17)),  rep(2015, length(S_Y18)),rep(2016, length(S_Y19)), rep(2017, length(S_Y20)),
                               rep(2018, length(S_Y21)),  rep(2019, length(S_Y22)),rep(2020, length(S_Y23))))



PArich = data.frame(RichH1=c(dat.iH1$S), RichH2=c(dat.iH2$S), RichH4=c(dat.iH4$S),
                    RichH31=c(dat.iH31$S), RichH32=c(dat.iH32$S), RichH33=c(dat.iH33$S),
                    Richall=c(dat.iALL$S), SiteID = rep(SiteFIn_un, 23),
                    Year = c(rep(1998, length(S_Y1)), rep(1999, length(S_Y2)), rep(2000, length(S_Y3)), rep(2001, length(S_Y4)),
                               rep(2002, length(S_Y5)),  rep(2003, length(S_Y6)),rep(2004, length(S_Y7)), rep(2005, length(S_Y8)),
                               rep(2006, length(S_Y9)),  rep(2007, length(S_Y10)),rep(2008, length(S_Y11)), rep(2009, length(S_Y12)),
                               rep(2010, length(S_Y13)),  rep(2011, length(S_Y14)),rep(2012, length(S_Y15)), rep(2013, length(S_Y16)),
                               rep(2014, length(S_Y17)),  rep(2015, length(S_Y18)),rep(2016, length(S_Y19)), rep(2017, length(S_Y20)),
                               rep(2018, length(S_Y21)),  rep(2019, length(S_Y22)),rep(2020, length(S_Y23))))


summary(PArich)
head(PArich)
PArich$Year = as.factor(PArich$Year)
PArich$Year = as.integer(PArich$Year)

Trend_site_FIN = Resid_site_FIN = Resid2_site_FIN = coef_trend_FIN = matrix(NA, nrow= 12340, ncol=11)
for (i in 1:length(SiteFIn_un)) {
  MapdataALL_i = PArich[PArich$SiteID == SiteFIn_un[i],]
  
  for (j in colnames(MapdataALL_i)[1:7]) {
    index_j = match(j, colnames(MapdataALL_i)[1:7])
    Fit_site_FIN = lm(reformulate('Year', j), data = MapdataALL_i)
    
    Resid_site_FIN[i, index_j] = var(resid(Fit_site_FIN))
    Resid2_site_FIN[i, index_j] = var(rstandard(Fit_site_FIN), na.rm=TRUE)
    coef_trend_FIN[i, index_j] = Fit_site_FIN$coefficients[2]
  }
}

Resid_site_FIN = as.data.frame(Resid_site_FIN)
colnames(Resid_site_FIN)[1:7] = colnames(MapdataALL_i)[1:7]
Resid_site_FIN$X = Fin_coord$X
Resid_site_FIN$Y = Fin_coord$Y
Resid_site_FIN$spatial = 'FINLAND'

Resid2_site_FIN = as.data.frame(Resid2_site_FIN)
colnames(Resid2_site_FIN)[1:7] = colnames(MapdataALL_i)[1:7]
Resid2_site_FIN$X = Fin_coord$X
Resid2_site_FIN$Y = Fin_coord$Y
Resid2_site_FIN$spatial = 'FINLAND'

coef_trend_FIN = as.data.frame(coef_trend_FIN)
colnames(coef_trend_FIN)[1:7] = colnames(MapdataALL_i)[1:7]
coef_trend_FIN$X = Fin_coord$X
coef_trend_FIN$Y = Fin_coord$Y
coef_trend_FIN$spatial = 'FINLAND'

## 

save(PArich, Resid_site_FIN, Resid2_site_FIN, coef_trend_FIN,
     file='PAtrendHP_FIN_2808.RDATA')


