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


getStatusFromGenstat <- function(
  assessment, type = c("region", "determinand", "region by determinand", "determinand by region"), 
  ACchoice, referenceValue = 0) {

  if (ACchoice == "CONC") {
    txt <- paste("status results", switch(
      type, region = "by region", determinand = "by determinand", "by region determinand"))
    data <- readFromGenstat(txt, assessment, z = 1.96)
    data$colour <- "black"
    return(data)
  }
  
  status <- sapply(assessment$AC, simplify = FALSE, FUN = function(AC) {
    txt <- paste("status results", AC)
    txt <- paste(txt, switch(
      type, region = "by region", determinand = "by determinand", "by region determinand"))
    data <- readFromGenstat(txt, assessment, z = 1.645)
    
    # omit combinations that don't exist
    
    if (type %in% c("region by determinand", "determinand by region")) {
      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)))
    
      ok <- with(data, paste(region, determinand)) %in% regDetExists
      data[!ok, c("est", "se", "upper", "lower")] <- NA
    } 
    
    data
  })

  mergeID <- intersect(c("region", "determinand"), names(status[[1]]))
    
  status <- mapply(function(data, id) {
    data <- within(data, colour <- switch(
      id, 
      BAC = ifelse(upper < referenceValue, "blue", "green"),
      ifelse(is.na(upper), "orange", ifelse(upper <  referenceValue, "green", "red"))
#      ifelse(upper < 0, "green", ifelse(est <  0, "orange", "red"))
    ))
    
    if (! id %in% ACchoice) 
      data <- data[c(mergeID, "colour")]

    if (id %in% "BAC")
      names(data)[names(data) %in% "colour"] <- "colourBAC"
    
    data
  }, status, names(status), SIMPLIFY = FALSE)
  
  status <- merge(status[[1]], status[[2]], by = mergeID, all.x = TRUE)
  status <- within(status, {
    colour <- ifelse(colourBAC %in% "blue", "blue", colour)
    rm(colourBAC)
  })
  status
}




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.4))
        with(data[subscripts, ], lsegments(lower, y, upper, y, col = "black"), lwd= 1.5)
        pch <- with(data[subscripts, ], ifelse(upper < 0, 25, ifelse(lower > 0, 24, 21)))
        lpoints(x, y, pch = pch, cex = 1.5, col = "black", fill = "white", lwd = 2)
      })
  }
  
  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, greenID = "green") {

  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, na.rm = TRUE)),
                   extendrange(range(data$upper, data$est, 0, na.rm = TRUE)))
    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, 
      drop.unused.levels = FALSE, 
      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))
        if (nrow(data[subscripts, ]) == 0) return()
        if (AC == "CONC")
          with(data[subscripts, ], lsegments(lower, y, upper, y, col = "black"))
        else
          with(data[subscripts, ], lsegments(x, y, upper, y, col = "black"))
        col <- data$col[subscripts]
        col[col %in% "green"] <- greenID
        col[col %in% "orange"] <- "goldenrod"
        lpoints(x, y, col = "black", fill = col, pch = 21, cex = 1.5)
      })
  }

  data <- getStatusFromGenstat(assessment, type, AC)

  if (type %in% c("region by determinand", "determinand by region"))
    data <- data[with(data, paste(region, determinand)) %in% regDetExists, ]
  
  switch(
    type, 
    region = plotEngine(data, as.formula(region ~ est), AC),
    determinand = plotEngine(data, as.formula(determinand ~ est), AC),
    "region by determinand" = plotEngine(data, as.formula(region ~ est | determinand), AC),
    "determinand by region" = plotEngine(data, as.formula(determinand ~ est | region), AC)
  )
}


plotTrend <- function(assessment) {
  require(lattice)
  path <- assessment$webpath
  if (is.null(assessment$trend)) return(invisible())
  data <- assessment$trend
  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) {
      pch <- data$shape[subscripts]
      pch[pch == "filled_circle"] <- "21"
      pch[pch == "upward_triangle"] <- "24"
      pch[pch == "downward_triangle"] <- "25"
      pch <- as.numeric(pch)
      panel.abline(v = 0, col = grey(0.4))
      panel.stripplot(x, y, jitter.data = TRUE, pch = pch, col = "black")
    })
  out
}


