R: Narysuj wielokąt z kolorem warunkowym

Question

Chcę pokolorować obszar pod krzywą. Obszar z y> 0 powinien być czerwony, obszar z y < 0 powinien być zielony.

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

Korzystanie ifelse() nie działa:

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

Co osiągnąłem do tej pory to:

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

Ale wtedy czerwony obszar jest zbyt duży. Czy masz jakieś pomysły na uzyskanie pożądanego rezultatu?

Answer 1

Jeśli chcesz mieć dwa różne kolory, potrzebujesz dwóch różnych wielokątów. Możesz wielokrotnie wywoływać wielokąty lub dodawać wartości NA do wektorów x i y, aby wskazać nowy wielokąt. R nie będzie automatycznie obliczać dla ciebie przecięcia. Musisz to zrobić sam. Oto, jak można to narysować w różnych kolorach.

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) 

W bardziej ogólnym przypadku, to może nie być tak łatwo podzielić wielokąt w różnych regionach. Wydaje się, że istnieje pewne wsparcie dla tego rodzaju operacji w pakietach GIS, gdzie tego typu rzeczy są bardziej powszechne. Złożyłem jednak nieco ogólny przypadek, który może zadziałać w przypadku prostych wielokątów.

Najpierw definiuję zamknięcie, które definiuje linię cięcia. Funkcja przyjmie nachylenie i punkt przecięcia z osią y dla linii i zwróci funkcje potrzebne do wycięcia wielokąta.

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

Teraz zdefiniujemy funkcję podziału wielokąta za pomocą rozdzielacza, który właśnie zdefiniowaliśmy.

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 
} 

Wejście jest tylko wartości x i y jak można przejść do polygon wraz z kutra. Zwróci obiekt data.frame z wartościami x i y, które mogą być używane z polygon.

Na przykład

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

To powinno działać dla prostych wielokątów wklęsłych, jak również ze skośnymi liniami. Oto kolejny przykład

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

Teraz wszystko może się nieco niechlujny gdy wierzchołek wielokąta spada bezpośrednio na linii rozłupywania. Mogę spróbować to poprawić w przyszłości.