In imarpe/imarpe: Tools for graphics and summary reports for IMARPE

\pagenumbering{gobble}

knitr::opts_chunk$set(echo = TRUE)

allDates <- x$allDates
allMarks <- x$allMarks

firstTable <- matrix(data = NA, nrow = ncol(x$catchByWeek), ncol = 6, 
                     dimnames = list(c(colnames(x$catchByWeek)), 
                                     c("$N (x10~^6)$", "Biomasa (t)", "$Juv N (x10~^6)$", "%juvN", "Juv Peso (t)", "%juv peso (t)")))
for(i in seq(ncol(x$catchByWeek))){

  tempValues <- as.numeric(x$catchByWeek[,i])
  juvIndex <- allMarks < 12

  if(i == 1){
    firstTable[i, 1:2] <- c(sum(tempValues), 
                            sum(x$a*allMarks^x$b*tempValues))
  }

  firstTable[i, 3:6] <- c(sum(tempValues[juvIndex]), 
                          sum(tempValues[juvIndex])/sum(tempValues)*100, 
                          sum(x$a*allMarks[allMarks < 12]^x$b*tempValues[juvIndex]), 
                          sum(x$a*allMarks[allMarks < 12]^x$b*tempValues[juvIndex])/sum(x$a*allMarks^x$b*tempValues)*100)
}
firstTable <- round(firstTable, 0)

pander::pander(firstTable, align = "c", keep.line.breaks = TRUE, split.table = 120, missing = "")

secondTable <- matrix(data = character(), nrow = 2, ncol = 6, 
                      dimnames = list(c("", paste("Al", format(x$endDate, format = "%d/%m/%y"))), 
                                      c("Pesca acum N", "Pesca acum P", "Juv N", "Juv Peso", "Población---", "---remanente")))

secondTable[1,] <- c("$N (x10~^6)$", "Peso (t)", "En N (%)", "En peso (%)", "$N (x10~^6)$", "B (millones t)")

sobrevivientes <- x$projByWeek[,ncol(x$projByWeek)]

catchInfo <- x$getInfo(x = rowSums(x$catchByDay), millionT = FALSE, allMarks = x$allMarks, a = x$a, b = x$b)

secondTable[2, 1:4] <- c(round(catchInfo[c(1, 2)], 0), round(catchInfo[c(5, 6)], 2))

secondTable[2, 5] <- as.character(round(sum(sobrevivientes), 0))
secondTable[2, 6] <- as.character(round(sum(x$a*x$allMarks^x$b*sobrevivientes)/1e6, 2))

pander::pander(secondTable, keep.line.breaks = TRUE, split.table = 250)

# Subplot 1-1
xlim <- c(x$startDate, x$endSeasonDate)
y1 <- colSums(sweep(x = x$projByDay[,-1], MARGIN = 1, STATS = x$a*x$allMarks^x$b, FUN = "*"))
y2 <- cumsum(colSums(sweep(x = x$catchByDay, MARGIN = 1, STATS = x$a*x$allMarks^x$b, FUN = "*")))
pretty_y1 <- pretty(c(0, 1.01 * y1))
pretty_y2 <- pretty(c(0, cuota))
by_y1 <- (range(pretty_y1)[2]-range(pretty_y1)[1])/(length(pretty_y1)-1)
by_y2 <- (range(pretty_y2)[2]-range(pretty_y2)[1])/(length(pretty_y2)-1)
# ylim1 <- axisPlot$ylim_p1$axis2_1
# ylim2 <- axisPlot$ylim_p1$axis4_1
ylim1 <- c(range(pretty_y1), by_y1)
ylim2 <- c(range(pretty_y2), by_y2)

# layout(matrix(c(1, 1, 2, 3), nrow = 2, byrow = TRUE))
par(mar = c(0, 0, 0, 0), oma = c(6, 4, 1, 4), yaxs = "i", xaxs = "i", mfrow = c(3, 1))

plot(1, 1, axes = FALSE, xlab = NA, ylab = NA, xlim = xlim, ylim = ylim1[1:2])

