Changin color and pch of Taylor diagram in R?

Question

I want to change the pch shape and color of my reference data in taylor.diagram(), from plotrix package. If anyone can show me how to do it, i am really appreciated. The reference's model is highlighted by dark circle at the bottom of x axis. Thanks in advance. It is my script using

library(plotrix)

# windows()
taylor.diagram(ref=a1, model=b1, pos.cor=T, col=cols[1], pch=17,add=FALSE,
               pcex=3,main="",cex=1.3)
taylor.diagram(ref=a1, model=c1, pos.cor=T, col=cols[2], pch=17,add=TRUE,
               pcex=3)
taylor.diagram(ref=a1, model=d1, pos.cor=T, col=cols[3], pch=17,add=TRUE,
               pcex=3)
taylor.diagram(ref=a1, model=e1, pos.cor=T, col=cols[4], pch=17,add=TRUE,
               pcex=3)
taylor.diagram(ref=a1, model=f1, pos.cor=T, col=cols[5], pch=17,add=TRUE,
               pcex=3)
taylor.diagram(ref=a1, model=g1, pos.cor=T, col=cols[6], pch=15,add=TRUE,
               pcex=3)
title(main="Spatial Taylor Diagram of Precipitation (CRU)",cex.main=2)
legend("topright",cex=1.8,pt.cex=2,title="GCM-RCM models",
       legend=c("CNRMr1-ALADIN63","ECEARTHr12-CCLM","HadGEMr1-HIRHAM",
                "MPIr1-RACMO","NorESMr1-REMO","EOBS"),
       pch=c(17,17,17,17,17,15),col=cols)

Some trial data:

a <- array(c(1,4,8,5,6,2,2,7,3), dim=c(3,3))
b <- array(c(9,5,3,6,5,8,1,2,7) , dim=c(3,3))
c <- array(c(1,6,8,5,8,5,8,9,3) , dim=c(3,3))
a2 <- as.numeric(a)
b2 <- as.numeric(b)
c2 <- as.numeric(c)

a2 is the reference, b2 and c2 are the models

Answer 1

It looks like it's not possible using the existing function. The piece of code that plots the reference point looks like this:

points(sd.r, 0, cex = pcex)

So, there is no way to change anything else about the point than its cex value. Since the change is pretty minor, we could modify the function to add refcol and refpch arguments for the color and symbol shape of the reference.

