library(fields)
library(sm)  
library(spatstat)
library(splancs)
library(spatial)
library(geoR)

## Simulate a Gaussian Random Field 
## with spherical variogram 

## simulate on a regular matrix, 10 by 10
c0 = .01
cs = .03
as = .3
mysim = grf(100, grid = "reg",cov.model = "spherical",cov.pars = c(cs, as), 
    nugget=c0, mean=3)
points(mysim) ## shows the data
plot(mysim) ## shows the semivariogram
par(mfrow=c(1,2))
image(mysim, zlim=range(mysim$data),col=grey(c(64:20)/64), main="gas prices") ## plots the simulation
zmin = min(mysim$data)
zmax = max(mysim$data)
zrng = zmax - zmin
zmid = zmin + zrng/2
plot(c(0,10),c(zmid-2*zrng/3,zmid+2*zrng/3),type="n",axes=F,xlab="",ylab="")
zgrid = seq(zmin,zmax,length=100)
## zgrid = vector of 100 equally-spaced numbers spanning range of the values.
image(c(-1:1),zgrid,matrix(rep(zgrid,2),ncol=100,byrow=T),add=T,col=gray((64:20)/64))
text(2.5,zmin,as.character(signif(zmin,2)),cex=1)
text(2.5,zmax,as.character(signif(zmax,2)),cex=1)
text(2.5,zmid,as.character(signif(zmid,2)),cex=1)
text(4.5,zmid,"$/gallon",srt=-90)

## simulate on an irregular matrix, 7 by 8
c0 = .01
cs = .03
as = .3
z1 = rep(seq(0,1,l=7),rep(8,7))
z2 = rep(seq(0,1,l=8),7)
mysim = grf(56, grid = cbind(z1,z2),cov.model = "spherical",cov.pars = c(cs, as), 
    nugget=c0, mean = 3)
plot(mysim) ## shows the semivariogram
par(mfrow=c(1,3))
points(mysim) ## shows the data
image(seq(0,1,l=7), seq(0,1,l=8), matrix(mysim$data,ncol=8,byrow=T), xlab="x", ylab="y",
    zlim=range(mysim$data),col=grey(c(64:20)/64), main="gas prices") ## plots the simulation
zmin = min(mysim$data)
zmax = max(mysim$data)
zrng = zmax - zmin
zmid = zmin + zrng/2
plot(c(0,10),c(zmid-2*zrng/3,zmid+2*zrng/3),type="n",axes=F,xlab="",ylab="")
zgrid = seq(zmin,zmax,length=100)
## zgrid = vector of 100 equally-spaced numbers spanning range of the values.
image(c(-1:1),zgrid,matrix(rep(zgrid,2),ncol=100,byrow=T),add=T,col=gray((64:20)/64))
text(2.5,zmin,as.character(signif(zmin,2)),cex=1)
text(2.5,zmax,as.character(signif(zmax,2)),cex=1)
text(2.5,zmid,as.character(signif(zmid,2)),cex=1)
text(4.5,zmid,"$/gallon",srt=-90)