Skip to content
Prev 324942 / 398503 Next

Conditional coloring of area between curves

That code has some problems if your time series do not cross at the time
points.  E.g., Duncan's code with some diagnostic plotting is:

  f0 <- function (site1, site2, x = time(site1), check = FALSE) 
  {
      stopifnot(length(x) == length(site1), length(x) == length(site2),  all(diff(x) > 0))
      if (check) {
          oldMfrow <- par(mfrow = c(2, 1))
          on.exit(par(oldMfrow))
          matplot(x, cbind(site1, site2), type = "b")
      }
      smaller <- pmin(site1, site2)
      plot(x, site1, ylim = range(c(site1, site2)), type = "n")
      polygon(c(x, rev(x)), c(smaller, rev(site1)), col = "red")
      polygon(c(x, rev(x)), c(smaller, rev(site2)), col = "blue")
  }
  # See how it messes the crossings in
  f0(c(0, 1, -1, 0), c(1, 0, -2, 1), check=TRUE)

The following is a quick and dirty way of adding all the crossing
points to the original data sets so the plot is right.
  addCrossingPoints <- function (y1, y2, x) 
  {
      stopifnot(length(y1)==length(x), length(y2)==length(x), all(diff(x)>0))
      # isCrossingSegment[i]==TRUE means that y1 and y2 cross between
      #    x[i] and x[i+1]  (strictly, not at x[i] or x[i+1]).
      i <- seq_len(length(y1) - 1)
      isCrossingSegment <- ((y1[i] < y2[i]) & (y1[i+1] > y2[i+1])) |
                           ((y1[i] > y2[i]) & (y1[i+1] < y2[i+1]))
      if (any(isCrossingSegment)) {
          i <- which(isCrossingSegment)
          dx <- x[i+1] - x[i]
          m1 <- (y1[i+1] - y1[i]) / dx
          m2 <- (y2[i+1] - y2[i]) / dx
          newDX <- (y1[i] - y2[i]) / (m2 - m1)
          newY <- y1[i] + newDX * m1 
          o <- order( x <- c(x, x[i] + newDX))
          x <- x[o]
          y1 <- c(y1, newY)[o]
          y2 <- c(y2, newY)[o]
      }
      list(y1=y1, y2=y2, x=x)
  }
  # Combine it with Duncan's code to get a better looking plot
  f1 <- function(site1, site2, x=time(site1), check = FALSE) {
      tmp <- addCrossingPoints(site1, site2, x)
      f0(tmp$y1, tmp$y2, tmp$x, check=check)
  }
  f1(c(0, 1, -1, 0), c(1, 0, -2, 1), check=TRUE)


Bill Dunlap
Spotfire, TIBCO Software
wdunlap tibco.com