R/plot.aep.R

In chgrl/bReeze: Functions for Wind Resource Assessment

plot.aep <-
function(x, show.total=TRUE, ...) {
###	plotting annual energy production rose
	
	aep <- x$aep
	dim.data <- dim(aep)
	num.sectors <- dim.data[1] - 1
	num.classes <- dim.data[2] - 3
	total <- aep$total[num.sectors+1]
	unit <- attr(aep$total, "unit")
	
	if(num.classes>1) {
		aep.cum <- aep[1:num.sectors,4:dim.data[2]]
		for(i in 2:num.classes) {
			aep.cum[,i] <- aep.cum[,i] + aep.cum[,i-1]
		}
	} else {
		aep.cum <- data.frame(aep[1:num.sectors,3])
		num.classes <- 1
	}
	
	# prepare plot
	sectors <- seq(0, 360-360/num.sectors, by=360/num.sectors)
	sectors <- sectors+90
	sector.edges <- sectors*pi/180
	sector.width <- sector.edges[2] - sector.edges[1]
	
	plot.param <- list(...)
	if(any(names(plot.param)=="col")) colset <- plot.param$col
	else {
		if(num.classes==1) colset <- c("#3182BD")
		else if(num.classes==2) colset <- c("#3182BD", "#9ECAE1")
		else if(num.classes>2 && num.classes<=11) {
	    	if(requireNamespace("RColorBrewer", quietly=TRUE)) colset <- rev(RColorBrewer::brewer.pal(num.classes, "Blues"))
	    	else colset <- rev(rainbow(num.classes, start=0.0, end=0.7))
	    } else {
	    	colset <- rev(rainbow(num.classes, start=0.0, end=0.7))
	    }
	}
	if(any(names(plot.param)=="cex")) cex <- plot.param$cex
	else cex <- 1
	if(any(names(plot.param)=="cex.axis")) cex.axis <- plot.param$cex.axis
	else cex.axis <- cex-0.2
	if(any(names(plot.param)=="cex.lab")) cex.lab <- plot.param$cex.lab
	else cex.lab <- cex
	if(any(names(plot.param)=="cex.leg")) cex.leg <- plot.param$cex.leg
	else cex.leg <- cex-0.2
	if(any(names(plot.param)=="width.leg")) width.leg <- plot.param$width.leg
	else width.leg <- 0.2
	if(any(names(plot.param)=="x.intersp")) x.intersp <- plot.param$x.intersp
	else x.intersp <- 0.4
	if(any(names(plot.param)=="y.intersp")) y.intersp <- plot.param$y.intersp
	else y.intersp <- 0.8
	if(any(names(plot.param)=="circles")) circles <- seq(plot.param$circles[1], plot.param$circles[2], by=plot.param$circles[3])
	else circles <- NULL
	if(any(names(plot.param)=="fg")) fg <- plot.param$fg
	else fg <- FALSE
	if(any(names(plot.param)=="pos.axis")) pos.axis <- plot.param$pos.axis
	else pos.axis <- 60
	if(any(names(plot.param)=="sec.space")) sec.space <- 1-plot.param$sec.space
	else sec.space <- 0.8
	if(any(names(plot.param)=="col.circle")) col.circle <- plot.param$col.circle
	else col.circle <- "gray45"
	if(any(names(plot.param)=="col.cross")) col.cross <- plot.param$col.cross
	else col.cross <- "gray45"
	if(any(names(plot.param)=="col.axis")) col.axis <- plot.param$col.axis
	else col.axis <- "gray45"
	if(any(names(plot.param)=="col.lab")) col.lab <- plot.param$col.lab
	else col.lab <- "black"
	if(any(names(plot.param)=="col.leg")) col.leg <- plot.param$col.leg
	else col.leg <- "black"
	if(any(names(plot.param)=="lwd.circle")) lwd.circle <- plot.param$lwd.circle
	else lwd.circle <- 0.7
	if(any(names(plot.param)=="lwd.cross")) lwd.cross <- plot.param$lwd.cross
	else lwd.cross <- 0.7
	if(any(names(plot.param)=="lty.circle")) lty.circle <- plot.param$lty.circle
	else lty.circle <- 2
	if(any(names(plot.param)=="lty.cross")) lty.cross <- plot.param$lty.cross
	else lty.cross <- 1
	if(any(names(plot.param)=="title.leg")) title.leg <- plot.param$title.leg
	else title.leg <- "Wind speed\n[m/s]"
	if(any(names(plot.param)=="border.leg")) border.leg <- plot.param$border.leg
	else border.leg <- colset[1:num.classes]
	if(any(names(plot.param)=="bty.leg")) bty.leg <- plot.param$bty.leg
	else bty.leg <- "n"
	if(any(names(plot.param)=="col.border")) col.border <- plot.param$col.border
	else col.border <- NULL
	if(any(names(plot.param)=="lwd.border")) lwd.border <- plot.param$lwd.border
	else lwd.border <- 0.5
	if(any(names(plot.param)=="type")) {
		type <- plot.param$type
		if(type=="n") type <- FALSE
		else type <- TRUE
	} else type <- TRUE
	
	if(type) {  # make a plot
		if(is.null(circles)) {
			aep.max <- max(aep.cum, na.rm=TRUE)
			mag <-length(strsplit(as.character(aep.max),"")[[1]])-1
			circ.max <- ceiling(aep.max/10^mag)*10^mag
			if(circ.max<=2*10^mag) circles <- tail(seq(0, 2*10^mag, length.out=5), -1)
			if(circ.max>2*10^mag && circ.max<=3*10^mag) circles <- tail(seq(0, 3*10^mag, length.out=4), -1)
			if(circ.max>3*10^mag && circ.max<=4*10^mag) circles <- tail(seq(0, 4*10^mag, length.out=5), -1)
			if(circ.max>4*10^mag && circ.max<=5*10^mag) circles <- tail(seq(0, 5*10^mag, length.out=5), -1)
			if(circ.max>5*10^mag && circ.max<=6*10^mag) circles <- tail(seq(0, 6*10^mag, length.out=5), -1)
			if(circ.max>6*10^mag && circ.max<=8*10^mag) circles <- tail(seq(0, 8*10^mag, length.out=5), -1)
			if(circ.max>8*10^mag) circles <- tail(seq(0, 10^(mag+1), length.out=5), -1)
		}
	
		if(is.numeric(width.leg)) if(width.leg>1) stop("'width.leg' must be a numeric value between 0 and 1")
		if(num.classes>1 && width.leg!=0) lo <- layout(matrix(1:2, 1, 2), widths=c(1, width.leg))
		
		# plot
		par(mar=c(1,1,1,1), las=1)
		plot.new()
		pin <- par("pin")
		xlim <- ylim <- c(-1, 1)
		if (pin[1] > pin[2]) {
			xlim <- (pin[1]/pin[2]) * xlim
		} else {
			ylim <- (pin[2]/pin[1]) * ylim
		}
		plot.window(xlim, ylim, "", asp=1)
			
		if(!fg) {
			for(c in num.classes:1) {
				plot.data <- c(tail(rev(as.vector(aep.cum[,c])), n=1), head(rev(as.vector(aep.cum[,c])), n=-1))
				
				for (i in 1:num.sectors) {
					arc.pts <- seq(sector.edges[i] - sector.width/2*sec.space, sector.edges[i] + sector.width/2*sec.space, length.out=trunc(360/num.sectors*sec.space))
					rad <- 0.9 * plot.data[i]/tail(circles, 1)
					xlist <- c(0, rad * cos(arc.pts), 0)
					ylist <- c(0, rad * sin(arc.pts), 0)
					polygon(xlist, ylist, col=colset[c], border=colset[c], lwd=0.01)
				}
			}
			if(!is.null(col.border)) {
				aep.cum[is.na(aep.cum)] <- 0
				plot.max <- apply(aep.cum, 1, max)
				plot.max <- c(tail(rev(plot.max), n=1), head(rev(plot.max), n=-1))
				for(i in 1:num.sectors) {
					arc.pts <- seq(sector.edges[i] - sector.width/2*sec.space, sector.edges[i] + sector.width/2*sec.space, length.out=trunc(360/num.sectors*sec.space))
					rad <- 0.9 * plot.max[i]/tail(circles, 1)
					xlist <- c(0, rad * cos(arc.pts), 0)
					ylist <- c(0, rad * sin(arc.pts), 0)
					lines(xlist, ylist, col=col.border, lwd=lwd.border)
				}
			}
		}
		
		circle.pts <- seq(0, 2*pi, length.out=360)
		pos.axis <- pi/2 - pos.axis*pi/180
		
		for(i in 1:length(circles)) {
			rad <- 0.9 * circles[i]/tail(circles, 1)
			circle.x <- cos(circle.pts)*rad
			circle.y <- sin(circle.pts)*rad
			lines(circle.x, circle.y, lty=lty.circle, lwd=lwd.circle, col=col.circle)
			text(cos(pos.axis)*rad, sin(pos.axis)*rad, circles[i], cex=cex.axis, col=col.axis)
		}
		     
		lines(c(-0.92, 0.92), c(0, 0), lty=lty.cross, lwd=lwd.cross, col=col.cross)
		lines(c(0, 0), c(0.92, -0.92), lty=lty.cross, lwd=lwd.cross, col=col.cross)
		text(0, -0.90, "S", pos=1, cex=cex.lab, col=col.lab)
		text(-0.90, 0, "W", pos=2, cex=cex.lab, col=col.lab)
		text(0, 0.90, "N", pos=3, cex=cex.lab, col=col.lab)
		text(0.90, 0, "E", pos=4, cex=cex.lab, col=col.lab)
		if(show.total) text(-1, 1, paste("Total:", total), pos=4, cex=cex.axis, col=col.axis)
		if(unit=="MWh/a") text(1, -1, "MWh/a", pos=2, cex=cex.axis, col=col.axis)
		else text(1, -1, unit, pos=2, cex=cex.axis, col=col.axis)
		
		if(fg) {
			for(c in num.classes:1) {
				plot.data <- c(tail(rev(as.vector(aep.cum[,c])), n=1), head(rev(as.vector(aep.cum[,c])), n=-1))
				
				for (i in 1:num.sectors) {
					arc.pts <- seq(sector.edges[i] - sector.width/2*sec.space, sector.edges[i] + sector.width/2*sec.space, length.out=trunc(360/num.sectors*sec.space))
					rad <- 0.9 * plot.data[i]/tail(circles, 1)
					xlist <- c(0, rad * cos(arc.pts), 0)
					ylist <- c(0, rad * sin(arc.pts), 0)
					polygon(xlist, ylist, col=colset[c], border=colset[c], lwd=0.01)
				}
			}
			if(!is.null(col.border)) {
				aep.cum[is.na(aep.cum)] <- 0
				plot.max <- apply(aep.cum, 1, max)
				plot.max <- c(tail(rev(plot.max), n=1), head(rev(plot.max), n=-1))
				for(i in 1:num.sectors) {
					arc.pts <- seq(sector.edges[i] - sector.width/2*sec.space, sector.edges[i] + sector.width/2*sec.space, length.out=trunc(360/num.sectors*sec.space))
					rad <- 0.9 * plot.max[i]/tail(circles, 1)
					xlist <- c(0, rad * cos(arc.pts), 0)
					ylist <- c(0, rad * sin(arc.pts), 0)
					lines(xlist, ylist, col=col.border, lwd=lwd.border)
				}
			}
		}
		
		if(num.classes>1 && width.leg!=0) {
			par(mar=c(0,0,0,0))
			plot(0, type="n", axes=FALSE, xlab="", ylab="")
			legend("left", legend=names(aep.cum), title=title.leg, fill=colset[1:num.classes], xjust=0, bty=bty.leg, border=border.leg, cex=cex.leg, x.intersp=x.intersp, y.intersp=y.intersp, text.col=col.leg)
		}
	}
	
	plt <- list(aep.cum=aep.cum, total=total, unit=unit, attr=attr(x, "call"))
	class(plt) <- "aep.plot"
	invisible(plt)
}