diff --git a/NAMESPACE b/NAMESPACE
index d74cc97a892ba0fa115f4135007ec2295501a75f..bfb1808d9adbdbaafa50e1d7cbb079f7b34af96a 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -52,6 +52,7 @@ export(PlotClim)
 export(PlotEquiMap)
 export(PlotLayout)
 export(PlotMatrix)
+export(PlotRobinson)
 export(PlotSection)
 export(PlotStereoMap)
 export(PlotVsLTime)
@@ -86,6 +87,8 @@ export(clim.palette)
 import(NbClust)
 import(SpecsVerification)
 import(bigmemory)
+import(cowplot)
+import(ggplot2)
 import(graphics)
 import(mapproj)
 import(maps)
@@ -94,6 +97,8 @@ import(multiApply)
 import(ncdf4)
 import(parallel)
 import(plyr)
+import(rnaturalearth)
+import(sf)
 import(stats)
 importFrom(ClimProjDiags,CombineIndices)
 importFrom(ClimProjDiags,Subset)
diff --git a/R/PlotRobinson.R b/R/PlotRobinson.R
new file mode 100644
index 0000000000000000000000000000000000000000..d8210258275d481c9e2724e0e860281fc2e6e048
--- /dev/null
+++ b/R/PlotRobinson.R
@@ -0,0 +1,548 @@
+#'Plot map in Robinson or other projections
+#'
+#'Transform a regular grid longitude-latitude data to a different projection and 
+#'plot the map. The target projection must be a valid CRS string, preferrably be 
+#'EPSG or ESRI code; check \link[sf]{st_crs} for more explanation. This function
+#'is mainly tested for Robinson projection (ESRI:54030), but it can work with 
+#'other projection types in theory.\n
+#'The map can be plotted by points or polygon. A legend can be plotted as either
+#'a color bar or a discrete ggplot legend. Dots can be drawn on top of the data,
+#'which can be used for significance test. A mask can be added to not plot the 
+#'data specified. A number of options is provided to adjust aesthetics, like 
+#'position, size, colors, etc. 
+#'
+#'@param data A numeric array with longitude and latitude dimensions. The grid
+#'  should be regular grid. It can contain NA values. 
+#'@param lon A numeric vector of longitude locations of the cell centers of the 
+#'  grid of 'data'. Expected to be regularly spaced, within the range of either
+#'  [-180, 180] or [0, 360]. 
+#'@param lat A numeric vector of latitude locations of the cell centers of the 
+#'  grid of 'data'. Expected to be regularly spaced, within the range [-90, 90]
+#'  of ascending or descending order.
+#'@param lon_dim A character string indicating the longitude dimension name in
+#'  'data'. If it is NULL, the function tries to find the name in 
+#'  \code{s2dv:::.KnownLonNames}. The default value is NULL.
+#'@param lat_dim A character string indicating the latitude dimension name in
+#'  'data'. If it is NULL, the function tries to find the name in 
+#'  \code{s2dv:::.KnownLatNames}. The default value is NULL.
+#'@param target_proj A character string indicating the target projection. It 
+#'  should be a valid crs string. The default projection is Robinson 
+#'  (ESRI:54030). Note that the character string may work differently depending
+#'  on PROJ and GDAL module version.
+#'@param legend A character string indicating the legend style. It can be 's2dv'
+#'  (color bar by \code{ColorBar()}), 'ggplot2' (discrete legend by ggplot2), or
+#'  NULL (no legend),
+#'@param style A character string indicating the plotting style. It can be 
+#'  'point' or 'polygon'. The default value is 'point'. Note that 'polygon' may
+#'  be time- and memory-consuming for global or high-resolution data.
+#'@param dots An array with the same dimensions as 'data' of [0, 1] or logical
+#'  indicating the grids to plot dots. The value 0 or FALSE is the point to be
+#'  dotted.
+#'@param mask An array with the same dimensions as 'data' of [0, 1] or logical
+#'  indicating the grids to not plot data. The value 0 or FALSE is the point not
+#'  to be plotted.
+#'@param brks,cols,bar_limits,triangle_ends Usually only providing 'brks' is 
+#'  enough to generate the desired color bar. These parameters allow to 
+#'  define n breaks that define n - 1 intervals to classify each of the values 
+#'  in 'data'. The corresponding grid cell of a given value in 'data' will be 
+#'  colored in function of the interval it belongs to. These parameters are 
+#'  sent to \code{ColorBar()} to generate the breaks and colours. Additional 
+#'  colours for values beyond the limits of the colour bar are also generated 
+#'  and applied to the plot if 'bar_limits' or 'brks' and 'triangle_ends' are 
+#'  properly provided to do so. See ?ColorBar for a full explanation.
+#'@param col_inf,col_sup,colNA Colour identifiers to color the values that 
+#'  excess the extremes of the color bar and to color NAs, respectively. 'colNA'
+#'  takes attr(cols, 'na_color') if available by default, where cols is the 
+#'  parameter 'cols' if provided or the vector of colors returned by 
+#'  'color_fun'. 'col_inf' and 'col_sup' will take the value of 'colNA' if not 
+#'  specified. See ?ColorBar for a full explanation.
+#'@param color_fun,bar_extra_margin Set of 
+#'  parameters to control the visual aspect of the drawn colour bar 
+#'  (1/3). See ?ColorBar for a full explanation.
+#'@param vertical A logical value indicating the direction of colorbar if 
+#' parameter 'legend' is 's2dv'. The default value is TRUE.
+#'@param toptitle A character string of the top title of the figure, scalable 
+#'  with parameter 'title_size'.
+#'@param caption A character string of the caption located at left-bottom of the
+#'  plot.
+#'@param units A character string of the data units, which is the title of the 
+#'  legend.
+#'@param crop_coastlines A named numeric vector [lonmin, lonmax, latmin, latmax] 
+#'  indicating the region to plot coastlines. Note that the longitude range 
+#'  cannot exceed 180 degrees.
+#'@param point_size A number of the size of the data points if "style = 'point'".
+#'  The default is 'auto' and the function tries to find the appropriate size.
+#'@param title_size A number of the size of the top title. The default is 16.
+#'@param dots_size A number of the size of the dots. The default is 0.5.
+#'@param dots_shape A number indicating the dot shape recognized by parameter 
+#'  'shape' in \code{geom_point()}.
+#'@param coastlines_width A number indicating the width of the coastlines.
+#'@param fileout A character string of the path to save the plot. If not 
+#'  specified (default), a graphic device will pop up. The extension should be  
+#'  accepted by \code{ggsave()}.
+#'@param width A number of the plot width, in the units specified in parameter
+#'  'size_units'. The default is 8.
+#'@param height A number of the plot height, in the units specified in parameter 
+#'  'size_units'. The default is 4.
+#'@param size_units A character string of the units of the size of the device
+#'  (file or window) to plot in. The default is 'in' (inches). See ?ggsave and
+#'  ?Devices for details of the corresponding device.
+#'@param res Resolution of the device (file or window) to plot in. The default 
+#'  value is 300. See ?ggsave 'dpi' and ?Devices for details of the 
+#'  corresponding device.
+#'
+#'@return A map plot with speficied projection, either in pop-up window or a 
+#'  saved file.
+#'
+#'@examples 
+#'data <- array(rep(seq(-10, 10, length.out = 181), 360) + rnorm(360), 
+#'              dim = c(lat = 181, lon = 360))
+#'dots <- data
+#'dots[which(dots < 4 & dots > -4)] <- 0
+#'dots[which(dots != 0)] <- 1
+#'PlotRobinson(data, lon = 0:359, lat = -90:90, dots = dots, 
+#'             brks = seq(-10, 10, length.out = 11), 
+#'             toptitle = 'synthetic example', vertical = F,
+#'             caption = 'Robinson Global\ns2dv::PlotRobinson example',
+#'             bar_extra_margin = c(0, 1, 0, 1), width = 8, height = 6)
+#'PlotRobinson(data, lon = 0:359, lat = -90:90, mask = dots, legend = 'ggplot2',
+#'             target_proj = "+proj=moll", brks = seq(-10, 10, length.out = 11),
+#'             color_fun = clim.palette("purpleorange"), colNA = 'green',
+#'             toptitle = 'synthetic example', 
+#'             caption = 'Mollweide Global\ns2dv::PlotRobinson example',
+#'             width = 8, height = 6)
+#'
+#'@import sf ggplot2 rnaturalearth cowplot
+#'@export
+PlotRobinson <- function(data, lon, lat, lon_dim = NULL, lat_dim = NULL, 
+                         target_proj = 54030, legend = 's2dv', style = 'point', 
+                         dots = NULL, mask = NULL, brks = NULL, cols = NULL, bar_limits = NULL,
+                         triangle_ends = NULL, col_inf = NULL, col_sup = NULL, colNA = NULL, 
+                         color_fun = clim.palette(), bar_extra_margin = rep(0, 4), vertical = TRUE,
+                         toptitle = NULL, caption = NULL, units = NULL, crop_coastlines = NULL,
+                         point_size = "auto", title_size = 16, dots_size = 0.5,
+                         dots_shape = 47, coastlines_width = 0.3,
+                         fileout = NULL, width = 8, height = 4, size_units = "in", 
+                         res = 300) {
+
+  # Sanity check
+  # data
+  data <- drop(data)
+  if (length(dim(data)) != 2) {
+    stop("Parameter 'data' must have two dimensions.")
+  }
+  dims <- dim(data)
+  # lon, lon_dim
+  if (is.null(lon_dim)) {
+    lon_dim <- names(dims)[names(dims) %in% .KnownLonNames()]
+    if (identical(lon_dim, character(0))) {
+      stop("Cannot find known longitude name in data dimension. Please define parameter 'lon_dim'.")
+    }
+  }
+  if (is.unsorted(lon)) {
+    .warning("Parameter 'lon' should be sorted to guarantee the correct result.")
+  }
+  # lat, lat_dim
+  if (is.null(lat_dim)) {
+    lat_dim <- names(dims)[names(dims) %in% .KnownLatNames()]
+    if (identical(lat_dim, character(0))) {
+      stop("Cannot find known latitude name in data dimension. Please define parameter 'lat_dim'.")
+    }
+  }
+  if (!all(names(dims) %in% c(lat_dim, lon_dim))) {
+    stop("Dimensions names in paramter 'data' should match 'lat_dim' and 'lon_dim.")
+  }
+  if (length(lon) != dims[lon_dim]) {
+    stop("Length of parameter 'lon' should match longitude dimension in 'data'.")
+  }
+  if (length(lat) != dims[lat_dim]) {
+    stop("Length of parameter 'lat' should match latitude dimension in 'data'.")
+  }
+  data <- s2dv::Reorder(data, c(lon_dim, lat_dim))
+  # Make lat always from 90 to -90
+  sort_lat <- FALSE
+  if (!is.unsorted(lat)) {
+    lat <- rev(lat)
+    data <- ClimProjDiags::Subset(data, along = lat_dim, indices = seq(length(lat), 1, -1))
+    sort_lat <- TRUE
+  }
+
+  # original_proj: it can only be regular grid now
+  original_proj <- st_crs(4326)
+  # tartget_proj
+  if (is.null(target_proj)) {
+    stop("Parameter 'target_proj' cannot be NULL.")
+  } else {
+    target_proj_tmp <- st_crs(target_proj)
+    if (is.na(target_proj_tmp)) {
+      .warning(paste0("Try ESRI code: ESRI:", target_proj))
+      target_proj <- st_crs(paste0("ESRI:", target_proj))
+    } else {
+      target_proj <- target_proj_tmp
+    }
+  }
+
+  # legend
+  if (!is.null(legend) && (!legend %in% c('s2dv', 'ggplot2'))) {
+    stop("Parameter 'legend' must be NULL, ggplot2 or s2dv.")
+  }
+  # style
+  if (!style %in% c('point', 'polygon') || length(style) != 1) {
+    stop("Parameter 'style' must be 'point' or 'polygon'.")
+  }
+  if (style == 'polygon') {
+    # polygon is slow for global map (and may be wrong) Confirm if users want to proceed
+    if ((abs(diff(range(lon))) > 350 & abs(diff(range(lat))) > 175) |
+        (prod(dim(data)) >= (180 * 360))) {
+      if (!isTRUE(utils::askYesNo("The region seems to be global and style 'polygon' is chosen. It may be time- and memory-consuming to plot the map. Are you sure that you want to continue?"))) {
+       return(invisible())
+      }
+    }
+  }
+  # dots
+  if (!is.null(dots)) {
+    dots <- drop(dots)
+    if (any(!names(dim(dots)) %in% c(lon_dim, lat_dim))) {
+      stop("Parameter 'dots' must have two dimensions named as longitude and latitude dimensions in 'data'.")
+    } else {
+      dots <- Reorder(dots, c(lon_dim, lat_dim))
+    }
+    if (!identical(dim(dots), dim(data))) {
+      stop("Parameter 'dots' must have the same dimensions as 'data'.")
+    } else if (is.numeric(dots)) {
+      if (all(dots %in% c(0, 1))) {
+        dots <- array(as.logical(dots), dim = dim(dots))
+      } else {
+        stop("Parameter 'dots' must have only TRUE/FALSE or 0/1.")
+      }
+    } else if (is.logical(dots)) {
+      if (!all(dots %in% c(T, F))) {
+        stop("Parameter 'dots' must have only TRUE/FALSE or 0/1.")
+      }
+    } else {
+      stop("Parameter 'dots' must be a logical or numerical array.")
+    }
+  }
+  # mask
+  if (!is.null(mask)) {
+    mask <- drop(mask)
+    if (any(!names(dim(mask)) %in% c(lon_dim, lat_dim))) {
+      stop("Parameter 'mask' must have two dimensions named as longitude and latitude dimensions in 'data'.")
+    } else {
+      mask <- Reorder(mask, c(lon_dim, lat_dim))
+    }
+    if (!identical(dim(mask), dim(data))) {
+      stop("Parameter 'mask' must have the same dimensions as 'data'.")
+    } else if (is.numeric(mask)) {
+      if (all(mask %in% c(0, 1))) {
+        mask <- array(as.logical(mask), dim = dim(mask))
+      } else {
+        stop("Parameter 'mask' must have only TRUE/FALSE or 0/1.")
+      }
+    } else if (is.logical(mask)) {
+      if (!all(mask %in% c(T, F))) {
+        stop("Parameter 'mask' must have only TRUE/FALSE or 0/1.")
+      }
+    } else {
+      stop("Parameter 'mask' must be a logical or numerical array.")
+    }
+  }
+
+  # Color bar
+  ## Check: brks, cols, bar_limits, color_fun, bar_extra_margin, units 
+  ## Build: brks, cols, bar_limits, col_inf, col_sup
+  var_limits <- c(min(data, na.rm = TRUE), max(data, na.rm = TRUE))
+  colorbar <- ColorBar(brks = brks, cols = cols, vertical = vertical, subsampleg = NULL,
+                       bar_limits = bar_limits, var_limits = var_limits, triangle_ends = triangle_ends, 
+                       col_inf = col_inf, col_sup = col_sup, color_fun = color_fun, 
+                       plot = FALSE, draw_ticks = TRUE, draw_separators = FALSE,
+                       triangle_ends_scale = 1, extra_labels = NULL,
+                       title = units, title_scale = 1, # units_scale
+                       label_scale = 1, tick_scale = 1, #bar_tick_scale
+                       extra_margin = bar_extra_margin, label_digits = 4)
+  brks <- colorbar$brks
+  cols <- colorbar$cols
+  col_inf <- colorbar$col_inf
+  col_sup <- colorbar$col_sup
+  bar_limits <- c(head(brks, 1), tail(brks, 1))
+  # colNA
+  if (is.null(colNA)) {
+    if ('na_color' %in% names(attributes(cols))) {
+      colNA <- attr(cols, 'na_color')
+      if (!.IsColor(colNA)) {
+        stop("The 'na_color' provided as attribute of the colour vector must be a valid colour identifier.")
+      }
+    } else {
+      colNA <- 'pink'
+    }
+  } else if (!.IsColor(colNA)) {
+    stop("Parameter 'colNA' must be a valid colour identifier.")
+  }
+  # toptitle
+  if (!is.null(toptitle) && !is.character(toptitle)) {
+    stop("Parameter 'toptitle' must be a character string.")
+  }
+  # caption
+  if (!is.null(caption) && !is.character(caption)) {
+    stop("Parameter 'caption' must be a character string.")
+  }
+  # crop_coastlines
+  if (!is.null(crop_coastlines)) {
+    # if crop_coastlines doesn't have name, [lonmin, lonmax, latmin, latmax]
+    if (is.null(names(crop_coastlines))) {
+      names(crop_coastlines) <- c("lonmin", "lonmax", "latmin", "latmax")
+    } else if (!identical(sort(names(crop_coastlines)), sort(c("latmax", "latmin", "lonmax", "lonmin")))) {
+      stop("Parameter 'crop_coastlines' needs to have names 'latmax', 'latmin', 'lonmax', 'lonmin'.")
+    }
+  }
+
+  # point_size
+  if (point_size == 'auto') {
+    # 360x181 with default setting, 0.05
+    point_size <- round(0.05 * (360 * 181) / (length(lon) * length(lat)), 2)
+  } else if (!is.numeric(point_size)) {
+    stop("Parameter 'point_size' must be a number.")
+  }
+  #
+
+#=================================================================
+
+  # Adapt s2dv ColorBar parameters to ggplot plot
+  # If legend is NULL, still tune with s2dv legend way
+  if (is.null(legend) || legend == 's2dv') {
+    # the colorbar triangle color. If it is NULL (no triangle plotted), use colNA
+    col_inf_image <- ifelse(is.null(col_inf), colNA, col_inf)
+    col_sup_image <- ifelse(is.null(col_sup), colNA, col_sup)
+    cols_ggplot <- c(col_inf_image, cols, col_sup_image)
+
+    # Add triangles to brks
+    brks_ggplot <- brks
+    if (max(data, na.rm = T) > tail(brks, 1)) {
+      brks_ggplot <- c(brks_ggplot, max(data, na.rm = T))
+    } else {
+      brks_ggplot <- c(brks_ggplot, tail(brks, 1) + diff(tail(brks, 2)))
+    }
+    if (min(data, na.rm = T) < brks[1]) {
+      brks_ggplot <- c(min(data, na.rm = T), brks_ggplot)
+    } else {
+      brks_ggplot <- c(brks[1] - diff(brks[1:2]), brks_ggplot)
+    }
+
+  } else {  # ggplot2 legend
+    brks_ggplot <- brks
+    cols_ggplot <- cols
+  }
+
+  # Build data dataframe
+  lonlat_df <- data.frame(lon = rep(as.vector(lon), length(lat)),
+                          lat = sort(rep(as.vector(lat), length(lon)), decreasing = TRUE))
+  data_df <- lonlat_df %>%
+    dplyr::mutate(dat = as.vector(data))
+
+  lonlat_df_ori <- NULL
+  # Remove the points where mask = FALSE
+  if (!is.null(mask)) {
+    # Save original lonlat_df to plot with expected region
+    lonlat_df_ori <- st_as_sf(lonlat_df, coords = c("lon", "lat"), crs = original_proj)
+    lonlat_df_ori <- st_transform(lonlat_df_ori, crs = target_proj)
+    lonlat_df_ori <- as.data.frame(st_coordinates(lonlat_df_ori))
+    names(lonlat_df_ori) <- c('long', 'lat')
+
+    if (sort_lat) {
+      mask <- ClimProjDiags::Subset(mask, along = lat_dim, indices = seq(length(lat), 1, -1))
+    }
+    mask_df <- data.frame(lon = rep(as.vector(lon), length(lat)),
+                          lat = sort(rep(as.vector(lat), length(lon)), decreasing = TRUE),
+                          mask = as.vector(mask))
+    data_df <- data_df[mask_df$mask == TRUE, ]
+    lonlat_df <- data_df[, 1:2]
+  }
+
+  #NOTE: if target_proj = "ESRI:54030", Nord3v2 has different behavior from hub and ws!!
+  data_df <- st_as_sf(data_df, coords = c("lon", "lat"), crs = original_proj)
+  data_df <- st_transform(data_df, crs = target_proj)
+  data_df <- data_df %>%
+    dplyr::mutate(long = st_coordinates(data_df)[, 1],
+                  lat = st_coordinates(data_df)[, 2])
+
+  # Re-project dots
+  if (!is.null(dots)) {
+    if (sort_lat) {
+      dots <- ClimProjDiags::Subset(dots, along = lat_dim, indices = seq(length(lat), 1, -1))
+    }
+    dots_df <- data.frame(lon = rep(as.vector(lon), length(lat)),
+                          lat = sort(rep(as.vector(lat), length(lon)), decreasing = TRUE),
+                          dot = as.vector(dots))
+
+    dots_df <- st_as_sf(dots_df, coords = c("lon", "lat"), crs = original_proj)
+    dots_df <- st_transform(dots_df, crs = target_proj)
+    dots_df <- dots_df %>%
+      dplyr::mutate(long = st_coordinates(dots_df)[, 1],
+                    lat = st_coordinates(dots_df)[, 2])
+    dots_df <- subset(dots_df, dot == FALSE)
+  }
+
+  # coastlines
+  coastlines <- rnaturalearth::ne_coastline(scale = "medium", returnclass = "sf")
+  ## crop the coastlines to the desired range
+  if (!is.null(crop_coastlines)) {
+    suppressWarnings({
+    coastlines <- st_crop(coastlines, 
+                          xmin = as.numeric(crop_coastlines['lonmin']), 
+                          xmax = as.numeric(crop_coastlines['lonmax']), 
+                          ymin = as.numeric(crop_coastlines['latmin']), 
+                          ymax = as.numeric(crop_coastlines['latmax']))
+    })
+  }
+  coastlines <- st_transform(coastlines, crs = target_proj) 
+
+  if (style == 'polygon') {
+    # Calculate polygon points from regular lat/lon
+    #NOTE: The original grid must be regular grid with same space
+    d_lon <- abs(lon[2] - lon[1]) / 2
+    d_lat <- abs(lat[2] - lat[1]) / 2
+    lon_poly <- lat_poly <- rep(NA, 4 * dim(lonlat_df)[1])
+    for (ii in 1:dim(lonlat_df)[1]) {
+      lon_poly[(ii*4-3):(ii*4)] <- c(lonlat_df$lon[ii] - d_lon, lonlat_df$lon[ii] + d_lon,
+                                     lonlat_df$lon[ii] + d_lon, lonlat_df$lon[ii] - d_lon)
+      lat_poly[(ii*4-3):(ii*4)] <- c(lonlat_df$lat[ii] - d_lat, lonlat_df$lat[ii] - d_lat,
+                                     lonlat_df$lat[ii] + d_lat, lonlat_df$lat[ii] + d_lat)
+    }
+    #  # To prevent out-of-global lon
+    #  lon_poly[which(lon_poly >= 180)] <- 179.9
+    #  lon_poly[which(lon_poly < -180)] <- -180
+    # To prevent out-of-global lat
+    lat_poly[which(lat_poly > 90)] <- 90
+    lat_poly[which(lat_poly < -90)] <- -90
+
+    lonlat_df <- data.frame(lon = lon_poly, lat = lat_poly)
+    # Transfer lon/lat to projection
+    proj_lonlat <- st_as_sf(lonlat_df, coords = c("lon", "lat"), crs = original_proj)
+    #NOTE: if target_proj = "ESRI:54030", on Nord3v2, st_transform has lon and lat swapped!
+    proj_lonlat <- st_transform(proj_lonlat, crs = target_proj)
+    lonlat_df_proj <- st_coordinates(proj_lonlat)
+
+    # Use id to create groups for each polygon
+    ids <- factor(paste0("id_", 1:dim(data_df)[1]))
+    values <- data.frame(id = ids,
+                         value = data_df$dat)
+    positions <- data.frame(id = rep(ids, each = 4),
+                            x = lonlat_df_proj[, 1],
+                            y = lonlat_df_proj[, 2])
+    datapoly <- merge(values, positions, by = "id")
+    datapoly <- st_as_sf(datapoly, coords = c("x", "y"), crs = target_proj)
+    datapoly <- datapoly %>%
+                dplyr::group_by(id) %>%
+                dplyr::summarise() %>%  #NOTE: VERY SLOW if plot global
+                dplyr::mutate(value = values[order(values$id), ]$value) %>%
+                st_cast("POLYGON") %>%
+                st_convex_hull()  # maintain outer polygen (no bowtie shape)
+  }
+
+  # Plots
+  if (style == 'polygon') {
+    res_p <- ggplot(data = data_df) +  #NOTE: must be data_df?
+             geom_sf(data = datapoly,
+                     aes(col = cut(value, breaks = brks_ggplot, include.lowest = T),
+                         fill = cut(value, breaks = brks_ggplot, include.lowest = T)))
+  } else if (style == 'point') {
+    res_p <- ggplot(data = data_df) +
+             geom_point(aes(x = long, y = lat,
+                            col = cut(dat, breaks = brks_ggplot, include.lowest = T)),
+                        #NOTE: These two lines make point size vary with lat
+                        #size = point_size /  (data_df$lat / min(data_df$lat))) +
+                        #size = (sort(rep(as.vector(lat), length(lon))) / max(lat)) * point_size) +
+                        size = point_size)
+  }
+
+  if (is.null(lonlat_df_ori)) {
+    coord_sf_lim <- c(range(data_df$long), range(data_df$lat))
+  } else {
+    coord_sf_lim <- c(range(lonlat_df_ori$long), range(lonlat_df_ori$lat))
+  }
+  res_p <- res_p +
+           geom_sf(data = coastlines, colour ='black', size = coastlines_width) +
+           # Remove background grid and lat/lon label; add white background
+           theme_void() + theme(plot.background = element_rect(fill = 'white', colour = 'white')) +
+           # crop the projection
+           coord_sf(xlim = coord_sf_lim[1:2], ylim = coord_sf_lim[3:4],
+                    expand = TRUE, datum = target_proj)
+
+  if (!is.null(dots)) {
+    res_p <- res_p + geom_point(data = dots_df, aes(x = long, y = lat),
+                                shape = dots_shape, size = dots_size)
+                                #NOTE: This line makes point size vary with lat
+                                #size = dots_size / (dots_df$lat / min(dots_df$lat)))
+  }
+
+  if (identical(legend, 'ggplot2')) {
+    if (style == 'polygon') {
+    res_p <- res_p + scale_colour_manual(values = cols_ggplot,
+                                         aesthetics = c("colour", "fill"),
+                                         drop = FALSE, na.value = colNA) +
+             guides(fill = guide_legend(title = units, override.aes = list(size = 1)),
+                    color = "none")
+   } else if (style == 'point') {
+    res_p <- res_p + scale_colour_manual(values = cols_ggplot,
+                                         drop = FALSE, na.value = colNA) +
+             guides(colour = guide_legend(title = units, override.aes = list(size = 1)))
+   }
+
+  } else {  # s2dv or NULL
+    if (style == 'polygon') {
+      res_p <- res_p + scale_colour_manual(values = cols_ggplot,
+                                           aesthetics = c("colour", "fill"),
+                                           drop = FALSE, na.value = colNA)
+    } else if (style == 'point') {
+      res_p <- res_p + scale_colour_manual(values = cols_ggplot,
+                                           drop = FALSE, na.value = colNA)
+    }
+    # Remove ggplot legend
+    res_p <- res_p + theme(legend.position = "none", plot.margin = margin(0.5, 0, 0, 0, 'cm'))
+  }
+
+  if (!is.null(toptitle)) {
+    res_p <- res_p + ggtitle(toptitle) +
+      theme(plot.title = element_text(size = title_size, hjust = 0.5, vjust = 3))
+  }
+  if (!is.null(caption)) {
+    res_p <- res_p + labs(caption = caption) +
+      theme(plot.caption = element_text(hjust = 0, vjust = 1, margin = margin(0, 0, 0, 0, 'cm')))
+  }
+
+  # s2dv legend fun to put in cowplot::plot_grid
+  if (identical(legend, 's2dv')) {
+    fun_legend <- function() {
+      if (vertical) {
+        par(mar = c(7.1, 2.2, 7.1, 3.1), mgp = c(3, 1, 0))
+      } else {
+        par(mar = c(1.1, 1.2, 0.1, 1.1), mgp = c(3, 1, 0))
+      }
+      ColorBar(brks = brks, cols = cols, vertical = vertical, subsampleg = NULL,
+               bar_limits = bar_limits, var_limits = var_limits, triangle_ends = triangle_ends,
+               col_inf = col_inf, col_sup = col_sup, color_fun = color_fun,
+               plot = TRUE, draw_ticks = TRUE, draw_separators = FALSE,
+               triangle_ends_scale = 1, extra_labels = NULL,
+               title = units, title_scale = 1, # units_scale
+               label_scale = 1, tick_scale = 1, #bar_tick_scale
+               extra_margin = bar_extra_margin, label_digits = 4)
+    }
+    if (vertical) {
+      res_p <- cowplot::plot_grid(res_p, fun_legend, rel_widths = c(6, 1))
+    } else {
+      res_p <- cowplot::plot_grid(res_p, fun_legend, rel_heights = c(5, 1), ncol = 1)
+    }
+    res_p <- res_p +  theme(plot.background = element_rect(fill = "white", colour = "white"))
+  }
+
+  if (!is.null(fileout)) {
+    ext <- regmatches(fileout, regexpr("[a-zA-Z0-9]*$", fileout))
+    ggsave(fileout, res_p, width = width, height = height, dpi = res, units = size_units, 
+           device = ext)
+  } else { # pop-up window
+    dev.new(units = size_units, res = res, width = width, height = height)
+    res_p 
+  } 
+
+}
+
diff --git a/R/PlotStereoMap.R b/R/PlotStereoMap.R
index 4b4fbd22b93110e15e51ffe63071fac24f3b1d41..92831459f16c1998707e4440b59fd15d7800235e 100644
--- a/R/PlotStereoMap.R
+++ b/R/PlotStereoMap.R
@@ -172,7 +172,7 @@
 #'y <- seq(from = -90, to = 90, length.out = 50)
 #'PlotStereoMap(data, x, y, latlims = c(60, 90), brks = 50,
 #'              toptitle = "This is the title")
