Newer
Older
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 = '')
}
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
1043
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)
}
1167
1168
1169
1170
1171
1172
1173
1174
1175
1176
1177
1178
1179
1180
1181
1182
1183
1184
1185
1186
1187
1188
1189
1190
1191
1192
1193
1194
1195
1196
1197
1198
1199
1200
1201
1202
1203
1204
1205
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')
1208
1209
1210
1211
1212
1213
1214
1215
1216
1217
1218
1219
1220
1221
1222
1223
1224
1225
1226
1227
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
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))
}