diff --git a/R/EOF.R b/R/EOF.R
index 6ed8c6f751bfba8c52ddfdfd4fbcda87d64bbdee..a1c0c7e3f548fd7f6df218b05c57420fa66a0da9 100644
--- a/R/EOF.R
+++ b/R/EOF.R
@@ -5,7 +5,9 @@
 #'\code{TRUE}.
 #'
 #'@param ano Array of anomalies with dimensions (number of timesteps, 
-#'  number of latitudes, number of longitudes).
+#'  number of latitudes, number of longitudes). 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 lon Vector of longitudes of \code{ano}.
 #'@param lat Vector of latitudes of \code{ano}.
 #'@param neofs Number of modes to be kept. Default = 15.
@@ -152,6 +154,17 @@ EOF <- function(ano, lon, lat, neofs = 15, corr = FALSE) {
     stop("Inconsistent number of longitudes and input field dimensions.")
   }
   
+  # 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.")
+    }
+  }
+
   # Buildup of the mask
   mask <- ano[1, , ]
   mask[!is.finite(mask)] <- NA
diff --git a/R/PlotEquiMap.R b/R/PlotEquiMap.R
index e0e02fbcf06644933f72bad95febd959de3c202a..3376a1436dfa010cfe05c81e8a002da04c94af71 100644
--- a/R/PlotEquiMap.R
+++ b/R/PlotEquiMap.R
@@ -66,7 +66,7 @@
 #'  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'.
+#'  'white'. For now, it is only functional if longitude range is [0, 360].
 #'@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 +116,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. 
@@ -235,6 +241,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, 
@@ -443,8 +450,8 @@ PlotEquiMap <- function(var, lon, lat, varu = NULL, varv = NULL,
       lake_color <- 'white'
     }
   } else {
-    if (!.IsColor(coast_color)) {
-      stop("Parameter 'coast_color' must be a valid colour identifier.")
+    if (!.IsColor(lake_color)) {
+      stop("Parameter 'lake_color' must be a valid colour identifier.")
     }
   }
 
@@ -544,6 +551,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 {
@@ -742,10 +759,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)
diff --git a/man/EOF.Rd b/man/EOF.Rd
index 25934f5237641092912a654c1374c3a35ac0a7bb..4ad3d80d452e7a5b7d4143925cecc181bb116861 100644
--- a/man/EOF.Rd
+++ b/man/EOF.Rd
@@ -8,7 +8,9 @@ EOF(ano, lon, lat, neofs = 15, corr = FALSE)
 }
 \arguments{
 \item{ano}{Array of anomalies with dimensions (number of timesteps, 
-number of latitudes, number of longitudes).}
+number of latitudes, number of longitudes). 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{lon}{Vector of longitudes of \code{ano}.}
 
diff --git a/man/PlotEquiMap.Rd b/man/PlotEquiMap.Rd
index 18456bed5819c5e53bac90042f070f119e43d8bc..ee2e56b672a753f0025fa6451b2c144ee2c76811 100644
--- a/man/PlotEquiMap.Rd
+++ b/man/PlotEquiMap.Rd
@@ -44,6 +44,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,
@@ -146,7 +148,7 @@ 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'.}
+'white'. For now, it is only functional if longitude range is [0, 360].}
 
 \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 +217,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.}