Skip to content

Conditional coloring of area between curves

4 messages · Duncan Murdoch, Roland Pape, William Dunlap

Duncan Murdoch
#
On 06/06/2013 10:41 AM, Roland Pape wrote:
Suppose the series are named "site1" and "site2", and they have common 
times in a variable "times".  Then the following should do what you want:

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")

If the times for the two series are different it's a little harder; 
first you need to give them common times, and that will depend on how 
you decide to evaluate the values between observations. Probably linear 
interpolation (using approx()) is fine, but it's up to you.

Duncan Murdoch
Roland Pape
# 
Hi Duncan,

that's exactly what I was looking for! The series have common times, so that's no problem. Thanks a lot!

Roland
William Dunlap
# 
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