MVKSA | Air quality
Mvksa Home   Daily Temp   World Bank  
Earth Quake   Stem Cell   Cardio  
Grevillea   Dumbbell   Air Quality  
 
Air quality in the Châtelet underground train station in the Paris metro.
library(ks)
library(colorspace)
library(RColorBrewer)

dd <- "~/Documents/research/www/mvstat.net/mvksa/air/" 
png.create <- function(file) { png(file=paste0(dd,file), res=150, height=4*1.175, width=4*1.175, units="in"); par(oma=c(0,0,0,0)+0.1, mgp=c(1.8,0.5,0), mar=c(2.9,2.9,0,0)+0.1, cex.axis=1.2, cex.lab=1.2) }

## colour functions
pt.col <- function(pos=1, alpha=1, ...) {return(paste0(brewer.pal(12, "Paired"), format(as.hexmode(round(alpha*255,0)), width=2))[pos])}
seq.col <- function(n, alpha=1, ...){c("transparent", rev(sequential_hcl(n=n-1, h=290, power=1.1, alpha=alpha, ...)))}

## plot parameters
par(oma=c(0,0,0,0)+0.1, mgp=c(1.8,0.5,0), mar=c(2.9,2.9,0,0)+0.1, cex.axis=1.2, cex.lab=1.2)
xlab <- expression(CO[2]*" g/"*m^3)
ylab <- expression(PM[10]*" (ppm)")
xlim <- c(350, 1050)
ylim <- c(30, 340)

## prepare data
data(air)
air <- air[, c("date", "time", "co2", "pm10")]
air2 <- reshape(air, idvar="date", timevar="time", direction="wide")
air <- as.matrix(na.omit(air2[,c("co2.20:00", "pm10.20:00")]))
## Fig 7.6a
## scatter plot
fhat.air <- kde(x=air)

plot(fhat.air, cont=0, xaxs="i", yaxs="i", xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab)
points(air, cex=1, pch=16, col=pt.col(4,alpha=0.3))

## Fig 7.6b
## classical kernel density estimate
plot(fhat.air, drawlabels=FALSE, display="filled.contour", lwd=1, col.fun=seq.col, cont=seq(20,80,by=20), xaxs="i", yaxs="i", xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab)

## Fig 7.6c
## deconvolution kernel density estimate
Sigma.air <- diag(c(var(air2[,"co2.19:00"] - air2["co2.21:00"], na.rm=TRUE), var(air2[,"pm10.19:00"] - air2[,"pm10.21:00"], na.rm=TRUE)))
fhat.air.dec <- kdcde(x=air, Sigma=Sigma.air, reg=0.00021)   ## kdcde is alias for dckde

plot(fhat.air.dec, drawlabels=FALSE, display="filled.contour", lwd=1, col.fun=seq.col, cont=seq(20,80,by=20), xaxs="i", yaxs="i", xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab)