readFromGenstat <- function(x, info, z = 1.96) {
  inFile <- file.path("meta analysis results",  paste0(info$ID, " ", x, ".csv"))
  out <- read.csv(inFile, na = "*", strip.white = TRUE)
  names(out) <- sub("region_", "", names(out), fixed = TRUE)
  
  out <- out[!is.na(out$est), ]
  
  within(out, {
    lower <- est - z * se
    upper <- est + z * se
    if ("region" %in% names(out)) region <- droplevels(factor(region, levels = info$regionID))
    if ("determinand" %in% names(out)) determinand <- droplevels(factor(determinand, levels = info$detID))
  })
}


plotTrendResults <- function(
  assessment, type = c("region", "determinand", "region by determinand", "determinand by region")) {
  
  require(lattice)
  
  if (is.null(assessment$regionalTrend)) return(invisible())

  type <- match.arg(type)
  
  # ensure plotting doesn't include combinations that don't exist
  
  regDetExists <- with(assessment$regionalTrend, unique(paste(region, determinand)))
  
  
  plotEngine <- function(data, formula) {
    xlim <- extendrange(range(data$upper, data$lower, 0))
    labels <- pretty(100 * (exp(xlim / 100) - 1), n = 5)
    at <- 100 * (log(as.numeric(labels) / 100 + 1))
    stripplot(
      formula,
      data = data, 
      xlab = "% yearly change in concentration",
      xlim = xlim, 
      scales = list(alternating = FALSE, x = list(at = at, labels = labels)),
      panel = function(x, y, subscripts) {
        panel.abline(v = 0, col = grey(0.7))
        with(data[subscripts, ], lsegments(lower, y, upper, y, col = "black"))
        col <- with(data[subscripts, ], ifelse(upper < 0, "green3", ifelse(lower > 0, "red", "goldenrod")))
        lpoints(x, y, col = col, pch = 16, cex = 1.5)
      })
  }
  
  switch(
    type, 
    region = {
      data <- readFromGenstat("trend results by region", assessment)
      plotEngine(data, as.formula(region ~ est))
    },
    determinand = {
      data <- readFromGenstat("trend results by determinand", assessment)
      plotEngine(data, as.formula(determinand ~ est))
    }, 
    "region by determinand" = {
      data <- readFromGenstat("trend results by region determinand", assessment)
      data <- subset(data, paste(region, determinand) %in% regDetExists)
      plotEngine(data, as.formula(region ~ est | determinand))
    },
    "determinand by region" = {
      data <- readFromGenstat("trend results by region determinand", assessment)
      data <- subset(data, paste(region, determinand) %in% regDetExists)
      plotEngine(data, as.formula(determinand ~ est | region))
    }
  )
}



plotStatusResults <- function(
  assessment, 
  type = c("region", "determinand", "region by determinand", "determinand by region"), 
  AC) {

  require(lattice)
  
  if (is.null(assessment$regionalStatus)) return(invisible())
  
  type <- match.arg(type)
  
  if (!(AC %in% c("CONC", assessment$AC))) stop("AC not recognised")

  # ensure plotting doesn't include combinations that don't exist
  
  originalData <- assessment$regionalStatus
  ok <- with(originalData, !is.na(switch(AC, "CONC" = status, get(paste0("status.", AC)))))
  originalData <- originalData[ok, ]
  regDetExists <- with(originalData, unique(paste(region, determinand)))
  

  plotEngine <- function(data, formula, AC) {
    xlim <- switch(AC, 
                   CONC = extendrange(range(data$upper, data$lower)),
                   extendrange(range(data$upper, data$lower, 0)))
    labels <- plot.scales(range(xlim, na.rm = TRUE), logData = TRUE, n = 5)
    xlab <- switch(AC, 
                   CONC = "mean concentration in final year", 
                   paste("mean concentration relative to", AC))
    stripplot(
      formula,
      data = data, 
      xlab = xlab,
      xlim = xlim, 
      scales = list(alternating = FALSE, x = list(at = log(labels), labels = labels)),
      panel = function(x, y, subscripts) {
        if (AC != "CONC") panel.abline(v = 0, col = grey(0.7))
        with(data[subscripts, ], lsegments(lower, y, upper, y, col = "black"))
        col <- switch(AC, 
                      CONC = "black", 
                      BAC = with(data[subscripts, ], ifelse(upper < 0, "blue", "green3")),
                      with(data[subscripts, ], ifelse(upper < 0, "green3", "red")))
        lpoints(x, y, col = col, pch = 16, cex = 1.5)
      })
  }

  txt <- paste0("status results", switch(AC, CONC = "", paste0(" ", AC)))
  z <- switch(AC, CONC = 1.96, 1.645)
    
  switch(
    type, 
    region = {
      txt2 <- paste(txt, "by region")
      data <- readFromGenstat(txt2, assessment, z = z)
      plotEngine(data, as.formula(region ~ est), AC)
    },  
    determinand = {
      txt2 <- paste(txt, "by determinand")
      data <- readFromGenstat(txt2, assessment, z = z)
      plotEngine(data, as.formula(determinand ~ est), AC)
    },  
    "region by determinand" = {
      txt2 <- paste(txt, "by region determinand")
      data <- readFromGenstat(txt2, assessment, z = z)
      data <- subset(data, paste(region, determinand) %in% regDetExists)
      plotEngine(data, as.formula(region ~ est | determinand), AC)
    },
    "determinand by region" = {
      txt2 <- paste(txt, "by region determinand")
      data <- readFromGenstat(txt2, assessment, z = z)
      data <- subset(data, paste(region, determinand) %in% regDetExists)
      plotEngine(data, as.formula(determinand ~ est | region), AC)
    }
  )
  
}



