Commit f7018404 authored by aho's avatar aho
Browse files

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

Add one check to ensure the grid point has all NAs or no NA along time dim. Refine the documentation.
parent 9d4099a1
......@@ -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) {
......
......@@ -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'.}
......
......@@ -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