lines(x$allDates, cumsum(colSums(sweep(x = x$catchByDay, MARGIN = 1, STATS = x$a*x$allMarks^x$b, FUN = "*")))*ylim1[2]/ylim2[2],
      col = "blue", lwd = 2)
lines(x$allDates, colSums(sweep(x = x$projByDay[,-1], MARGIN = 1, STATS = x$a*x$allMarks^x$b, FUN = "*")), 
      col = "red", lwd = 2)

mtext(text = "Biomasa (millones t)", side = 2, line = 2.5)
mtext(text = "Captura (millones t)", side = 4, line = 3)

legend("topright", legend = c("Biomasa", "Captura acum"), bty = "n", col = c("red", "blue"), lwd = 2, cex = 1.2)

yLabs <- seq(ylim1[1], ylim1[2], ylim1[3])
axis(side = 2, at = yLabs, labels = yLabs/1e6, las = 2)

yLabs2 <- seq(ylim2[1], ylim2[2], ylim2[3])
yLabs2[1] <- NA
axis(side = 4, at = seq(ylim1[1], ylim1[2], length.out = length(yLabs2)), labels = yLabs2/1e6, las = 2)
box()

# Subplot 1-2
juvN <- apply(x$catchByDay, 2, function(z) sum(z[x$allMarks < 12])/sum(z))*100

juvB <- sweep(x = x$catchByDay, MARGIN = 1, STATS = x$a*x$allMarks^x$b, FUN = "*")
juvB <- apply(juvB, 2, function(z) sum(z[x$allMarks < 12])/sum(z))*100

xValues <- as.Date(colnames(x$catchByDay))
xValues2 <- c(min(xValues) - 1e3, xValues, max(xValues) + 1e3)

xlim <- c(x$startDate, x$endSeasonDate)
ylim1 <- axisPlot$ylim_p1$axis4_2
ylim2 <- axisPlot$ylim_p1$axis2_3

plot(1, 1, axes = FALSE, xlab = NA, ylab = NA, xlim = xlim, ylim = ylim1[1:2])
if(!is.null(juvLimits$number)){
  polygon(x = c(xValues2, rev(xValues2)), 
          y = rep(juvLimits$number[1] + juvLimits$number[2]*c(-1, 1), each = length(xValues2)), 
          border = FALSE, col = "gray80")
  abline(h = juvLimits$number[1], lty = "dotted")
} 
abline(v = x$startDate, lty = "dashed", col = "gray50")

lines(xValues, juvN, col = "blue", lty = "solid", lwd = 2, type = "o", pch = 19, cex = 0.8)

mtext(text = "% juveniles (número)", side = 3, line = -2, adj = 0.99)
axis(side = 4, at = seq(ylim1[1], ylim1[2], ylim1[3]), labels = paste(seq(ylim1[1], ylim1[2], ylim1[3]), "%"), las = 2)
box()

# Subplot 1-3
plot(1, 1, axes = FALSE, xlab = NA, ylab = NA, xlim = xlim, ylim = ylim2[1:2])
if(!is.null(juvLimits$weight)){
  polygon(x = c(xValues2, rev(xValues2)), 
          y = rep(juvLimits$weight[1] + juvLimits$weight[2]*c(-1, 1), each = length(xValues2)), 
          border = FALSE, col = "gray80")
  abline(h = juvLimits$weight[1], lty = "dotted")
} 
abline(v = x$startDate, lty = "dashed", col = "gray50")
lines(xValues, juvB, col = "blue", lty = "solid", lwd = 2, type = "o", pch = 19, cex = 0.8)
mtext(text = "% juveniles (peso)", side = 3, line = -2, adj = 0.99)
axis(side = 2, at = seq(ylim2[1], ylim2[2], ylim2[3]), labels = paste(seq(ylim2[1], ylim2[2], ylim2[3]), "%"), las = 2)
box()

xValues <- seq(xlim[1], xlim[2], by = "day")
axis(side = 1, at = xValues, labels = NA, las = 2, tcl = -0.25)
xValues2 <- xValues[is.element(day(xValues), seq(0, 30, 5))]
axis(side = 1, at = xValues2, labels = as.character(format(xValues2, format = "%d/%m/%Y")), las = 2)

\pagebreak

Captura y proyecciones