getData <- function(x, dropByMstat = FALSE) {

  AC <- x$AC
  assessment <- x$assessment
  timeSeries <- x$timeSeries


  # get status stations
  
  ok <- with(timeSeries, paramFit)
  if (!any(ok)) return(x)
  
  timeSeries <- droplevels(timeSeries[ok, ])
  assessment <- assessment[row.names(timeSeries)]
  
  Status <- sapply(assessment, USE.NAMES = TRUE, simplify = FALSE, FUN = function(assessment) 
    cbind(tail(assessment$pred, 1), assessment$summary[AC])
  )

  Status <- do.call(rbind, Status)
  
  id <- intersect(c("determinand", "country", "station", "species", "OSPARregion", "MSFDregion", 
                    "region", "offshore", "MSTAT", "latitude", "longitude"), names(timeSeries))
  Status[id] <- timeSeries[row.names(Status), id]
  
  if (!is.null(AC)) Status[paste("status", AC, sep = ".")] <- Status$fit - log(Status[AC])

  Status <- within(Status, {
    status <- fit
    se.status <- se
    station <- factor(station)
    determinand <- factor(as.character(determinand), levels = x$detID)
  })

  if (dropByMstat) 
    Status <- Status[Status$MSTAT %in% c("B", "RH"), ]

  x$status <- droplevels(Status)
  
  
  # now look at only time series with trend information (given by nypos) 
  
  ok <- with(timeSeries, nypos >= 5)
  if (!any(ok)) return(x)
  
  timeSeries <- droplevels(timeSeries[ok, ])
  assessment <- assessment[row.names(timeSeries)]
  
  Trend <- sapply(assessment, USE.NAMES = TRUE, simplify = FALSE, function(assessment)
    assessment$contrasts["recent", ]
  )
  Trend <- do.call(rbind, Trend)
    
  Trend[id] <- timeSeries[row.names(Trend), id]

  Trend <- within(Trend, {
    trend <- 100 * estimate / (end - start)
    se.trend <- 100 * se / (end - start)
    station <- factor(station)
    determinand <- factor(as.character(determinand), levels = x$detID)
  })

  if (dropByMstat) 
    Trend <- Trend[Trend$MSTAT %in% "RH", ]

  x$trend <- droplevels(Trend)
  
  x
}


regionalData <- function(assessment) {
  
  isMetals <- strsplit(assessment$ID, " ")[[1]][2] == "Metals"
  
  assessment[c("regionalTrend", "regionalStatus")] <- lapply(c("trend", "status"), function(id) {
    
    if (is.null(assessment[[id]])) return(NULL)
    
    data <- assessment[[id]]
    
    data <- droplevels(subset(data, MSTAT %in% switch(id, trend = "RH", status = c("B", "RH"))))
  
    # must have at least three unique sampling locations in each region
    # NB can't use station because (conceivably) different countries could monitor at the same spot
    
    coverage <- with(unique(data[c("latitude", "longitude", "region")]), table(region))
    ok <- names(coverage)[coverage >= 3]
    data <- droplevels(data[data$region %in% ok, ])

    # must have at least three unique sampling locations for each determinand over all regions
        
    coverage <- with(unique(data[c("latitude", "longitude", "determinand")]), table(determinand))
    ok <- names(coverage)[coverage >= 3]
    data <- droplevels(data[data$determinand %in% ok, ])
    
    stopifnot(with(data, table(region) >= 3), with(data, table(determinand) >= 3))  

    # extra reduction for metals, because we don't analyse them as a random effect (assume 
    # they arise from different processes)
        
    if (isMetals) {
      det.reg <- with(unique(data[c("latitude", "longitude", "determinand", "region")]), 
                      paste(determinand, region))
      coverage <- table(det.reg)
      ok <- names(coverage)[coverage >= 3]
      data <- droplevels(data[with(data, paste(determinand, region)) %in% ok, ])
    }

    data
  })
  
  assessment
}


finaliseRegionalData <- function(assessment) {
  
  assessment[c("regionalTrend", "regionalStatus")] <- 
    lapply(c("regionalTrend", "regionalStatus"), function(id) {
    
      if (is.null(assessment[[id]])) return(NULL)
      
      data <- assessment[[id]]
    
      data <- data[data$MSTAT %in% switch(id, regionalTrend = "RH", regionalStatus = c("B", "RH")), ]
      data <- droplevels(data)
      
      coverage <- with(unique(data[c("station", "region")]), table(region))
      ok <- names(coverage)[coverage >= 3]
      data <- droplevels(subset(data, region %in% ok))
      
      coverage <- with(data, table(determinand))
      ok <- names(coverage)[coverage >= 3]
      data <- droplevels(subset(data, determinand %in% ok))
      
      stopifnot(with(data, table(region) >= 3), with(data, table(determinand) >= 3))  
      
      data
    })
  
  assessment
}



plotTrend <- function(assessment) {
  require(lattice)
  path <- assessment$webpath
  if (is.null(assessment$trend)) return(invisible())
  data <- assessment$trend
#  png(file.path(path, paste("trend stripplot.png")), height = 16, width = 16, units = 'cm', res = 150)
  labels <- range(data$trend, na.rm = TRUE)
  labels <- pretty(100 * (exp(labels / 100) - 1), n = 3)
  at <- 100 * (log(as.numeric(labels) / 100 + 1))
  out <- stripplot(
    region ~ trend | determinand,
    data = data, 
    xlab = "% yearly change in concentration",
    scales = list(alternating = FALSE, x = list(at = at, labels = labels)),
    panel = function(x, y, subscripts) {
      col <- as.character(data$MSTAT[subscripts])
      col <- as.character(factor(col, levels = c("B", "RH", "IH"), labels = c("blue", "green3", "red")))
      panel.abline(v = 0, col = grey(0.7))
      panel.stripplot(x, y, jitter.data = TRUE, pch = 16, col = col)
    })
  out
#  dev.off()
#  invisible()
}



plotStatus <- function(assessment, AC) {

  require(lattice)

  if (!(AC %in% c("CONC", assessment$AC))) stop("AC not recognised")

#  path <- assessment$webpath
  data <- assessment$status
  data$status.CONC <- data$status
  data$response <- data[[paste0("status.", AC)]]

  labels <- plot.scales(range(data$response, na.rm = TRUE), logData = TRUE, n = 3)

#  fileName <- switch(AC, 
#                     CONC = "status stripplot mean concentration in final year.png", 
#                     paste0("status stripplot relative to ", AC, ".png"))

  ylab <- switch(AC, 
                 CONC = "mean concentration in final year", 
                 paste("mean concentration relative to", AC))

  #    png(file.path(path, fileName), height = 16, width = 16, units = 'cm', res = 150)

  out <- stripplot(
    region ~ response | determinand,  
    data = data,
    scales = list(alternating = FALSE, x = list(at = log(labels), labels = labels)),
    xlab = ylab,
    panel = function(x, y, subscripts) {
      col <- as.character(data$MSTAT[subscripts])
      col <- as.character(factor(col, levels = c("B", "RH", "IH"), labels = c("blue", "green3", "red")))
      panel.abline(v = 0, col = grey(0.7))
      panel.stripplot(x, y, jitter.data = TRUE, pch = 16, col = col)
    })
  out
  #    dev.off()
#  invisible()
}


