CodHump lexis plot functions

plottingFunctions.R

# Author: tim
###############################################################################
#get packages sorted out
if( !"colorspace" %in% rownames(installed.packages())){
	install.packages("colorspace")
}
# this is installed for certain grayscale tools in it.
# too bad we end up with this huge package just for that,
# but really the functions are convenient.
if( !"spatstat" %in% rownames(installed.packages())){
	install.packages("spatstat")
}

library(colorspace)
library(spatstat)

# darken by changing luminance
darkenhex <- function(hexcols,fac=.3){
	# from @Roland:
	#https://stackoverflow.com/questions/30219738/is-there-a-way-to-programmatically-darken-the-color-given-rgb-values
	require(colorspace)
	zz   <- textConnection(hexcols)
	cols <- readhex(file = zz,	class = "RGB")
	close(zz)
# transform to hue/lightness/saturation colorspace
	cols <- as(cols, "HLS")
	cols@coords[, "L"] <- pmax(0, cols@coords[, "L"] - fac)
	cols <- as(cols, "RGB")
	hex(cols)
}
# a gives a number from 0 to 1. 0 is black, 1 is white. vectorized
howdark <- function(col){
	colorspace::hex2RGB(spatstat::to.grey(col))@coords[, 1]
}
# darkens a color to a specific grayscale target darkness. not vectorized
darken.to <- function(hexcol,target=.2){
	require(colortools)
	require(spatstat)

	# a function to optimize the ideal darkening factor
	mintarg <- function(fac,col,target){
		hexcol <- darkenhex(hexcol,fac)
		abs(target - howdark(hexcol))
	}
    # get darkening factor
	fac    <- optimize(mintarg,c(-.8,.8),col=hexcol,target=target)$minimum
	# now darken by that much...
	colout <- darkenhex(hexcol,fac)
	colout
}

# will use this for surface ramps.
my.palette <- function(col,n,from=.2,to=.9){
	fromcol <- darken.to(col,from)
	tocol   <- darken.to(col,to)
	colorRampPalette(c(fromcol,tocol),space="Lab")(n)
}

# this function used to draw both brute and hump surfaces.
draw.surf <- function(Mbrj, Msmj, col, axes = c(1, 2)){
	
	breaks <- pretty(Msmj, n = 10)
	
	colsj  <- rev(my.palette(col,n=length(breaks)-1,from=.3,to=.95))
	
	ages   <- as.integer(rownames(Mbrj))
	years  <- as.integer(colnames(Mbrj))
	
	# test, preserve 0s in plot
	Mbrj[Mbrj < 1e-6] <- NA
	image(
			xlim = range(years) + c(0,1),
			ylim = range(ages) + c(0,1),
			axes = FALSE, 
			asp = 1,
			x = years + .5, 
			y = ages + .5, 
			z = t(Mbrj), 
			xlab = "", 
			ylab = "",
			main = "",
			col = colsj,
			breaks = breaks)
	LexRefN(ages,years,col="#1A1A1A30",N=10)
	
	# label every second contour.
	breaks0 <- breaks[1:length(breaks) %% 2 == 0]
	breaks1 <- breaks[1:length(breaks) %% 2 == 1]
	contour(years + .5, 
			ages + .5, 
			Msmj, 
			levels = breaks1, 
			add = TRUE,
			lwd = 0.7, 
			col = gray(0),
			labcex=.8)
	contour(years + .5, 
			ages + .5, 
			Msmj, 
			levels = breaks0, 
			add = TRUE,
			lwd = 0.5, 
			col = gray(.2),
			drawlabels = FALSE)
	
	if (length(axes)>0){
		if (any(axes == 1)){
			# year bottom
			axis(1,at=years[years%%20==0],cex.axis=1.5,pos=min(ages),lwd=0,lwd.ticks=1)
		} else {
			axis(1,at=years[years%%20==0],pos=min(ages), labels = FALSE,lwd=0,lwd.ticks=1)
		}
		if (any(axes == 3)){
			# year top
			axis(3,at=years[years%%20==0],cex.axis=1.5,pos=max(ages)+1,lwd=0,lwd.ticks=1)
		} else {
			axis(3,at=years[years%%20==0],pos=max(ages)+1, labels = FALSE,lwd=0,lwd.ticks=1)
		}
		if (any(axes == 2)){
			# age bottom
			axis(2,at=ages[ages%%20==0],cex.axis=1.5,pos=min(years),las=2,lwd=0,lwd.ticks=1)
		} else {
			axis(2,at=ages[ages%%20==0],pos=min(years), labels = FALSE,lwd=0,lwd.ticks=1)
		}
		if (any(axes == 4)){
			# age top
			axis(4,at=ages[ages%%20==0],cex.axis=1.5,pos=max(years)+1,las=2,lwd=0,lwd.ticks=1)
		} else {
			axis(4,at=ages[ages%%20==0],pos = max(years)+1, labels = FALSE,lwd=0,lwd.ticks=1)
		}
	}
	box()
}


LexRefN <- function(ages, years, col = "#DDDDDD50", N=5,...){
	# vertical
	segments(years[years %% N == 0],
			min(ages),
			years[years %% N == 0],
			max(ages) + 1,
			col = col,
			...
	)
	# horiz lexis lines
	segments(min(years),
			ages[ages %% N == 0],
			max(years) + 1,
			ages[ages %% N == 0],
			col = col,
			...
	)
	# diag cohort references, bit more tricky:
	l.years <- (min(years) - ages) 
	coh.ext <- range(c(l.years[l.years %% N == 0], years[years %% N == 0]))
	cohs  <- seq(coh.ext[1], coh.ext[2], by = N)
	
	# bottom, left:
	xl  <- cohs + min(ages)
	yb  <- rep(min(ages), length(cohs))
	yb[xl < min(years)] <- yb[xl < min(years)] + min(years) - xl[xl < min(years)]
	xl[xl < min(years)] <- min(years)
	
	# top, right:
	xr  <- cohs + max(ages) + 1
	yt  <- rep(max(ages) + 1, length(cohs))
	yt[xr > (max(years) + 1)] <- yt[xr > (max(years) + 1)] - xr[xr > (max(years) + 1)] + (max(years) + 1)
	xr[xr > (max(years) + 1)] <- max(years) + 1
	
	# cut down one last time:
	xr <- xr[yt >= min(ages)]
	xl <- xl[yt >= min(ages)]
	yb <- yb[yt >= min(ages)]
	yt <- yt[yt >= min(ages)]
	
	# draw cohort refs:
	segments(xl, yb, xr, yt, col = col, ...)
}