# library(spatial)
# library(MASS)

# from Ripley and Venables
# Modern Applied Statistics with S, chapter 15


topo.ls = surf.ls(2,topo)
trsurf = trmat( topo.ls, 0, 6.5, 0, 6.5, 30)
eqscplot(trsurf, ,xlab="", ylab="", type = "n")
contour( trsurf, levels = seq(600, 1000, 25), add = T)
points(topo)




#  fit a model with exponential covariance


topo.kr = surf.gls(2, expcov, topo, d = 0.7)
correlogram( topo.kr, 25)




#  bottom of p 420 
topo.ls = surf.ls(4,topo)
trsurf = trmat( topo.ls, 0, 6.5, 0, 6.5, 30)
trs = trsurf
trs[ c("x", "y")] = expand.grid( x = trs$x, y = trs$y)
plt1 = levelplot(z ~ x*y, trs, aspect = 1, at = seq(650, 1000, 10))
plt2 = wireframe(  z~x*y, trs , aspect = c(1,0.5), screen = list( z = -30, x = -60))
print( plt1, position = c(0,0, .5, 1), more = T)
print(plt2, position = c(.45, 0, 1, 1))


# Kriging 


topo.ls = surf.ls(2,topo)
trsurf = trmat( topo.ls, 0, 6.5, 0, 6.5, 30)
eqscplot(trsurf, ,xlab="", ylab="", type = "n")
contour( trsurf, levels = seq(600, 1000, 25), add = T)
points(topo)
title(""LS Surface degree 2")


topo.kr = surf.ls(2,topo)
correlogram(topo.kr, 25)
d = seq( 0, 7, .1)
lines( d, expcov(d, .7))
lines( d, gaucov(d, 1, .3), lty = 3)



topo.kr = surf.gls( 2, gaucov, topo, d = 1, alpha = .3)
prsurf = prmat( topo.kr, 0, 6.5, 0, 6.5, 50)
eqscplot(prsurf, ,xlab="", ylab="", type = "n")
contour( prsurf, levels = seq(600, 1000, 25), add = T)
points(topo)


# Standard error surface 
sesurf = semat( topo.kr, 0, 6.5, 0, 6.5, 25)
hist(c(sesurf$z))  # we can see the range of SE
eqscplot(sesurf, ,xlab="", ylab="", type = "n")
contour( sesurf, levels = c(15, 33, 5), add = T)
points(topo)



# fit surface of degree 2
topo.kr = surf.gls( 2, gaucov, topo, d = 2, alpha = .3)
prsurf = prmat( topo.kr, 0, 6.5, 0, 6.5, 50)
eqscplot(prsurf, ,xlab="", ylab="", type = "n")
contour( prsurf, levels = seq(600, 1000, 25), add = T)
points(topo)


# Standard error surface 
sesurf = semat( topo.kr, 0, 6.5, 0, 6.5, 25)
hist(c(sesurf$z))  # we can see the range of SE
eqscplot(sesurf, ,xlab="", ylab="", type = "n")
contour( sesurf, levels = c(20, 25, 30), add = T)
points(topo)