plot.scales <- function(x, n = 5, min.n = 3, logData = FALSE, f = 0.05) {
  
  # x gives data on log (base e) scale as used in e.g. cstm
  # n is desired number of ticks
  # adapted from axTicks
  
  x <- x[!is.na(x)]
  if (logData) x <- exp(x)
  rng <- range(x)
  
  small <- .Machine$double.eps^0.5
  if (diff(rng) < small) 
  {
    rng <- 
      if (abs(rng[1]) < small) c(-f, f)
    else rng + c(-f, f) * abs(rng[1])
  }
  
  
  lin.calc <- !logData | (rng[2] / rng[1] < 10)     # linear scale
  
  if (lin.calc)
  {
    scales <- mean(rng)
    wk.n <- n
    while (length(scales) < min.n)
    {
      scales <- pretty(rng, wk.n, min.n)
      scales <- scales[scales >= rng[1] & scales <= rng[2]]
      wk.n <- wk.n + 1
    }    
    if (n == 3 & length(scales) %in% c(5, 6)) scales <- scales[c(1,3,5)]
    scales.lin <- scales
  }
  
  if (!logData) return(scales.lin)
  
  ii <- c(ceiling(log10(rng[1])), floor(log10(rng[2])))
  scales.log <- lapply(1:3, function(j)
  {
    x10 <- 10^((ii[1] - (j >= 2)):ii[2])
    scales <- switch(j, x10, c(outer(c(1, 3), x10))[-1], c(outer(c(1, 2, 5), x10))[-1])
    scales[scales >= rng[1] & scales <= rng[2]]  
  })
  
  n.choice <- which.min(abs(sapply(scales.log, length) - n))
  if (length(scales.log[[n.choice]]) < min.n & n.choice < 3) n.choice <- n.choice + 1
  scales.log <- scales.log[[n.choice]]
  if (n == 3 & length(scales.log) %in% c(5, 6)) scales.log <- scales.log[c(1,3,5)]
  
  if (lin.calc && (length(scales.lin) < length(scales.log) | length(scales.log) < n)) 
    scales.lin 
  else 
    scales.log
}