# Tabla diaria
index <- match(x$endDate, x$allDates)
output <- x$getInfo(x = x$catchByDay[,index], allMarks = x$allMarks, a = x$a, b = x$b, nDecimalsBiomass = 3)
dailyTable <- matrix(data = output, ncol = 1, 
                     dimnames = list(names(output), format(x$endDate, "%d/%m/%Y")))

# Tabla de datos acumulados
output <- cbind(rowSums(x$catchByDay), x$surveyVector, x$projByWeek[,ncol(x$projByWeek)])

weeklyTable <- apply(output, 2, x$getInfo, allMarks = x$allMarks, a = x$a, b = x$b)

# Bind tables
outputTable <- cbind.data.frame(dailyTable, weeklyTable)
colnames(outputTable) <- c(paste0("Pesca día\n", format(x$endDate, "%d/%m/%Y")), 
                           "Acumulado\nPesca", "Población inicio\nCrucero", "\nSobrevivientes")

outputTable[is.na(outputTable)] <- " "

pander::pander(outputTable, justify = "center", keep.line.breaks = TRUE, split.table = 140)

index <- match(x$endDate, x$allDates)
output <- cbind(x$allMarks, x$catchByDay[,index], x$catchByDay[,index]*x$a*x$allMarks^x$b)
colnames(output) <- c("Talla", "Abundancia", "Biomasa")

# Get axis limits
xlim <- range(x$allMarks)
y1 <- output[,2]
y2 <- output[,3]
pretty_y1 <- pretty(c(0, 1.01 * y1))
pretty_y2 <- pretty(c(0, 1.01 * y2))
by_y1 <- (range(pretty_y1)[2]-range(pretty_y1)[1])/(length(pretty_y1)-1)
by_y2 <- (range(pretty_y2)[2]-range(pretty_y2)[1])/(length(pretty_y2)-1)
# ylim1 <- axisPlot$ylim_p2$axis2
# ylim2 <- axisPlot$ylim_p2$axis4
ylim1 <- c(range(pretty_y1), by_y1)
ylim2 <- c(range(pretty_y2), by_y2)

# Make empty plot
par(mar = c(3, 5, 1, 4), xaxs = "i", yaxs = "i")
plot(1, 1, pch = NA, axes = FALSE, xlim = xlim, ylim = ylim1[1:2], xlab = NA, ylab = NA)

if(sum(output[,2],  na.rm = TRUE)){
  # Make plot and lines Abundance
  tempLines <- output[,2]
  tempLines[tempLines < 0.001] <- NA
  lines(x$allMarks, tempLines, lwd = 2, col = "blue4")

  # Make plot and lines Biomass
  tempLines <- output[,3]/ylim2[2]*ylim1[2]
  tempLines[tempLines < 0.001] <- NA
  lines(x$allMarks, tempLines, type = "o", lwd = 2, col = "red4", pch = 16)
}

# Add left axis
axis(side = 2, at = seq(ylim1[1], ylim1[2], ylim1[3]), las = 2, cex.axis = 0.8)

# Add right axis
yLabs <- seq(ylim2[1], ylim2[2], ylim2[3])
axis(side = 4, at = seq(ylim1[1], ylim1[2], length.out = length(yLabs)), labels = yLabs, las = 2, cex.axis = 0.8)

# Add juvenile limit
abline(v = 12, lty = "dotted", col = "red", lwd = 2)

# Add legend
legend("topleft", legend = c("Abundancia", "Peso"), col = c("blue4", "red4"), 
       pch = c(NA, 16), lty = rep("solid", 2), lwd = 2)

# Add date text
mtext(text = format(x$endDate, "%d/%m/%Y"), side = 3, line = -2, adj = .99)

# Add axis labels 
mtext(text = "Peso (t)", side = 4, line = 3)
mtext(text = "Abundancia (millones ind)", side = 2, line = 3)
mtext(text = "Longitud total (cm)", side = 1, line = 2, outer = TRUE)

# Add axis X label
axis(side = 1, at = x$allMarks, labels = NA, tcl = -0.25)
axis(side = 1, at = seq(xlim[1], xlim[2]), cex.axis = 0.8)

