library(ks)
library(colorspace)
library(RColorBrewer)
library(maps)

## colour function
pt.col <- function(pos=1, alpha=1, ...) {return(paste0(brewer.pal(12, "Paired"), format(as.hexmode(round(alpha*255,0)), width=2))[pos])}

## Pacific map
pacific.map <- function(xlim, ylim, col=grey(0, 1), mar=c(0,0,0,0), interior=FALSE, lty=3, ...)
{
    if (missing(xlim)) xlim <- c(85, 305)
    if (missing(ylim)) ylim <- c(-70, 70)
    map("world2", xlim=xlim, ylim=ylim, col=col, mar=mar, interior=interior, lty=lty, ...)
    box()
}

## plot parameters
par(oma=c(0,0,0,0), mgp=c(0,0,0), mar=c(0,0,0,0), cex.axis=1.2, cex.lab=1.2)
xlim <- c(85, 305)
ylim <- c(-70, 70)
		
## prepare data
data(quake)
## Restrict to approximate Pacific Ring of Fire
quake <- quake[quake$prof==1,]
## translate to [-180,180] -> [0,360] longitude
quake$long[quake$long<0] <- quake$long[quake$long<0] + 360
quake <- quake[, c("long", "lat")]

## tectonic plate boundaries
data(plate)
plate <- plate[plate$long < -20 | plate$long > 20,]
plate$long[plate$long<0 & !is.na(plate$long)] <- plate$long[plate$long<0 & !is.na(plate$long)] + 360

## Fig 1.3
## scatter plot
pacific.map()
lines(plate[,1:2], col=pt.col(2, alpha=0.5), lwd=2)
points(quake, cex=0.7, pch=16, col=pt.col(4,alpha=0.3))

## Fig 6.9
## first principal components + rug plot
pca.quake <- princomp(quake)
U <- eigen(var(quake))$vectors
lambda <- eigen(var(quake))$values
pca.quake.scores <- as.matrix(quake) %*% U
pca.quake.scores <- sweep(pca.quake.scores, 2, apply(pca.quake.scores, 2, min))
pca.quake.scores <- sweep(pca.quake.scores, 2, apply(pca.quake.scores, 2, max), FUN="/")

pacific.map()
lines(plate[,1:2], col=pt.col(2, alpha=0.5), lwd=2)
points(quake, cex=0.7, pch=16, col=pt.col(4,alpha=0.3))
points(pca.quake$center[1], pca.quake$center[2], col=1, pch=16)
arrows(x0=pca.quake$center[1], y0=pca.quake$center[2], x1=pca.quake$center[1]-0.015*lambda[1]*U[1,1], y1=pca.quake$center[2]-0.015*lambda[1]*U[2,1], length=0.1, lwd=2)
arrows(x0=pca.quake$center[1], y0=pca.quake$center[2], x1=pca.quake$center[1]-0.015*lambda[2]*U[1,2], y1=pca.quake$center[2]-0.015*lambda[2]*U[2,2], length=0.1, lwd=2)
rug(pca.quake.scores[,1]*0.9*diff(range(quake$long)) + min(quake$long), lwd=2, col=pt.col(10, alpha=0.1))

## Fig 6.10
## kernel density ridge
dr.quake <- kdr(x=quake, xmin=c(70,-70), xmax=c(310, 80), verbose=TRUE)

pacific.map()
lines(plate[,1:2], col=pt.col(2, alpha=0.5), lwd=2)
points(dr.quake$end.points, cex=0.5, pch=16, col=pt.col(10,alpha=1))