#' @title
#' Plot maps with areal results
#'
#' @description
#' \code{PlotMap_Fn} is a hidden function to plot a map and fill in regions with colors to represent intensity in an areal-interpretion of model results
#'
#' @inheritParams plot_maps
#' @param plot_legend_fig Boolean, whether to plot a separate figure for the heatmap legend or not
#' @param land_color color for filling in land (use \code{land_color=rgb(0,0,0,alpha=0)} for transparent land)
#' @param ... arguments passed to \code{par}
#'
#' @details
#' This function was necessary to build because \code{mapproj::mapproject} as used in \code{maps::map} has difficulties with both rotations (for most projections) and
#' truncating the cocuntry boundaries within the plotting region (which \code{mapproj::mapproject} appears to do prior to projection,
#' so that the post-projection is often missing boundaries that are within the plotting rectangle). I use rectangular projections by default, but Lamberts or Albers conformal
#' projections would also be useful for many cases.
#' @return A data frame with location information and predicted values, for the years given, for each grid point in the underlying extrapolation grid
#' @export
PlotMap_Fn_JRW <- function (MappingDetails, Mat, PlotDF, MapSizeRatio = c(`Width(in)` = 4,
`Height(in)` = 4), Xlim, Ylim, FileName = paste0(getwd(),
"/"), Year_Set, Rescale = FALSE, Rotate = 0, Format = "png",
Res = 200, zone = NA, Cex = 0.01, textmargin = "", add = FALSE,
pch = 15, outermargintext = c("Eastings", "Northings"), zlim = NULL,
Col = NULL, Legend = list(use = FALSE, x = c(10, 30), y = c(10,
30)), mfrow = c(1, 1), plot_legend_fig = FALSE, land_color = "grey",
ignore.na = FALSE, ...)
{
require(maps)
require(mapdata)
# on.exit(detach("package:mapdata"))
# on.exit(detach("package:maps"), add = TRUE)
# Mat <<- Mat
# PlotDF <<- PlotDF
Mat = Mat[PlotDF[, "x2i"], , drop = FALSE]
Which = which(PlotDF[, "Include"] > 0)
if (Rescale)
Mat = Mat/outer(rep(Rescale, nrow(Mat)), colMeans(Mat[Which, ]))
f = function(Num, zlim = NULL) {
if (is.null(zlim))
Return = ((Num) - min(Num, na.rm = TRUE))/max(diff(range(Num,
na.rm = TRUE)), 0.01)
if (!is.null(zlim))
Return = ((Num) - zlim[1])/max(diff(zlim), 0.01)
return(Return)
}
if (is.null(Col))
Col = colorRampPalette(colors = c("darkblue", "blue",
"lightblue", "lightgreen", "yellow", "orange", "red"))
Par = list(mfrow = mfrow, ...)
if (!add) {
if (Format == "png") {
png(file = paste0(FileName, ".png"), width = Par$mfrow[2] *
MapSizeRatio["Width(in)"], height = Par$mfrow[1] *
MapSizeRatio["Height(in)"], res = Res, units = "in")
}
if (Format == "jpg") {
jpeg(file = paste0(FileName, ".jpg"), width = Par$mfrow[2] *
MapSizeRatio["Width(in)"], height = Par$mfrow[1] *
MapSizeRatio["Height(in)"], res = Res, units = "in")
}
if (Format %in% c("tif", "tiff")) {
tiff(file = paste0(FileName, ".tif"), width = Par$mfrow[2] *
MapSizeRatio["Width(in)"], height = Par$mfrow[1] *
MapSizeRatio["Height(in)"], res = Res, units = "in")
}
par(Par)
}
for (tI in 1:length(Year_Set)) {
if (is.null(MappingDetails)) {
plot(1, type = "n", ylim = Ylim, xlim = Xlim, main = "",
xlab = "", ylab = "")
points(x = PlotDF[Which, "Lon"], y = PlotDF[Which,
"Lat"], col = Col(n = 50)[ceiling(f(Mat[Which,
], zlim = zlim)[, t] * 49) + 1], cex = 0.01)
}
else {
boundary_around_limits = 3
Map = maps::map(MappingDetails[[1]], MappingDetails[[2]],
plot = FALSE, ylim = mean(Ylim) + boundary_around_limits *
c(-0.5, 0.5) * diff(Ylim), xlim = mean(Xlim) +
boundary_around_limits * c(-0.5, 0.5) * diff(Xlim),
fill = TRUE)
Tmp1 <- na.omit(cbind(PID = cumsum(is.na(Map$x)),
POS = 1:length(Map$x), X = Map$x, Y = Map$y,
matrix(0, ncol = length(Year_Set), nrow = length(Map$x),
dimnames = list(NULL, Year_Set))))
TmpLL <- rbind(Tmp1, cbind(PID = max(Tmp1[, 1]) +
1, POS = 1:length(Which) + max(Tmp1[, 2]), X = PlotDF[Which,
"Lon"], Y = PlotDF[Which, "Lat"], Mat[Which,
]))
tmpUTM = TmpLL
tmpUTM[, c("X", "Y")] = as.matrix(FishStatsUtils:::Convert_LL_to_UTM_Fn(Lon = TmpLL[,
"X"], Lat = TmpLL[, "Y"], zone = zone, flip_around_dateline = ifelse(MappingDetails[[1]] %in%
c("world2", "world2Hires"), TRUE, FALSE))[, c("X",
"Y")])
tmpUTM <- data.frame(tmpUTM)
sp::coordinates(tmpUTM) = c("X", "Y")
tmpUTM_rotated <- maptools::elide(tmpUTM, rotate = Rotate)
plot(1, type = "n", xlim = range(tmpUTM_rotated@coords[-c(1:nrow(Tmp1)),
"x"]), ylim = range(tmpUTM_rotated@coords[-c(1:nrow(Tmp1)),
"y"]), xaxt = "n", yaxt = "n")
Col_Bin = ceiling(f(tmpUTM_rotated@data[-c(1:nrow(Tmp1)),
-c(1:2), drop = FALSE], zlim = zlim)[, tI] *
49) + 1
if (ignore.na == FALSE && any(Col_Bin < 1 | Col_Bin >
50))
stop("zlim doesn't span the range of the variable")
points(x = tmpUTM_rotated@coords[-c(1:nrow(Tmp1)),
"x"], y = tmpUTM_rotated@coords[-c(1:nrow(Tmp1)),
"y"], col = Col(n = 50)[Col_Bin], cex = Cex,
pch = pch)
lev = levels(as.factor(tmpUTM_rotated@data$PID))
for (levI in 1:(length(lev) - 1)) {
indx = which(tmpUTM$PID == lev[levI])
if (var(sign(TmpLL[indx, "Y"])) == 0) {
polygon(x = tmpUTM_rotated@coords[indx, "x"],
y = tmpUTM_rotated@coords[indx, "y"], col = land_color)
}
else {
warning("Skipping map polygons that straddle equation, because PBSmapping::convUL doesn't work for these cases")
}
}
}
title(Year_Set[tI], line = 0.1, cex.main = ifelse(is.null(Par$cex.main),
1.8, Par$cex.main), cex = ifelse(is.null(Par$cex.main),
1.8, Par$cex.main))
box()
}
if (Legend$use) {
FishStatsUtils:::smallPlot(FishStatsUtils:::Heatmap_Legend(colvec = Col(50),
heatrange = list(range(Mat[Which, ], na.rm = TRUE),
zlim)[[ifelse(is.null(zlim), 1, 2)]], dopar = FALSE),
x = Legend$x, y = Legend$y, mar = c(0, 0, 0, 0),
mgp = c(2, 0.5, 0), tck = -0.2, font = 2)
}
if (!add)
mtext(side = 1, outer = TRUE, outermargintext[1], cex = 1.75,
line = par()$oma[1]/2)
if (!add)
mtext(side = 2, outer = TRUE, outermargintext[2], cex = 1.75,
line = par()$oma[2]/2)
if (Format %in% c("png", "jpg", "tif", "tiff"))
# dev.off()
if (plot_legend_fig) {
if (Format == "png") {
png(file = paste0(FileName, "_Legend.png", sep = ""),
width = 1, height = 2 * MapSizeRatio["Height(in)"],
res = Res, units = "in")
}
if (Format == "jpg") {
jpeg(file = paste0(FileName, "_Legend.jpg", sep = ""),
width = 1, height = 2 * MapSizeRatio["Height(in)"],
res = Res, units = "in")
}
if (Format %in% c("tif", "tiff")) {
tiff(file = paste0(FileName, "_Legend.tif", sep = ""),
width = 1, height = 2 * MapSizeRatio["Height(in)"],
res = Res, units = "in")
}
if (Format %in% c("png", "jpg", "tif", "tiff")) {
FishStatsUtils:::Heatmap_Legend(colvec = Col(n = 50),
heatrange = list(range(Mat, na.rm = TRUE), zlim)[[ifelse(is.null(zlim),
1, 2)]], textmargin = textmargin)
# dev.off()
}
}
invisible(data.frame(tmpUTM_rotated@coords[-(1:nrow(Tmp1)), ], TmpLL[-(1:nrow(Tmp1)), -(1:2)]))
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.