getMetaData <- function(compartmentID, detGroups, byGroup = FALSE) {

  # get assessment object 

  assessmentOb <- local({
    load(file.path("..", "RData backup", paste(compartmentID, "web.RData")))
    get(paste0(compartmentID, ".web"))$assessment
  })
  
  purpose <- assessmentOb$info$purpose
  timeSeries <- assessmentOb$timeSeries
  assessment <- assessmentOb$assessment
  data <- assessmentOb$data
  rm(assessmentOb)

  
  # only retain relevant determinand groups
  
  timeSeries <- timeSeries[timeSeries$detGroup %in% detGroups, ]

  
  # drop estuarine stations - need to provide info within assessmentOb$stations for CSEMP and 
  # need to code for OSPAR
  
  ok <- switch(
    purpose, 
    OSPAR = {
      warning("need to code deletion of estuarine stations", call. = FALSE)
      rep(TRUE, nrow(timeSeries))
    },
    CSEMP = {
      warning("need to pass estuarine information directly to assessment object", call. = FALSE)
      infile = file.path("..", "data", "stndict_merman_041016.xlsx")
      stations <- xlsx::read.xlsx(infile, sheetName = "Fixed Stations", check.names = FALSE)
      stations <- data.frame(stations[c("New Station Code", "ICES Station Type")], 
                             row.names = "New Station Code", check.names = FALSE)
      !(stations[as.character(timeSeries$station), "ICES Station Type"] %in% "Estuary")
    }
  )
      
  timeSeries <- droplevels(timeSeries[ok, ])
  
  
  # sort out region information and ensure ordered sensibly geographically
  
  if (purpose == "OSPAR") {
  
    # combine Dogger Bank with Southern North Sea
    # order regions from North to South

    regionID <- c("Gulf of Cadiz", "Iberian Sea", "Northern Bay of Biscay", "Celtic Sea", "Irish Sea", 
                  "Irish and Scottish West Coast", "Channel", "Southern North Sea", "Skagerrak", 
                  "Northern North Sea", "Norwegian Trench", "Norwegian Sea", "Greenland-Scotland ridge", 
                  "Barents Sea")    
    
    timeSeries <- within(timeSeries, {
      region[region %in% "Dogger Bank"] <- "Southern North Sea"
      region <- droplevels(factor(region, levels = regionID))
    })
    
  } else if (purpose == "CSEMP") {
    
    # drop stations outside reporting areas
    
    timeSeries <- timeSeries[!timeSeries$region %in% c("AtlanticNW", "ENorthSea", "SIrishSea"), ]

    regionID <- c("W Channel & Celtic Sea", "Irish Sea", "Minches & W Scotland", 
                  "Scottish Continental Shelf", "Atlantic & NW Approaches", 
                  "E Channel", "Southern North Sea", "Northern North Sea")

    timeSeries <- within(timeSeries, {
      region[region %in% c("Anglia", "HumWash")] <- "Southern North Sea"
      region[region %in% c("Bailey", "FaroeShetC", "Rockall")] <- "NW Approaches"
      region[region %in% c("CardBay", "Clyde", "IrishSea")] <- "Irish Sea"
      region[region %in% c("EastChan")] <- "E Channel"
      region[region %in% c("EScotland", "EShetland", "Fladen", "Forth", "Forties", "MorayF", 
                           "TyneTees")] <- "Northern North Sea"
      region[region %in% c("Hebrides", "NScotland", "WShetland")] <- "Scottish Continental Shelf"
      region[region %in% c("MinchMalin")] <- "Minches & W Scotland"
      region[region %in% c("Severn", "WestChan")] <- "W Channel & Celtic Sea"
      region <- droplevels(factor(region, levels = regionID))
    })  
        
    if (any(is.na(timeSeries$region)))
      stop("not all CSEMP regions mapped to CP2 regions")

    # create OSPARregion as temporary way of getting maps
        
    timeSeries <- within(timeSeries, {
      MSFDregion <- region
      levels(MSFDregion) <- list(
        "Celtic Seas" = c("W Channel & Celtic Sea", "Irish Sea", "Minches & W Scotland", 
                "Scottish Continental Shelf", "Atlantic & NW Approaches"), 
        "Greater North Sea" = c("E Channel", "Southern North Sea", "Northern North Sea")
      )
    })
  }    

  timeSeries <- droplevels(timeSeries) 
  
   
  # tidy up assessment and data structures
  
  assessment <- assessment[row.names(timeSeries)]
  data <- data[data$seriesID %in% row.names(timeSeries), ]
  data <- droplevels(data)

        
  # sort out MSTAT values
  
  cat("Summary of pre-corrected MSTAT values:\n\n")
  MSTAT.id <- intersect(c("country", "MSTAT"), names(timeSeries))
  MSTAT.table <- do.call(table, c(timeSeries[MSTAT.id], list(useNA = "ifany")))
  print(MSTAT.table)

  timeSeries <- within(timeSeries, {
    MSTAT[grepl("RH", MSTAT)] <- "RH"
    MSTAT[grepl("IH", MSTAT)] <- "IH"
    MSTAT[is.na(MSTAT) | MSTAT %in% ""] <- "RH"
    MSTAT <- MSTAT[, drop = TRUE]
  })  
  
  timeSeries <- switch(
    purpose, 
    OSPAR = within(timeSeries, {
      MSTAT[country %in% c("Belgium", "Denmark") & !(grepl("IH", MSTAT))] <- "RH"
      if (compartmentID %in% "biota") {
        MSTAT[country %in% "United Kingdom" & station %in% "MinchMalin_BroadfordBay_sh01"] <- "RH"
        MSTAT[country %in% "Norway" & grepl("43B2 Trom", station)] <- "IH"
        MSTAT[country %in% "Norway" & station %in% "13B Kristiansand harbour"] <- "IH"
        MSTAT[country %in% "Norway" & MSTAT %in% "B"] <- "RH"
      }
      if ("Imposex" %in% detGroups) {
        MSTAT[country %in% "Sweden"] <- "RH"
        MSTAT[country %in% "Ireland" & MSTAT %in% "IH"] <- "RH"
      }
    }), 
    CSEMP = within(timeSeries, {
      MSTAT[station %in% "MinchMalin_BroadfordBay_sh01"] <- "RH"
    })
  )
  
  cat("\n\nSummary of corrected MSTAT values:\n\n")
  MSTAT.table <- do.call(table, c(timeSeries[MSTAT.id], list(useNA = "ifany")))
  print(MSTAT.table)
  

  # identify timeSeries for which a parametric model has been fitted

  warning("ad hoc approach to testing for trend and status information for imposex: ", 
          "need to rationalise", call. = FALSE)

  paramFit <- sapply(assessment, function(i) !is.null(i$pred))
  timeSeries$paramFit <- paramFit[row.names(timeSeries)]

  
  # get nypos to identify whether there is trend information, or just status

  if (length(detGroups) == 1 && detGroups %in% "Imposex")
    nypos <- sapply(assessment, function(i) if (!is.na(i$summary$rtrend)) 5 else i$summary$nypos)
  else   
    nypos <- sapply(assessment, function(i) i$summary$nypos)

  timeSeries$nypos <- nypos[row.names(timeSeries)]
 


  if (!byGroup) 
    return(
      list(assessment = assessment, timeSeries = timeSeries, data = data, 
           regionID = levels(timeSeries$region), detID = levels(timeSeries$determinand))
    )
  

  # now extract relevant information for each determinand group
  
  sapply(detGroups, simplify = FALSE, FUN = function(ID) {
  
    webpath <- paste(compartmentID, ID)
  
    AC <- switch(webpath, 
                 "sediment Chlorobiphenyls" = c("EAC", "BAC"),
                 "sediment PAH (parent)" = c("ERL", "BAC"),
                 "sediment Metals" = c("ERL", "BAC"),
                 "biota Chlorobiphenyls" = c("EAC", "BAC"),
                 "biota PAH (parent)" = c("BAC", "EAC"),
                 "biota Metals" = c("BAC", "EC"),
                 "biota Imposex" = c("BAC", "EAC"),
                 NULL)
    
    detID <- switch(
      ID, 
      Chlorobiphenyls = c("CB28", "CB52", "CB101", "CB118", "CB138", "CB153", "CB180"), 
      "PAH (parent)" = c("PA", "ANT", "FLU", "PYR", "BAA", "CHR", "BAP", "BGHIP", "ICDP"), 
      "Organo-bromines" = switch(
        compartmentID, 
        sediment = c("BDE28", "BDE47", "BDE99", "BD100", "BD153", "BD154", "BD209"), 
        biota = c("BDE28", "BDE47", "BDE99", "BD100", "BD153", "BD154")), 
      "Organo-metals" = c("TBTIN", "DBTIN", "MBTIN", "TPTIN", "DPTIN", "MPTIN"),
      Metals = c("HG", "CD", "PB"), 
      Imposex = "VDS")
    

    # reduce to relevant determinand group 
    
    timeSeries <- droplevels(subset(timeSeries, determinand %in% detID))
    assessment <- assessment[row.names(timeSeries)]
    data <- droplevels(subset(data, seriesID %in% row.names(timeSeries)))

    timeSeries <- within(timeSeries, determinand <- droplevels(factor(determinand, levels = detID)))

    list(ID = webpath, assessment = assessment, timeSeries = timeSeries, data = data, AC = AC,
         regionID = levels(timeSeries$region), detID = levels(timeSeries$determinand))
  })
}





