#####  Inputting geostatistical data, sampled at n=34 fixed (non-random) locations.
#####  Range of the values will be from 2.3 to 21.
#####  X-coordinates will range from 0 to 1.  Same for y-coordinates.

library(sm)
library(spatstat)
library(spatial)
library(splancs)


n = 34
par(mfrow=c(1,1))
plot(c(0,1),c(0,1),type="n",xlab="x-coordinate",ylab="y-coordinate")
x = rep(0,n)
y = rep(0,n)
z = rep(0,n)
for(i in 1:n){
    a = locator(1)
    x[i] = a$x
    y[i] = a$y
    z[i] = ask("% Edwards")
    text(x[i],y[i],as.character(z[i]))
}

zmin = 2.3
zmax = 20.2

### the range to plot (should be a bit wider than the data's range)
zmin1 = -5
zmax1 = 30

######  Save the data into a file:
sink("myzzfile")
cat(n,x,y,z)
sink()

#####  To read the data back in:
z1 = scan("myzzfile")
n = z1[1]
x = z1[2:(n+1)]
y = z1[(n+2):(2*n+1)]
z = z1[(2*n+2):(3*n+1)]

    

##### Plot the numerical values again:
plot(c(0,1),c(0,1),type="n",xlab="x-coordinate",ylab="y-coordinate")
text(x,y,as.character(z))



par(mfrow=c(2,3))
data2 = data.frame(x=x,y=y,z=z)
## data2 has 3 long vectors: the x-coordinates, the y's, and the z's.
my.ls = surf.ls(4,data2)
## my.ls is the least squares fit of a degree 4 polynomial surface to the data
my.matrix.ls = trmat(my.ls,0,1,0,1,17)  ## you can use whatever you like here, for 17
### it's going to output the results on a 17 x 17 grid
image(my.matrix.ls,zlim=c(zmin1,zmax1),col=grey(c(64:20)/64),main="LS fit, degree 4")
## the next 9 lines are just to make a legend:
zrng = zmax1 - zmin1
zmid = zmin1 + zrng/2
plot(c(0,10),c(zmid-2*zrng/3,zmid+2*zrng/3),type="n",axes=F,xlab="",ylab="")
zgrid = seq(zmin1,zmax1,length=100)
image(c(-1:1),zgrid,matrix(rep(zgrid,2),ncol=100,byrow=T),add=T,col=gray((64:20)/64))
text(2.5,zmin1,as.character(signif(zmin1,2)),cex=1)
text(2.5,zmax1,as.character(signif(zmax1,2)),cex=1)
text(2.5,zmid,as.character(signif(zmid,2)),cex=1)
text(4.5,zmid,"% for Edwards",srt=-90)
## GLS fit, degree 4:
my.gls = surf.gls(4,expcov,data2,d=0.1)
## my.gls is the fit of a degree 4 surface to the data by GLS,
## estimating the variogram exponentially with range d=0.1. (change the 0.1 if you'd like)
my.matrix.gls = trmat(my.gls,0,1,0,1,17)  ## again, feel free to change 17 
image(my.matrix.gls, zlim=c(zmin1,zmax1), col=grey(c(64:20)/64),
    main="GLS fit, degree 4\n with expcov(0.1)")
my.krig = prmat(my.gls,0,1,0,1,17)  ## Kriging estimates
image(my.krig,zlim=c(zmin1,zmax1), col=grey(c(64:20)/64),
    main="Kriging of deg. 4 fit\n with expcov(5)")
my.ses = semat(my.gls,0,1,0,1,17) ## SEs for Kriging estimates
image(my.ses,zlim=range(my.ses$z), col=grey(c(64:20)/64),main="Kriging SEs")
## the next 11 lines are just to make a legend for the kriging estimates
semax = max(my.ses$z)
semin = min(my.ses$z)
serng = semax - semin
semid = semin + serng/2
plot(c(0,10),c(semid-2*serng/3,semid+2*serng/3),type="n",axes=F,xlab="",ylab="")
segrid = seq(semin,semax,length=100)
image(c(-1:1),segrid,matrix(rep(segrid,2),ncol=100,byrow=T),add=T,col=gray((64:20)/64))
text(2.5,semin,as.character(signif(semin,2)),cex=1)
text(2.5,semax,as.character(signif(semax,2)),cex=1)
text(2.5,semid,as.character(signif(semid,2)),cex=1)
text(4.5,semid,"SE (%)",srt=-90)


