code for panel graphics

profilePlot.R

### code to make decent-looking, simple panel graphics

profilePlot <- function(xdata = dated$OM.pct, ydata = dated$midpoint, siteList = dated$site, 
                        xlabel = "Organic matter content (decimal fraction)",
                        ylabel = "depth (cm)",
                        rowNos = 4, colNos = 4,
                        xTiers = 1,
                        filename = "OM_profiles.png",
                        figylims  = c(80, 0),
                        figxlims = c(0, 0.8)
) {

  
  
  text.size   <- 0.85
  cex.size    <- 0.65
  p           <- 0.5    # point marker size
  t1          <- 0.3    # tck, tick mark direction
  t2          <- 0.15   # tcl, tick mark length
  
  
  left_buffer <- 0.081
  bot_buffer  <- 0.11
  xgap        <- 0.005
  ygap        <- 0.005
  r           <- rowNos # number of rows
  cols        <- colNos # number of columns
  extra.bit.for.polygon <- abs(diff(figylims) * 0.13) # formerly a flat 15
  
  xlabs <- pretty(figxlims, n = 4)
    xlabs <- xlabs[xlabs < max(figxlims)]
  
  # need to be able to accomodate depth profiles (inverted y axis without negative signs)
  # while also printing other plots 
  ifelse(figylims[2] > figylims[1],
         figylims_plot  <- c(figylims[1], figylims[2] + extra.bit.for.polygon),
         figylims_plot  <- c(figylims[1], figylims[2] - extra.bit.for.polygon)  
        )
  
  # link site names, numbers
  site_names <- list(
    'AND' = "Oyster R.",
    'CIC' = "Canfield",
    'EBD' = "Banca",
    'GF' = "Greens Farms",
    'GPE' = "Gulf Pond East",
    'HV1' = "Harborview 1",
    'HV2' = "Harborview 2",
    'HV4' = "Harborview 3",
    'MC' = "Mar. Cons.",
    'OC' = "Otter Creek", 
    'PBA' = "Pelham Bay",
    'PBB' = "Hutchinson R.",
    'SMC' = "Saw Mill Creek",
    'SPR' = "Spring Creek",
    'SWD' = "Sherwood",
    'VCR' = "Village Creek"
  )
  
  site_nos <- list(
    'GF' = "11",
    'GPE' = "12",
    'AND' = "13",
    'EBD' = "14",
    'MC' = "6",
    'VCR' = "7",
    'CIC' = "9",
    'SWD' = "10",
    'OC' = "5", 
    'HV1' = "8A",
    'HV2' = "8B",
    'HV4' = "8C",
    'SMC' = "1",
    'SPR' = "2",
    'PBB' = "3",
    'PBA' = "4"
    )
  
  
  # store fig coordinates, incorporating buffers & gaps:
  coords <- data.frame(x1 = c(left_buffer, left_buffer + 1:(cols - 1)*((1-left_buffer)/cols)), 
                       x2 = c(left_buffer + 1:cols*((1-left_buffer)/cols) - xgap), 
                       y1 = c(rep(bot_buffer, times = cols), 
                              rep(bot_buffer + ((1-bot_buffer)/r), times = cols), 
                              rep(bot_buffer + 2*((1-bot_buffer)/r), times = cols), 
                              rep(bot_buffer + 3*((1-bot_buffer)/r), times = cols)), 
                       y2 = c(rep(bot_buffer + ((1-bot_buffer)/r) - ygap, times = cols),  
                              rep(bot_buffer + 2*((1-bot_buffer)/r) - ygap, times = cols), 
                              rep(bot_buffer + 3*((1-bot_buffer)/r) - ygap, times = cols), 
                              rep(0.995, times = cols))
  )
  
  plotTier1  <- function(site){
    output <- match.call(plot, call("plot", x = xdata[siteList == site], 
                                    y = ydata[siteList == site], 
                                    ylab = '' , xlab = '' , type = 'p', xlim = figxlims, ylim = figylims_plot, 
                                    pch = 19, yaxt = 'n', xaxt = 'n', xaxs = 'i', yaxs = 'i', cex = p,
                                    cex.lab = cex.size, cex.axis = cex.size, tck = t1, tcl = t2)
    )
  }
  
  
  png(filename = filename, width = 140, height = 180, units = "mm", res = 300)
  
  for(i in 1:16) {
    par(fig = c(coords$x1[as.numeric(i)], coords$x2[as.numeric(i)], coords$y1[as.numeric(i)], coords$y2[as.numeric(i)]), mar = c(0, 0, 0, 0))
    
    site <- names(site_nos[i])
    # make the plot
    eval(plotTier1(site))
    # add annotations
    axis(1, at = xlabs, 
         labels = F, cex.axis = cex.size, las = 1, tck = t1, tcl = t2)
    if(figylims[2] > figylims[1]){
           polygon(c(figxlims, figxlims[2], figxlims[1]), y = c(max(figylims_plot), max(figylims_plot), max(figylims), max(figylims)), 
              density = NA, col = "lightgray", border = "black")
           text(mean(figxlims), figylims_plot[2] - 0.5*(extra.bit.for.polygon), 
              as.character(site_nos[names(site_nos) == site]), cex = text.size)
    } else {
           polygon(c(figxlims, figxlims[2], figxlims[1]), y = c(min(figylims), min(figylims), min(figylims_plot), min(figylims_plot)), 
                   density = NA, col = "lightgray", border = "black")
           text(mean(figxlims), figylims_plot[2] + 0.5*(extra.bit.for.polygon), 
                as.character(site_nos[names(site_nos) == site]), cex = text.size)
    }
    
    axis(2, labels = F, cex.axis = cex.size, las = 1, tck = t1, tcl = t2)
    axis(4, labels = F, cex.axis = cex.size, las = 1, tck = t1, tcl = t2)
    
    if(i %in% c(1, 5, 9, 13)){ # add y axis to rightmost plots
       if(figylims[2] > figylims[1]){
            mtext(text = axTicks(2)[axTicks(2) <= max(figylims)], at = axTicks(2)[axTicks(2) <= max(figylims)], 
                  side = 2, cex = cex.size, las = 1, line = 0.2)
            } else {
              mtext(text = axTicks(2)[axTicks(2) >= min(figylims)], at = axTicks(2)[axTicks(2) >= min(figylims)], 
                  side = 2, cex = cex.size, las = 1, line = 0.2)
            }
    }
    if(i %in% c(1, 3)){ # add x axis to bottom row, alternating plots
      mtext(text = xlabs, 
            at = xlabs, 
            side = 1, cex = cex.size, las = 1, line = 0.2)
    }
    par(new = T) # so that the subsequent plot will be superimposed
  }
  
  
  
  # now add centered axis labels
  par(fig = c(0, 1, bot_buffer, 1))
  mtext(ylabel, side = 2, line = -1)
  par(new = T)
  par(fig = c(left_buffer, 1, 0, 1))
  mtext(xlabel, side = 1, line = -2)
  dev.off()
  
  print(paste0(filename, " saved to ", getwd()))
}