R：用条件颜色绘制多边形

Question

我想给曲线下的区域着色。 y > 0 的区域应该是红色的，y

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

使用ifelse() 不起作用：

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

到目前为止，我取得的成果如下：

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

但是红色区域太大了。你有什么想法如何得到想要的结果吗？

Answer 1

一个更快但不是很准确的解决方案是根据分组变量（例如，上面=红色和下面=蓝色）将数据框拆分为列表。对于相当大的（我会说> 100个元素）数据集，这是一个非常好的解决方法。对于较小的块，可能会出现一些不连续性：

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

正如我所说，对于较小的数据集，如果正面和负面区域之间发生急剧变化，则可能会出现一些不连续性，但如果水平线区分这两者，或者绘制更多元素，则可以忽略这种影响：

Answer 2

如果你想要两种不同的颜色，你需要两个不同的多边形。您可以多次调用多边形，也可以在x 和y 向量中添加NA 值以指示新多边形。 R 不会自动为您计算交点。你必须自己做。以下是您可以使用不同颜色绘制的方法。

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)

在更一般的情况下，将多边形分割成不同的区域可能并不容易。 GIS 包中似乎对这种类型的操作有一些支持，其中这种类型的东西更常见。不过，我整理了一个可能适用于简单多边形的一般情况。

首先，我定义了一个将定义切割线的闭包。该函数将获取直线的斜率和 y 轴截距，并返回切割多边形所需的函数。

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

现在我们将定义一个函数，使用我们刚刚定义的分割器实际分割多边形。

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
}

输入只是x 和y 的值，因为您将与切割器一起传递给polygon。它将返回一个带有 x 和 y 值的 data.frame，可用于polygon。

例如

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

这应该适用于简单的凹多边形以及斜线。这是另一个例子

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

现在，当多边形的顶点直接落在分割线上时，事情会变得有些混乱。以后我可能会尝试纠正它。