plotTrendSE <- function(assessment) {
  require(lattice)
  if (is.null(assessment$trend)) return(invisible())
  data <- assessment$trend
  data <- within(data, se.trend <- se.trend * sign(trend))
  out <- stripplot(
    region ~ se.trend | determinand,
    data = data, 
    xlab = "standard error on trend",
    scales = list(alternating = FALSE),
    panel = function(x, y, subscripts) {
      pch <- data$newShape[subscripts]
      pch[pch == "circle"] <- "21"
      pch[pch == "upward_triangle"] <- "24"
      pch[pch == "downward_triangle"] <- "25"
      pch <- as.numeric(pch)
      panel.abline(v = 0, col = grey(0.4))
      panel.stripplot(x, y, jitter.data = TRUE, pch = pch, col = "black")
    })
  out
}


plotStatus <- function(assessment, AC, greenID = "green") {

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

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

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

  out <- stripplot(
    region ~ response | determinand,  
    data = data,
    scales = list(alternating = FALSE, x = list(at = log(labels), labels = labels)),
    xlab = xlab,
    panel = function(x, y, subscripts) {
      col <- data$colour[subscripts]
      col[col == "orange"] <- "goldenrod"
      if (! AC %in% "CONC") 
        panel.abline(v = 0, col = grey(0.4))
      panel.stripplot(x, y, jitter.data = TRUE, pch = 16, col = col)
    })
  out
}


plotStatusSE <- function(assessment, AC) {
  
  require(lattice)
  
  if (!(AC %in% c("CONC", assessment$AC))) stop("AC not recognised")
  
  data <- assessment$status
  data <- within(data, se.status <- 100 * se.status)

  xlab <- "cv on status estimates (%)" 

  out <- stripplot(
    region ~ se.status | determinand,  
    data = data,
    scales = list(alternating = FALSE),
    xlab = xlab,
    panel = function(x, y, subscripts) {
      panel.stripplot(x, y, jitter.data = TRUE, pch = 16, col = "black")
    })
  out
}



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
}


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"), add = TRUE)
    
  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))
}


plotStations <- function(assessment, regionID, purpose = c("OSPAR", "CSEMP"), includeAll = TRUE) {
  
  purpose <- match.arg(purpose)
  
  require(rgeos)
  require(plyr)
  
  if (purpose == "OSPAR") {
    
    attach(
      file.path("..", "..", "mapping functions", "OSPAR mapping info.RData"), name = "OSPAR shape files")
    on.exit(detach("OSPAR shape files"))

    if (missing(regionID))
      regionID <- as.character(1:4)
    
    plotOSPARregions(regionID)
  }
  
  
  # get stations 
  # includeAll can get from timeSeries
  # ! includeAll need to use trend and status data sets (or have to recode selection and check correct)

  varID <- c("OSPARregion", "latitude", "longitude")

  if (includeAll) {
    
    stations <- assessment$timeSeries
    stations <- stations[! stations$shape %in% "open_circle", ]
    stations <- within(stations, type <- ifelse(nypos >= 5, "trend", "status"))

    # order by nypos, so trend stations come first
    
    stations <- stations[order(stations$nypos, decreasing = TRUE), ]
    stations <- stations[c(varID, "type")]

  } else {
  
    trendStations <- unique(assessment$trend[varID])
    statusStations <- unique(assessment$status[varID])
  
    trendStations$type <- "trend"
    statusStations$type <- "status"
  
    # ensure trend stations come first
    
    stations <- rbind(trendStations, statusStations)
    
  }
    
  # trend stations come first, so a station will only have type status if it doesn't have trend info
  
  stations <- stations[!duplicated(stations[varID]), ]
  
  
  stations <- stations[stations$OSPARregion %in% regionID, ]
    

  # get appropriate colour scheme 
  
  stations <- within(stations, {
    type <- factor(type, levels = c("status", "trend"))
    bg <- ifelse(type %in% "trend", "magenta", "cyan")
  })
  
  stations <- stations[order(stations$type), ]

  # plot stations
  
  if (purpose == "OSPAR")
    stations <- transformCoordinates(stations)
  
  cex <- switch(purpose, OSPAR = 1, CSEMP = 1.2)
  
  with(stations, points(longitude, latitude, pch = 21, cex = cex, col = "black", bg = bg))
  
}