-#'@import mapproj
+#'@import mapproj maps
 #'@importFrom grDevices dev.cur dev.new dev.off gray
 #'@importFrom stats median
 #'@export
diff --git a/man/PlotRobinson.Rd b/man/PlotRobinson.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..559edb403fdf7f5487521cfb4fb406dd3f9aadf5
--- /dev/null
+++ b/man/PlotRobinson.Rd
@@ -0,0 +1,186 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/PlotRobinson.R
+\name{PlotRobinson}
+\alias{PlotRobinson}
+\title{Plot map in Robinson or other projections}
+\usage{
+PlotRobinson(
+  data,
+  lon,
+  lat,
+  lon_dim = NULL,
+  lat_dim = NULL,
+  target_proj = 54030,
+  legend = "s2dv",
+  style = "point",
+  dots = NULL,
+  mask = NULL,
+  brks = NULL,
+  cols = NULL,
+  bar_limits = NULL,
+  triangle_ends = NULL,
+  col_inf = NULL,
+  col_sup = NULL,
+  colNA = NULL,
+  color_fun = clim.palette(),
+  bar_extra_margin = rep(0, 4),
+  vertical = TRUE,
+  toptitle = NULL,
+  caption = NULL,
+  units = NULL,
+  crop_coastlines = NULL,
+  point_size = "auto",
+  title_size = 16,
+  dots_size = 0.5,
+  dots_shape = 47,
+  coastlines_width = 0.3,
+  fileout = NULL,
+  width = 8,
+  height = 4,
+  size_units = "in",
+  res = 300
+)
+}
+\arguments{
+\item{data}{A numeric array with longitude and latitude dimensions. The grid
+should be regular grid. It can contain NA values.}
+
+\item{lon}{A numeric vector of longitude locations of the cell centers of the 
+grid of 'data'. Expected to be regularly spaced, within the range of either
+[-180, 180] or [0, 360].}
+
+\item{lat}{A numeric vector of latitude locations of the cell centers of the 
+grid of 'data'. Expected to be regularly spaced, within the range [-90, 90]
+of ascending or descending order.}
+
+\item{lon_dim}{A character string indicating the longitude dimension name in
+'data'. If it is NULL, the function tries to find the name in 
+\code{s2dv:::.KnownLonNames}. The default value is NULL.}
+
+\item{lat_dim}{A character string indicating the latitude dimension name in
+'data'. If it is NULL, the function tries to find the name in 
+\code{s2dv:::.KnownLatNames}. The default value is NULL.}
+
+\item{target_proj}{A character string indicating the target projection. It 
+should be a valid crs string. The default projection is Robinson 
+(ESRI:54030). Note that the character string may work differently depending
+on PROJ and GDAL module version.}
+
+\item{legend}{A character string indicating the legend style. It can be 's2dv'
+(color bar by \code{ColorBar()}), 'ggplot2' (discrete legend by ggplot2), or
+NULL (no legend),}
+
+\item{style}{A character string indicating the plotting style. It can be 
+'point' or 'polygon'. The default value is 'point'. Note that 'polygon' may
+be time- and memory-consuming for global or high-resolution data.}
+
+\item{dots}{An array with the same dimensions as 'data' of [0, 1] or logical
+indicating the grids to plot dots. The value 0 or FALSE is the point to be
+dotted.}
+
+\item{mask}{An array with the same dimensions as 'data' of [0, 1] or logical
+indicating the grids to not plot data. The value 0 or FALSE is the point not
+to be plotted.}
+
+\item{brks, cols, bar_limits, triangle_ends}{Usually only providing 'brks' is 
+enough to generate the desired color bar. These parameters allow to 
+define n breaks that define n - 1 intervals to classify each of the values 
+in 'data'. The corresponding grid cell of a given value in 'data' will be 
+colored in function of the interval it belongs to. These parameters are 
+sent to \code{ColorBar()} to generate the breaks and colours. Additional 
+colours for values beyond the limits of the colour bar are also generated 
+and applied to the plot if 'bar_limits' or 'brks' and 'triangle_ends' are 
+properly provided to do so. See ?ColorBar for a full explanation.}
+
+\item{col_inf, col_sup, colNA}{Colour identifiers to color the values that 
+excess the extremes of the color bar and to color NAs, respectively. 'colNA'
+takes attr(cols, 'na_color') if available by default, where cols is the 
+parameter 'cols' if provided or the vector of colors returned by 
+'color_fun'. 'col_inf' and 'col_sup' will take the value of 'colNA' if not 
+specified. See ?ColorBar for a full explanation.}
+
+\item{color_fun, bar_extra_margin}{Set of 
+parameters to control the visual aspect of the drawn colour bar 
+(1/3). See ?ColorBar for a full explanation.}
+
+\item{vertical}{A logical value indicating the direction of colorbar if 
+parameter 'legend' is 's2dv'. The default value is TRUE.}
+
+\item{toptitle}{A character string of the top title of the figure, scalable 
+with parameter 'title_size'.}
+
+\item{caption}{A character string of the caption located at left-bottom of the
+plot.}
+
+\item{units}{A character string of the data units, which is the title of the 
+legend.}
+
+\item{crop_coastlines}{A named numeric vector [lonmin, lonmax, latmin, latmax] 
+indicating the region to plot coastlines. Note that the longitude range 
+cannot exceed 180 degrees.}
+
+\item{point_size}{A number of the size of the data points if "style = 'point'".
+The default is 'auto' and the function tries to find the appropriate size.}
+
+\item{title_size}{A number of the size of the top title. The default is 16.}
+
+\item{dots_size}{A number of the size of the dots. The default is 0.5.}
+
+\item{dots_shape}{A number indicating the dot shape recognized by parameter 
+'shape' in \code{geom_point()}.}
+
+\item{coastlines_width}{A number indicating the width of the coastlines.}
+
+\item{fileout}{A character string of the path to save the plot. If not 
+specified (default), a graphic device will pop up. The extension should be  
+accepted by \code{ggsave()}.}
+
+\item{width}{A number of the plot width, in the units specified in parameter
+'size_units'. The default is 8.}
+
+\item{height}{A number of the plot height, in the units specified in parameter 
+'size_units'. The default is 4.}
+
+\item{size_units}{A character string of the units of the size of the device
+(file or window) to plot in. The default is 'in' (inches). See ?ggsave and
+?Devices for details of the corresponding device.}
+
+\item{res}{Resolution of the device (file or window) to plot in. The default 
+value is 300. See ?ggsave 'dpi' and ?Devices for details of the 
+corresponding device.}
+}
+\value{
+A map plot with speficied projection, either in pop-up window or a 
+ saved file.
+}
+\description{
+Transform a regular grid longitude-latitude data to a different projection and 
+plot the map. The target projection must be a valid CRS string, preferrably be 
+EPSG or ESRI code; check \link[sf]{st_crs} for more explanation. This function
+is mainly tested for Robinson projection (ESRI:54030), but it can work with 
+other projection types in theory.\n
+The map can be plotted by points or polygon. A legend can be plotted as either
+a color bar or a discrete ggplot legend. Dots can be drawn on top of the data,
+which can be used for significance test. A mask can be added to not plot the 
+data specified. A number of options is provided to adjust aesthetics, like 
+position, size, colors, etc.
+}
+\examples{
+data <- array(rep(seq(-10, 10, length.out = 181), 360) + rnorm(360), 
+             dim = c(lat = 181, lon = 360))
+dots <- data
+dots[which(dots < 4 & dots > -4)] <- 0
+dots[which(dots != 0)] <- 1
+PlotRobinson(data, lon = 0:359, lat = -90:90, dots = dots, 
+            brks = seq(-10, 10, length.out = 11), 
+            toptitle = 'synthetic example', vertical = F,
+            caption = 'Robinson Global\ns2dv::PlotRobinson example',
+            bar_extra_margin = c(0, 1, 0, 1), width = 8, height = 6)
+PlotRobinson(data, lon = 0:359, lat = -90:90, mask = dots, legend = 'ggplot2',
+            target_proj = "+proj=moll", brks = seq(-10, 10, length.out = 11),
+            color_fun = clim.palette("purpleorange"), colNA = 'green',
+            toptitle = 'synthetic example', 
+            caption = 'Mollweide Global\ns2dv::PlotRobinson example',
+            width = 8, height = 6)
+
+}
diff --git a/tests/testthat/test-PlotRobinson.R b/tests/testthat/test-PlotRobinson.R
new file mode 100644
index 0000000000000000000000000000000000000000..32c5653ba1a335be20f07aaa6dcd18d5395fb552
--- /dev/null
+++ b/tests/testthat/test-PlotRobinson.R
@@ -0,0 +1,148 @@
+# data1: global
+set.seed(19)
+data1 <- array(rep(seq(-10, 10, length.out = 181), 360) + rnorm(360), 
+               dim = c(sdate = 1, lat = 181, lon = 360))
+mask1 <- drop(data1)
+mask1[which(mask1 > 9 & mask1 < -9)] <- 0
+mask1[which(mask1 != 0)] <- 1
+dots1 <- mask1
+lon1_1 <- -180:179
+lon1_2 <- 0:359
+lat1 <- -90:90
+
+data1_NA <- data1; data1_NA[10000:15000] <- NA
+
+# data2: Europe
+data2 <- data1[,21:61, 161:221]
+lon2 <- lon1_1[161:221]
+lat2 <- 70:30
+dots2 <- array(t(dots1[21:61, 161:221]), dim = c(lon = length(lon2), lat = length(lat2)))
+
+#--------------------------------------------------------------------
+
+test_that("1. Input checks", {
+
+# data
+expect_error(
+  PlotRobinson(array(data1*2, dim = c(sdate = 2, lat = 181, lon = 360))),
+  "Parameter 'data' must have two dimensions."
+)  
+# lon, lon_dim
+tmp <- data1; dim(tmp) <- c(sdate = 1, latt = 181, lonn = 360)
+expect_error(
+  PlotRobinson(tmp), 
+  "Cannot find known longitude name in data dimension. Please define parameter 'lon_dim'."
+)
+expect_error(
+  PlotRobinson(data1, lon = 1:10),
+  "Length of parameter 'lon' should match longitude dimension in 'data'."
+)
+# target_proj
+expect_error(
+  PlotRobinson(data1, lon = lon1_1, lat = lat1, target_proj = NULL),
+  "Parameter 'target_proj' cannot be NULL."
+)
+# legend
+expect_error(
+  PlotRobinson(data1, lon = lon1_1, lat = lat1, legend = 1),
+  "Parameter 'legend' must be NULL, ggplot2 or s2dv."
+)
+# style
+expect_error(
+  PlotRobinson(data1, lon = lon1_1, lat = lat1, style = 'line'),
+  "Parameter 'style' must be 'point' or 'polygon'."
+)
+# dots
+expect_error(
+  PlotRobinson(data1, lon = lon1_1, lat = lat1, dots = c(1, 0)),
+  "Parameter 'dots' must have two dimensions named as longitude and latitude dimensions in 'data'."
+)
+tmp <- dots1; dim(tmp) <- c(sdate = 1, lat = 360, lon = 181)
+expect_error(
+  PlotRobinson(data1, lon = lon1_1, lat = lat1, dots = tmp),
+  "Parameter 'dots' must have the same dimensions as 'data'."
+)
+tmp <- dots1; tmp[1] <- NA
+expect_error(
+  PlotRobinson(data1, lon = lon1_1, lat = lat1, dots = tmp),
+  "Parameter 'dots' must have only TRUE/FALSE or 0/1."
+)
+# colNA
+expect_error(
+  PlotRobinson(data1, lon = lon1_1, lat = lat1, colNA = 'fair'),
+  "Parameter 'colNA' must be a valid colour identifier."
+)
+# toptitle
+expect_error(
+  PlotRobinson(data1, lon = lon1_1, lat = lat1, toptitle = 1:10),
+  "Parameter 'toptitle' must be a character string."
+)
+# caption
+expect_error(
+  PlotRobinson(data1, lon = lon1_1, lat = lat1, caption = 1:10),
+  "Parameter 'caption' must be a character string."
+)
+# crop_coastlines
+expect_error(
+  PlotRobinson(data1, lon = lon1_1, lat = lat1, crop_coastlines = c(x1 = 1, x2 = 10, y1 = -10, y2 = 10)),
+  "Parameter 'crop_coastlines' needs to have names 'latmax', 'latmin', 'lonmax', 'lonmin'."
+)
+# point_size
+expect_error(
+  PlotRobinson(data1, lon = lon1_1, lat = lat1, point_size = 'small'),
+  "Parameter 'point_size' must be a number."
+)
+
+})
+
+#--------------------------------------------------------
+# PlotRobinson
+save_PlotRobinson <- function(...) {
+  path <- tempfile(fileext = ".png")
+  do.call(PlotRobinson, list(..., fileout = path))
+  path
+}
+
+test_that("2. data1: Global", {
+
+# up triangle end only
+expect_snapshot_file(
+  save_PlotRobinson(data1, lon = lon1_1, lat = lat1, dots = dots1, brks = seq(-12.6, 12.6, 1), vertical = F, toptitle = 'PlotRobinson', caption = 'unit test: 1a'),
+  name = 'PlotRobinson_1a.png' 
+)
+# NAs
+expect_snapshot_file(
+  save_PlotRobinson(data1_NA, lon = lon1_1, lat = lat1, dots = dots1, brks = seq(-13, 13, 2), toptitle = 'PlotRobinson', caption = 'unit test: 1b', bar_extra_margin = c(5, 0, 5, 0)),
+  name = 'PlotRobinson_1b.png'
+)
+
+expect_snapshot_file(
+  save_PlotRobinson(data1, lon = lon1_1, lat = lat1, mask = mask1, brks = seq(-11, 11, 2), color_fun = clim.palette('purpleorange'), colNA = 'red', legend = 'ggplot2', toptitle = 'PlotRobinson', caption = 'unit test: 1c'),
+  name = 'PlotRobinson_1c.png'
+)
+
+
+})
+
+
+test_that("3. data2: Europe", {
+
+expect_snapshot_file(
+  save_PlotRobinson(data2, lon = lon2, lat = lat2, dots = dots2, target_proj = 102014,
+    brks = seq(-10, 0, 2), cols = c("#FDD49E", "#FDBB84", "#FC8D59", "#E34A33", "#B30000"), col_inf = "#FEF0D9",
+    toptitle = 'PlotRobinson', caption = 'unit test: 1a\n projection: Lambert Europe (ESRI:102014)', point_size = 1.1, title_size = 12,
+    width = 6, vertical = F, bar_extra_margin = c(0.5, 4, 0.5, 4), crop_coastlines = c(latmin = 27, latmax = 70, lonmin = -20, lonmax = 40)),
+  name = 'PlotRobinson_2a.png'
+)
+
+#polygon
+expect_snapshot_file(
+  save_PlotRobinson(data2, lon = lon2, lat = lat2, dots = dots2, target_proj = 102014,
+    brks = seq(-10, 0, 2), cols = c("#FDD49E", "#FDBB84", "#FC8D59", "#E34A33", "#B30000"), col_inf = "#FEF0D9",
+    toptitle = 'PlotRobinson', caption = 'unit test: 1a\n projection: Lambert Europe (ESRI:102014)', point_size = 1.1, title_size = 12,
+    width = 6, vertical = F, bar_extra_margin = c(0.5, 4, 0.5, 4), crop_coastlines = c(latmin = 27, latmax = 70, lonmin = -20, lonmax = 40), style = 'polygon'),
+  name = 'PlotRobinson_2a.png'
+)
+
+
+})
diff --git a/vignettes/vis_proj_eu_1.png b/vignettes/vis_proj_eu_1.png
new file mode 100644
index 0000000000000000000000000000000000000000..7275ea88d10278c12702bb3419e4aa7802839af4
Binary files /dev/null and b/vignettes/vis_proj_eu_1.png differ
diff --git a/vignettes/vis_proj_eu_2.png b/vignettes/vis_proj_eu_2.png
new file mode 100644
index 0000000000000000000000000000000000000000..dae2c4d5b6d8380365e623bc39c18867e92c129b
Binary files /dev/null and b/vignettes/vis_proj_eu_2.png differ
diff --git a/vignettes/vis_proj_global_1.png b/vignettes/vis_proj_global_1.png
new file mode 100644
index 0000000000000000000000000000000000000000..dbd18c15d35392f1993e20fe9cee8c053f1cc6cf
Binary files /dev/null and b/vignettes/vis_proj_global_1.png differ
diff --git a/vignettes/vis_proj_global_2.png b/vignettes/vis_proj_global_2.png
new file mode 100644
index 0000000000000000000000000000000000000000..e82903bc2c2b136b0c9c49618d8c35e1162f3294
Binary files /dev/null and b/vignettes/vis_proj_global_2.png differ
diff --git a/vignettes/vis_proj_np_1.png b/vignettes/vis_proj_np_1.png
new file mode 100644
index 0000000000000000000000000000000000000000..a3f4a1995290c4813729d4a9f152aff49a5d5286
Binary files /dev/null and b/vignettes/vis_proj_np_1.png differ
diff --git a/vignettes/vis_proj_np_2.png b/vignettes/vis_proj_np_2.png
new file mode 100644
index 0000000000000000000000000000000000000000..b2a7e49ed4726fd2278e8e6345289aca39b32d9a
Binary files /dev/null and b/vignettes/vis_proj_np_2.png differ
diff --git a/vignettes/visualization_projection.md b/vignettes/visualization_projection.md
new file mode 100644
index 0000000000000000000000000000000000000000..404f86cb3440ac3e1ac9148eed29d81c2cc21e94
--- /dev/null
+++ b/vignettes/visualization_projection.md
@@ -0,0 +1,236 @@
+---
+author: "An-Chi Ho"
+date: "2023-05-03"
+output: html_document
+vignette: >
+  %\VignetteEngine{R.rsp::md}
+  %\VignetteIndexEntry{Visualisation}
+  %\usepackage[utf8]{inputenc}
+---
+
+# Visualization of different projections
+
+This vignette shows how to use function `PlotRobinson()`, which can plot maps with the projection of your choice. Since packages `sf` and `rnaturalearth` are required by the function, make sure you have modules `PROJ`, `GEOS`, `GDAL` loaded.
+The projection types and IDs can be found on https://epsg.io/.
+
+```r
+## Load packages
+library(sf)
+library(ggplot2)
+library(rnaturalearth)
+library(RColorBrewer)
+library(cowplot)
+library(startR)
+library(s2dv)
+
+## Source function
+#URL: https://earth.bsc.es/gitlab/es/s2dv/-/raw/develop-PlotRobinson/R/PlotRobinson.R
+source("/esarchive/scratch/aho/tmp/PlotRobinson.R")
+## Source unexported s2dv functions used in PlotRobinson
+.KnownLonNames <- s2dv:::.KnownLonNames
+.KnownLatNames <- s2dv:::.KnownLatNames
+.IsColor <- s2dv:::.IsColor
+.warning <- s2dv:::.warning
+```
+Specify where you want to save the figures later.
+```r
+fileout_path <- "/esarchive/scratch/<userid>"
+```
+
+## Global map
+
+### 1. Load data
+
+We use `startR::Start` to load system5c3s monthly temperature data from esarchive. To focus on visualization, we simply load one file with one ensemble and one time step.
+
+```r
+repos <- "/esarchive/exp/ecmwf/system5c3s/monthly_mean/$var$_f6h/$var$_$sdate$.nc"
+data <- Start(dat = list(list(name = 'system5c3s', path = repos)),
+              var = 'tas',
+              sdate = '20170101',
+              ensemble = indices(1),
+              time = indices(1),
+              lon = 'all',
+              lon_reorder = CircularSort(-180, 180),
+              lat = 'all',
+              lat_reorder = Sort(decreasing = T),
+              synonims = list(lat = c('lat', 'latitude'),
+                              lon = c('lon', 'longitude')),
+              return_vars = list(time = 'sdate', longitude = 'dat', latitude = 'dat'),
+              retrieve = T)
+lon <- attributes(data)$Variable$system5c3s$lon
+lat <- attributes(data)$Variable$system5c3s$lat
+data <- data - 273.15 # K to degC
+```
+
+### 2. Create mask
+
+We can create a logical or [0, 1] array to serve as mask in the plot.
+```r
+mask <- drop(data)
+mask[which(mask > -20 & mask < -10)] <- 0
+mask[which(mask != 0)] <- 1
+```
+
+### 3. Create colorbar breaks and colors
+
+Create the adequete colorbar elements for the plot
+```r
+## have exceeded values
+brks1 <- seq(-40, 40, 10)
+cols1 <- brewer.pal(9, 'PuOr')[-9]
+## cover all values
+#brks1 <- seq(-ceiling(max(abs(range(data)))), ceiling(max(abs(range(data)))), length.out = 11)
+```
+### 4. Plot
+
+You can find more projection types through the URL above. `PlotRobinson()` is mainly designed for Robinson projection, but it is possible to use other projections as well.
+E.g., Mollweide projection is "ESRI:54009" or "ESRI:53009".
+Note that different PROJ and GDAL module versions have different CRS codes. If on Nord3v2, please use "ESRI:53030" or "+proj=robin" to get the correct Robinson projection.
+
+There are two types of legend, 'ggplot2' or 's2dv'. You can also set `legend` as `NULL` to not plot legend.
+
+```r
+## legend = 'ggplot2'
+PlotRobinson(data, lon, lat, brks = brks1, cols = NULL, legend = 'ggplot2',
+             target_proj = 54030, style = 'point',
+             toptitle = "system5c3s monthly tas Jan 2017", caption = "Projection: Robinson",
+             dots = mask, units = 'degC', colNA = 'red',
+             width = 10, height = 6)
+
+## legend = 's2dv'
+PlotRobinson(data, lon, lat, brks = brks1, cols = cols1, legend = 's2dv',
+             target_proj = 54030, style = 'point',
+             toptitle = "system5c3s monthly tas Jan 2017", caption = "Projection: Robinson",
+             dots = mask, units = "degC", col_inf = 'brown',
+             width = 10, height = 5,
+             bar_extra_margin = c(0, 5, 0, 5),
+             vertical = F)
+```
+
+<img src="vis_proj_global_1.png" width="600" />
+
+<img src="vis_proj_global_2.png" width="600" />
+
+
+## European region
+
+### 1. Subset data
+
+Crop the european area from the global data we loaded above.
+
+```r
+data_eu <- ClimProjDiags::Subset(data, c('lat', 'lon'), list(21:61, 161:221))
+lon_eu <- lon[161:221]
+lat_eu <- lat[21:61]
+```
+### 2. Crop mask
+
+```r
+mask_eu <- ClimProjDiags::Subset(mask, c('lat', 'lon'), list(21:61, 161:221))
+```
+
+### 3. Create colorbar breaks
+
+```r
+## have exceeded values
+brks_eu <- seq(-18, 18, length.out = 11)
+## cover all values
+#brks_eu <- seq(-ceiling(max(abs(range(data_eu)))), ceiling(max(abs(range(data_eu)))), length.out = 11)
+```
+
+### 4. Plot
+In the previous plots, we use `dots` to mark the masked area. This time, we use parameter `mask` to remove the unwanted data. The removed points are left as blank.
+
+The parameter `crop_coastlines` helps to adjust the map coastline plotting. You can test and find the best range. Note that this parameter is not suitable for the case that longitude range exceeds half of the earth (i.e., 180 degrees.)
+
+We can also choose to plot the map by polygon instead of point. Note that polygon is relatively **slow** so it takes a long time for global map, and it **does not work on our workstation** for some technical reason.
+
+Note that `target_proj = 102014` does not work well on Nord3v2. Use PROJ.4 string `target_proj = "+proj=lcc +lat_0=30 +lon_0=10 +lat_1=43 +lat_2=62 +x_0=0 +y_0=0 +ellps=intl +units=m +no_defs +type=crs"` instead.
+
+```r
+## style = point
+PlotRobinson(data_eu, lon_eu, lat_eu, brks = brks_eu, cols = NULL, legend = 's2dv',
+             target_proj = 102014, style = 'point',
+             toptitle = "system5c3s monthly tas Jan 2017",
+             caption = paste0("Projection: Lambert Europe\nESRI:102014\n", "Europe"),
+             mask = mask_eu, units = "degC",
+             width = 8, height = 6, point_size = 1.3,
+             crop_coastlines = c(lonmin = -20, lonmax = 40, latmin = 27, latmax = 70),
+             triangle_ends = c(T, T), bar_extra_margin = c(0.4, 4, 0.4, 4),
+             vertical = F)
+
+## style = polygon
+PlotRobinson(data_eu, lon_eu, lat_eu, brks = brks_eu, cols = NULL, legend = 's2dv',
+             target_proj = 102014, style = 'polygon',
+             toptitle = "system5c3s monthly tas Jan 2017",
+             caption = paste0("Projection: Lambert Europe\nESRI:102014\n", "Europe"),
+             mask = mask_eu, units = "degC",
+             width = 8, height = 6,
+             crop_coastlines = c(lonmin = -20, lonmax = 40, latmin = 27, latmax = 70),
+             triangle_ends = c(T, T), bar_extra_margin = c(0.4, 4, 0.4, 4),
+             vertical = F)
+```
+
+<img src="vis_proj_eu_1.png" width="450" />
+
+<img src="vis_proj_eu_2.png" width="450" />
+
+
+## North Pole region
+
+### 1. Subset data
+
+Crop the north pole area from the global data we loaded above.
+```r
+data_np <- ClimProjDiags::Subset(data, c('lat'), list(1:31))
+lon_np <- lon
+lat_np <- lat[1:31]
+```
+### 2. Crop mask
+
+```r
+mask_np <- ClimProjDiags::Subset(mask, c('lat'), list(1:31))
+```
+
+### 3. Create colorbar breaks
+
+```r
+## have exceeded values
+brks_np <- seq(-40, 0, length.out = 11)
+## cover all values
+#brks_np <- seq(-ceiling(max(abs(range(data_np)))), ceiling(max(abs(range(data_np)))), length.out = 11)
+```
+
+### 4. Plot
+
+We compare `PlotRobinson()` with `PlotStereoMap()`, which is specifically for the stereographic projection and it uses polygon to plot. 
+Note that we cannot use `crop_coastlines` here in `PlotRobinson()` because the longitude range is global.
+It is recommended saving the figures for these two plots (especially `PlotStereoMap()`) because polygon plotting is time-consuming when it is in the pop-up window.
+
+The projection used here is Arctic Polar Stereographic (EPSG:3995).
+```r
+#NOTE: A bit slow, be patience.
+PlotRobinson(data_np, lon_np, lat_np, brks = brks_np, cols = NULL, legend = 's2dv',
+             target_proj = 3995, style = 'polygon',
+             toptitle = "system5c3s monthly tas Jan 2017",
+             caption = paste0("Projection: WGS 84 / Arctic Polar Stereographic\nEPSG:3995\n", "North Pole"),
+             dots = mask_np, dots_size = 1.5,
+             units = "degC", color_fun = clim.palette("orangepurple"),
+             width = 8, height = 8,
+             triangle_ends = c(T, T), bar_extra_margin = c(0.4, 4, 0.4, 4),
+             vertical = F, fileout = paste0(fileout_path, '/PlotRobinson_np.png'))
+
+# Compare with PlotStereoMap
+#NOTE: Save the figure instead of plotting in pop-up window because it is slow.
+PlotStereoMap(data_np, lon = lon_np, lat = lat_np, brks = brks_np, cols = NULL,
+              color_fun = clim.palette("orangepurple"),
+              toptitle = "system5c3s monthly tas Jan 2017 (PlotStereoMap)",
+              width = 8, height = 8, filled.continent = F,
+              dots = !mask_np, units = "degC", title_scale = 0.8,
+              fileout = paste0(fileout_path, '/PlotStereoMap_np.png'))
+```
+
+<img src="vis_proj_np_1.png" width="400" />   
+
+<img src="vis_proj_np_2.png" width="400" />