Commit 36b9fcd5 authored by aho's avatar aho
Browse files

Merge branch 'develop-remain_issues' into 'master'

Add one check to ensure the grid point has all NAs or no NA along time dim....

See merge request !70
parents d2cb9328 abeb0962
Pipeline #5792 passed with stage
in 4 minutes and 21 seconds
......@@ -5,7 +5,9 @@
#'set to TRUE.
#'
#'@param ano A numerical array of anomalies with named dimensions to calculate
#' EOF. The dimensions must have at least 'time_dim' and 'space_dim'.
#' EOF. The dimensions must have at least 'time_dim' and 'space_dim'. NAs
#' could exist but it should be consistent along time_dim. That is, if one grid
#' point has NAs, all the time steps at this point should be NAs.
#'@param lat A vector of the latitudes of 'ano'.
#'@param lon A vector of the longitudes of 'ano'.
#'@param time_dim A character string indicating the name of the time dimension
......@@ -196,11 +198,22 @@ EOF <- function(ano, lat, lon, time_dim = 'sdate', space_dim = c('lat', 'lon'),
ny <- dim(ano)[2]
nx <- dim(ano)[3]
# Check if all the time steps at one grid point are NA-consistent.
# The grid point should have all NAs or no NA along time dim.
if (any(is.na(ano))) {
ano_latlon <- array(ano, dim = c(nt, ny * nx)) # [time, lat*lon]
na_ind <- which(is.na(ano_latlon), arr.ind = T)
if (dim(na_ind)[1] != nt * length(unique(na_ind[, 2]))) {
stop("Detect certain grid points have NAs but not consistent across time ",
"dimension. If the grid point is NA, it should have NA at all time step.")
}
}
# Build the mask
mask <- ano[1, , ]
mask[!is.finite(mask)] <- NA
mask[is.finite(mask)] <- 1
dim(mask) <- dim(ano)[c(2, 3)]
dim(mask) <- c(ny, nx)
# Replace mask of NAs with 0s for EOF analysis.
ano[!is.finite(ano)] <- 0
......@@ -263,23 +276,12 @@ EOF <- function(ano, lat, lon, time_dim = 'sdate', space_dim = c('lat', 'lon'),
var.eof <- 100 * pca$d[1:neofs]^2/tot.var
for (e in 1:neofs) {
# Factor to normalize the EOF.
eof.patt.nn <- EOF[e, , ] * mask
eof.patt.ms <- sum(eof.patt.nn^2, na.rm = TRUE)
# Normalize the EOF
eof.patt <- eof.patt.nn/eof.patt.ms
# PC is multiplied by the normalization factor and the
# weights, then the reconstruction is only EOF * PC (we have
# multiplied ano by weight)
eof.pc <- PC[, e] * eof.patt.ms * W[e]
eof.patt <- eof.patt/wght
EOF[e, , ] <- eof.patt
PC[, e] <- eof.pc
# Set all masked grid points to NA in the EOFs
# Divide patterns by area weights so that EOF * PC gives unweigthed (original) data
EOF[e, , ] <- EOF[e, , ] * mask / wght
# PC is multiplied by the explained variance,
# so that the reconstruction is only EOF * PC
PC[, e] <- PC[, e] * W[e]
}
if (neofs == 1) {
......
......@@ -6,8 +6,9 @@
#'A colour bar (legend) can be plotted and adjusted. It is possible to draw
#'superimposed arrows, dots, symbols, contour lines and boxes. A number of
#'options is provided to adjust the position, size and colour of the
#'components. This plot function is compatible with figure layouts if colour
#'bar is disabled.
#'components. Some parameters are provided to add and adjust the masks that
#'include continents, oceans, and lakes. This plot function is compatible with
#'figure layouts if colour bar is disabled.
#'
#'@param var Array with the values at each cell of a grid on a regular
#' rectangular or gaussian grid. The array is expected to have two
......@@ -60,13 +61,15 @@
#'@param filled.continents Colour to fill in drawn projected continents.
#' Takes the value gray(0.5) by default or, if 'square = FALSE', takes the
#' value FALSE. If set to FALSE, continents are not filled in.
#'@param filled.oceans A logical value or the color name to fill in drawn
#' projected oceans. The default value is FALSE. If it is TRUE, the default
#' colour is "light blue".
#'@param coast_color Colour of the coast line of the drawn projected continents.
#' Takes the value gray(0.5) by default.
#'@param coast_width Line width of the coast line of the drawn projected
#' continents. Takes the value 1 by default.
#'@param lake_color Colour of the lake or other water body inside continents.
#' It is only functional when 'filled.continents = TRUE'. The default value is
#' 'white'.
#' The default value is NULL.
#'@param contours Array of same dimensions as 'var' to be added to the plot
#' and displayed with contours. Parameter 'brks2' is required to define the
#' magnitude breaks for each contour curve. Disregarded if 'square = FALSE'.
......@@ -116,6 +119,12 @@
#' TRUE by default.
#'@param labW Whether to label the longitude axis with a 'W' instead of minus
#' for negative values. Defaults to FALSE.
#'@param lab_dist_x A numeric of the distance of the longitude labels to the
#' box borders. The default value is NULL and is automatically adjusted by
#' the function.
#'@param lab_dist_y A numeric of the distance of the latitude labels to the
#' box borders. The default value is NULL and is automatically adjusted by
#' the function.
#'@param intylat Interval between latitude ticks on y-axis, in degrees.
#' Defaults to 20.
#'@param intxlon Interval between latitude ticks on x-axis, in degrees.
......@@ -218,6 +227,7 @@ PlotEquiMap <- function(var, lon, lat, varu = NULL, varv = NULL,
triangle_ends = NULL, col_inf = NULL, col_sup = NULL,
colNA = NULL, color_fun = clim.palette(),
square = TRUE, filled.continents = NULL,
filled.oceans = FALSE,
coast_color = NULL, coast_width = 1, lake_color = NULL,
contours = NULL, brks2 = NULL, contour_lwd = 0.5,
contour_color = 'black', contour_lty = 1,
......@@ -227,6 +237,7 @@ PlotEquiMap <- function(var, lon, lat, varu = NULL, varv = NULL,
arr_ref_len = 15, arr_units = "m/s",
arr_scale_shaft = 1, arr_scale_shaft_angle = 1,
axelab = TRUE, labW = FALSE,
lab_dist_x = NULL, lab_dist_y = NULL,
intylat = 20, intxlon = 20,
axes_tick_scale = 1, axes_label_scale = 1,
drawleg = TRUE, subsampleg = NULL,
......@@ -408,10 +419,20 @@ PlotEquiMap <- function(var, lon, lat, varu = NULL, varv = NULL,
} else if (!is.logical(filled.continents)) {
continent_color <- filled.continents
filled.continents <- TRUE
} else if (filled.continents) {
} else {
continent_color <- gray(0.5)
}
# Check filled.oceans
if (!.IsColor(filled.oceans) & !is.logical(filled.oceans)) {
stop("Parameter 'filled.oceans' must be logical or a colour identifier.")
} else if (!is.logical(filled.oceans)) {
ocean_color <- filled.oceans
filled.oceans <- TRUE
} else if (filled.oceans) {
ocean_color <- "light blue"
}
# Check coast_color
if (is.null(coast_color)) {
if (filled.continents) {
......@@ -430,13 +451,13 @@ PlotEquiMap <- function(var, lon, lat, varu = NULL, varv = NULL,
}
# Check lake_color
if (is.null(lake_color)) {
if (filled.continents) {
lake_color <- 'white'
}
} else {
if (!.IsColor(coast_color)) {
stop("Parameter 'coast_color' must be a valid colour identifier.")
if (!is.null(lake_color)) {
# if (filled.continents) {
# lake_color <- 'white'
# }
# } else {
if (!.IsColor(lake_color)) {
stop("Parameter 'lake_color' must be a valid colour identifier.")
}
}
......@@ -536,6 +557,16 @@ PlotEquiMap <- function(var, lon, lat, varu = NULL, varv = NULL,
if (!is.logical(labW)) {
stop("Parameter 'labW' must be logical.")
}
if (!is.null(lab_dist_x)) {
if (!is.numeric(lab_dist_x)) {
stop("Parameter 'lab_dist_x' must be numeric.")
}
}
if (!is.null(lab_dist_y)) {
if (!is.numeric(lab_dist_y)) {
stop("Parameter 'lab_dist_y' must be numeric.")
}
}
if (!is.numeric(intylat)) {
stop("Parameter 'intylat' must be numeric.")
} else {
......@@ -734,10 +765,13 @@ PlotEquiMap <- function(var, lon, lat, varu = NULL, varv = NULL,
}
if (axelab) {
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,
mgp = c(0, spaceticklab + 0.2, 0))
mgp = c(0, lab_distance_y, 0))
axis(1, at = xpos, labels = xlabs, cex.axis = cex_axes_labels, tcl = cex_axes_ticks,
mgp = c(0, spaceticklab + cex_axes_labels / 2 - 0.3, 0))
mgp = c(0, lab_distance_x, 0))
}
title(toptitle, cex.main = cex_title)
rect(par("usr")[1], par("usr")[3], par("usr")[2], par("usr")[4], col = colNA)
......@@ -792,30 +826,67 @@ PlotEquiMap <- function(var, lon, lat, varu = NULL, varv = NULL,
} else { # [0, 360]
xlim_conti <- c(0.01, 359.99)
}
coast <- map(continents, interior = FALSE, wrap = TRUE,
xlim = xlim_conti, ylim = c(-89.99, 89.99),
fill = filled.continents, add = TRUE, plot = FALSE)
if (filled.continents) {
old_lwd <- par('lwd')
par(lwd = coast_width)
if (min(lon) >= 0) {
ylat <- latmin:latmax
xlon <- lonmin:lonmax
proj <- GEOmap::setPROJ(1, LON0 = mean(xlon),
LAT0 = mean(ylat), LATS = ylat, LONS = xlon)
lakes <- which(coastmap$STROKES$col == "blue")
old_lwd <- par('lwd')
par(lwd = coast_width)
# If [0, 360], use GEOmap; if [-180, 180], use maps
if (min(lon) >= 0) {
ylat <- latmin:latmax
xlon <- lonmin:lonmax
proj <- GEOmap::setPROJ(1, LON0 = mean(xlon),
LAT0 = mean(ylat), LATS = ylat, LONS = xlon)
lakes <- which(coastmap$STROKES$col == "blue")
par(new = TRUE)
if (filled.continents) {
coastmap$STROKES$col[which(coastmap$STROKES$col != "blue")] <- continent_color
coastmap$STROKES$col[lakes] <- lake_color #"white"
par(new = TRUE)
if (is.null(lake_color)) {
coastmap$STROKES$col[lakes] <- continent_color
} else {
coastmap$STROKES$col[lakes] <- lake_color #"white"
}
GEOmap::plotGEOmap(coastmap, PROJ = proj, border = coast_color,
add = TRUE, lwd = coast_width)
} else {
polygon(coast, col = continent_color, border = coast_color, lwd = coast_width)
coastmap$STROKES$col[which(coastmap$STROKES$col != "blue")] <- coast_color
if (is.null(lake_color)) {
coastmap$STROKES$col[lakes] <- coast_color
} else {
coastmap$STROKES$col[lakes] <- lake_color #"white"
}
GEOmap::plotGEOmap(coastmap, PROJ = proj, #MAPcol = coast_color,
add = TRUE, lwd = coast_width, MAPstyle = 2)
}
par(lwd = old_lwd)
} else {
lines(coast, col = coast_color, lwd = coast_width)
# [-180, 180]
coast <- map(continents, interior = FALSE, wrap = TRUE,
xlim = xlim_conti, ylim = c(-89.99, 89.99),
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)) {
map('lakes', add = TRUE, fill = filled.continents, col = lake_color)
}
}
par(lwd = old_lwd)
# filled.oceans
if (filled.oceans) {
old_lwd <- par('lwd')
par(lwd = coast_width)
outline <- map(continents, fill = T, plot = FALSE) # must be fill = T
xbox <- xlim_conti + 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)
}
box()
# Draw rectangle on the map
if (!is.null(boxlim)) {
......
......@@ -17,7 +17,9 @@ EOF(
}
\arguments{
\item{ano}{A numerical array of anomalies with named dimensions to calculate
EOF. The dimensions must have at least 'time_dim' and 'space_dim'.}
EOF. The dimensions must have at least 'time_dim' and 'space_dim'. NAs
could exist but it should be consistent along time_dim. That is, if one grid
point has NAs, all the time steps at this point should be NAs.}
\item{lat}{A vector of the latitudes of 'ano'.}
......
......@@ -23,6 +23,7 @@ PlotEquiMap(
color_fun = clim.palette(),
square = TRUE,
filled.continents = NULL,
filled.oceans = FALSE,
coast_color = NULL,
coast_width = 1,
lake_color = NULL,
......@@ -44,6 +45,8 @@ PlotEquiMap(
arr_scale_shaft_angle = 1,
axelab = TRUE,
labW = FALSE,
lab_dist_x = NULL,
lab_dist_y = NULL,
intylat = 20,
intxlon = 20,
axes_tick_scale = 1,
......@@ -138,6 +141,10 @@ the spaces in between with colours (FALSE). In the latter case,
Takes the value gray(0.5) by default or, if 'square = FALSE', takes the
value FALSE. If set to FALSE, continents are not filled in.}
\item{filled.oceans}{A logical value or the color name to fill in drawn
projected oceans. The default value is FALSE. If it is TRUE, the default
colour is "light blue".}
\item{coast_color}{Colour of the coast line of the drawn projected continents.
Takes the value gray(0.5) by default.}
......@@ -145,8 +152,7 @@ Takes the value gray(0.5) by default.}
continents. Takes the value 1 by default.}
\item{lake_color}{Colour of the lake or other water body inside continents.
It is only functional when 'filled.continents = TRUE'. The default value is
'white'.}
The default value is NULL.}
\item{contours}{Array of same dimensions as 'var' to be added to the plot
and displayed with contours. Parameter 'brks2' is required to define the
......@@ -215,6 +221,14 @@ TRUE by default.}
\item{labW}{Whether to label the longitude axis with a 'W' instead of minus
for negative values. Defaults to FALSE.}
\item{lab_dist_x}{A numeric of the distance of the longitude labels to the
box borders. The default value is NULL and is automatically adjusted by
the function.}
\item{lab_dist_y}{A numeric of the distance of the latitude labels to the
box borders. The default value is NULL and is automatically adjusted by
the function.}
\item{intylat}{Interval between latitude ticks on y-axis, in degrees.
Defaults to 20.}
......@@ -300,8 +314,9 @@ grid cells. Only the region for which data has been provided is displayed.
A colour bar (legend) can be plotted and adjusted. It is possible to draw
superimposed arrows, dots, symbols, contour lines and boxes. A number of
options is provided to adjust the position, size and colour of the
components. This plot function is compatible with figure layouts if colour
bar is disabled.
components. Some parameters are provided to add and adjust the masks that
include continents, oceans, and lakes. This plot function is compatible with
figure layouts if colour bar is disabled.
}
\examples{
# See examples on Load() to understand the first lines in this example
......
......@@ -96,71 +96,111 @@ test_that("1. Input checks", {
})
##############################################
test_that("2. dat1", {
res1 <- EOF(dat1, lon = lon1, lat = lat1, neofs = 10)
expect_equal(
names(EOF(dat1, lon = lon1, lat = lat1)),
names(res1),
c("EOFs", "PCs", "var", "tot_var", "mask", "wght")
)
expect_equal(
dim(EOF(dat1, lon = lon1, lat = lat1)$EOFs),
dim(res1$EOFs),
c(mode = 10, lat = 6, lon = 2)
)
expect_equal(
dim(EOF(dat1, lon = lon1, lat = lat1)$PCs),
dim(res1$PCs),
c(sdate = 10, mode = 10)
)
expect_equal(
dim(EOF(dat1, lon = lon1, lat = lat1)$var),
dim(res1$var),
c(mode = 10)
)
expect_equal(
dim(EOF(dat1, lon = lon1, lat = lat1)$mask),
dim(res1$mask),
c(lat = 6, lon = 2)
)
expect_equal(
dim(EOF(dat1, lon = lon1, lat = lat1)$wght),
dim(res1$wght),
c(lat = 6, lon = 2)
)
expect_equal(
EOF(dat1, lon = lon1, lat = lat1)$EOFs[1:5],
res1$EOFs[1:5],
c(-0.2888168, 0.2792765, 0.1028387, 0.1883640, -0.2896943),
tolerance = 0.0001
)
expect_equal(
mean(EOF(dat1, lon = lon1, lat = lat1)$EOFs),
mean(res1$EOFs),
0.01792716,
tolerance = 0.0001
)
expect_equal(
EOF(dat1, lon = lon1, lat = lat1)$PCs[1:5],
res1$PCs[1:5],
c(3.2061306, -0.1692669, 0.5420990, -2.5324441, -0.6143680),
tolerance = 0.0001
)
expect_equal(
mean(EOF(dat1, lon = lon1, lat = lat1)$PCs),
mean(res1$PCs),
0.08980279,
tolerance = 0.0001
)
expect_equal(
EOF(dat1, lon = lon1, lat = lat1)$var[1:5],
res1$var[1:5],
array(c(29.247073, 25.364840, 13.247046, 11.121006, 8.662517), dim = c(mode = 5)),
tolerance = 0.0001
)
expect_equal(
sum(EOF(dat1, lon = lon1, lat = lat1)$mask),
sum(res1$mask),
12
)
expect_equal(
EOF(dat1, lon = lon1, lat = lat1)$wght[1:5],
res1$wght[1:5],
c(0.9923748, 0.9850359, 0.9752213, 0.9629039, 0.9480475),
tolerance = 0.0001
)
expect_equal(
EOF(dat1, lon = lon1, lat = lat1)$tot_var,
res1$tot_var,
88.20996,
tolerance = 0.0001
)
# rebuild the field
latlon_eof <- array(res1$EOFs, dim = c(mode = 10, latlon = 12))
field <- res1$PCs %*% latlon_eof
latlon_dat1<- array(dat1, dim = c(sdate = 10, laton = 12))
expect_equal(
as.vector(latlon_dat1),
as.vector(field)
)
dat1_1 <- dat1
dat1_1[, 2, 1] <- NA
res1_1 <- EOF(dat1_1, lon = lon1, lat = lat1, neofs = 10)
expect_equal(
mean(res1_1$EOFs, na.rm = T),
0.02270081,
tolerance = 0.0001
)
expect_equal(
mean(res1_1$PCs, na.rm = T),
0.1092327,
tolerance = 0.0001
)
# rebuild the field
latlon_eof <- array(res1_1$EOFs, dim = c(mode = 10, latlon = 12))
field <- res1_1$PCs %*% latlon_eof
latlon_dat1<- array(dat1_1, dim = c(sdate = 10, laton = 12))
expect_equal(
as.vector(latlon_dat1),
as.vector(field)
)
dat1_2 <- dat1
dat1_2[2:5, 2, 1] <- NA
expect_error(
EOF(dat1_2, lon = lon1, lat = lat1, neofs = 10),
"Detect certain grid points have NAs but not consistent across time dimension. If the grid point is NA, it should have NA at all time step."
)
})
##############################################
test_that("3. dat2", {
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment