# sample program to map kriging output
# in R, first load "maps" and "fields" 
#
x1<-matrix(scan('krigout'),ncol=4,byrow=T)
# define boundaries of x and y coordinates
x<-0.56*(1:100)-124
y<-0.24*(1:100)+25
z<-matrix(NA,ncol=100,nrow=100)
# rewrite third column of x1 into main matrix z,
# ignoring NA values
for(k in 1:length(x1[,1])){
i<-round((x1[k,1]+124)/0.56)
j<-round((x1[k,2]-25)/0.24)
z[i,j]<-x1[k,3]}
# compute range: can change this if desired
zr<-range(z,na.rm=T)
# draw map
map('state',xlim=c(-130,-60),ylim=c(25,54))
image.plot(x,y,z,add=T,zlim=zr)
map('state',xlim=c(-130,-60),ylim=c(25,54),add=T)