Newer
Older
} else {
letters[ypos < 0] <- paste(intToUtf8(176), 'S')
letters[ypos > 0] <- paste(intToUtf8(176), 'N')
}
ylabs <- paste(as.character(abs(ypos)), letters, sep = '')
# X axis label
if (!is.null(xlabels)) {
xpos <- seq(lonmin, lonmax, intxlon) + xlonshft
if (length(xpos) != length(xlabels)) {
stop(paste0("Parameter 'xlabels' must have the same length as the longitude ",
"vector spaced by 'intxlon' (length = ", length(xpos), ")."))
}
xlabs <- xlabels
xpos <- seq(lonmin, lonmax, intxlon) + xlonshft
letters <- array('', length(xpos))
if (labW) {
xpos2 <- xpos
xpos2[xpos2 > 180] <- 360 - xpos2[xpos2 > 180]
}
letters[xpos < 0] <- 'W'
letters[xpos > 0] <- 'E'
letters[xpos < 0] <- paste(intToUtf8(176), 'W')
letters[xpos > 0] <- paste(intToUtf8(176), 'E')
}
if (labW) {
letters[xpos == 0] <- ' '
letters[xpos == 180] <- ' '
if (degree_sym == FALSE) {
letters[xpos > 180] <- 'W'
} else {
letters[xpos > 180] <- paste(intToUtf8(176), 'W')
}
xlabs <- paste(as.character(abs(xpos2)), letters, sep = '')
} else {
xlabs <- paste(as.character(abs(xpos)), letters, sep = '')
}
1043
1044
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062
1063
spaceticklab <- max(-cex_axes_ticks, 0)
margins[1] <- margins[1] + 1.2 * cex_axes_labels + spaceticklab
margins[2] <- margins[2] + 1.2 * cex_axes_labels + spaceticklab
}
bar_extra_margin[2] <- bar_extra_margin[2] + margins[2]
bar_extra_margin[4] <- bar_extra_margin[4] + margins[4]
if (toptitle != '') {
margins[3] <- margins[3] + cex_title + 1
}
if (!is.null(varu)) {
margins[1] <- margins[1] + 2.2 * units_scale
}
if (drawleg) {
layout(matrix(1:2, ncol = 1, nrow = 2), heights = c(5, 1))
}
plot.new()
# Load the user parameters
par(userArgs)
par(mar = margins, cex.main = cex_title, cex.axis = cex_axes_labels,
mgp = c(0, spaceticklab, 0), las = 0)
#NOTE: Here creates the window for later plot. If 'usr' for par() is not specified,
# use the lat/lon as the borders. If 'usr' is specified, use the assigned values.
if (is.null(userArgs$usr)) {
#NOTE: The grids are assumed to be equally spaced
xlim_cal <- c(lonb$x[1] - (lonb$x[2] - lonb$x[1]) / 2,
lonb$x[length(lonb$x)] + (lonb$x[2] - lonb$x[1]) / 2)
ylim_cal <- c(latb$x[1] - (latb$x[2] - latb$x[1]) / 2,
latb$x[length(latb$x)] + (latb$x[2] - latb$x[1]) / 2)
plot.window(xlim = xlim_cal, ylim = ylim_cal, xaxs = 'i', yaxs = 'i')
# Below is Old code. The border grids are only half plotted.
# plot.window(xlim = range(lonb$x, finite = TRUE), ylim = range(latb$x, finite = TRUE),
# xaxs = 'i', yaxs = 'i')
} else {
plot.window(xlim = par("usr")[1:2], ylim = par("usr")[3:4], xaxs = 'i', yaxs = 'i')
}
lab_distance_y <- ifelse(is.null(lab_dist_y), spaceticklab + 0.2, lab_dist_y)
lab_distance_x <- ifelse(is.null(lab_dist_x), spaceticklab + cex_axes_labels / 2 - 0.3, lab_dist_x)
axis(2, at = ypos, labels = ylabs, cex.axis = cex_axes_labels, tcl = cex_axes_ticks,
axis(1, at = xpos, labels = xlabs, cex.axis = cex_axes_labels, tcl = cex_axes_ticks,
}
title(toptitle, cex.main = cex_title)
rect(par("usr")[1], par("usr")[3], par("usr")[2], par("usr")[4], col = colNA)
col_inf_image <- ifelse(is.null(col_inf), colNA, col_inf)
col_sup_image <- ifelse(is.null(col_sup), colNA, col_sup)
if (square) {
# If lat and lon are both regular-spaced, "useRaster = TRUE" can avoid
# artifact white lines on the figure. If not, useRaster has to be FALSE (default)
tryCatch({
image(lonb$x, latb$x, var[lonb$ix, latb$ix],
col = c(col_inf_image, cols, col_sup_image),
breaks = c(-.Machine$double.xmax, brks, .Machine$double.xmax),
axes = FALSE, xlab = "", ylab = "", add = TRUE, useRaster = TRUE)
}, error = function(x) {
image(lonb$x, latb$x, var[lonb$ix, latb$ix],
col = c(col_inf_image, cols, col_sup_image),
breaks = c(-.Machine$double.xmax, brks, .Machine$double.xmax),
axes = FALSE, xlab = "", ylab = "", add = TRUE)
})
} else {
.filled.contour(lonb$x, latb$x, var[lonb$ix, latb$ix],
levels = c(.Machine$double.xmin, brks, .Machine$double.xmax),
col = c(col_inf_image, cols, col_sup_image))
}
if (!is.null(contours)) {
#NOTE: 'labcex' is the absolute size of contour labels. Parameter 'contour_label_scale'
# is provided in PlotEquiMap() but it was not used. Here, 'cex_axes_labels' was used
# and it was calculated from 'axes_label_scale', the size of lat/lon axis label.
# It is changed to use contour_label_scale*par('cex').
contour(lonb$x, latb$x, contours[lonb$ix, latb$ix], levels = brks2,
method = "edge", add = TRUE,
# labcex = cex_axes_labels,
labcex = contour_label_scale * par("cex"),
lwd = contour_lwd, lty = contour_lty,
aho
committed
col = contour_color, drawlabels = contour_draw_label)
}
#
# Adding black dots or symbols
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#
if (!is.null(dots)) {
data_avail <- !is.na(var)
for (counter in 1:(dim(dots)[1])) {
points <- which(dots[counter, , ] & data_avail, arr.ind = TRUE)
points(lon[points[, 1]], lat[points[, 2]],
pch = dot_symbol[counter],
cex = dot_size[counter] * 3 / sqrt(sqrt(length(var))),
lwd = dot_size[counter] * 3 / sqrt(sqrt(length(var))))
}
}
#
# Plotting continents
# ~~~~~~~~~~~~~~~~~~~~~
wrap_vec <- c(lonb$x[1], lonb$x[1] + 360)
# If [0, 360], use GEOmap; if [-180, 180], use maps::map
# UPDATE: Use maps::map for both cases. The difference between GEOmap and
# maps is trivial. The only thing we can see for now is that
# GEOmap has better lakes.
coast <- maps::map(interior = country.borders, wrap = wrap_vec,
fill = filled.continents, add = TRUE, plot = FALSE)
if (filled.continents) {
polygon(coast, col = continent_color, border = coast_color, lwd = coast_width)
} else {
lines(coast, col = coast_color, lwd = coast_width)
}
if (!is.null(lake_color)) {
maps::map('lakes', add = TRUE, wrap = wrap_vec, fill = filled.continents, col = lake_color)
# filled.oceans
if (filled.oceans) {
old_lwd <- par('lwd')
par(lwd = coast_width)
outline <- maps::map(wrap = wrap_vec, fill = T, plot = FALSE) # must be fill = T
xbox <- wrap_vec + c(-2, 2)
ybox <- c(-92, 92)
outline$x <- c(outline$x, NA, c(xbox, rev(xbox), xbox[1]))
outline$y <- c(outline$y, NA, rep(ybox, each = 2), ybox[1])
polypath(outline, col = ocean_color, rule = 'evenodd', border = NA)
par(lwd = old_lwd)
}
# Plot shapefile
#NOTE: the longitude range cannot cut shapefile range, or not all the shapefile will be plotted.
if (!is.null(shapefile)) {
maps::map(shape, interior = country.borders, #wrap = wrap_vec,
fill = filled.continents, add = TRUE, plot = TRUE,
lwd = shapefile_lwd, col = shapefile_color)
}
1187
1188
1189
1190
1191
1192
1193
1194
1195
1196
1197
1198
1199
1200
1201
1202
1203
1204
1205
1206
1207
1208
1209
1210
1211
1212
1213
1214
1215
1216
1217
1218
1219
1220
1221
1222
1223
1224
1225
box()
# Draw rectangle on the map
if (!is.null(boxlim)) {
counter <- 1
for (box in boxlim) {
if (box[1] > box[3]) {
box[1] <- box[1] - 360
}
if (length(box) != 4) {
stop(paste("The", counter, "st box defined in the parameter 'boxlim' is ill defined."))
} else if (box[2] < latmin || box[4] > latmax ||
box[1] < lonmin || box[3] > lonmax) {
stop(paste("The limits of the", counter, "st box defined in the parameter 'boxlim' are invalid."))
} else if (box[1] < 0 && box[3] > 0) {
#segments south
segments(box[1], box[2], 0, box[2], col = boxcol[counter], lwd = boxlwd[counter])
segments(0, box[2], box[3], box[2], col = boxcol[counter], lwd = boxlwd[counter])
#segments north
segments(box[1], box[4], 0, box[4], col = boxcol[counter], lwd = boxlwd[counter])
segments(0, box[4], box[3], box[4], col = boxcol[counter], lwd = boxlwd[counter])
#segments west
segments(box[1], box[2], box[1], box[4], col = boxcol[counter],
lwd = boxlwd[counter])
#segments est
segments(box[3], box[2], box[3],box[4], col = boxcol[counter],
lwd = boxlwd[counter])
} else {
rect(box[1], box[2], box[3], box[4], border = boxcol[counter], col = NULL,
lwd = boxlwd[counter], lty = 'solid')
}
counter <- counter + 1
}
}
#
# PlotWind
# ~~~~~~~~~~
#
if (!is.null(varu) && !is.null(varv)) {
# Create a two dimention array of longitude and latitude
lontab <- InsertDim(lonb$x, 2, length(latb$x), name = 'lat')
lattab <- InsertDim(latb$x, 1, length(lonb$x), name = 'lon')
1228
1229
1230
1231
1232
1233
1234
1235
1236
1237
1238
1239
1240
1241
1242
1243
1244
1245
1246
1247
1248
1249
1250
1251
1252
1253
1254
1255
1256
1257
1258
1259
1260
1261
1262
1263
1264
1265
1266
1267
1268
1269
1270
1271
1272
1273
1274
1275
1276
1277
1278
1279
1280
1281
1282
1283
1284
1285
1286
varplotu <- varu[lonb$ix, latb$ix]
varplotv <- varv[lonb$ix, latb$ix]
# Select a subsample af the points to an arrow
#for each "subsample" grid point
sublon <- seq(1,length(lon), arr_subsamp)
sublat <- seq(1,length(lat), arr_subsamp)
uaux <- lontab[sublon, sublat] + varplotu[sublon, sublat] * 0.5 * arr_scale
vaux <- lattab[sublon, sublat] + varplotv[sublon, sublat] * 0.5 * arr_scale
lenshaft <- 0.18 * arr_scale * arr_scale_shaft
angleshaft <- 12 * arr_scale_shaft_angle
# Plot Wind
arrows(lontab[sublon, sublat], lattab[sublon, sublat],
uaux, vaux,
angle = angleshaft,
length = lenshaft)
# Plotting an arrow at the bottom of the plot for the legend
posarlon <- lonb$x[1] + (lonmax - lonmin) * 0.1
posarlat <- latmin - ((latmax - latmin) + 1) / par('pin')[2] *
(spaceticklab + 0.2 + cex_axes_labels + 0.6 * units_scale) * par('csi')
arrows(posarlon, posarlat,
posarlon + 0.5 * arr_scale * arr_ref_len, posarlat,
length = lenshaft, angle = angleshaft,
xpd = TRUE)
#save the parameter value
xpdsave <- par('xpd')
#desactivate xpd to be able to plot in margen
par(xpd = NA)
#plot text
mtext(paste(as.character(arr_ref_len), arr_units, sep = ""),
line = spaceticklab + 0.2 + cex_axes_labels + 1.2 * units_scale, side = 1,
at = posarlon + (0.5 * arr_scale * arr_ref_len) / 2,
cex = units_scale)
#come back to the previous xpd value
par(xpd = xpdsave)
}
#
# Colorbar
# ~~~~~~~~~~
#
if (drawleg) {
ColorBar(brks, cols, FALSE, subsampleg, bar_limits, var_limits,
triangle_ends, col_inf = col_inf, col_sup = col_sup,
extra_labels = bar_extra_labels, draw_ticks = draw_bar_ticks,
draw_separators = draw_separators, title = units,
title_scale = units_scale, triangle_ends_scale = triangle_ends_scale,
label_scale = bar_label_scale, tick_scale = bar_tick_scale,
extra_margin = bar_extra_margin, label_digits = bar_label_digits)
}
# If the graphic was saved to file, close the connection with the device
if (!is.null(fileout)) dev.off()
invisible(list(brks = brks, cols = cols, col_inf = col_inf, col_sup = col_sup))
}