plotTimeSeries.CSEMP <- function(assessment, detID, purpose = c("OSPAR", "CSEMP")) {
  
  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"), add = TRUE)
  
  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("determinand", "latitude", "longitude", "newColour", "newShape")
  stations <- assessment$timeSeries[varID]
  ok <- with(stations, !is.na(newColour) & determinand %in% detID)
  stations <- stations[ok, ]


  # get appropriate colour scheme 
  
  stations <- within(stations, {
    bg <- factor(newColour, levels = c("blue", "green", "orange", "red"), 
                 labels = c("blue", "green", "goldenrod", "red"))
    pch <- factor(newShape, levels = c("circle", "downward_triangle", "upward_triangle"), 
                  labels = c("21", "25", "24"))

    ord <- order(bg, pch)

    bg <- as.character(bg)
    pch <- as.numeric(as.character(pch))
  })
    
  # plot stations
  
  with(stations, points(longitude, latitude, pch = pch, cex = 1.2, bg = bg))
}


plotTimeSeries <- function(assessment, detID, regionID, purpose = c("OSPAR", "CSEMP")) {
  
  purpose <- match.arg(purpose)

  require(rgeos)
  require(plyr)
  
  if (purpose == "OSPAR") {

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

    if (missing(regionID))
      regionID <- as.character(1:4)
    
    plotOSPARregions(regionID)
  }
   
  
  # get appropriate subset of stations
  
  varID <- c("determinand", "latitude", "longitude", "colour", "shape")
  stations <- assessment$timeSeries[varID]
  
  ok <- with(stations, ! shape %in% "open_circle" &  determinand %in% detID)
  stations <- stations[ok, ]
  
  
  # get appropriate colour scheme 
  
  stations <- within(stations, {
    bg <- factor(colour, levels = c("black", "blue", "green", "orange", "red"), 
             labels = c("black", "blue", "green", "goldenrod", "red"))

    pch <- factor(shape, levels = c("filled_circle", "downward_triangle", "upward_triangle"), 
                  labels = c("21", "25", "24"))
    
    ord <- order(bg, pch)
    
    bg <- as.character(bg)
    pch <- as.numeric(as.character(pch))
  })
  
  # plot stations
  
  stations <- stations[stations$ord, ]
  
  if (purpose == "OSPAR")
    stations <- transformCoordinates(stations)
  
  cex <- switch(purpose, OSPAR = 1, CSEMP = 1.2)
  
  with(stations, points(longitude, latitude, pch = pch, cex = cex, col = "black", bg = bg))
}


transformCoordinates <- function(data) {
  stopifnot(c("latitude", "longitude") %in% names(data))
  coor_trans <- ctsm.projection(data$latitude, data$longitude)
  data[c("latitude", "longitude")] <- coor_trans[c("latitude", "longitude")]
  data
}


plotOSPARregions <- function(regionID) {
  
  regionID <- as.numeric(as.character(regionID))
  
  if (length(regionID) == 1) {

    plot(latitude ~ longitude, data = regionLines[[regionID]], type = "l", xlab = "", ylab = "")
  
  } else {

    regionLines <- regionLines[regionID]
        
    xlim <- extendrange(do.call(range, c(lapply(regionLines, "[[", "longitude"), na.rm = TRUE)), 
                        f = 0.01)
    ylim <- extendrange(do.call(range, c(lapply(regionLines, "[[", "latitude"), na.rm = TRUE)), 
                        f = 0.01)
    
    plot(latitude ~ longitude, data = regionLines[[1]], type = "n", xlab = "", ylab = "", xlim = xlim, 
         ylim = ylim)
    
    l_ply(regionLines, function(x) lines(x$longitude, x$latitude))
    
  }  
    
  cols_bg <- c("grey90", "white")
  l_ply(regionPolygons[regionID], 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))
}


