diff --git a/DESCRIPTION b/DESCRIPTION
index 2c87e1472aeed292a148f0e89adbc82e9fc0453b..98dabd6ad06c4afc0d409099f234725053cba7f8 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -22,6 +22,7 @@ Authors@R: c(
     person("Raul", "Marcos", , "raul.marcos@bsc.es", role = "ctb"),
     person("Niti", "Mishra", , "niti.mishra@bsc.es", role = "ctb"),
     person("Jesus", "Peña", , "jesus.pena@bsc.es", role = "ctb"),
+    person("Francesc", "Roura-Adserias", , "francesc.roura@bsc.es", role = "ctb"),
     person("Silvia", "Terzago", , "s.terzago@isac.cnr.it", role = "ctb"),
     person("Danila", "Volpi", , "d.volpi@isac.cnr.it", role = "ctb"),
     person("BSC-CNS", role = c("cph")))
diff --git a/NEWS.md b/NEWS.md
index 86f68cf453151690ec56fb6640686a3e4ee95db6..04f986fb49d99cddcb1a5d19837fbbf69f24b556 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -2,7 +2,8 @@
 **Submission date to CRAN: XX-06-2020**
 
 - Fixes
-    + CST_Anomaly handles exp, obs or both.
+    + CST_Anomaly handles exp, obs or both
+    + PlotForecastPDF vignette displays figures correctly.
 
 
 ### CSTools 3.0.0
diff --git a/R/PlotForecastPDF.R b/R/PlotForecastPDF.R
index c40cc32c3b57211e9d14f2555c67541ec18b523d..14e6cd3aadca4cda3f6d040e526b8bcc5a0b206a 100644
--- a/R/PlotForecastPDF.R
+++ b/R/PlotForecastPDF.R
@@ -220,7 +220,7 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits = NULL, obs = N
         })
         attr <- attr(hatch.ls, "split_labels")
         for (i in 1:length(hatch.ls)) {
-            hatch.ls[[i]] <- cbind(hatch.ls[[i]], attr[i, ])
+            hatch.ls[[i]] <- cbind(hatch.ls[[i]], attr[i, ], row.names = NULL)
         }
         hatch.df <- do.call("rbind", hatch.ls)
         # Compute max y for each extreme category
@@ -229,7 +229,7 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits = NULL, obs = N
         })
         attr <- attr(max.ls, "split_labels")
         for (i in 1:length(max.ls)) {
-            max.ls[[i]] <- cbind(max.ls[[i]], attr[i, ])
+            max.ls[[i]] <- cbind(max.ls[[i]], attr[i, ], row.names = NULL)
         }
         max.df <- do.call("rbind", max.ls)
     }
@@ -323,7 +323,8 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits = NULL, obs = N
         pct2$lab.pos <-  as.vector(apply(extreme.limits, 1, function(x) {c(x[1], NA, x[2])}))
         pct2 <- merge(pct2, max.df, by = c("init", "extremes"))
         # include potentially missing groups
