From 9bbd1ed347b44d01ca8a1c80fc160df4e6f3a735 Mon Sep 17 00:00:00 2001
From: aho <an.ho@bsc.es>
Date: Tue, 8 Jun 2021 14:17:10 +0200
Subject: [PATCH 1/5] Add arguments lab_dist_x and lab_dist_y

---
 R/PlotEquiMap.R    | 24 ++++++++++++++++++++++--
 man/PlotEquiMap.Rd | 10 ++++++++++
 2 files changed, 32 insertions(+), 2 deletions(-)

diff --git a/R/PlotEquiMap.R b/R/PlotEquiMap.R
index 8892985..8f59ae2 100644
--- a/R/PlotEquiMap.R
+++ b/R/PlotEquiMap.R
@@ -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. 
@@ -227,6 +233,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, 
@@ -536,6 +543,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 +751,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/PlotEquiMap.Rd b/man/PlotEquiMap.Rd
index 2673353..674e770 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,
@@ -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.}
 
-- 
GitLab


From 9d4099a14625332503224c0f6b32c9f961179d2e Mon Sep 17 00:00:00 2001
From: aho <an.ho@bsc.es>
Date: Tue, 8 Jun 2021 17:13:18 +0200
Subject: [PATCH 2/5] Add a sentence saying that 'lake_color' only works if lon
 = [0, 360].

---
 R/PlotEquiMap.R    | 2 +-
 man/PlotEquiMap.Rd | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/R/PlotEquiMap.R b/R/PlotEquiMap.R
index 8f59ae2..9bc36ea 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'.
diff --git a/man/PlotEquiMap.Rd b/man/PlotEquiMap.Rd
index 674e770..b5c3e03 100644
--- a/man/PlotEquiMap.Rd
+++ b/man/PlotEquiMap.Rd
@@ -148,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 
-- 
GitLab


From f7018404c6977637b3d9fc45d1678432b9fe7400 Mon Sep 17 00:00:00 2001
From: aho <an.ho@bsc.es>
Date: Wed, 9 Jun 2021 16:49:28 +0200
Subject: [PATCH 3/5] Add one check to ensure the grid point has all NAs or no
 NA along time dim. Refine the documentation.

---
 R/EOF.R                   | 40 ++++++++++++-----------
 man/EOF.Rd                |  4 ++-
 tests/testthat/test-EOF.R | 68 +++++++++++++++++++++++++++++++--------
 3 files changed, 78 insertions(+), 34 deletions(-)

diff --git a/R/EOF.R b/R/EOF.R
index 4221ccf..8f8d640 100644
--- a/R/EOF.R
+++ b/R/EOF.R
@@ -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) {
diff --git a/man/EOF.Rd b/man/EOF.Rd
index ae84c55..3d67308 100644
--- a/man/EOF.Rd
+++ b/man/EOF.Rd
@@ -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'.}
 
diff --git a/tests/testthat/test-EOF.R b/tests/testthat/test-EOF.R
index c8c6d30..7966518 100644
--- a/tests/testthat/test-EOF.R
+++ b/tests/testthat/test-EOF.R
@@ -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", {
-- 
GitLab


From cc33fb19be7b05a539456b926ea355ddc8aa1c8c Mon Sep 17 00:00:00 2001
From: aho <an.ho@bsc.es>
Date: Thu, 10 Jun 2021 15:47:27 +0200
Subject: [PATCH 4/5] Add filled.oceans feature

---
 R/PlotEquiMap.R    | 34 ++++++++++++++++++++++++++++++++--
 man/PlotEquiMap.Rd |  5 +++++
 2 files changed, 37 insertions(+), 2 deletions(-)

diff --git a/R/PlotEquiMap.R b/R/PlotEquiMap.R
index 9bc36ea..8988a42 100644
--- a/R/PlotEquiMap.R
+++ b/R/PlotEquiMap.R
@@ -60,6 +60,9 @@
 #'@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 
@@ -224,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,
@@ -419,6 +423,16 @@ PlotEquiMap <- function(var, lon, lat, varu = NULL, varv = NULL,
     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) {
@@ -442,8 +456,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.")
     }
   }
 
@@ -836,6 +850,22 @@ PlotEquiMap <- function(var, lon, lat, varu = NULL, varv = NULL,
   } else {
     lines(coast, col = coast_color, lwd = coast_width)
   }
+
+  # 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)) {
diff --git a/man/PlotEquiMap.Rd b/man/PlotEquiMap.Rd
index b5c3e03..159bd57 100644
--- a/man/PlotEquiMap.Rd
+++ b/man/PlotEquiMap.Rd
@@ -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,
@@ -140,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.}
 
-- 
GitLab


From abeb096201df6937b9e88933b7e7ef311a408e33 Mon Sep 17 00:00:00 2001
From: aho <an.ho@bsc.es>
Date: Thu, 10 Jun 2021 20:37:50 +0200
Subject: [PATCH 5/5] Refine the code for masks

---
 R/PlotEquiMap.R    | 75 +++++++++++++++++++++++++++++-----------------
 man/PlotEquiMap.Rd |  8 ++---
 2 files changed, 52 insertions(+), 31 deletions(-)

diff --git a/R/PlotEquiMap.R b/R/PlotEquiMap.R
index 8988a42..9c78873 100644
--- a/R/PlotEquiMap.R
+++ b/R/PlotEquiMap.R
@@ -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 
@@ -68,8 +69,7 @@
 #'@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'. For now, it is only functional if longitude range is [0, 360].
+#'  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'.
@@ -419,7 +419,7 @@ 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)
   }
 
@@ -451,11 +451,11 @@ 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 (!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.")
     }
@@ -826,30 +826,51 @@ 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) {
diff --git a/man/PlotEquiMap.Rd b/man/PlotEquiMap.Rd
index 159bd57..c025afd 100644
--- a/man/PlotEquiMap.Rd
+++ b/man/PlotEquiMap.Rd
@@ -152,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'. For now, it is only functional if longitude range is [0, 360].}
+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 
@@ -315,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
-- 
GitLab