plotRegionalSummary.CSEMP <- function(assessment, detID) {
  
  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"), add = TRUE)
  
  
  # get region names and appropriate status colour
  
  status <- CP2.borders @ data ["REGION"]
  names(status) <- "region"
  status$colour <- "white"

  status <- within(status, {
    region <- as.character(region)
    region[region == "Western Channel & Celtic Sea"] <- "W Channel & Celtic Sea"
    region[region == "Minches & Western Scotland"] <- "Minches & W Scotland"
    region[region == "Atlantic North-West Approaches, Rockall Trough and"] <- 
      "Atlantic & NW Approaches"
    region[region == "Eastern Channel"] <- "E Channel"
  })
  
  row.names(status) <- status$region
  
  

  # first adapt for where there is some data

  statusData <- assessment$status
  if (!missing(detID))
    statusData <- statusData[statusData$determinand %in% detID, ]
  statusData <- as.character(unique(statusData$region))
  
  regionID <- with(status, region[region %in% statusData])
  status[regionID, "colour"] <- "lightgray"


  # and now for where there was a regional assessment
  
  type <- if (!missing(detID)) "region by determinand" else "region"
  statusData <- getStatusFromGenstat(assessment, type, setdiff("assessment$AC", "BAC"))
  if (!missing(detID))
    statusData <- statusData[statusData$determinand %in% detID, ]
  row.names(statusData) <- as.character(statusData$region)
  
  regionID <- with(status, region[region %in% row.names(statusData)])
  status[regionID, "colour"] <- statusData[regionID, "colour"]

  
  # now adapt colour so can see trend lines
  
  status <- within(status, {
    colour[colour == "green"] <- "green3"
    colour[colour == "blue"] <- "deepskyblue"
  })
  
  
  # now do trends

  type <- if (!missing(detID)) "region determinand" else "region"
  trendData <- readFromGenstat(paste("trend results by", type), assessment)
  if (!missing(detID))
    trendData <- trendData[trendData$determinand %in% detID, ]
  trendData <- within(trendData, region <- as.character(region))
  row.names(trendData) <- trendData$region
  
  trendData <- within(trendData, {
    trend <- ifelse(upper < 0, "down", ifelse(lower > 0, "up", "stable"))
  })

  regionID <- with(trendData, region[region %in% row.names(trendData)])
  status$trend <- NA
  status[regionID, "trend"] <- trendData[regionID, "trend"]
  
  
  require(rgeos)
  plot(CP2.borders, axes=TRUE, col = status$colour, xlim = c(-10, 4), ylim = c(49, 62), bg = "lightcyan")
  box()

  
  status <- status[!is.na(status$trend), ]
  
  lapply(row.names(status), function(regionID) {
    pos <- switch(
      regionID, 
      "Northern North Sea" = c(0, 57), 
      "Southern North Sea" = c(1.75, 53.5), 
      "Irish Sea" = c(-3.7, 53.75), 
      "Minches & W Scotland" = c(-7.1, 56.2),
      stop("option not coded yet")
    )

    x <- pos[1]
    y <- pos[2]
    z <- switch(status[regionID, "trend"], stable = 0, up = pi/6, down = -pi/6)
    adj <- 0.5 * c(-cos(z), -sin(z), cos(z), sin(z))
    arrows(x + adj[1], y + adj[2], x + adj[3], y + adj[4], col = "black", length = 0.1, lwd = 2)
  })
  
  box(col = "black")
  
  plyr::l_ply(CP2.borders @ polygons, function(y) 
    plyr::l_ply(y @ Polygons, function(x) 
      lines(x @ coords[, 1], x @ coords[, 2], lwd = 1)))

  plyr::l_ply(MSFD.borders @ polygons, function(y) 
    plyr::l_ply(y @ Polygons, function(x) 
      lines(x @ coords[, 1], x @ coords[, 2], lwd = 2)))
}


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