-        pct2 <- pct2[CJ(levels(pct2$init), factor(c("Below P10", "Normal", "Above P90"), 
+        pct2 <- pct2[CJ(factor(levels(pct2$init), levels = levels(pct2$init)), 
+            factor(c("Below P10", "Normal", "Above P90"), 
             levels = c("Below P10", "Normal", "Above P90"))), ]
     }
     #------------------------
diff --git a/vignettes/Figures/PlotForecastPDF_ex1.png b/vignettes/Figures/PlotForecastPDF_ex1.png
new file mode 100644
index 0000000000000000000000000000000000000000..d7067133054902e6b661f828f8a89b0bb285aa68
Binary files /dev/null and b/vignettes/Figures/PlotForecastPDF_ex1.png differ
diff --git a/vignettes/Figures/PlotForecastPDF_ex2.png b/vignettes/Figures/PlotForecastPDF_ex2.png
new file mode 100644
index 0000000000000000000000000000000000000000..67b4af1e1b523c564047fb545baecc11a4dbac2c
Binary files /dev/null and b/vignettes/Figures/PlotForecastPDF_ex2.png differ
diff --git a/vignettes/Figures/PlotForecastPDF_ex3.png b/vignettes/Figures/PlotForecastPDF_ex3.png
new file mode 100644
index 0000000000000000000000000000000000000000..10b1e63b4766e8aff8d10419b849d8677f197c58
Binary files /dev/null and b/vignettes/Figures/PlotForecastPDF_ex3.png differ
diff --git a/vignettes/Figures/PlotForecastPDF_ex4.png b/vignettes/Figures/PlotForecastPDF_ex4.png
new file mode 100644
index 0000000000000000000000000000000000000000..34df30d73e73f1d1d8e5ff4c34e57c2d487f4c2f
Binary files /dev/null and b/vignettes/Figures/PlotForecastPDF_ex4.png differ
diff --git a/vignettes/PlotForecastPDF.Rmd b/vignettes/PlotForecastPDF.Rmd
index c6d8961828e70e8a4b17b4bb6a4144258402d140..98b2ae144003e7d876689a827325cdf0eadf9fdc 100644
--- a/vignettes/PlotForecastPDF.Rmd
+++ b/vignettes/PlotForecastPDF.Rmd
@@ -1,5 +1,5 @@
 ---
-author: "Francesc Roura and Llorenç Lledó"
+author: "Francesc Roura-Adserias and Llorenç Lledó"
 date: "`r Sys.Date()`"
 output: rmarkdown::html_vignette
 vignette: >
@@ -23,8 +23,9 @@ library(CSTools)
 The first step is to put your forecasts in an appropriate format. For this vignette we generate some random values from two normal distributions. The PlotForecastPDF by default will plot the ensemble members, the estimated density distributions and the tercile probabilities.
 
 ```{r,fig.show = 'hide',warning=F}
-fcst <- data.frame(fcst1=rnorm(mean=25,sd=3,n=30),fcst2=rnorm(mean=23,sd=4.5,n=30))
-PlotForecastPDF(fcst,tercile.limits=c(20,26))
+fcst <- data.frame(fcst1 = rnorm(mean = 25, sd = 3, n = 30), 
+	fcst2 = rnorm(mean = 23, sd = 4.5, n = 30))
+PlotForecastPDF(fcst, tercile.limits = c(20, 26))
 ```
 
 ![Example 1](./Figures/PlotForecastPDF_ex1.png)
@@ -33,15 +34,47 @@ PlotForecastPDF(fcst,tercile.limits=c(20,26))
 Changing the title, the forecast labels or the units will be needed in most cases.
 
 ```{r,fig.show = 'hide',warning=F}
-fcst <- data.frame(fcst1=rnorm(mean=25,sd=3,n=30),fcst2=rnorm(mean=23,sd=4.5,n=30))
-PlotForecastPDF(fcst,tercile.limits=c(20,26),var.name="Temperature (ºC)",title="Forecasts valid on 2019-01-01 at Sunny Hills",fcst.names = c("model a","model b"))
+fcst <- data.frame(fcst1 = rnorm(mean = 25, sd = 3, n = 30),
+	fcst2 = rnorm(mean = 23, sd = 4.5, n = 30))
+PlotForecastPDF(fcst, tercile.limits = c(20, 26), var.name = "Temperature (ºC)",
+	title = "Forecasts valid on 2019-01-01 at Sunny Hills", 
+	fcst.names = c("model a", "model b"))
 ```
 ![Example 2](./Figures/PlotForecastPDF_ex2.png)
 
 ### 3.- Adding extremes and observed values
 We can add the probability of extreme values and the observed values. The tercile and extreme limits can be specified for each panel separately, as well as the observed values.
+
+```{r,fig.show = 'hide',warning=F}
+fcst <- data.frame(fcst1 = rnorm(mean = 25, sd = 3, n = 30), 
+	fcst2 = rnorm(mean = 28, sd = 4.5, n = 30), fcst3 = rnorm(mean = 17, sd = 3, n = 30))
+PlotForecastPDF(fcst, tercile.limits = rbind(c(20, 26), c(22, 28), c(15, 22)), 
+	var.name = "Temperature (ºC)", title = "Forecasts at Sunny Hills", 
+	fcst.names = c("January", "February", "March"), obs = c(21, 24, 17), 
+	extreme.limits = rbind(c(18, 28), c(20, 30), c(12, 24)))
+```
+
+The same example using a forecast in array format is provided.
 ```{r,fig.show = 'hide',warning=F}
-fcst <- data.frame(fcst1=rnorm(mean=25,sd=3,n=30),fcst2=rnorm(mean=28,sd=4.5,n=30),fcst3=rnorm(mean=17,sd=3,n=30))
-PlotForecastPDF(fcst,tercile.limits=rbind(c(20,26),c(22,28),c(15,22)),var.name="Temperature (ºC)",title="Forecasts at Sunny Hills",fcst.names = c("January","February","March"),obs=c(21,24,17),extreme.limits = rbind(c(18,28),c(20,30),c(12,24)))
+fcst <- array(cbind(cbind(rnorm(mean = 25, sd = 3, n = 30),
+	rnorm(mean = 23, sd = 4.5, n = 30)), rnorm(mean = 17, sd = 3, n = 30)), 
+	dim = c(members = 30, 3))
+PlotForecastPDF(fcst, tercile.limits = rbind(c(20, 26), c(22, 28), c(15, 22)), 
+	var.name = "Temperature (ºC)", title = "Forecasts at Sunny Hills", 
+	fcst.names = c("January", "February", "March"), obs = c(21, 24, 17), 
+	extreme.limits = rbind(c(18, 28), c(20, 30), c(12, 24)))
 ```
 ![Example 3](./Figures/PlotForecastPDF_ex3.png)
+
+### 4.- Example using lonlat_data
+An example using the lonlat_data from CSTools is provided.
+
+```{r,fig.show = 'hide',warning=F}
+fcst <- data.frame(fcst1 = lonlat_data$exp$data[1,,1,1,1,1] - 273.15,
+                   fcst2 = lonlat_data$exp$data[1,,1,2,1,1] - 273.15)
+PlotForecastPDF(fcst, tercile.limits = c(5, 7), extreme.limits = c(4, 8), 
+  var.name = "Temperature (ºC)",
+  title = "Forecasts valid on 2000-11 at sample mediterranean region", 
+  fcst.names = c("November", "December"))
+```
+![Example 4](./Figures/PlotForecastPDF_ex4.png)