plotStations <- function(assessment) {
  
  require(plyr)
  
  timeSeries <- assessment$timeSeries
  
  by(timeSeries, timeSeries$OSPARregion, function(data) {
    
    regionID <- as.numeric(as.character(data$OSPARregion))[1]
    
    # status stations

    stations <- unique(data[c("OSPARregion", "latitude", "longitude", "station", "MSTAT")])    

    png(file.path(assessment$webpath, paste("stations region", regionID, "status.png")), 
        height = 16, width = 16, units = 'cm', res = 150)
    
    # plot map
    plot(latitude ~ longitude, data = regionLines[[regionID]], type = "l", xlab = "", ylab = "")
    
    cols_bg <- c("paleturquoise1", "white")
    l_ply(regionPolygons[[regionID]], function(x) 
      polygon(x$longitude, x$latitude, col = cols_bg[x$hole + 1]))
    
    l_ply(subRegionLines, function(x) lines(x$longitude, x$latitude))
    
    # Plot stations
    cols <- factor(stations$MSTAT, levels = c("B", "RH", "IH"), labels = c("blue", "green", "red"))
    cols <- as.character(cols)
    coor_trans <- ctsm.projection(stations$latitude, stations$longitude)
    points(coor_trans$longitude, coor_trans$latitude, pch=21, cex=1.2, col="black", bg=cols)
    
    dev.off()
    
    
    # and now trend stations
    
    ok <- with(data, nypos >= 5)
    if (!any(ok)) return(invisible())
    
    data <- droplevels(data[ok, ])
    stations <- unique(data[c("OSPARregion", "latitude", "longitude", "station", "MSTAT")])    

    png(file.path(assessment$webpath, paste("stations region", regionID, "trend.png")), 
        height = 16, width = 16, units = 'cm', res = 150)
    
    # plot map
    plot(latitude ~ longitude, data = regionLines[[regionID]], type = "l", xlab = "", ylab = "")
    
    cols_bg <- c("paleturquoise1", "white")
    l_ply(regionPolygons[[regionID]], function(x) 
      polygon(x$longitude, x$latitude, col = cols_bg[x$hole + 1]))
    
    l_ply(subRegionLines, function(x) lines(x$longitude, x$latitude))
    
    # Plot stations
    cols <- factor(stations$MSTAT, levels = c("B", "RH", "IH"), labels = c("blue", "green", "red"))
    cols <- as.character(cols)
    coor_trans <- ctsm.projection(stations$latitude, stations$longitude)
    points(coor_trans$longitude, coor_trans$latitude, pch=21, cex=1.2, col="black", bg=cols)
    
    dev.off()
  })
}


plotStations.OSPAR <- function(assessment, subset = c("trend", "status"), OSPARregion) {
  
  subset <- match.arg(subset)
  
  attach(
    file.path("..", "..", "mapping functions", "OSPAR mapping info.RData"), name = "OSPAR shape files")
  on.exit(detach("OSPAR shape files"))
  
  require(plyr)
  require(rgeos)

  # plot map
  
  regionID <- as.numeric(as.character(OSPARregion))

  plot(latitude ~ longitude, data = regionLines[[regionID]], type = "l", xlab = "", ylab = "")
  
  cols_bg <- c("paleturquoise1", "white")
  l_ply(regionPolygons[[regionID]], function(x) 
    polygon(x$longitude, x$latitude, col = cols_bg[x$hole + 1]))
  
  l_ply(subRegionLines, function(x) lines(x$longitude, x$latitude))
  

  # get stations
  
  stations <- assessment$timeSeries
  stations <- stations[stations$OSPARregion %in% OSPARregion, ]
    
  ok <- switch(
    subset, 
    status = with(stations, paramFit), 
    trend = with(stations, paramFit & nypos >= 5)
  )
  stations <- stations[ok, ]
  
  stations <- unique(stations[c("latitude", "longitude", "MSTAT")])

  stations <- within(stations, {
    bg <- factor(MSTAT, levels = c("RH", "B", "IH"), labels = c("green3", "blue", "red"))
    ord <- order(bg)
    bg <- as.character(bg)
  })
  stations <- stations[stations$ord, ]
  
    
  # Plot stations

  coor_trans <- ctsm.projection(stations$latitude, stations$longitude)
  points(coor_trans$longitude, coor_trans$latitude, pch = 21, cex = 1.2, col = "black", bg = stations$bg)
}




plotAllStations.OSPAR <- function(assessment) {

  attach(
    file.path("..", "..", "mapping functions", "OSPAR mapping info.RData"), name = "OSPAR shape files")
  on.exit(detach("OSPAR shape files"))
  
  require(plyr)
  require(rgeos)


  # get stations
  
  id <- c("OSPARregion", "latitude", "longitude", "station", "paramFit", "nypos")
  timeSeries <- assessment$timeSeries[id]
  
  timeSeries <- timeSeries[timeSeries$paramFit, ]
  timeSeries$paramFit <- NULL
  
  
  # categorise as status or trend station
  
  maxNypos <- aggregate(timeSeries["nypos"], by = timeSeries["station"], max)
  maxNypos <- within(maxNypos, type <- ifelse(nypos >= 5, "trend", "status"))
  
  stations <- merge(unique(timeSeries[c("OSPARregion", "latitude", "longitude", "station")]), maxNypos, 
                    all.x = TRUE)
  
  
    
  # plot map

  regionID <- unique(as.numeric(as.character(stations$OSPARregion)))
  xlim <- extendrange(do.call(range, c(lapply(regionLines[regionID], "[[", "longitude"), na.rm = TRUE)), 
                      f = 0.01)
  ylim <- extendrange(do.call(range, c(lapply(regionLines[regionID], "[[", "latitude"), na.rm = TRUE)), 
                      f = 0.01)

  plot(latitude ~ longitude, data = regionLines[[regionID[1]]], type = "n", xlab = "", ylab = "", 
       xlim = xlim, ylim = ylim)
  
  l_ply(regionLines, function(x) lines(x$longitude, x$latitude))
  
  cols_bg <- c("paleturquoise1", "white")
  l_ply(regionPolygons, function(y)
    l_ply(y, function(x) 
      polygon(x$longitude, x$latitude, col = cols_bg[x$hole + 1])))
  
  l_ply(subRegionLines, function(x) lines(x$longitude, x$latitude))
  
  
  # plot stations
  
  cols <- ifelse(stations$type == "trend", "black", "white")
  coor_trans <- ctsm.projection(stations$latitude, stations$longitude)
  points(coor_trans$longitude, coor_trans$latitude, pch=21, cex=1.2, col="black", bg=cols)
  
}


plotAllStations.CSEMP <- function(assessment) {
 
  attach(file.path("..", "..", "mapping functions", "CP2 borders.RData"), name = "CP2 shape files")
  on.exit(detach("CP2 shape files"))
    
  require(rgeos)
  plot(CP2.borders, axes=TRUE, col = "white", xlim = c(-10, 4), ylim = c(49, 62), bg = grey(0.9))
  box()
  
  stations <- unique(assessment$timeSeries[c("latitude", "longitude", "station", "MSTAT")])
  stations <- within(stations, {
    bg <- factor(MSTAT, levels = c("RH", "B", "IH"), labels = c("green3", "blue", "red"))
    ord <- order(bg)
    bg <- as.character(bg)
  })
  stations <- stations[stations$ord, ]
  
  with(stations, points(longitude, latitude, pch = 21, bg = bg))
}
  

