## Display of spatial data loc=cbind(rep(1:4, 3), rep(1:3, each=4)) plot(loc[,1], loc[,2], xlim=c(0, 4), ylim=c(0,3), xlab="x-coordinate", ylab="y-coordinate") for(i in 1:nrow(loc)){ text(loc[i,1]-0.2, loc[i,2], i) } #postscript(file="AirText.eps", horizontal=FALSE, onefile=FALSE, # height=4, width=5, pointsize=10) plot(paws\$Lat, paws\$Long, xlab="Latitude", ylab="Longtitude", type="n") for(i in 1:nrow(paws)){ text(paws\$Lat, paws\$Long, paws\$Air, cex=0.5) } #dev.off() #### Generate spatial normal data ## Define the Matern covariance function cov.matern=function(x, nu = 2, alpha = 1, vars=1) { if(nu == 0.5) return(vars*exp( - x * alpha)) ismatrix <- is.matrix(x) if(ismatrix){nr=nrow(x); nl=ncol(x)} x <- c(alpha * x) output <- rep(1, length(x)) n <- sum(x > 0) if(n > 0) { x1 <- x[x > 0] output[x > 0] <- (1/((2^(nu - 1)) * gamma(nu))) * (x1^nu) * besselK( x1, nu) } if(ismatrix){ output <- matrix(output, nr, nl) } vars*output } for(nu in seq(0.5, 4, by=0.25)){ temp=cov.matern(2.8*sqrt(nu), nu) cat(temp, " ", nu, "\n") } ## This function defines the distance matrix distmat=function(x) { as.matrix(dist(x)) } d=distmat(loc) set.seed(10) library(MASS) y=mvrnorm(mu=rep(0, 12), Sigma=10*exp(-d)) plot(loc[,1], loc[,2], xlim=c(0.5, 4.5), ylim=c(0,3.5), xlab="x-coordinate", ylab="y-coordinate") for(i in 1:nrow(loc)){ text(loc[i,1], loc[i,2]+0.2, i) text(loc[i,1]-0.1, loc[i,2]-0.2, round(y[i],3) ) } yma=matrix(y, ncol=3) persp(x=1:4, y=1:3, z=yma, xlab="x", ylab="y", zlab=" z", theta=45, phi=35, r=5, expand=0.6, axes=T, ticktype="detailed") V=cov.matern(distmat(cbind(rep(0:20, 21)/20, rep(0:20, each=21)/20)),nu=3, alpha=7*sqrt(3)) set.seed(20) z=mvrnorm(mu=rep(0, 21^2), Sigma=V) z=matrix(z, ncol=21) postscript(file="SimuProc4.eps", horizontal=FALSE, onefile=FALSE, height=5, width=5, pointsize=10) persp(x=(0:20)/21, y=(0:20)/21, z, theta=45, phi=35, r=5, expand=0.6, axes=T, ticktype="detailed", xlab="x", ylab="y", zlab="z") dev.off() postscript(file="SimuProc5.eps", horizontal=FALSE, onefile=FALSE, height=5, width=5, pointsize=10) filled.contour(x=0:20, y=0:20, z, color.palette=gray.colors) dev.off()