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

#####  Inputting points, using the mouse:
n = 23
plot(c(0,1),c(0,1),type="n",xlab="longitude",ylab="latitude",main="Star Homes")
x1 = rep(0,n)
y1 = rep(0,n)
for(i in 1:n){
z1 = locator(1)
x1[i] = z1$x
y1[i] = z1$y
points(x1[i],y1[i])
}

##### Plot the points
plot(c(0,1),c(0,1),type="n",xlab="x-coordinate",ylab="y-coordinate",
main="Star homes")
points(x1,y1,pch=1)

##### Image example
data(lennon)
image(lennon,col=grey(c(20:64)/64))

## Data as a list
mydata1 = 
data.frame(x=rep(seq(0,1,l=50),times=rep(50,50)),y=rep(seq(0,1,l=50),50),
    z=c(t(lennon[121:170,131:180])))

## Data as a matrix
mydata2 = lennon[121:170,131:180]
image(mydata2,zlim=range(mydata2),col=grey(c(64:20)/64))

#####  Inputting spatial data, on a 7 x 8 grid
#####  Range of the values will be from 2.80 to 3.40
#####  X-coordinates will range from 0 to 1.  Same for y-coordinates.

ask = function(message = "Type in datum"){
    eval(parse(prompt = paste(message, ": ", sep="")))
}

n1 = 7
n2 = 8
par(mfrow=c(1,2))
plot(c(0,1),c(0,1),type="n",xlab="x-coordinate",ylab="y-coordinate")
zmin = 2.70
zmax = 3.50
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)
zz = matrix(0, nrow=n1,ncol=n2)
for(i in 1:n1){
    for(j in 1:n2){
	zz[i,j] = ask(paste("row ",as.character(i),", column ",as.character(j)))
	par(mfg=c(1,1,1,2))
	plot(c(0,1),c(0,1),type="n",xlab="x-coordinate",ylab="y-coordinate")
	image(zz,zlim=c(zmin,zmax),col=grey(c(64:20)/64),add=T)
    }
}

##### change a particular value
zz[3,3] = 3.1
	
##### plot the final result
par(mfrow=c(1,2))
plot(c(0,1),c(0,1),type="n",xlab="x-coordinate",ylab="y-coordinate",main="gas prices")
image(zz,zlim=c(zmin,zmax),col=grey(c(64:20)/64),add=T)
zmin = 2.80
zmax = 3.40
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)

##### you might also try these:
filled.contour(zz)
persp(zz)
contour(zz)

######  Save the data into a file:
sink("mygas")
cat(zz)
sink()

#####  To read the data back in:
zz = matrix(scan("mygas"),ncol=n2)

######## Pocket plot
## e = unit vector pointing north
y = matrix(0,nrow=n1,ncol=n1)
for(i in c(1:(n1-1))){
for(j in c((i+1):n1)){
y[i,j] = mean(sqrt(abs(zz[i,] - zz[j,])))
y[j,i] = y[i,j]
}
}

z = rep(0,(n1-1))  ### z[k] will represent y-double-bar[k]
for(k in c(1:(n1-2))){
cat(k,"\n")
temp = y[1,(1+k)]
for(i in c(2:(n1-k))){
temp = c(temp,y[i,(i+k)])
}
z[k] = mean(temp)
}
z[n1-1] = y[1,n1]


p = matrix(rep(0,(n1*n1)),ncol=n1)
for(i in c(1:(n1-1))){
for(j in c((i+1):n1)){
k = j-i
p[i,j] = y[i,j] - z[k]
p[j,i] = p[i,j]
}
}

### This is just if you want to only plot boxplots of a subset of rows, since 
### otherwise you have too many rows to plot:
#p1 = p[5,c(1:4,6:n1)]
#for(i in seq(5,n1-5,by=10)){
#p1 = cbind(p1,p[i,c(1:n1)[-i]])
#}
## Then you'd replace p by p1 in what's below.

par(mfrow=c(1,1))
plot(c(1,7),range(p),xlab="Row number",ylab="P_jk",type="n",xaxt="n")
boxplot(p[,1],p[,2],p[,3],p[,4],p[,5],p[,6],p[,7],
    ### I don't know any way to do this other than to list all the columns individually.
names=as.character(seq(1,n1)),add=T)

### Convert lattice data zz to geostatistical data mydata1
mydata1 = 
data.frame(x=rep(seq(0,1,l=7),times=rep(8,7)),y=rep(seq(0,1,l=8),7),
    z=c(t(zz)))

## zz = mydata2, to estimate the variogram of the Lennon image
####### SEMI-VARIOGRAM:
par(mfrow=c(1,1))
sill1 = var(c(zz))
my.ls0 = surf.ls(0,mydata1)
c1 = correlogram(my.ls0,100,plot=F)  
m1 = max((1:length(c1$x))[c1$x<.5]) ## to only go up to distance .5
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") ## semi-variogram
lines(c2$x[1:m1], c2$y[1:m1]+.5*1.96*sqrt(2/c2$cnt[1:m1])*2*c2$y[1:m1],lty=2)
lines(c2$x[1:m1], c2$y[1:m1]-.5*1.96*sqrt(2/c2$cnt[1:m1])*2*c2$y[1:m1],lty=2)
r = c1$x[1:m1]
b1 = 2*sill1 - 2*c1$y[1:m1] * sill1