plotStations.CSEMP <- function(assessment) {
  
  attach(file.path("..", "..", "mapping functions", "CP2 borders.RData"), name = "CP2 shape files")
  on.exit(detach("CP2 shape files"))

  attach(file.path("..", "..", "mapping functions", "MSFD borders.RData"), name = "MSFD shape files")
  on.exit(detach("MSFD shape files"))
    
  require(rgeos)
  plot(CP2.borders, axes=TRUE, col = "white", xlim = c(-10, 4), ylim = c(49, 62), bg = grey(0.9))
  box()
  
  plyr::l_ply(MSFD.borders @ polygons, function(y) 
    plyr::l_ply(y @ Polygons, function(x) 
      lines(x @ coords[, 1], x @ coords[, 2], lwd = 2)))

    
  # get appropriate subset of stations
  
  varID <- c("station", "latitude", "longitude")
  trendStations <- unique(assessment$trend[varID])
  statusStations <- unique(assessment$status[varID])

  trendStations$type <- "trend"
  statusStations$type <- "status"
  
  # since trend stations are first, a station will only have type status if it doesn't
  # have any trend information
  
  stations <- rbind(trendStations, statusStations)
  stations <- stations[!duplicated(stations[varID]), ]
  
  # get appropriate colour scheme 

  stations <- unique(stations[c("latitude", "longitude", "type")])
  
  stations <- within(stations, {
    # pch <- ifelse(type %in% "trend", 16, 1)
    bg <- ifelse(type %in% "trend", "magenta", "cyan")
  })

  # plot stations
  
  with(stations, points(longitude, latitude, pch = 21, cex = 1.2, bg = bg))
}




writeData <- function(assessment) {
  path <- "data for assessors"
  ID <- assessment$ID
  if (!is.null(assessment$trend)) {
    fileName <- paste(ID, "trend all.csv")
    write.csv(assessment$trend, file.path(path, fileName), row.names = FALSE, na = "")
  }
  if (!is.null(assessment$status)) {
    fileName <- paste(ID, "status all.csv")
    write.csv(assessment$status, file.path(path, fileName), row.names = FALSE, na = "")
  }
  if (!is.null(assessment$regionalTrend)) {
    fileName <- paste(ID, "trend meta analysis.csv")
    write.csv(assessment$regionalTrend, file.path(path, fileName), row.names = FALSE, na = "")
  }
  if (!is.null(assessment$regionalStatus)) {
    fileName <- paste(ID, "status meta analysis.csv")
    write.csv(assessment$regionalStatus, file.path(path, fileName), row.names = FALSE, na = "")
  }
}



# functions to calculate percentage compliance
# get appropriate data set

getSampleData <- function(assessment, yearRange, detGroup, detOrder, ACtype = c("EAC", "BAC")) {
  
  ACtype <- match.arg(ACtype)
  ACtype <- switch(ACtype, EAC = c("EAC", "ERL", "EC"), BAC = "BAC")
  
  data <- assessment$data

  # reduce data to relevant detGroup
  
  if (!missing(detGroup))
    data <- droplevels(data[data$group %in% detGroup, ])

  
  # merge region information with data and also pick up latitude and longitude
  
  id <- intersect(c("seriesID", "species", "year", "determinand", "concentration"), names(data))
  data <- data[id]
  data <- within(data, seriesID <- as.character(seriesID))
  data <- data[data$year %in% yearRange & !is.na(data$concentration), ]
  data <- cbind(data, assessment$timeSeries[data$seriesID, c("region", "MSFDregion", "latitude", "longitude")])
  
  # get relevant EAC values
  
  AC <- lapply(assessment$assessment, function(x) {
    out <- x$summary
    out <- out[intersect(ACtype, names(out))]
    out <- as.data.frame(t(out))
    names(out) <- "AC"
    out$ACid <- row.names(out)
    if (all(is.na(out$AC))) return (data.frame(AC = NA, ACid = NA))
    out <- na.omit(out)
    if (nrow(out) > 1) stop("multiple EAC found")
    out
  })
  
  AC <- do.call(rbind, AC)
  
  # merge with data, drop data with no AC, drop data with missing concentrations
  # (either due to poor uncertainty or no normalising information) and order determinand
  
  data <- cbind(data, AC[data$seriesID, ])
  data <- data[!is.na(data$AC), ]
  data <- data[!is.na(data$concentration), ]
  
  if (!missing(detOrder))
    data <- within(data, determinand <- factor(determinand, levels = detOrder))

  data <- droplevels(data)
  data
}


# percentage compliance

tableCompliance <- function(data, regionID = c("region", "MSFDregion")) {
  regionID <- match.arg(regionID)
  data$regionID <- data[[regionID]]
  with(data, {
    diff <- AC - concentration
    tapply(diff, list(regionID, determinand), function(x) 
      floor(100 * sum(x > 0) / length(x)))
  })
}  


# number of samples

tableSamples <- function(data, regionID = c("region", "MSFDregion")) {
  regionID <- match.arg(regionID)
  data$regionID <- data[[regionID]]
  out <- with(data, table(regionID, determinand))
  names(dimnames(out)) <- NULL
  out
}


# number of stations from sampleData structure

tableStations <- function(data, regionID = c("region", "MSFDregion")) {
  
  regionID <- match.arg(regionID)
  data$regionID <- data[[regionID]]
  data <- within(data, regionID <- factor(regionID, levels = rev(levels(regionID))))
  
  # use latitude and longitude, because there could be multiple species
  # at a particular station, or multiple countries at a particular place (OSPAR)
  
  data <- unique(data[c("latitude", "longitude", "regionID", "determinand")])
  out <- with(data, table(regionID, determinand))
  names(dimnames(out)) <- NULL
  out
}  


# AC summary