box()

output <- cbind(x$surveyVector, rowSums(x$catchByDay), x$projByWeek[,ncol(x$projByWeek)])
colnames(output) <- c("Crucero", "Capturas", "Sobrevivientes")

output[output < 1e-5] <- NA

allMarks <- x$allMarks

cols <- c("red3", "blue3", "green4")

par(mar = rep(0, 4), oma = c(4, 6.5, 1, 5), xaxs = "i", yaxs = "i", mfrow = c(3, 1))

for(i in seq(ncol(output))){

  xlim <- range(allMarks)
  if(i == 2){
    y1 <- output[,i]
    y2 <- output[,i]*x$a*x$allMarks^x$b
  }else{
    y1 <- c(output[,1], output[,ncol(output)])
    y2 <- c(output[,1]*x$a*x$allMarks^x$b, output[,ncol(output)]*x$a*x$allMarks^x$b)
  }
  pretty_y1 <- pretty(c(0, 1.01 * y1))
  pretty_y2 <- pretty(c(0, 1.01 * y2))
  by_y1 <- (range(pretty_y1)[2]-range(pretty_y1)[1])/(length(pretty_y1)-1)
  by_y2 <- (range(pretty_y2)[2]-range(pretty_y2)[1])/(length(pretty_y2)-1)
  # ylim1 <- axisPlot$ylim_p3[[paste0("axis2_", i)]]
  # ylim2 <- axisPlot$ylim_p3[[paste0("axis4_", i)]]
  ylim1 <- c(range(pretty_y1), by_y1)
  ylim2 <- c(range(pretty_y2), by_y2)

  plot(1, 1, pch = NA, axes = FALSE, xlim = xlim, ylim = ylim1[1:2], xlab = NA, ylab = NA)

  lines(allMarks, output[,i], col = cols[i], lwd = 2)
  lines(allMarks, output[,i]*x$a*x$allMarks^x$b*(ylim1[2]/ylim2[2]), col = cols[i], lwd = 2, lty = "dashed")

  abline(v = 12, col = "red", lwd = 2, lty = "dotted")

  mtext(text = colnames(output)[i], side = 3, line = -1.5, adj = 0.01)

  index <- allMarks < 12

  juvValue <- round(sum(output[index, i], na.rm = TRUE)/sum(output[,i], na.rm = TRUE)*100, 0)
  mtext(text = paste("JuvN =", juvValue, "%"), side = 3, line = -1.5, adj = 0.99)

  juvValue <- round(sum((output[,i]*x$a*x$allMarks^x$b)[index], na.rm = TRUE)/sum(output[,i]*x$a*x$allMarks^x$b, na.rm = TRUE)*100, 0)
  mtext(text = paste("JuvB =", juvValue, "%"), side = 3, line = -3, adj = 0.99)

  if(i == ncol(output)){
    axis(side = 1, at = allMarks, labels = NA, tcl = -0.25)
    axis(side = 1, at = seq(xlim[1], xlim[2]))
  }else if(i == 2){
    legend("bottomright", legend = c("Abundancia", "Biomasa"), col = cols[i], lty = c("solid", "dashed"), 
           lwd = 2, pch = NA, bty = "n")
  }

  axis(side = 2, at = seq(ylim1[ifelse(i == ncol(output), 1, 3)], ylim1[2], ylim1[3]), las = 2)

  yLabs2 <- seq(ylim2[1], ylim2[2], ylim2[3])
  if(i != 2){
    yLabs2 <- paste(yLabs2/1e6, "m")
  }

  if(i != ncol(output)){
    yLabs2[1] <- NA
  }

  axis(side = 4, at = seq(ylim1[1], ylim1[2], length.out = length(yLabs2)), 
       labels = yLabs2, las = 2)

  box()
}

mtext(text = "Longitud total (cm)", side = 1, outer = TRUE, line = 3)
mtext(text = "m: millones", side = 1, outer = TRUE, line = 3, adj = .99, cex = 0.8)
mtext(text = "Abundancia (millones ind)", side = 2, outer = TRUE, line = 5)
mtext(text = "Biomasa (t)", side = 4, outer = TRUE, line = 4)