## This fits just the middle part of the spherical variogram model
#par(mfrow=c(1,1))
#plot(r, b1, xlab="distance, h", ylab="2gamma(h)",type="l",
#    lty=2,main=" spherical variogram") 
#g2 = data.frame(h = r, b1 = b1)
#b2 = nls(b1 ~ c0 + cs*(1.5*abs(h)/as - .5*(abs(h)/as)^3),
#    start = list(c0 = .001, cs = .01, as = .3),data=g2,trace=TRUE)
#summary(b2)
#c0 = 0.0015284
#cs = -0.0198162
#as = -0.3643008
#lines(r, (c0 + cs*(1.5*abs(r)/as - .5*(abs(r)/as)^3))*(r <= as)+(c0+cs)*(r > as))

## This fits the whole spherical variogram model better
plot(r, b1, xlab="distance, h", ylab="2gamma(h)",type="l",
    lty=1,main=" spherical variogram") 
points(r,b1)
lines(r, b1+1.96*sqrt(2/c1$cnt[1:m1])*b1,lty=2)
lines(r, b1-1.96*sqrt(2/c1$cnt[1:m1])*b1,lty=2)
myg = function(p){
    c0 = p[1]
    cs = p[2]
    as = p[3]
    if((c0 < 0) || (cs < 0) || (as < 0)) return(9999999e+99)
    aa = r
    for(i in 1:length(aa)){
	if(r[i] == 0) aa[i] = 0
	else if(r[i] <= as) aa[i] = c0 + cs*(1.5*abs(r[i])/as - .5*(abs(r[i])/as)^3)
	if(r[i] > as) aa[i] = c0 + cs
    }
    sum((b1-aa)^2)
}
p = c(3,5,1)
g3 = optim(p,myg,hessian=T)
g3$par
c0 = g3$par[1]
cs = g3$par[2]
as = g3$par[3]
aa = r
for(i in 1:length(aa)){
	if(r[i] == 0) aa[i] = 0
	else if(r[i] <= as) aa[i] = c0 + cs*(1.5*abs(r[i])/as - .5*(abs(r[i])/as)^3)
	if(r[i] > as) aa[i] = c0 + cs
    }
lines(r, aa,lty=3)
p = g3$par
g4 = optim(p,myg,hessian=T)
c0 = g4$par[1]
cs = g4$par[2]
as = g4$par[3]
aa = r
for(i in 1:length(aa)){
	if(r[i] == 0) aa[i] = 0
	else if(r[i] <= as) aa[i] = c0 + cs*(1.5*abs(r[i])/as - .5*(abs(r[i])/as)^3)
	if(r[i] > as) aa[i] = c0 + cs
    }
lines(r, aa,lty=4)

plot(r, b1, xlab="distance, h", ylab="2gamma(h)",type="l",
    lty=1,main="rational quadratic variogram") 
points(r,b1)
lines(r, b1+1.96*sqrt(2/c1$cnt[1:m1])*b1,lty=2)
lines(r, b1-1.96*sqrt(2/c1$cnt[1:m1])*b1,lty=2)
myg = function(p){
    c0 = p[1]
    cr = p[2]
    ar = p[3]
    if((c0 < 0) || (cr < 0) || (ar < 0)) return(9999999e+99)
    aa = r
    for(i in 1:length(aa)){
	if(r[i] == 0) aa[i] = 0
	else aa[i] = c0 + cr*r[i]^2/(1+r[i]^2/ar)
    }
    sum((b1-aa)^2)
}
p = c(3,10000,.1)
g3 = optim(p,myg,hessian=T)
g3$par
c0 = g3$par[1]
cr = g3$par[2]
ar = g3$par[3]
aa = r
for(i in 1:length(aa)){
	if(r[i] == 0) aa[i] = 0
	else aa[i] = c0 + cr*r[i]^2/(1+r[i]^2/ar)
    }
lines(r, aa,lty=3)

p = g3$par
g4 = optim(p,myg,hessian=T)
c0 = g4$par[1]
cr = g4$par[2]
ar = g4$par[3]
aa = r
for(i in 1:length(aa)){
	if(r[i] == 0) aa[i] = 0
	else aa[i] = c0 + cr*r[i]^2/(1+r[i]^2/ar)
    }
lines(r, aa,lty=4)

## Fitting surfaces and kriging
n1 = n2 = 50
x = seq(0,1,length=n1)
y = seq(0,1,length=n2)
zz = mydata2
data1 = list(x=x,y=y,z=zz)
data2 = data.frame(x=rep(x,times=rep(n2,n1)),y=rep(y,n1),z=c(t(zz)))
## 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
par(mfrow=c(2,2))
zmin = min(c(zz))
zmax = max(c(zz))
image(my.matrix.ls,zlim=c(.8*zmin,1.2*zmax),col=grey(c(64:20)/64),
main="degree 4 polynomial smoothing")
image(zz,zlim=c(.8*zmin,1.2*zmax),col=grey(c(64:20)/64),
main="raw data, 50 x 50") 
my.gls = surf.gls(2,expcov,data2,d=.1)
## my.gls is the fit of a quadratic surface to the data by GLS,
## estimating the variogram exponentially with range d=5. (change the 5 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(.8*zmin,1.2*zmax), col=grey(c(64:20)/64),
main="degree 2 polynomial smoothing \n using GLS and expcov(d=5)")
my.krig = prmat(my.gls,0,1,0,1,100)  ## Kriging estimates
image(my.krig,zlim=c(.8*zmin,1.2*zmax), col=grey(c(64:20)/64),
main="Kriging of raw data \n onto a 100 x 100 grid \n using GLS and expcov(d=5)")