Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

R: Draw a polygon with conditional colour

Tags:

plot

r

polygon

I want to colour the area under a curve. The area with y > 0 should be red, the area with y < 0 should be green.

x <- c(1:4)
y <- c(0,1,-1,2,rep(0,4))
plot(y[1:4],type="l")
abline(h=0)

Using ifelse() does not work:

polygon(c(x,rev(x)),y,col=ifelse(y>0,"red","green"))

What I achieved so far is the following:

polygon(c(x,rev(x)),y,col="green")
polygon(c(x,rev(x)),ifelse(y>0,y,0),col="red")

But then the red area is too large. Do you have any ideas how to get the desired result?

like image 661
mrub Avatar asked Jul 02 '14 15:07

mrub


2 Answers

If you want two different colors, you need two different polygons. You can either call polygon multiple times, or you can add NA values in your x and y vectors to indicate a new polygon. R will not automatically calculate the intersection for you. You must do that yourself. Here's how you could draw that with different colors.

x <- c(1,2,2.5,NA,2.5,3,4)
y <- c(0,1,0,NA,0,-1,0)

#calculate color based on most extreme y value
g <- cumsum(is.na(x))
gc <- ifelse(tapply(y, g, 
    function(x) x[which.max(abs(x))])>0, 
    "red","green")

plot(c(1, 4),c(-1,1), type = "n")
polygon(x, y, col = gc)
abline(h=0)

enter image description here

In the more general case, it might not be as easy to split a polygon into different regions. There seems to be some support for this type of operation in GIS packages, where this type of thing is more common. However, I've put together a somewhat general case that may work for simple polygons.

First, I define a closure that will define a cutting line. The function will take a slope and y-intercept for a line and will return the functions we need to cut a polygon.

getSplitLine <- function(m=1, b=0) {
    force(m); force(b)
    classify <- function(x,y) {
        y >= m*x + b
    }
    intercepts <- function(x,y, class=classify(x,y)) {
        w <- which(diff(class)!=0)
        m2 <- (y[w+1]-y[w])/(x[w+1]-x[w])
        b2 <- y[w] - m2*x[w]

        ix <- (b2-b)/(m-m2)
        iy <- ix*m + b
        data.frame(x=ix,y=iy,idx=w+.5, dir=((rank(ix, ties="first")+1) %/% 2) %% 2 +1)
    }
    plot <- function(...) {
    abline(b,m,...)
    }
    list(
        intercepts=intercepts,
        classify=classify,
        plot=plot
    )
}

Now we will define a function to actually split a polygon using the splitter we've just defined.

splitPolygon <- function(x, y, splitter) {
    addnullrow <- function(x) if (!all(is.na(x[nrow(x),]))) rbind(x, NA) else x
    rollup <- function(x,i=1) rbind(x[(i+1):nrow(x),], x[1:i,])
    idx <- cumsum(is.na(x) | is.na(y))
    polys <- split(data.frame(x=x,y=y)[!is.na(x),], idx[!is.na(x)])
    r <- lapply(polys, function(P) {
        x <- P$x; y<-P$y
        side <- splitter$classify(x, y)
        if(side[1] != side[length(side)]) {
            ints <- splitter$intercepts(c(x,x[1]), c(y, y[1]), c(side, side[1]))
        } else {
            ints <- splitter$intercepts(x, y, side)
        }
        sideps <- lapply(unique(side), function(ss) {
            pts <- data.frame(x=x[side==ss], y=y[side==ss], 
                idx=seq_along(x)[side==ss], dir=0)
            mm <- rbind(pts, ints)
            mm <- mm[order(mm$idx), ]
            br <- cumsum(mm$dir!=0 & c(0,head(mm$dir,-1))!=0 & 
                c(0,diff(mm$idx))>1)
            if (length(unique(br))>1) {
                mm<-rollup(mm, sum(br==br[1]))
            }
            br <- cumsum(c(FALSE,abs(diff(mm$dir*mm$dir))==3))
            do.call(rbind, lapply(split(mm, br), addnullrow))
        })
        pss<-rep(unique(side), sapply(sideps, nrow))
        ps<-do.call(rbind, lapply(sideps, addnullrow))[,c("x","y")]
        attr(ps, "side")<-pss
        ps
   })
   pss<-unname(unlist(lapply(r, attr, "side")))
   src <- rep(seq_along(r), sapply(r, nrow))
   r <- do.call(rbind, r)
   attr(r, "source")<-src
   attr(r, "side")<-pss
   r
}

The input is just the values of x and y as you would pass to polygon along with the cutter. It will return a data.frame with x and y values that can be used with polygon.

For example

x <- c(1,2,2.5,NA,2.5,3,4)
y <- c(1,-2,2,NA,-1,2,-2)

sl<-getSplitLine(0,0)

plot(range(x, na.rm=T),range(y, na.rm=T), type = "n")
p <- splitPolygon(x,y,sl)
g <- cumsum(c(F, is.na(head(p$y,-1))))
gc <- ifelse(attr(p,"side")[is.na(p$y)],  
    "red","green")
polygon(p, col=gc)
sl$plot(lty=2, col="grey")

enter image description here

This should work for simple concave polygons as well with sloped lines. Here's another example

x <- c(1,2,3,4,5,4,3,2)
y <- c(-2,2,1,2,-2,.5,-.5,.5)

sl<-getSplitLine(.5,-1.25)

plot(range(x, na.rm=T),range(y, na.rm=T), type = "n")
p <- splitPolygon(x,y,sl)
g <- cumsum(c(F, is.na(head(p$y,-1))))
gc <- ifelse(attr(p,"side")[is.na(p$y)],  
    "red","green")
polygon(p, col=gc)
sl$plot(lty=2, col="grey")

enter image description here

Right now things can get a bit messy when the the vertex of the polygon falls directly on the splitting line. I may try to correct that in the future.

like image 161
MrFlick Avatar answered Oct 30 '22 18:10

MrFlick


A faster, but not very accurate solution is to split data frame to list according to grouping variable (e.g. above=red and below=blue). This is a pretty nice workaround for rather big (I would say > 100 elements) datasets. For smaller chunks some discontinuity may be visible:

x <- 1:100
y1 <- sin(1:100/10)*0.8
y2 <- sin(1:100/10)*1.2
plot(x, y2, type='l')
lines(x, y1, col='red')

df <- data.frame(x=x, y1=y1, y2=y2)

df$pos_neg <- ifelse(df$y2-df$y1>0,1,-1) # above (1) or below (-1) average

# create the number for chunks to be split into lists:
df$chunk <- c(1,cumsum(abs(diff(df$pos_neg)))/2+1) # first element needs to be added`
df$colors <- ifelse(df$pos_neg>0, "red","blue") # colors to be used for filling the polygons
# create lists to be plotted:
l <- split(df, df$chunk) # we should get 4 sub-lists
lapply(l, function(x) polygon(c(x$x,rev(x$x)),c(x$y2,rev(x$y1)),col=x$colors))

As I said, for smaller dataset some discontinuity may be visible if sharp changes occur between positive and negative areas, but if horizontal line distinguishes between those two, or more elements are plotted then this effect is neglected:

results

like image 35
speleo Avatar answered Oct 30 '22 18:10

speleo