library(fields,lib.loc='/Users/ssain/myR2') load('nc2.save') source('multiKrig.r') s <- nc.winter[,c("lon","lat")] Y <- nc.winter[,c("prcp","temp")] Z <- matrix(nc.winter[,"ele"],ncol=1) dimnames(Z) <- list(NULL,"ele") cat('some initial plots...\n') par(mfrow=c(1,2),mar=c(2,2,5,2)) quilt.plot(nc.winter$lon,nc.winter$lat,nc.winter$prcp,main='prcp') US(add=TRUE) quilt.plot(nc.winter$lon,nc.winter$lat,nc.winter$temp,main='temp') US(add=TRUE) cat('parameter estimation...\n') obj.est <- multi.est(pms=list(rho=c(1,1),sp=c(1,0.5,2)), s=s,Y=Y,Z=Z) cat('model fitting...\n') obj.fit <- multiKrig(s=s,Y=Y,Z=Z, cov.function='multi.cov', cov.args=list(theta=2.5,rho=obj.est$rho[2]/obj.est$rho[1],smoothness=1.5), wght.function='multi.wght', wght.args=list(sp=obj.est$sp[2:3]/obj.est$sp[1])) cat('output...\n') print(obj.fit) cat('some plots...\n') par(mfrow=c(2,3),ask=TRUE) plot(obj.fit) cat('prediction....\n') obj.pre <- predict(obj.fit,s=nc.dem[nc.dem[,4]==1,1:2],Z=nc.dem[nc.dem[,4]==1,3]) obj.pre <- matrix(obj.pre,ncol=2) pr <- rep(NA,nrow(nc.dem)) pr[nc.dem[,4]==1] <- obj.pre[,1] tm <- rep(NA,nrow(nc.dem)) tm[nc.dem[,4]==1] <- obj.pre[,2] x <- sort(unique(nc.dem[,1])) y <- sort(unique(nc.dem[,2])) # some plots par(mfrow=c(1,2),mar=c(2,2,5,2)) image.plot(x,y,matrix(pr,nrow=length(x),ncol=length(y)),zlim=range(c(pr,Y[,1]),na.rm=T), main='prcp: predicted values') US(add=TRUE) image.plot(x,y,matrix(tm,nrow=length(x),ncol=length(y)),zlim=range(c(tm,Y[,2]),na.rm=T), main='temp: predicted values') US(add=TRUE) #quilt.plot(nc.winter$lon,nc.winter$lat,nc.winter$temp,zlim=range(c(tm,Y[,2]),na.rm=T),add.legend=FALSE,add=TRUE)