### quadrat means, on a 10x10 grid
x2 = matrix(mean(z),ncol=10,nrow=10)
for(i in 1:10){
for(j in 1:10){
a1 = c(1:n)[(x<i/10) & (x >= (i-1)/10) & (y < j/10) & (y >= (j-1)/10)]
if(length(a1)>0) x2[i,j] = mean(z[a1])
}}
mean(x2) ## should be same as mean(z)
mean(z)

##### Plot the quadrat means
par(mfrow=c(1,2))  ## makes a 1x2 grid of plots on the graphic screen
plot(c(0,1),c(0,1),type="n",xlab="x-coordinate",ylab="y-coordinate",
main="Quadrat means of \n Edwards' %")
image(x=c(0:9)/10+.05,y=c(0:9)/10+.05,z=x2,zlim=c(zmin1,zmax1),col=grey(c(64:20)/64),add=T)
######### LEGEND:
plot(c(0,10),c(zmid-2*zrng/3,zmid+2*zrng/3),type="n",axes=F,xlab="",ylab="")
zgrid = seq(zmin1,zmax1,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,
    zlim=c(zmin1,zmax1),col=gray((64:20)/64))
text(2.5,zmin1,as.character(signif(zmin1,2)),cex=1)
text(2.5,zmax1,as.character(signif(zmax1,2)),cex=1)
text(2.5,zmid,as.character(signif(zmid,2)),cex=1)
text(4.5,zmid,"%",srt=-90)

## compare with:
par(mfrow=c(1,1))  ## back to a 1x1 grid of plots
filled.contour(x2)



####### VARIOGRAMS AND COVARIOGRAMS:
sill1 = var(c(z))
my.ls0 = surf.ls(0,data2)
c1 = correlogram(my.ls0,100,plot=F)  ## 100 lags -- you can change this
plot(c1$x, c1$y * sill1, xlab="distance", ylab="covariance",type="l") ## covariogram
m1 = max((1:length(c1$x))[c1$x<.5]) ## to only go up to distance 0.5
plot(c1$x[1:m1], c1$y[1:m1] * sill1, xlab="distance", ylab="covariance",type="l") 
plot(c1$x[1:m1], c1$y[1:m1], xlab="distance", ylab="correlation",type="l") ## correlogram
c2 = variogram(my.ls0,100,plot=F)  ## again, 100 lags -- you may change this.
plot(c2$x[1:m1], c2$y[1:m1], xlab="distance, h", ylab="gamma(h)",type="l",
main="semi-variogram of raw data") 
c3 = variogram(my.ls,100,plot=F)  ## fits semi-variogram to LS residuals
plot(c3$x[1:m1], c3$y[1:m1], xlab="distance, h", ylab="gamma(h)",type="l", 
    main="semi-variogram of LS residuals")
c4 = variogram(my.gls,100,plot=F) ## fits semi-variogram to GLS residuals
plot(c4$x[1:m1], c4$y[1:m1], xlab="distance, h", ylab="gamma(h)",type="l",
    main="semi-variogram of GLS residuals")



### VARIOGRAM MODELS:
par(mfrow=c(2,2))
a1 = seq(0.01,0.5,length=50)
plot(c1$x[1:m1], 2*sill1 - 2*c1$y[1:m1] * sill1, 
    xlab="distance, h", ylab="2gamma(h)",type="l",
    lty=2,main="exponential, d=5") 
