diff --git a/.Rbuildignore b/.Rbuildignore
index 91114bf2f2bba5e0c5252e75018da19b869776f1..bcf15358754c52b155221f2bb7f1caeabfb81988 100644
--- a/.Rbuildignore
+++ b/.Rbuildignore
@@ -1,2 +1,5 @@
 ^.*\.Rproj$
 ^\.Rproj\.user$
+^data-raw$
+^README\.Rmd$
+^README-.*\.png$
diff --git a/DESCRIPTION b/DESCRIPTION
index 27e268813f9e581af0a4ba039fde37a0dfd1a647..df1ccdef6ed243827892a514dd3b3f6ab30efb8b 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -14,6 +14,10 @@ URL: https://gitlab.rrz.uni-hamburg.de/saskiaotto/IEAtools
 BugReports: https://gitlab.rrz.uni-hamburg.de/saskiaotto/IEAtools/issues
 Depends: 
     R(>= 4.0.0)
+Imports:
+    stats,
+    graphics,
+    magrittr
 Encoding: UTF-8
 LazyData: true
 RoxygenNote: 7.1.1
diff --git a/NAMESPACE b/NAMESPACE
index d75f824ec6278db24891505b14ab3d915514dba7..b8ed91efa77e50253361407d52624b2f79da78f8 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -1 +1,5 @@
-exportPattern("^[[:alpha:]]+")
+# Generated by roxygen2: do not edit by hand
+
+export("%>%")
+export(trafficlight)
+importFrom(magrittr,"%>%")
diff --git a/R/hello.R b/R/hello.R
deleted file mode 100644
index 3d348f227d8a6bb6214cb759621bcf33e2fe4147..0000000000000000000000000000000000000000
--- a/R/hello.R
+++ /dev/null
@@ -1,18 +0,0 @@
-# Hello, world!
-#
-# This is an example function named 'hello' 
-# which prints 'Hello, world!'.
-#
-# You can learn more about package authoring with RStudio at:
-#
-#   http://r-pkgs.had.co.nz/
-#
-# Some useful keyboard shortcuts for package authoring:
-#
-#   Install Package:           'Cmd + Shift + B'
-#   Check Package:             'Cmd + Shift + E'
-#   Test Package:              'Cmd + Shift + T'
-
-hello <- function() {
-  print("Hello, world!")
-}
diff --git a/R/trafficlight.R b/R/trafficlight.R
new file mode 100644
index 0000000000000000000000000000000000000000..50bd99d576a8fc4dea338a6618de5286f1f2ee4f
--- /dev/null
+++ b/R/trafficlight.R
@@ -0,0 +1,235 @@
+#' Create a traffic light plot
+#'
+#' This function creates for multiple time series an image plot where the color
+#' code is based on selected quantiles or evenly spaced intervals.
+#'
+#' @name trafficlight
+#' @param x matrix or data frame containing time series of multiple variables.
+#' @param time vector of time units that will be used for the x axis.
+#' @param sort_5yrmean logical; should the variables be sorted by the first 5yr mean? Default
+#'   is set to TRUE.
+#' @param sort_vec If specific order of variable is desired the sorting index should be
+#'   provided here.
+#' @param main A title.
+#' @param xlab The x axis title. Default is none.
+#' @param ylab The y axis title. Default is none.
+#' @param adj_xlab integer; vertical adjustment of the x axis title.
+#' @param adj_ylab integer; horizontal adjustment of the y axis title.
+#' @param cols A character vector with colors for each quantile or interval.
+#' @param method character; the type of methods to create the color code. Choose
+#'   between using "quantiles" (default) or "intervals".
+#' @param probs a vector of probabilities to calculate the quantiles using the function
+#'   \code{quantile}. Should include the probabilities 0 and 1.
+#'   numeric vector of probabilities with values in [0,1]. (Values up to
+#'   2e-14 outside that range are accepted and moved to the nearby endpoint.)
+#' @param quantile_type an integer between 1 and 9 selecting one of the nine quantile
+#'   algorithms detailed below to be used. Default is 7 (see also \code{\link{quantile}}).
+#' @param intervals logical; number of evenly-spaced intervals. Default is 5.
+#' @param hadj_x double; horizontal adjustment of the x labels.
+#' @param vadj_x double; vertical adjustment of the x labels.
+#' @param hadj_y double; horizontal adjustment of the x labels.
+#' @param vadj_y double; vertical adjustment of the x labels.
+#' @param madj_bottom double; adjustment of bottom margin relative to the
+#'   current setting.
+#' @param madj_left double; adjustment of left margin relative to the
+#'   current setting.
+#' @param madj_top double; adjustment of top margin relative to the
+#'   current setting.
+#' @param madj_right double; adjustment of right margin relative to the
+#'   current setting.
+#' @param legend_pos character; the legend is shown on the "top" (default),
+#'   "center", or "bottom" right of the plot.
+#' @param legend_intersp_x double; character interspacing factor for horizontal (x) spacing in legend.
+#' @param legend_intersp_y double; interspacing factor for vertical (y) line distances in legend.
+#' @param tick_x logical; set to TRUE if ticks on x-axis should be displayed.
+#' @param tick_y logical; set to TRUE if ticks on x-axis should be displayed.
+#' @param cex_x double; the magnification to be used for the x labels relative to the
+#'   current setting of cex.
+#' @param cex_y double; the magnification to be used for the y labels relative to the
+#'   current setting of cex.
+#' @param cex_xlab double; the magnification to be used for the x axis title relative
+#'   to the current setting of cex.
+#' @param cex_ylab double; the magnification to be used for the y axis title  relative
+#'   to the current setting of cex.
+#' @param cex_legend double; the magnification to be used for the legend relative to
+#'   the current setting of cex.
+#' @param cex_main double; the magnification to be used for the title relative to the
+#'   current setting of cex.
+#' @param respect logical; this argument controls whether a unit column-width is the
+#'   same physical measurement on the device as a unit row-height (default is TRUE).
+#'
+#'
+#' @export
+#' @examples
+#' df <- data.frame(aaaaaaaa = 1:20, bbbbbbbbbb= rnorm(20, 100, 30), c = 1:20 + rnorm(20))
+#' trafficlight(x = df, time = 1981:2000)
+#' df <- matrix(rnorm(100), ncol = 5)
+#' colnames(df) <- letters[1:5]
+#' trafficlight(x = df, time = 1981:2000, legend_pos = "bottom", method = "intervals")
+
+trafficlight <- function(x, time = NULL, sort_5yrmean = TRUE, sort_vec = NULL,
+  main = "", xlab = "", ylab = "", adj_xlab = NULL, adj_ylab = NULL,
+  cols = c("green", "yellowgreen", "yellow", "gold", "red"),
+  method = "quantiles",
+  probs = seq(0, 1, 0.2), quantile_type = 7, intervals = 5,
+  hadj_x = 0.5, vadj_x = -0.7, hadj_y = -0.7, vadj_y = 0.3,
+  madj_bottom = 0, madj_left = 0, madj_top = 0, madj_right = 0,
+  legend_pos = "top", legend_intersp_x = 1, legend_intersp_y = 1,
+  tick_x = FALSE, tick_y = FALSE,
+  cex_x = 0.8, cex_y = 0.8, cex_xlab = 1, cex_ylab = 1,
+  cex_legend = 0.8, cex_main = 1, respect = TRUE) {
+
+
+  ### Data input validation
+  if (!method %in% c("quantiles", "intervals")) {
+    stop("Choose as method either 'quantiles' or 'intervals'.")
+  }
+  if (method == "intervals") {
+    if (length(cols) != intervals) {
+      stop("The number of colours need to match your chosen interval number.")
+    }
+  }
+  # --------------------------
+
+  z <- x
+
+  # Save the dimensions (number of rows[1] and columns[2] ) in a vector
+  n <- dim(z)
+  # Get names for the y-axis
+  ylabel <- colnames(z)
+  if (!is.null(time)) {
+    xlabel <- time
+  } else {
+    xlabel <- 1:n[1]
+  }
+
+  # Converting the original values into quantiles or even intervals
+  convert2quantiles <- function(v, probs, type, var) {
+    br <- stats::quantile(v, probs = probs, na.rm = TRUE, type = type)
+    if (any(diff(br) == 0)) {
+      sel <- which(diff(br) == 0)
+      br[sel+1] <- br[sel+1] + br[sel]/10000
+      print(paste0("For variable '", var, "' the ", probs[sel+1]*100,
+        "%-quantile is the same as the ", probs[sel]*100,
+        "%-quantile. All values are grouped under the lower quantile!"))
+    }
+    qv <- cut(v, breaks = br, include.lowest = TRUE, labels = FALSE)
+    return(qv)
+  }
+
+  convert2intervals <- function(v, intervals) {
+    qv <- cut(v, breaks = intervals, include.lowest = TRUE, labels = FALSE)
+    return(qv)
+  }
+
+  zc <- z
+  if (method == "quantiles") {
+    for (i in 1:n[2]) {
+      zc[,i] <- convert2quantiles(zc[,i], probs, type = quantile_type, var = names(z)[i])
+      legend_txt <- paste0("< ", probs[-1]*100, "%")
+      nl <- length(legend_txt)
+      legend_txt[nl] <- paste0("> ", probs[nl]*100, "%")
+      legend_txt <- c(legend_txt, "missing value")
+    }
+  } else {
+    if (intervals %% 1 == 0) { # is full number?
+      for (i in 1:n[2]) {
+        zc[,i] <- convert2intervals(zc[,i], intervals)
+        legend_txt <- as.character(1:intervals)
+        legend_txt[1] <- paste0(legend_txt[1], " (low)")
+        nl <- length(legend_txt)
+        legend_txt[nl] <- paste0(legend_txt[nl], " (high)")
+        legend_txt <- c(legend_txt, "missing value")
+      }
+    } else {
+      stop(" If you want to use evenly-spaced intervals, provide an integer for the number of intervals!")
+    }
+  }
+
+
+  ### Sort variables according to settings
+  if (isTRUE(sort_5yrmean)) {
+    m5 <- vector(length = n[2])
+    # Then fill the vector with the standardised variable averages over the
+    # first five data points (mean of first 5 years - mean of full time series/
+    # standard deviation of full time series)
+    for (i in 1:(n[2])) {
+      m5[i] <- (mean(x[c(1:5), i], na.rm = TRUE) -
+          mean(x[, i], na.rm = TRUE)) /
+        stats::sd(x[, i], na.rm = TRUE) }
+    # Finally, order the variable and create an index vector
+    ordvar <- order(m5)
+  } else {
+    if (is.null(sort_vec)) {
+      ordvar <- n[2]:1
+    } else {
+      ordvar <- rev(sort_vec)
+    }
+  }
+  zc_sort <- as.matrix(zc[ ,ordvar])
+
+
+  ### Plot settings
+  x <- 1:n[1]
+  y <- 1:n[2]
+  # Position of legend:
+  if (legend_pos == "top") {
+    xleg <- max(x)+1
+    yleg <- max(y)+.5
+    yjustleg <- 1
+  } else {
+    if (legend_pos == "center") {
+      xleg <- max(x)+1
+      yleg <- max(y)/2+.5
+      yjustleg <- 0.5
+    } else {
+      if (legend_pos == "bottom") {
+        xleg <- max(x)+1
+        yleg <- min(y)-.5
+        yjustleg <- 0
+      } else {
+        stop("You need to choose as legend position 'bottom', 'center' or 'top'.")
+      }
+    }
+  }
+
+  mar <- c(2+madj_bottom, 5+madj_left, 1+madj_top, 8+madj_right)
+  if (nchar(xlab) > 0) mar[1] <- mar[1] + 1
+  if (nchar(ylab) > 0) mar[2] <- mar[2] + 1
+  if (nchar(main) > 0) mar[3] <- mar[3] + 2
+  if (is.null(adj_xlab)) {
+    adj_xlab <- mar[1]-1
+  } else {
+    adj_xlab <- adj_xlab
+  }
+  if (is.null(adj_ylab)) {
+    adj_ylab <- mar[2]-1
+  } else {
+    adj_ylab <- adj_ylab
+  }
+
+  ### Plot
+  graphics::layout(matrix(c(1,1,1,1), ncol = 2), widths = c(3.5,3.5), heights = c(2,2),
+    respect = respect)
+  graphics::par(mar = c(mar[1], mar[2], mar[3], mar[4]), oma = c(0.5,.5,0,0), xpd = TRUE)
+  graphics::image(x, y, z = zc_sort, zlim = c(1,5), col = cols,
+    axes = FALSE, xlab = "", ylab = "")
+  if (isTRUE(tick_x)) graphics::axis(1, at = x, tick = -.015, labels = NA)
+  if (isTRUE(tick_y)) graphics::axis(2, at = y, tick = -.015, labels = NA)
+  graphics::axis(1, at = x, tick = FALSE, labels = xlabel, cex.axis = cex_x, las = 3,
+    line = vadj_x, padj = hadj_x)
+  graphics::axis(2, at = y, tick = FALSE, labels = ylabel[ordvar], cex.axis = cex_y,
+    las = 1, line = hadj_y, padj = vadj_y)
+  graphics::box()
+  graphics::legend(x = xleg, y = yleg, legend = legend_txt, fill = c(cols, "white"),
+    cex = cex_legend, bty = "n", xjust = .1, yjust = yjustleg,
+    x.intersp = legend_intersp_x, y.intersp = legend_intersp_y, )
+  graphics::title(main, cex.main = cex_main)
+  if (nchar(xlab) > 0) {
+    graphics::mtext(text = xlab, side = 1, line = adj_xlab, cex = cex_xlab)
+  }
+  if (nchar(ylab) > 0) {
+    graphics::mtext(text = ylab, side = 2, line = adj_ylab, cex = cex_ylab)
+  }
+
+}
diff --git a/R/utils-pipe.R b/R/utils-pipe.R
new file mode 100644
index 0000000000000000000000000000000000000000..e79f3d80856576d6db32c9fb05cdbc02b28f4ff8
--- /dev/null
+++ b/R/utils-pipe.R
@@ -0,0 +1,11 @@
+#' Pipe operator
+#'
+#' See \code{magrittr::\link[magrittr:pipe]{\%>\%}} for details.
+#'
+#' @name %>%
+#' @rdname pipe
+#' @keywords internal
+#' @export
+#' @importFrom magrittr %>%
+#' @usage lhs \%>\% rhs
+NULL
diff --git a/man/hello.Rd b/man/hello.Rd
deleted file mode 100644
index 0fa7c4b8817591c2dff2b3997d2566320ac6d9fc..0000000000000000000000000000000000000000
--- a/man/hello.Rd
+++ /dev/null
@@ -1,12 +0,0 @@
-\name{hello}
-\alias{hello}
-\title{Hello, World!}
-\usage{
-hello()
-}
-\description{
-Prints 'Hello, world!'.
-}
-\examples{
-hello()
-}
diff --git a/man/pipe.Rd b/man/pipe.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..0eec752616534f156cac398f07bef8cc3f527935
--- /dev/null
+++ b/man/pipe.Rd
@@ -0,0 +1,12 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/utils-pipe.R
+\name{\%>\%}
+\alias{\%>\%}
+\title{Pipe operator}
+\usage{
+lhs \%>\% rhs
+}
+\description{
+See \code{magrittr::\link[magrittr:pipe]{\%>\%}} for details.
+}
+\keyword{internal}
diff --git a/man/trafficlight.Rd b/man/trafficlight.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..316f733988f862d14068de548effccab1e20400a
--- /dev/null
+++ b/man/trafficlight.Rd
@@ -0,0 +1,142 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/trafficlight.R
+\name{trafficlight}
+\alias{trafficlight}
+\title{Create a traffic light plot}
+\usage{
+trafficlight(
+  x,
+  time = NULL,
+  sort_5yrmean = TRUE,
+  sort_vec = NULL,
+  main = "",
+  xlab = "",
+  ylab = "",
+  adj_xlab = NULL,
+  adj_ylab = NULL,
+  cols = c("green", "yellowgreen", "yellow", "gold", "red"),
+  method = "quantiles",
+  probs = seq(0, 1, 0.2),
+  quantile_type = 7,
+  intervals = 5,
+  hadj_x = 0.5,
+  vadj_x = -0.7,
+  hadj_y = -0.7,
+  vadj_y = 0.3,
+  madj_bottom = 0,
+  madj_left = 0,
+  madj_top = 0,
+  madj_right = 0,
+  legend_pos = "top",
+  legend_intersp_x = 1,
+  legend_intersp_y = 1,
+  tick_x = FALSE,
+  tick_y = FALSE,
+  cex_x = 0.8,
+  cex_y = 0.8,
+  cex_xlab = 1,
+  cex_ylab = 1,
+  cex_legend = 0.8,
+  cex_main = 1,
+  respect = TRUE
+)
+}
+\arguments{
+\item{x}{matrix or data frame containing time series of multiple variables.}
+
+\item{time}{vector of time units that will be used for the x axis.}
+
+\item{sort_5yrmean}{logical; should the variables be sorted by the first 5yr mean? Default
+is set to TRUE.}
+
+\item{sort_vec}{If specific order of variable is desired the sorting index should be
+provided here.}
+
+\item{main}{A title.}
+
+\item{xlab}{The x axis title. Default is none.}
+
+\item{ylab}{The y axis title. Default is none.}
+
+\item{adj_xlab}{integer; vertical adjustment of the x axis title.}
+
+\item{adj_ylab}{integer; horizontal adjustment of the y axis title.}
+
+\item{cols}{A character vector with colors for each quantile or interval.}
+
+\item{method}{character; the type of methods to create the color code. Choose
+between using "quantiles" (default) or "intervals".}
+
+\item{probs}{a vector of probabilities to calculate the quantiles using the function
+\code{quantile}. Should include the probabilities 0 and 1.
+numeric vector of probabilities with values in [0,1]. (Values up to
+2e-14 outside that range are accepted and moved to the nearby endpoint.)}
+
+\item{quantile_type}{an integer between 1 and 9 selecting one of the nine quantile
+algorithms detailed below to be used. Default is 7 (see also \code{\link{quantile}}).}
+
+\item{intervals}{logical; number of evenly-spaced intervals. Default is 5.}
+
+\item{hadj_x}{double; horizontal adjustment of the x labels.}
+
+\item{vadj_x}{double; vertical adjustment of the x labels.}
+
+\item{hadj_y}{double; horizontal adjustment of the x labels.}
+
+\item{vadj_y}{double; vertical adjustment of the x labels.}
+
+\item{madj_bottom}{double; adjustment of bottom margin relative to the
+current setting.}
+
+\item{madj_left}{double; adjustment of left margin relative to the
+current setting.}
+
+\item{madj_top}{double; adjustment of top margin relative to the
+current setting.}
+
+\item{madj_right}{double; adjustment of right margin relative to the
+current setting.}
+
+\item{legend_pos}{character; the legend is shown on the "top" (default),
+"center", or "bottom" right of the plot.}
+
+\item{legend_intersp_x}{double; character interspacing factor for horizontal (x) spacing in legend.}
+
+\item{legend_intersp_y}{double; interspacing factor for vertical (y) line distances in legend.}
+
+\item{tick_x}{logical; set to TRUE if ticks on x-axis should be displayed.}
+
+\item{tick_y}{logical; set to TRUE if ticks on x-axis should be displayed.}
+
+\item{cex_x}{double; the magnification to be used for the x labels relative to the
+current setting of cex.}
+
+\item{cex_y}{double; the magnification to be used for the y labels relative to the
+current setting of cex.}
+
+\item{cex_xlab}{double; the magnification to be used for the x axis title relative
+to the current setting of cex.}
+
+\item{cex_ylab}{double; the magnification to be used for the y axis title  relative
+to the current setting of cex.}
+
+\item{cex_legend}{double; the magnification to be used for the legend relative to
+the current setting of cex.}
+
+\item{cex_main}{double; the magnification to be used for the title relative to the
+current setting of cex.}
+
+\item{respect}{logical; this argument controls whether a unit column-width is the
+same physical measurement on the device as a unit row-height (default is TRUE).}
+}
+\description{
+This function creates for multiple time series an image plot where the color
+code is based on selected quantiles or evenly spaced intervals.
+}
+\examples{
+df <- data.frame(aaaaaaaa = 1:20, bbbbbbbbbb= rnorm(20, 100, 30), c = 1:20 + rnorm(20))
+trafficlight(x = df, time = 1981:2000)
+df <- matrix(rnorm(100), ncol = 5)
+colnames(df) <- letters[1:5]
+trafficlight(x = df, time = 1981:2000, legend_pos = "bottom", method = "intervals")
+}