tableTimeSeries <- function(assessment, type = c("trend", "status"), asNumbers = TRUE, 
                            purpose = c("OSPAR", "CSEMP")) {
  
  purpose <- match.arg(purpose)
  type <- match.arg(type)

  timeSeries <- assessment$timeSeries
  
  ok <- with(timeSeries, switch(
    type, 
    status = ! colour %in% "black" & ! shape %in% "open_circle", 
    trend = nypos >= 5
  ))
  timeSeries <- timeSeries[ok, ]
  timeSeries <- droplevels(timeSeries)
  
  if (type == "status") {
    colLevels <- c("blue", "green", "orange", "red")
    colLevels <- colLevels[colLevels %in% c("blue", timeSeries$colour, "red")]
  }
  
  timeSeries <- within(timeSeries, {
    
    response <- switch(
      type, 
      status = factor(colour, levels = colLevels), 
      trend = factor(shape, levels = c("upward_triangle", "filled_circle", "downward_triangle"), 
                     labels = c("upward trend", "no trend", "downward trend"))
    )
    
    stopifnot(!is.na(response)) 
    
    metaRegion <- switch(
      purpose, 
      OSPAR = ordered(OSPARregion, levels = as.character(1:4)),
      CSEMP = ordered(MSFDregion, levels = c("Greater North Sea", "Celtic Seas"))
    )
  })

  timeSeries <- split(timeSeries, timeSeries$metaRegion, drop = TRUE)

  out <- tableTimeSeriesCalc(timeSeries, asNumbers)

  out <- mapply(function(x, y) {
    x$metaRegion <- y
    x
  }, out, names(out), SIMPLIFY = FALSE)

  out <- do.call(rbind, out)
      
  out <- within(out, metaRegion[duplicated(metaRegion)] <- NA)
  out <- out[c("metaRegion", setdiff(names(out), "metaRegion"))]
  
  names(out)[1:2] <- switch(
    purpose,
    OSPAR = c("OSPAR region", "MIME region"),
    CSEMP = c("MSFD sub-region", "biogeographic region")
  )

  row.names(out) <- NULL
  out[-nrow(out), ]
}  


tableTimeSeriesCalc <- function(timeSeries, asNumbers) {

  timeSeries <- lapply(timeSeries, function(x) within(x, region <- droplevels(region)))

  detLevels <- levels(timeSeries[[1]]$determinand)
  nDet <- length(detLevels)
  
    
  tab.region <- sapply(timeSeries, simplify = FALSE, FUN = function(data) {
    
    tab <- with(data, table(region, response, determinand))
    
    # now convert to data frame for printing
    
    tab <- sapply(rev(dimnames(tab)$region), simplify = FALSE, FUN = function(x) {
      
      if (nDet > 1) {
        out <- as.data.frame.matrix(tab[x, , ])
        out$total <- as.integer(rowSums(out))
      } else {
        out <- as.data.frame(tab[x, , ])
        names(out) <- detLevels
      }

      if (!asNumbers)
        out[] <- lapply(out[], function(x) {
          if (sum(x) == 0) return(rep(NA, length(x)))
          prop <- as.integer(round(100 * x / sum(x)))
          if (sum(prop) == 101) {
            id <- which.max(prop * sum(x) / 100 - x)
            prop[id] <- prop[id] - 1L
          }
          if (sum(prop) == 99) {
            id <- which.min(prop * sum(x) / 100 - x)
            prop[id] <- prop[id] + 1L
          }
          stopifnot(sum(prop) == 100)
          prop
        })
      
      out <- within(out, {
        status <- row.names(out)
        region <- x
      })
      
      row.names(out) <- NULL
      out[nrow(out)+1, ] <- NA
      
      varID <- intersect(c("region", "status", dimnames(tab)$determinand, "total"), names(out))
      out[varID]
    })
    
    tab <- do.call(rbind, tab)
    tab <- within(tab, region[duplicated(region)] <- NA)
    tab
  })
  
  tab.metaRegion <- sapply(timeSeries, simplify = FALSE, FUN = function(data) {
    
    tab <- with(data, table(response, determinand))
    
    if (nDet > 1) {
      out <- as.data.frame.matrix(tab)
      out$total <- as.integer(rowSums(out))
    } else {
      out <- as.data.frame(tab[ , ])
      names(out) <- detLevels
    }

    if (!asNumbers)
      out[] <- lapply(out[], function(x) {
        if (sum(x) == 0) return(rep(NA, length(x)))
        prop <- as.integer(round(100 * x / sum(x)))
        if (sum(prop) == 101) {
          id <- which.max(prop * sum(x) / 100 - x)
          prop[id] <- prop[id] - 1L
        }
        if (sum(prop) == 99) {
          id <- which.min(prop * sum(x) / 100 - x)
          prop[id] <- prop[id] + 1L
        }
        stopifnot(sum(prop) == 100)
        prop
      })
    
    out <- within(out, {
      status <- row.names(out)
      region <- "total"
    })
    row.names(out) <- NULL
    out[nrow(out)+1, ] <- NA
    
    varID <- intersect(c("region", "status", dimnames(tab)$determinand, "total"), names(out))
    tab <- out[varID]
    
    tab <- within(tab, region[duplicated(region)] <- NA)
    tab
  })

  out <- mapply(rbind, tab.region, tab.metaRegion, SIMPLIFY = FALSE)

  out
}
  