taylor.diagram2 <- function (ref, model, add = FALSE, col = "red", pch = 19, pos.cor = TRUE, 
          xlab = "Standard deviation", ylab = "", main = "Taylor Diagram", 
          show.gamma = TRUE, ngamma = 3, gamma.col = 8, sd.arcs = 0, 
          ref.sd = FALSE, sd.method = "sample", grad.corr.lines = c(0.2, 
                                                                    0.4, 0.6, 0.8, 0.9), pcex = 1, cex.axis = 1, normalize = FALSE, 
          mar = c(4, 3, 4, 3), refcol = "black", refpch = 16, ...) {
  require("plotrix")
  grad.corr.full <- c(0, 0.2, 0.4, 0.6, 0.8, 0.9, 0.95, 0.99, 
                      1)
  R <- cor(ref, model, use = "pairwise")
  if (is.list(ref)) 
    ref <- unlist(ref)
  if (is.list(model)) 
    ref <- unlist(model)
  SD <- function(x, subn) {
    meanx <- mean(x, na.rm = TRUE)
    devx <- x - meanx
    ssd <- sqrt(sum(devx * devx, na.rm = TRUE)/(length(x[!is.na(x)]) - 
                                                  subn))
    return(ssd)
  }
  subn <- sd.method != "sample"
  sd.r <- SD(ref, subn)
  sd.f <- SD(model, subn)
  if (normalize) {
    sd.f <- sd.f/sd.r
    sd.r <- 1
  }
  maxsd <- 1.5 * max(sd.f, sd.r)
  oldpar <- par("mar", "xpd", "xaxs", "yaxs")
  if (!add) {
    par(mar = mar)
    if (pos.cor) {
      if (nchar(ylab) == 0) 
        ylab = "Standard deviation"
      plot(0, xlim = c(0, maxsd * 1.1), ylim = c(0, maxsd * 
                                                   1.1), xaxs = "i", yaxs = "i", axes = FALSE, main = main, 
           xlab = "", ylab = ylab, type = "n", cex = cex.axis, 
           ...)
      mtext(xlab, side = 1, line = 2.3)
      if (grad.corr.lines[1]) {
        for (gcl in grad.corr.lines) lines(c(0, maxsd * 
                                               gcl), c(0, maxsd * sqrt(1 - gcl^2)), lty = 3)
      }
      segments(c(0, 0), c(0, 0), c(0, maxsd), c(maxsd, 
                                                0))
      axis.ticks <- pretty(c(0, maxsd))
      axis.ticks <- axis.ticks[axis.ticks <= maxsd]
      axis(1, at = axis.ticks, cex.axis = cex.axis)
      axis(2, at = axis.ticks, cex.axis = cex.axis)
      if (sd.arcs[1]) {
        if (length(sd.arcs) == 1) 
          sd.arcs <- axis.ticks
        for (sdarc in sd.arcs) {
          xcurve <- cos(seq(0, pi/2, by = 0.03)) * sdarc
          ycurve <- sin(seq(0, pi/2, by = 0.03)) * sdarc
          lines(xcurve, ycurve, col = "blue", lty = 3)
        }
      }
      if (show.gamma[1]) {
        if (length(show.gamma) > 1) 
          gamma <- show.gamma
        else gamma <- pretty(c(0, maxsd), n = ngamma)[-1]
        if (gamma[length(gamma)] > maxsd) 
          gamma <- gamma[-length(gamma)]
        labelpos <- seq(45, 70, length.out = length(gamma))
        for (gindex in 1:length(gamma)) {
          xcurve <- cos(seq(0, pi, by = 0.03)) * gamma[gindex] + 
            sd.r
          endcurve <- which(xcurve < 0)
          endcurve <- ifelse(length(endcurve), min(endcurve) - 
                               1, 105)
          ycurve <- sin(seq(0, pi, by = 0.03)) * gamma[gindex]
          maxcurve <- xcurve * xcurve + ycurve * ycurve
          startcurve <- which(maxcurve > maxsd * maxsd)
          startcurve <- ifelse(length(startcurve), max(startcurve) + 
                                 1, 0)
          lines(xcurve[startcurve:endcurve], ycurve[startcurve:endcurve], 
                col = gamma.col)
          if (xcurve[labelpos[gindex]] > 0) 
            plotrix::boxed.labels(xcurve[labelpos[gindex]], ycurve[labelpos[gindex]], 
                         gamma[gindex], border = FALSE)
        }
      }
      xcurve <- cos(seq(0, pi/2, by = 0.01)) * maxsd
      ycurve <- sin(seq(0, pi/2, by = 0.01)) * maxsd
      lines(xcurve, ycurve)
      bigtickangles <- acos(seq(0.1, 0.9, by = 0.1))
      medtickangles <- acos(seq(0.05, 0.95, by = 0.1))
      smltickangles <- acos(seq(0.91, 0.99, by = 0.01))
      segments(cos(bigtickangles) * maxsd, sin(bigtickangles) * 
                 maxsd, cos(bigtickangles) * 0.97 * maxsd, sin(bigtickangles) * 
                 0.97 * maxsd)
      par(xpd = TRUE)
      if (ref.sd) {
        xcurve <- cos(seq(0, pi/2, by = 0.01)) * sd.r
        ycurve <- sin(seq(0, pi/2, by = 0.01)) * sd.r
        lines(xcurve, ycurve)
      }
      points(sd.r, 0, cex = pcex, col=refcol, pch = refpch)
      text(cos(c(bigtickangles, acos(c(0.95, 0.99)))) * 
             1.05 * maxsd, sin(c(bigtickangles, acos(c(0.95, 
                                                       0.99)))) * 1.05 * maxsd, c(seq(0.1, 0.9, by = 0.1), 
                                                                                  0.95, 0.99), cex = cex.axis)
      text(maxsd * 0.8, maxsd * 0.8, "Correlation", srt = 315, 
           cex = cex.axis)
      segments(cos(medtickangles) * maxsd, sin(medtickangles) * 
                 maxsd, cos(medtickangles) * 0.98 * maxsd, sin(medtickangles) * 
                 0.98 * maxsd)
      segments(cos(smltickangles) * maxsd, sin(smltickangles) * 
                 maxsd, cos(smltickangles) * 0.99 * maxsd, sin(smltickangles) * 
                 0.99 * maxsd)
    }
    else {
      x <- ref
      y <- model
      R <- cor(x, y, use = "pairwise.complete.obs")
      E <- mean(x, na.rm = TRUE) - mean(y, na.rm = TRUE)
      xprime <- x - mean(x, na.rm = TRUE)
      yprime <- y - mean(y, na.rm = TRUE)
      sumofsquares <- (xprime - yprime)^2
      Eprime <- sqrt(sum(sumofsquares)/length(complete.cases(x)))
      E2 <- E^2 + Eprime^2
      if (add == FALSE) {
        maxray <- 1.5 * max(sd.f, sd.r)
        plot(c(-maxray, maxray), c(0, maxray), type = "n", 
             asp = 1, bty = "n", xaxt = "n", yaxt = "n", 
             xlim = c(-1.1 * maxray, 1.1 * maxray), xlab = xlab, 
             ylab = ylab, main = main, cex = cex.axis)
        discrete <- seq(180, 0, by = -1)
        listepoints <- NULL
        for (i in discrete) {
          listepoints <- cbind(listepoints, maxray * 
                                 cos(i * pi/180), maxray * sin(i * pi/180))
        }
        listepoints <- matrix(listepoints, 2, length(listepoints)/2)
        listepoints <- t(listepoints)
        lines(listepoints[, 1], listepoints[, 2])
        lines(c(-maxray, maxray), c(0, 0))
        lines(c(0, 0), c(0, maxray))
        for (i in grad.corr.lines) {
          lines(c(0, maxray * i), c(0, maxray * sqrt(1 - 
                                                       i^2)), lty = 3)
          lines(c(0, -maxray * i), c(0, maxray * sqrt(1 - 
                                                        i^2)), lty = 3)
        }
        for (i in grad.corr.full) {
          text(1.05 * maxray * i, 1.05 * maxray * sqrt(1 - 
                                                         i^2), i, cex = cex.axis, adj = cos(i)/2)
          text(-1.05 * maxray * i, 1.05 * maxray * sqrt(1 - 
                                                          i^2), -i, cex = cex.axis, adj = 1 - cos(i)/2)
        }
        seq.sd <- seq.int(0, 2 * maxray, by = (maxray/10))[-1]
        for (i in seq.sd) {
          xcircle <- sd.r + (cos(discrete * pi/180) * 
                               i)
          ycircle <- sin(discrete * pi/180) * i
          for (j in 1:length(xcircle)) {
            if ((xcircle[j]^2 + ycircle[j]^2) < (maxray^2)) {
              points(xcircle[j], ycircle[j], col = "darkgreen", 
                     pch = ".")
              if (j == 10) 
                text(xcircle[j], ycircle[j], signif(i, 
                                                    2), cex = cex.axis, col = "darkgreen", 
                     srt = 90)
            }
          }
        }
        seq.sd <- seq.int(0, maxray, length.out = 5)
        for (i in seq.sd) {
          xcircle <- cos(discrete * pi/180) * i
          ycircle <- sin(discrete * pi/180) * i
          if (i) 
            lines(xcircle, ycircle, lty = 3, col = "blue")
          text(min(xcircle), -0.06 * maxray, signif(i, 
                                                    2), cex = cex.axis, col = "blue")
          text(max(xcircle), -0.06 * maxray, signif(i, 
                                                    2), cex = cex.axis, col = "blue")
        }
        text(0, -0.14 * maxray, "Standard Deviation", 
             cex = cex.axis, col = "blue")
        text(0, -0.22 * maxray, "Centered RMS Difference", 
             cex = cex.axis, col = "darkgreen")
        points(sd.r, 0, pch = 22, bg = "darkgreen", cex = pcex)
        text(0, 1.2 * maxray, "Correlation Coefficient", 
             cex = cex.axis)
      }
      S <- (2 * (1 + R))/(sd.f + (1/sd.f))^2
    }
  }
  points(sd.f * R, sd.f * sin(acos(R)), pch = pch, col = col, 
         cex = pcex)
  invisible(oldpar)
}

a <- array(c(1,4,8,5,6,2,2,7,3), dim=c(3,3))
b <- array(c(9,5,3,6,5,8,1,2,7) , dim=c(3,3))
c <- array(c(1,6,8,5,8,5,8,9,3) , dim=c(3,3))
a2 <- as.numeric(a)
b2 <- as.numeric(b)
c2 <- as.numeric(c)
cols <- c("red", "blue")
taylor.diagram2(ref=a2, model=c2, pos.cor=T, col=cols[2], pch=17,add=FALSE,
                pcex=3, refcol="purple", refpch=15)

^{Created on 2023-01-20 by the reprex package (v2.0.1)}