setwd("E:/Helsinki/ResultsAug23/RDATA/3.ModelRun")
load('test10908cent_yearPA.RDATA')
modelyear = list()
modelyear[[1]] = models[[1]]
load('test10908cent_yearAB.RDATA')
modelyear[[2]] = models[[2]]

#load('test1910_yearAB.RDATA')


{
  library(Hmsc)
  library(tidyr)
  library(dplyr)
  library(RColorBrewer)
  library(viridis)
  library(abind)
  library(grid)
  library(gridExtra)
  library(ggplot2)
  library(ggpubr)
  library(scales)
}




## Observe Habitat information from the model
XDataOld = modelyear[[1]]$XData
head(XDataOld)
dim(XDataOld)

modelyear[[1]]$XFormula

#### change in GDDwinter
ClimpropOLD = XDataOld[,1:2]

dim(ClimpropOLD)


# Moth Site info
#setwd("C:/Users/guilbaul/OneDrive/DataREC/Moth")

# Moth data to get site name
setwd("E:/Helsinki/ResultsAug23/RDATA/2.ObsInfo")
Moth <- read.csv("Copy_mothMonitoring.csv", header=T)
head(Moth)
dim(Moth)

M2 = Moth[Moth$Species %in% modelyear[[1]]$spNames,]
dim(M2)
# if start after a certain date
M3 = Moth[Moth$YearStart >= 1998,]
dim(M3)

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

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

Site = Moth.xy$SiteID

Datacurrent = data.frame(XDataOld, modelyear[[1]]$studyDesign)
head(Datacurrent)


# current value
setwd("E:/Helsinki/ResultsAug23/RDATA/1.EnvironmentalInfo")

## prediction over all Finland

load('CLimFingrid1408.RDATA')  # large Finland
head(Fin_dat_ALL)
dim(Fin_dat_ALL) ## all sites all years

load('HabFingrid_1008.RDATA') # Large FIN
head(Fin_dat_ALL)
dim(Fin_dat_ALL) ## all sites one year


load('ClimSmallFIN_1408.RDATA')  # small Finland
head(mapData2)
dim(mapData2) ## all sites all years

mapData2 = mapData2[mapData2$Year < 2021,]

load('HabSmallFIN1408.RDATA') # small Fin
head(mapData2hab)
dim(mapData2hab) ## all sites one year


dim(mapData2)  # clim 283820
dim(mapData2hab)  # hab  283820


save(mapData2, mapData2hab, file = 'CurrentSmallFIN_1408.RDATA')
head(mapData2)
head(mapData2hab)

Fin_all_small = cbind(mapData2[, c(2, 6,15, 10:11)], mapData2hab[, c(1:8)])
head(Fin_all_small)

Fin_all_small_sd = mapData2[, c(1, 16,17)]
dim(Fin_all_small_sd)

#######
# DEFINING THE NEW STUDY DESIGN

sample.id = paste(Fin_all_small_sd$X, '_', Fin_all_small_sd$Y)
year = Fin_all_small_sd$Year

load('Region_grid_BR.RDATA')

region = as.factor(rep(Region_corr, length(unique(year)))) # 283820

setwd("E:/Helsinki/ResultsAug23/RDATA/4.Predictions")
save(modelyear, Fin_all_small, Fin_all_small_sd, sample.id, year, region, file = 'Predsmall2408.RDATA')


