pop pyramid over time figure for first paadataviz flyer

PAAdatavizflyergraphic1.R

# Author: tim
###############################################################################

library(HMDHFDplus) # install from CRAN

# set 'us' and 'pw' as strings in the console equal to
# your HMD username and password
Pop <- readHMDweb("SWE", "Population", username = us, password = pw)

library(RColorBrewer)

range(Pop$Year)
ramp      <- colorRampPalette(brewer.pal(9,"YlGnBu"),space="Lab")

# single years
years.all <- unique(Pop$Year)
# years ending in 0 since 1870
years     <- years.all[years.all %% 10 == 0 & years.all >= 1870]

# color by mean age:
meana <- c()
for (i in 1:length(years)){
	x <- Pop$Male1[Pop$Year == years[i]]
	meana[i] <- sum(x*.5:110.5) / sum(x)
}
ints <- unique(floor(meana))

cols <- ramp(length(ints)+6)


pdf("/home/tim/workspace/Other/PAAdatavizgraphic.pdf",height=3,width=9)
par(mai=c(.1,0,0,.1), xaxs="i",yaxs="i",xpd=TRUE)
plot(NULL,xlim=c(-.028,.155),ylim=c(0,110),axes=FALSE,xlab="",ylab="")
x.at <- 0
for (i in 1:length(years)){
	x <- Pop$Male1[Pop$Year == years[i]]
	
	# this ugly bit is to cut off the top outline of the pyramid if
	# the highest ages have 0 pop. Ugly but swift to implement. Not
	# worth trying to understand what's happening in this little chunk
	max.x <- 111 - max(cumsum(diff(cumsum(cumsum(rev(x))==0)))) 
	x     <- x[1:max.x]
	
	# scale to proportion!
	x     <- x / sum(x)
	color <- cols[which(ints == floor(meana[i]))+2]

	# the pyramid slice (just males here)
	polygon(c(0, rep(-x, each = 2),0) + x.at, rep(0:max.x, each = 2), col = color, 
			border = "white", lwd = .5)
	x.at <- x.at + .01 # shift right
}
# legend
rect(x.at - .007, seq(0, 100, by = 10) * .7 + .1, x.at, seq(10, 110, by = 10) * .7 + .1, col = colors, border = "white")
text(x.at, seq(0, 100, by = 10) * .7 + .1, ints, pos = 4, cex = .8)
text(x.at - .003, 110 * .7 + .1, bquote(bar(x)), pos = 3)

dev.off()