lines(a1, 2*sill1 - 2*expcov(a1,5)*sill1)
plot(c1$x[1:m1], 2*sill1 - 2*c1$y[1:m1] * sill1, xlab="distance, h", ylab="2gamma(h)",type="l",
    lty=2,main="exponential, d=0.5") 
lines(a1, 2*sill1 - 2*expcov(a1,0.5)*sill1)
plot(c1$x[1:m1], 2*sill1 - 2*c1$y[1:m1] * sill1, xlab="distance, h", ylab="2gamma(h)",type="l",
    lty=2,main="exponential, d=0.1") 
lines(a1, 2*sill1 - 2*expcov(a1,0.1)*sill1)
plot(c1$x[1:m1], 2*sill1 - 2*c1$y[1:m1] * sill1, xlab="distance, h", ylab="2gamma(h)",type="l",
    lty=2,main="exponential, d=0.01") 
lines(a1, 2*sill1 - 2*expcov(a1,0.01)*sill1)

par(mfrow=c(2,2))
plot(c1$x[1:m1], 2*sill1 - 2*c1$y[1:m1] * sill1, xlab="distance, h", ylab="2gamma(h)",type="l",
    lty=2,main="Gaussian, d=5") 
lines(a1, 2*sill1 - 2*gaucov(a1,5)*sill1)
plot(c1$x[1:m1], 2*sill1 - 2*c1$y[1:m1] * sill1, xlab="distance, h", ylab="2gamma(h)",type="l",
    lty=2,main="Gaussian, d=0.5") 
lines(a1, 2*sill1 - 2*gaucov(a1,0.5)*sill1)
plot(c1$x[1:m1], 2*sill1 - 2*c1$y[1:m1] * sill1, xlab="distance, h", ylab="2gamma(h)",type="l",
    lty=2,main="Gaussian, d=0.1") 
lines(a1, 2*sill1 - 2*gaucov(a1,0.1)*sill1)
plot(c1$x[1:m1], 2*sill1 - 2*c1$y[1:m1] * sill1, xlab="distance, h", ylab="2gamma(h)",type="l",
    lty=2,main="Gaussian, d=0.01") 
lines(a1, 2*sill1 - 2*gaucov(a1,0.01)*sill1)

par(mfrow=c(2,2))
plot(c1$x[1:m1], 2*sill1 - 2*c1$y[1:m1] * sill1, xlab="distance, h", ylab="2gamma(h)",type="l",
    lty=2,main="spherical, d=5") 
lines(a1, 2*sill1 - 2*sphercov(a1,5)*sill1)
plot(c1$x[1:m1], 2*sill1 - 2*c1$y[1:m1] * sill1, xlab="distance, h", ylab="2gamma(h)",type="l",
    lty=2,main="spherical, d=0.5") 
lines(a1, 2*sill1 - 2*sphercov(a1,0.5)*sill1)
plot(c1$x[1:m1], 2*sill1 - 2*c1$y[1:m1] * sill1, xlab="distance, h", ylab="2gamma(h)",type="l",
    lty=2,main="spherical, d=0.1") 
lines(a1, 2*sill1 - 2*sphercov(a1,0.1)*sill1)
plot(c1$x[1:m1], 2*sill1 - 2*c1$y[1:m1] * sill1, xlab="distance, h", ylab="2gamma(h)",type="l",
    lty=2,main="spherical, d=0.01") 
lines(a1, 2*sill1 - 2*sphercov(a1,0.01)*sill1)

## Compare all 3 with d = 0.1:
par(mfrow=c(1,1))
plot(c1$x[1:m1], 2*sill1 - 2*c1$y[1:m1] * sill1, xlab="distance, h", ylab="2gamma(h)",type="l",
    lty=2,main="covariograms, d=0.1") 
lines(a1, 2*sill1 - 2*expcov(a1,0.1)*sill1,col=3)
lines(a1, 2*sill1 - 2*gaucov(a1,0.1)*sill1,col=4)
lines(a1, 2*sill1 - 2*sphercov(a1,0.1)*sill1,col=5)
legend(0.3,12,c("nonparametric","exponential", "Gaussian", "spherical"),
    col=c(1,3,4,5),lty=c(2,1,1,1))