tableTrend <- function(assessment, purpose = c("OSPAR", "CSEMP"), 
                       type = c("region", "determinand", "region determinand"), detLabel) {
  
  purpose <- match.arg(purpose)
  type <- match.arg(type)
  
  infile <- paste("trend results by", type)
  
  data <- readFromGenstat(infile, assessment)
  
  data <- within(data, {
    change <- 100 * (exp(est / 100) - 1)
    change.upper <- 100 * (exp(upper / 100) - 1)
    change.lower <- 100 * (exp(lower / 100) - 1)
  })
  
  data <- switch(
    type, 
    determinand = data[order(data$determinand), ],
    region = data[order(data$region, decreasing = TRUE), ], 
    "region determinand" = {
      data <- data[order(data$determinand), ]
      data[order(data$region, decreasing = TRUE), ]
    }
  )
  
  if (type == "region determinand")  {
    data <- split(data, data$region, drop = TRUE)
    
    data <- sapply(rev(data), simplify = FALSE, FUN = function(x) {
      x[nrow(x)+1, ] <- NA
      x
    })
    
    data <- do.call(rbind, data)
    data <- within(data, region[duplicated(region)] <- NA)  
    data <- data[-nrow(data), ]
  }

  id <- c("region", "determinand", "est", "se", "lower", "upper", "change", "change.lower", "change.upper")
  data <- data[intersect(id, names(data))]
  
  if ("region" %in% names(data)) {
    id <- match("region", names(data))
    names(data)[id] <- switch(purpose, OSPAR = "MIME region", CSEMP = "biogeographic region")
  }
  
  if ("determinand" %in% names(data)) {
    id <- match("determinand", names(data))
    names(data)[id] <- detLabel
  }

  id <- switch(type, "region determinand" = 3:9, 2:8)
  names(data)[id] <- c("trend", "se", "lower", "upper", "% yearly change", 
                      "%yc lower", "%yc upper")
  
  row.names(data) <- NULL
  
  
  data
}  


tableStatus <- function(assessment, AC, purpose = c("OSPAR", "CSEMP"), 
                        type = c("region", "determinand", "region determinand"), detLabel) {
  
  purpose <- match.arg(purpose)
  type <- match.arg(type)

  infile <- switch(type, "region determinand" = "region by determinand", type)

  data <- getStatusFromGenstat(assessment, infile, AC)

  
  data <- within(data, {
    conc <- exp(est)
    conc.upper <- exp(upper)
  })

  
  data <- switch(
    type, 
    determinand = data[order(data$determinand), ],
    region = data[order(data$region, decreasing = TRUE), ], 
    "region determinand" = {
      data <- data[order(data$determinand), ]
      data[order(data$region, decreasing = TRUE), ]
    }
  )
  

  if (type == "region determinand")  {
    data <- split(data, data$region, drop = TRUE)
    
    data <- sapply(rev(data), simplify = FALSE, FUN = function(x) {
      x[nrow(x)+1, ] <- NA
      x
    })
  
    data <- do.call(rbind, data)
    data <- within(data, region[duplicated(region)] <- NA)  
    data <- data[-nrow(data), ]
  }
  
  id <- c("region", "determinand", "est", "se", "upper", "conc", "conc.upper")
  data <- data[intersect(id, names(data))]
  
  if ("region" %in% names(data)) {
    id <- match("region", names(data))
    names(data)[id] <- switch(purpose, OSPAR = "MIME region", CSEMP = "biogeographic region")
  }
  
  if ("determinand" %in% names(data)) {
    id <- match("determinand", names(data))
    names(data)[id] <- detLabel
  }

  id <- switch(type, "region determinand" = 3:7, 2:6)
  names(data)[id] <- c("status", "se", "upper", "concentration", "conc upper")
  
  row.names(data) <- NULL
  
  data
}  




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


