# get relevant functions and information files

source("data work up functions.r")

ctsm.env <- attach(NULL, name = "ctsm")
function.path <- file.path("..", "..", "functions")
sys.source(file.path(function.path, "xml functions 1_9.r"), ctsm.env)


# can't read info files directly, because some of them require functions in base R to run

wk.obj <- objects()
source(file.path(function.path, "information functions 1_9_3.r"))
ctsm.env$get.info <- get.info
ctsm.env$info.determinand <- info.determinand
ctsm.env$info.species <- info.species
ctsm.env$get.basis <- get.basis
ctsm.env$convert.basis <- convert.basis
remove(list = setdiff(objects(), wk.obj))
rm(function.path)



# get meta data by group

sediment <- getMetaData(
  "sediment", 
  c("Metals", "PAH (parent)", "Chlorobiphenyls", "Organo-bromines", "Organo-metals"), 
  byGroup = TRUE)


biota <- getMetaData(
  "biota", 
  c("Metals", "PAH (parent)", "Chlorobiphenyls", "Organo-bromines", "Imposex"), 
  byGroup = TRUE)


# get relevant status and trend data - strip off timeseries with insufficient years and, for trend, 
# impacted and baseline stations

sediment <- lapply(sediment, getData)
biota <- lapply(biota, getData)


# reduce to regions with at least three stations 

sediment <- lapply(sediment, regionalData, "initialise")
biota <- lapply(biota, regionalData, "initialise")


# make manual adjustments 

# sediment Iberian Sea - all impacted and grouped in one small area, so insufficient geographic coverage 

sediment <- lapply(sediment, function(x) {
  within(x, regionalStatus <- droplevels(subset(regionalStatus, ! region %in% "Iberian Sea")))
})

# biota Norwegian Sea - sometimes four stations, but clustered in two pairs of two

id <- c("Metals", "PAH (parent)", "Chlorobiphenyls", "Organo-bromines")
biota[id] <- lapply(biota[id], function(x) {
  within(x, {
    regionalStatus <- droplevels(subset(regionalStatus, ! region %in% "Norwegian Sea"))
    regionalTrend <- droplevels(subset(regionalTrend, ! region %in% "Norwegian Sea"))
  })
})


# sediment$Chlorobiphenyls <- within(sediment$Chlorobiphenyls, {
#   regionalStatus <- droplevels(subset(regionalStatus, ! region %in% "Channel"))
#   regionalTrend <- droplevels(subset(regionalTrend, ! region %in% c("Channel", "Gulf of Cadiz")))
# })
# 
# sediment$PAHs <- within(sediment$PAHs, regionalTrend <- 
#                           droplevels(subset(regionalTrend, ! region %in% "Channel")))
# 
# sediment$BDEs <- within(sediment$BDEs, {
#   regionalStatus <- droplevels(subset(regionalStatus, ! region %in% "Irish Sea"))
#   rm(regionalTrend)
# })


# and strip out any further determinands that have less than three timeseries

sediment <- lapply(sediment, regionalData, "finalise")
biota <- lapply(biota, regionalData, "finalise")

save(sediment, biota, file = file.path("meta analysis objects", "regional data.RData"))
# load(file.path("meta analysis objects", "regional data.RData"))


# write files to respective folders for data assessors

lapply(sediment, writeData)
lapply(biota, writeData)



require(rmarkdown)

lapply(names(sediment)[1:4], function(x) {
  outfile <- paste0("regional_assessment_sediment_", x, ".html")
  outfile <- gsub(" ", "_", outfile, fixed = TRUE)
  render("meta analysis test.Rmd", output_file = outfile, output_dir = "output", 
         params = list(compartment = "sediment", detGroup = x)
  )
})

render("meta analysis tins.Rmd", 
       output_file = "regional_assessment_sediment_Organo-metals.html",  
       output_dir = "output", 
       params = list(compartment = "sediment", detGroup = "Organo-metals"))



lapply(names(biota)[1], function(x) {
  render("meta analysis test.Rmd", 
         output_file = paste0("regional_assessment_biota_", x, ".html"),  
         output_dir = "output", 
         params = list(compartment = "biota", detGroup = x)
  )
})

render("meta analysis imposex.Rmd", 
       output_file = "regional_assessment_biota_Imposex.html",  
       output_dir = "output", 
       params = list(compartment = "biota", detGroup = "Imposex"))