tableAC <- function(data) {
  
  ACid <- unique(data$ACid)
  if (length(ACid) > 1) stop("multiple AC types")
  
  
  # check if multiple AC - usually because for biota expressed on a non wet weight basis
  
  multipleAC <- with(unique(data[c("determinand", "AC")]), any(duplicated(determinand)))
  if (multipleAC & !("species" %in% names(data))) stop("multiple AC for a determinand")
  
  
  if (!multipleAC) {
    data <- unique(data[c("determinand", "AC")])
    data <- data[order(data$determinand), ]
    names(data)[1:2] <- c("", ACid)
    return(print(data, row.names = FALSE))
  }
  
  
  # now tabulate by species
  
  data <- unique(data[c("determinand", "AC", "species")])
  
  speciesLevels <- c("Limanda limanda", "Platichthys flesus", "Pleuronectes platessa",  
                     "Merlangius merlangus", "Mytilus edulis", "Crassostrea gigas")
  speciesLabels = c("dab", "flounder", "plaice", "whiting", "mussel", "oyster")
  
  data <- within(data, species <- factor(species, levels = speciesLevels, labels = speciesLabels))
  if (any(is.na(data$species))) stop("some common species names not coded")

  data <- reshape(data, direction = "wide", timevar = "species", idvar = "determinand")
  names(data) <- gsub("AC.", "", names(data))
  
  
  # order data set
  
  speciesLabels <- intersect(speciesLabels, names(data))
  
  data <- data[order(data$determinand), c("determinand", speciesLabels)]
  names(data)[1] <- ""
  print(data, row.names = FALSE)
}
  
  
# gets number of stations by region and determinand from trend or status data set

tableStationsTS <- function(data, AC) {
  
  # if status data, remove any time series with missing AC (might have a mean concentration and an EAC
  # but no BAC to go with it)
  
  if (!missing(AC) && AC != "CONC") {
    ok <- !is.na(data[[paste0("status.", AC)]])
    data <- data[ok, ]
  }
  
  data <- within(data, region <- factor(region, levels = rev(levels(region))))
  
  # get unique stations (might be multiple time series due to different species)
  # use lat and long here because haven't 
  
  data <- unique(data[c("latitude", "longitude", "region", "determinand")])
  out <- with(data, table(region, determinand))
  names(dimnames(out)) <- NULL
  out
}  



getData.Imposex <- function(x) {
  
  AC <- x$AC
  assessment <- x$assessment
  timeSeries <- x$timeSeries
  
  
  # get status stations
  
  ok <- with(timeSeries, paramFit)
  if (!any(ok)) return(x)
  
  timeSeries <- droplevels(timeSeries[ok, ])
  assessment <- assessment[row.names(timeSeries)]

  Status <- sapply(assessment, USE.NAMES = TRUE, simplify = FALSE, FUN = function(assessment) {
    out <- tail(assessment$pred, 1)
    out$index <- NULL # only in some time series!
    cbind(out, assessment$summary[c("meanLY", "clLY", AC)])
  })
  
  Status <- do.call(rbind, Status)
  
  id <- intersect(c("determinand", "country", "station", "species", "OSPARregion", "MSFDregion", 
                    "region", "offshore", "MSTAT", "latitude", "longitude"), names(timeSeries))
  Status[id] <- timeSeries[row.names(Status), id]
  
  warning("ad-hoc fix when meanLY = 0", call. = FALSE)
  
  Status <- within(Status, {
    fit <- sqrt(meanLY)
    se <- (sqrt(clLY) - sqrt(meanLY)) / 1.645
  })
  
  
  if (!is.null(AC)) Status[paste("status", AC, sep = ".")] <- Status$fit / sqrt(Status[AC])
  
  # se.status is based on an AC of 1 - will need to adjust in Genstat (too complicated to do that here!)
  # also ad-hoc adjustment to very small standard errors (usually when trend diving to minus infinity)
  # made as larges as the smallest standard error, when estimated values is above the BAC
  
  warning("ad-hoc adjustment to small standard errors", call. = FALSE)
  
  Status <- within(Status, {
    status <- fit
    se.status <- pmax(se, min(se[status > BAC], na.rm = TRUE))
    station <- factor(station)
    determinand <- factor(as.character(determinand), levels = x$detID)
  })
  
  x$status <- droplevels(Status)
  
  
  # now look at only time series with trend information (given by nypos) 
  
  ok <- with(timeSeries, nypos >= 5)
  if (!any(ok)) return(x)
  
  timeSeries <- droplevels(timeSeries[ok, ])
  assessment <- assessment[row.names(timeSeries)]
  
  Trend <- sapply(assessment, USE.NAMES = TRUE, simplify = FALSE, function(assessment) {
    if ("contrasts" %in% names(assessment)) {
      out <- assessment$contrasts["recent", ]
      out <- within(out, {
        trend <- estimate / (end - start)
        se.trend <- se / (end - start)
      })      
      out <- out[c("trend", "se.trend")]
    } 
    else {
      out <- assessment$coefficients["year", c("Estimate", "Std. Error")]
      names(out) <- c("trend", "se.trend")
      out
    }
  })
  
  Trend <- do.call(rbind, Trend)
  
  Trend[id] <- timeSeries[row.names(Trend), id]

  Trend <- within(Trend, {

    # look at on exponential scale to avoid extreme low value having an undue effect (perhaps)
    # but some extreme trends then have very small standard errors, so have to make an
    # ad-hoc adjustment to get predictions
    
    trend <- exp(trend)
    se.trend <- se.trend * trend
    se.trend <- pmax(se.trend, min(se.trend[trend > 0.1], na.rm = TRUE))
    
    station <- factor(station)
    determinand <- factor(as.character(determinand), levels = x$detID)
  })
  
  x$trend <- droplevels(Trend)
  
  x
}



plotTrend.Imposex <- function(assessment) {
  require(lattice)
  path <- assessment$webpath
  if (is.null(assessment$trend)) return(invisible())
  data <- assessment$trend
  #  png(file.path(path, paste("trend stripplot.png")), height = 16, width = 16, units = 'cm', res = 150)
  at <- c(0, 0.5, 1, 1.5, 2)
  labels <- at
  out <- stripplot(
    region ~ trend | determinand,
    data = data, 
    xlab = "VDS trend: odds of exceeding EAC this year relative to last",
    scales = list(alternating = FALSE, x = list(at = at, labels = labels)),
    panel = function(x, y, subscripts) {
      col <- as.character(data$MSTAT[subscripts])
      col <- as.character(factor(col, levels = c("B", "RH", "IH"), labels = c("blue", "green3", "red")))
      panel.abline(v = 1, col = grey(0.7))
      panel.stripplot(x, y, jitter.data = TRUE, pch = 16, col = col)
    })
  out
  #  dev.off()
  #  invisible()
}