### COVARIOGRAM MODELS:
par(mfrow=c(2,2))
a1 = seq(0.01,0.5,length=50)
plot(c1$x[1:m1], c1$y[1:m1] * sill1, xlab="distance, h", ylab="covariance",type="l",
    lty=2,main="exponential, d=5") 
lines(a1, expcov(a1,5)*sill1)
plot(c1$x[1:m1], c1$y[1:m1] * sill1, xlab="distance, h", ylab="covariance",type="l",
    lty=2,main="exponential, d=0.5") 
lines(a1, expcov(a1,0.5)*sill1)
plot(c1$x[1:m1], c1$y[1:m1] * sill1, xlab="distance, h", ylab="covariance",type="l",
    lty=2,main="exponential, d=0.1") 
lines(a1, expcov(a1,0.1)*sill1)
plot(c1$x[1:m1], c1$y[1:m1] * sill1, xlab="distance, h", ylab="covariance",type="l",
    lty=2,main="exponential, d=0.01") 
lines(a1, expcov(a1,0.01)*sill1)

par(mfrow=c(2,2))
plot(c1$x[1:m1], c1$y[1:m1] * sill1, xlab="distance, h", ylab="covariance",type="l",
    lty=2,main="Gaussian, d=5") 
lines(a1, gaucov(a1,5)*sill1)
plot(c1$x[1:m1], c1$y[1:m1] * sill1, xlab="distance, h", ylab="covariance",type="l",
    lty=2,main="Gaussian, d=0.5") 
lines(a1, gaucov(a1,0.5)*sill1)
plot(c1$x[1:m1], c1$y[1:m1] * sill1, xlab="distance, h", ylab="covariance",type="l",
    lty=2,main="Gaussian, d=0.1") 
lines(a1, gaucov(a1,0.1)*sill1)
plot(c1$x[1:m1], c1$y[1:m1] * sill1, xlab="distance, h", ylab="covariance",type="l",
    lty=2,main="Gaussian, d=0.01") 
lines(a1, gaucov(a1,0.01)*sill1)

par(mfrow=c(2,2))
plot(c1$x[1:m1], c1$y[1:m1] * sill1, xlab="distance, h", ylab="covariance",type="l",
    lty=2,main="spherical, d=5") 
lines(a1, sphercov(a1,5)*sill1)
plot(c1$x[1:m1], c1$y[1:m1] * sill1, xlab="distance, h", ylab="covariance",type="l",
    lty=2,main="spherical, d=0.5") 
lines(a1, sphercov(a1,0.5)*sill1)
plot(c1$x[1:m1], c1$y[1:m1] * sill1, xlab="distance, h", ylab="covariance",type="l",
    lty=2,main="spherical, d=0.1") 
lines(a1, sphercov(a1,0.1)*sill1)
plot(c1$x[1:m1], c1$y[1:m1] * sill1, xlab="distance, h", ylab="covariance",type="l",
    lty=2,main="spherical, d=0.01") 
lines(a1, sphercov(a1,0.01)*sill1)

## Compare all 3 with d = 0.1:
par(mfrow=c(1,1))
plot(c1$x[1:m1], c1$y[1:m1] * sill1, xlab="distance, h", ylab="covariance",type="l",
    lty=2,main="covariograms, d=0.1") 
lines(a1, expcov(a1,0.1)*sill1,col=3)
lines(a1, gaucov(a1,0.1)*sill1,col=4)
lines(a1, sphercov(a1,0.1)*sill1,col=5)
legend(0.3,18,c("nonparametric","exponential", "Gaussian", "spherical"),
    col=c(1,3,4,5),lty=c(2,1,1,1))


##  Variogram models:
vgm1 <- variogram(log(zinc)~1, ~x+y, meuse)
fit.variogram(vgm1, vgm(1,"Sph",300,1))