diff --git a/inst/doc/FAQ.md b/inst/doc/FAQ.md
index a8e49db9024f6e063a2fbe59db9274a4efd20b18..885e49c5250b5434fdc5b9587636495b99b2218a 100644
--- a/inst/doc/FAQ.md
+++ b/inst/doc/FAQ.md
@@ -5,6 +5,7 @@ This document intends to be the first reference for any doubts that you may have
 ## Index
 1. **How to**
    1. [Global Map with non-standard longitudinal boundaries](#1-global-map-with-non-standard-longitudinal-boundaries)
+   2. [Plot a region boundary from shape file with PlotEquiMap](#2-plot-a-region-boundary-from-shape-file-with-plotequimap)
 
 2. **Something goes wrong...**  
    1. [CDORemap() returns errors or warnings with specific module versions](#1-cdoremap-returns-errors-or-warnings-with-specific-module-versions)
@@ -45,6 +46,49 @@ Note: You can adjust many parameters to visualize the plot, here we are just sho
 If you want to add other information to the plot (e.g.: hatching, points, countours, ...), you can add it just before ColorBar() function.
 
 
+### 2. Plot a region boundary from shape file with PlotEquiMap()
+
+PlotEquiMap() provides the option to plot the country boundaries by argument
+`country.borders = TRUE`. If you want to plot specific area, you can save 
+the shapefile data as an RDS file or read the shapefile as a list, and use 
+argument `shapefile` to plot it. 
+
+If you have more than one shapefile to plot, or you want to add additional 
+shapefile with `country.borders = TRUE`, you have to plot it on top of the 
+figure after PlotEquiMap(). 
+The following script shows an example of plotting the Douro region on top of
+default country boundaries. The data in this example is synthetic.
+To use "rgdal" package, `module load GDAL` first (see wiki [R tips 3. How to load dependencies of R package rgdal](https://earth.bsc.es/wiki/doku.php?id=tools:Rtools&s%5B%5D=Rtools#r_tips) for details.)
+
+```r
+library(rgdal)
+library(sp)
+library(maps)
+library(s2dv)
+
+# Read the Douro shapefile and transform:
+shp1 <- readOGR(dsn = '/esarchive/scratch/cchou/MEDGOLD/others/Shapefiles/', layer = 'RDD_ETRS89')
+s1 <- spTransform(shp1, CRS("+proj=longlat"))
+douro <- SpatialPolygons2map(s1)
+str(douro)
+#List of 4
+# $ x    : num [1:37826] -7.77 -7.77 -7.77 -7.77 -7.77 ...
+# $ y    : num [1:37826] 41.1 41.1 41.1 41.1 41.1 ...
+# $ names: chr [1:10] "0:1" "0:2" "0:3" "0:4" ...
+# $ range: num [1:4] -7.92 -6.75 40.92 41.53
+# - attr(*, "class")= chr "map"
+
+# Synthetic data
+dat <- 1 : (21 * 10)
+dim(dat) <- c(lon = 21, lat = 10)
+
+# Plot Douro region on top of figure that has country boundaries plotted
+PlotEquiMap(dat, lon = c(1:10, 350:360), lat = 36 : 45, drawleg = FALSE, filled.continents = FALSE, country.border = TRUE)
+map(douro, interior = FALSE, add = TRUE, lwd = 1, col = "red")
+
+# Plot Douro region solely, using argument 'shapefile'
+PlotEquiMap(dat, lon = c(1:10, 350:360), lat = 36 : 45, drawleg = FALSE, filled.continents = FALSE, shapefile = douro)
+```
 
 ## 2. Something goes wrong...  
 
@@ -65,6 +109,3 @@ R/3.6.1-foss-2015a-bare
 HDF5/1.10.5-foss-2015a  
 
 
-
-
-