plotStatus.Imposex <- function(assessment, AC) {
  
  require(lattice)
  
  if (!(AC %in% c("CONC", assessment$AC))) stop("AC not recognised")
  
  data <- assessment$status
  data$status.CONC <- data$status
  data$response <- data[[paste0("status.", AC)]]

  data <- data[!is.na(data$response), ]
    
  if (sqrt(max(data$response)) < 1.5)
    at <- c(0, 0.5, 0.8, 1, 1.2, 1.4)
  else
    at <- c(0, 0.7, 1, 1.5, 2.0, 3.0)
  
  labels <- as.character(at)
  

  #  fileName <- switch(AC, 
  #                     CONC = "status stripplot mean concentration in final year.png", 
  #                     paste0("status stripplot relative to ", AC, ".png"))
  
  ylab <- switch(AC, 
                 CONC = "mean VDS in final year", 
                 paste("mean VDS relative to", AC))
  
  #    png(file.path(path, fileName), height = 16, width = 16, units = 'cm', res = 150)
  
  out <- stripplot(
    region ~ response | determinand,  
    data = data,
    scales = list(alternating = FALSE, x = list(at = at ^ 2, labels = labels)),
    xlab = ylab,
    panel = function(x, y, subscripts) {
      col <- as.character(data$MSTAT[subscripts])
      col <- as.character(factor(col, levels = c("B", "RH", "IH"), labels = c("blue", "green3", "red")))
      panel.abline(v = 1, col = grey(0.7))
      panel.stripplot(x, y, jitter.data = TRUE, pch = 16, col = col)
    })
  out
  #    dev.off()
  #  invisible()
}



plotStatusResults.Imposex <- function(
  assessment, 
  type = c("region", "determinand", "region by determinand", "determinand by region"), 
  AC) {
  
  require(lattice)
  
  if (is.null(assessment$regionalStatus)) return(invisible())
  
  type <- match.arg(type)
  
  if (!(AC %in% c("CONC", assessment$AC))) stop("AC not recognised")
  
  # ensure plotting doesn't include combinations that don't exist
  
  originalData <- assessment$regionalStatus
  ok <- with(originalData, !is.na(switch(AC, "CONC" = status, get(paste0("status.", AC)))))
  originalData <- originalData[ok, ]
  regDetExists <- with(originalData, unique(paste(region, determinand)))
  
  
  plotEngine <- function(data, formula, AC) {
    xlim <- switch(AC, 
                   CONC = extendrange(range(data$upper, data$lower)),
                   extendrange(range(data$upper, data$lower, 0)))
    if (sqrt(xlim[2]) < 1.5)
      at <- c(0, 0.5, 0.8, 1, 1.2, 1.4)
    else
      at <- c(0, 0.7, 1, 1.2, 1.5, 2.0, 3.0)
    labels <- as.character(at)
    xlab <- switch(AC, 
                   CONC = "mean VDS in final year", 
                   paste("mean VDS relative to", AC))
    stripplot(
      formula,
      data = data, 
      xlab = xlab,
      xlim = xlim, 
      scales = list(alternating = FALSE, x = list(at = at ^ 2, labels = labels)),
      panel = function(x, y, subscripts) {
        if (AC != "CONC") panel.abline(v = 1, col = grey(0.7))
        with(data[subscripts, ], lsegments(lower, y, upper, y, col = "black"))
        col <- switch(AC, 
                      CONC = "black", 
                      BAC = with(data[subscripts, ], ifelse(upper < 1, "blue", "green3")),
                      with(data[subscripts, ], ifelse(upper < 1, "green3", "red")))
        lpoints(x, y, col = col, pch = 16, cex = 1.5)
      })
  }
  
  txt <- paste0("status results", switch(AC, CONC = "", paste0(" ", AC)))
  z <- switch(AC, CONC = 1.96, 1.645)
  
  switch(
    type, 
    region = {
      txt2 <- paste(txt, "by region")
      data <- readFromGenstat(txt2, assessment, z = z)
      plotEngine(data, as.formula(region ~ est), AC)
    },  
    determinand = {
      txt2 <- paste(txt, "by determinand")
      data <- readFromGenstat(txt2, assessment, z = z)
      plotEngine(data, as.formula(determinand ~ est), AC)
    },  
    "region by determinand" = {
      txt2 <- paste(txt, "by region determinand")
      data <- readFromGenstat(txt2, assessment, z = z)
      data <- subset(data, paste(region, determinand) %in% regDetExists)
      plotEngine(data, as.formula(region ~ est | determinand), AC)
    },
    "determinand by region" = {
      txt2 <- paste(txt, "by region determinand")
      data <- readFromGenstat(txt2, assessment, z = z)
      data <- subset(data, paste(region, determinand) %in% regDetExists)
      plotEngine(data, as.formula(determinand ~ est | region), AC)
    }
  )
  
}


plotTrendResults.Imposex <- function(
  assessment, type = c("region", "determinand", "region by determinand", "determinand by region")) {
  
  require(lattice)
  
  if (is.null(assessment$regionalTrend)) return(invisible())
  
  type <- match.arg(type)
  
  # ensure plotting doesn't include combinations that don't exist
  
  regDetExists <- with(assessment$regionalTrend, unique(paste(region, determinand)))
  
  
  plotEngine <- function(data, formula) {
    xlim <- extendrange(range(data$upper, data$lower, 1))
    stripplot(
      formula,
      data = data, 
      xlab = "VDS trend: odds of exceeding EAC this year relative to last",
      xlim = xlim, 
      scales = list(alternating = FALSE), 
      panel = function(x, y, subscripts) {
        panel.abline(v = 1, col = grey(0.7))
        with(data[subscripts, ], lsegments(lower, y, upper, y, col = "black"))
        col <- with(data[subscripts, ], 
                    ifelse(upper < 1, "green3", ifelse(lower > 1, "red", "goldenrod")))
        lpoints(x, y, col = col, pch = 16, cex = 1.5)
      })
  }
  
  switch(
    type, 
    region = {
      data <- readFromGenstat("trend results by region", assessment)
      plotEngine(data, as.formula(region ~ est))
    },
    determinand = {
      data <- readFromGenstat("trend results by determinand", assessment)
      plotEngine(data, as.formula(determinand ~ est))
    }, 
    "region by determinand" = {
      data <- readFromGenstat("trend results by region determinand", assessment)
      data <- subset(data, paste(region, determinand) %in% regDetExists)
      plotEngine(data, as.formula(region ~ est | determinand))
    },
    "determinand by region" = {
      data <- readFromGenstat("trend results by region determinand", assessment)
      data <- subset(data, paste(region, determinand) %in% regDetExists)
      plotEngine(data, as.formula(determinand ~ est | region))
    }
  )
}