plotTrendImposex <- function(assessment) {
  require(lattice)
  path <- assessment$webpath
  if (is.null(assessment$trend)) return(invisible())
  data <- assessment$trend
  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) {
      pch <- data$shape[subscripts]
      pch[pch == "filled_circle"] <- "21"
      pch[pch == "upward_triangle"] <- "24"
      pch[pch == "downward_triangle"] <- "25"
      pch <- as.numeric(pch)
      panel.abline(v = 1, col = grey(0.4))
      panel.stripplot(x, y, jitter.data = TRUE, pch = pch, col = "black")
    })
  out
}


plotStatusImposex <- function(assessment, AC, greenID = "green") {
  
  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)]]
  
  if (max(data$response, na.rm = TRUE) < 2)
    at <- c(0, 0.5, 0.8, 1, 1.2, 1.5)
  else
    at <- c(0, 0.2, 0.5, 1, 2, 5)

  labels <- as.character(at)
  
  xlab <- switch(AC, 
                 CONC = "mean VDS in final year", 
                 paste("mean VDS relative to", AC))
  
  out <- stripplot(
    region ~ response | determinand,  
    data = data,
    scales = list(alternating = FALSE, x = list(at = sqrt(at), labels = labels)),
    xlab = xlab,
    panel = function(x, y, subscripts) {
      col <- data$colour[subscripts]
      col[col == "orange"] <- "goldenrod"
      if (! AC %in% "CONC") 
        panel.abline(v = 1, col = grey(0.4))
      panel.stripplot(x, y, jitter.data = TRUE, pch = 16, col = col)
    })
  out
}


plotStatusResultsImposex <- function(
  assessment, 
  type = c("region", "determinand", "region by determinand", "determinand by region"), 
  AC, greenID = "green") {
  
  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$est, 1)))
    if (xlim[2] < 2)
      at <- c(0, 0.5, 0.8, 1, 1.2, 1.5)
    else
      at <- c(0, 0.2, 0.5, 1, 2, 5)
    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, 
      drop.unused.levels = FALSE, 
      scales = list(alternating = FALSE, x = list(at = sqrt(at), labels = labels)),
      panel = function(x, y, subscripts) {
        if (AC != "CONC") panel.abline(v = 1, col = grey(0.7))
        if (nrow(data[subscripts, ]) == 0) return()
        if (AC == "CONC")
          with(data[subscripts, ], lsegments(lower, y, upper, y, col = "black"))
        else
          with(data[subscripts, ], lsegments(x, y, upper, y, col = "black"))
        col <- data$col[subscripts]
        col[col %in% "orange"] <- "goldenrod"
        lpoints(x, y, col = "black", fill = col, pch = 21, cex = 1.5)
      })
  }

  data <- getStatusFromGenstat(assessment, type, AC, referenceValue = 1)
  
  if (type %in% c("region by determinand", "determinand by region"))
    data <- data[with(data, paste(region, determinand)) %in% regDetExists, ]
  
  switch(
    type, 
    region = plotEngine(data, as.formula(region ~ est), AC),
    determinand = plotEngine(data, as.formula(determinand ~ est), AC),
    "region by determinand" = plotEngine(data, as.formula(region ~ est | determinand), AC),
    "determinand by region" = plotEngine(data, as.formula(determinand ~ est | region), AC)
  )
}


