# 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_1.r"))
ctsm.env$get.info <- get.info
ctsm.env$info.determinand <- info.determinand
remove(list = setdiff(objects(), wk.obj))
rm(function.path)


# get meta data by compartment to look at MSTAT values

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

biota <- getMetaData(
  "biota", 
  c("Metals", "PAH (parent)", "Chlorobiphenyls", "Organo-bromines"))

save(sediment, file = file.path("meta analysis objects", "MSTAT investigation.RData"))
# load(file.path("meta analysis objects", "MSTAT investigation.RData"))

require(rmarkdown)

lapply(c("sediment", "biota"), function(x) {
  render("MSTAT summary.Rmd", 
         output_file = paste0("MSTAT_screening_", x, ".html"),  
         output_dir = "output", 
         params = list(compartment = x)
  )
})

lapply(c("sediment", "biota"), function(x) {
  data <- get(x) 
  data <- subset(data$timeSeries, paramFit)[c(
    "OSPARregion", "region", "country", "station", "latitude", "longitude", "MSTATsubmitted", "MSTATedited", 
    "MSTAT")]
  data <- unique(data) 
  row.names(data) <- NULL
  data <- data[with(data, order(OSPARregion, region, country, MSTAT)), ]
  write.csv(data, file = file.path("output", paste0("MSTAT ", x, ".csv")), row.names = FALSE, na = "")
})

rm(sediment, biota)


# 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"), 
  byGroup = TRUE)


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

sediment <- lapply(sediment, getData, dropByMstat = TRUE)
biota <- lapply(biota, getData, dropByMstat = TRUE)


# now reduce to regions with at least three stations 

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


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


mapply(biota, names(biota), MoreArgs = list(compartment = "biota"), FUN = function(data, id, compartment) {
  trend <- data$trend
  trend$file <- do.call(paste, c(data$timeSeries[row.names(trend), c("filePrefix", "fileSuffix")], sep = "_"))
  trend <- within(trend, {
    assessmentPlot <- paste0("http://dome.ices.dk/OSPARMIME2016/graphics/", file, "da.jpg")
    auxiliaryPlot <- paste0("http://dome.ices.dk/OSPARMIME2016/graphics/", file, "au.jpg")
    rm("file")
  })
  xlsx::write.xlsx(trend, file.path("output", "data screening", paste(compartment, id, "trend.xlsx")))

  status <- data$status
  status$file <- do.call(
    paste, c(data$timeSeries[row.names(status), c("filePrefix", "fileSuffix")], sep = "_"))
  status <- within(status, {
    assessmentPlot <- paste0("http://dome.ices.dk/OSPARMIME2016/graphics/", file, "da.jpg")
    auxiliaryPlot <- paste0("http://dome.ices.dk/OSPARMIME2016/graphics/", file, "au.jpg")
    rm("file")
  })
  xlsx::write.xlsx(status, file.path("output", "data screening", paste(compartment, id, "status.xlsx")))
})


require(rmarkdown)

lapply(c("sediment", "biota"), function(x) {
  render("individual time series screening.Rmd", 
         output_file = paste0("time_series_screening_", x, ".html"),  
         output_dir = "output", 
         params = list(compartment = x)
  )
})




# and make manual adjustments

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)
})


# strip out any determinands that have less than three timeseries

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



# write files to respective folders for data assessors

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



require(rmarkdown)

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

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


wk <- sediment$Chlorobiphenyls$regionalStatus
with(wk, table(region, determinand))