plotTrendResultsImposex <- 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.4))
        with(data[subscripts, ], lsegments(lower, y, upper, y, col = "black"), lwd= 1.5)
        pch <- with(data[subscripts, ], ifelse(upper < 1, 25, ifelse(lower > 1, 24, 21)))
        lpoints(x, y, pch = pch, cex = 1.5, col = "black", fill = "white", lwd = 2)
      })
  }

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



plotMSTAT.OSPAR <- function(assessment, OSPARregion, subset = c("status", "trend"), MSTATid = "MSTAT") {
  
  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", MSTATid)])
  
  stations <- within(stations, {
    response <- get(MSTATid)
    bg <- factor(response, 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)
}



tableTrendImposex <- function(assessment, purpose = c("OSPAR", "CSEMP"), 
                       type = c("region", "determinand", "region determinand"), detLabel) {
  
  purpose <- match.arg(purpose)
  type <- match.arg(type)
  
  infile <- paste("trend results by", type)
  
  data <- readFromGenstat(infile, assessment)
  
  data <- switch(
    type, 
    determinand = data[order(data$determinand), ],
    region = data[order(data$region, decreasing = TRUE), ], 
    "region determinand" = {
      data <- data[order(data$determinand), ]
      data[order(data$region, decreasing = TRUE), ]
    }
  )
  
  if (type == "region determinand")  {
    data <- split(data, data$region, drop = TRUE)
    
    data <- sapply(rev(data), simplify = FALSE, FUN = function(x) {
      x[nrow(x)+1, ] <- NA
      x
    })
    
    data <- do.call(rbind, data)
    data <- within(data, region[duplicated(region)] <- NA)  
    data <- data[-nrow(data), ]
  }
  
  id <- c("region", "determinand", "est", "se", "lower", "upper", "change", "change.lower", "change.upper")
  data <- data[intersect(id, names(data))]
  
  if ("region" %in% names(data)) {
    id <- match("region", names(data))
    names(data)[id] <- switch(purpose, OSPAR = "MIME region", CSEMP = "biogeographic region")
  }
  
  if ("determinand" %in% names(data)) {
    id <- match("determinand", names(data))
    names(data)[id] <- detLabel
  }
  
  id <- switch(type, "region determinand" = 3:6, 2:5)
  names(data)[id] <- c("trend", "se", "lower", "upper")
  
  row.names(data) <- NULL
  
  
  data
}  


tableStatusImposex <- function(assessment, AC, purpose = c("OSPAR", "CSEMP"), 
                               type = c("region", "determinand", "region determinand"), detLabel) {
  
  purpose <- match.arg(purpose)
  type <- match.arg(type)
  
  infile <- switch(type, "region determinand" = "region by determinand", type)
  
  data <- getStatusFromGenstat(assessment, infile, AC)
  
  
  data <- within(data, {
    conc <- est ^ 2
    conc.upper <- upper ^ 2
  })
  
  
  data <- switch(
    type, 
    determinand = data[order(data$determinand), ],
    region = data[order(data$region, decreasing = TRUE), ], 
    "region determinand" = {
      data <- data[order(data$determinand), ]
      data[order(data$region, decreasing = TRUE), ]
    }
  )
  
  
  if (type == "region determinand")  {
    data <- split(data, data$region, drop = TRUE)
    
    data <- sapply(rev(data), simplify = FALSE, FUN = function(x) {
      x[nrow(x)+1, ] <- NA
      x
    })
    
    data <- do.call(rbind, data)
    data <- within(data, region[duplicated(region)] <- NA)  
    data <- data[-nrow(data), ]
  }
  
  id <- c("region", "determinand", "est", "se", "upper", "conc", "conc.upper")
  data <- data[intersect(id, names(data))]
  
  if ("region" %in% names(data)) {
    id <- match("region", names(data))
    names(data)[id] <- switch(purpose, OSPAR = "MIME region", CSEMP = "biogeographic region")
  }
  
  if ("determinand" %in% names(data)) {
    id <- match("determinand", names(data))
    names(data)[id] <- detLabel
  }
  
  id <- switch(type, "region determinand" = 3:7, 2:6)
  names(data)[id] <- c("status", "se", "upper", "VDS", "VDS upper")
  
  row.names(data) <- NULL
  
  data
}  
