### Make sure you can install and load the following packages:
### sm, spatstat, splancs, spatial
## install.packages("sm")  ## to download this package
## install.packages("spatstat") 
## install.packages("splancs")
## install.packages("spatial") 
## install.packages("fields") 
## This package is available at http://cran.r-project.org/ ,
## click on "contributed packages" which is at the top of the page,
## under "ALL PLATFORMS", then look alphabetically.
library(sm)  ## to load the package into R.  You have to do this
                    ## every time you quit and re-start R, if you want to
		   ## use these functions
library(spatstat)
library(splancs)
library(spatial)
library(fields)
## Also if you want you can do
## library(help=spatstat)      ### briefly describes functions in spatstat.
## library()		   ### lists available installed libraries.
################################################

#### Simulate a Neyman-Scott model
nclust =  function(x0, y0, radius, n) {
                                  return(runifdisc(n, radius, x0, y0))
                                }
z = rNeymanScott(10, 0.2, nclust, radius=0.05, n=20)
par(mfrow=c(1,1))
plot(z,pch=".")
z1 = cbind(z$x,z$y)

### K-function and L-function for Neyman-Scott, with confidence bounds:
par(mfrow=c(2,1)) ## makes a 2x1 grid of plots
s = seq(.001,.3,length=50)
bdry = matrix(c(0,0,1,0,1,1,0,1,0,0),ncol=2,byrow=T)
k4 = khat(z1,bdry,s)
L4 = sqrt(k4/pi)-s
k4conf = Kenv.csr(npts(z1), bdry, 100, s)
plot(c(0,max(s)),c(0,max(k4conf$upper,k4)),
	type="n",xlab="distance",ylab="K4-hat")
points(s,k4,pch="*")
lines(s,k4)
lines(s,k4conf$upper,lty=2)
lines(s,k4conf$lower,lty=2)
L4upper = sqrt(k4conf$upper/pi) - s
L4lower = sqrt(k4conf$lower/pi) - s
plot(c(0,max(s)),c(min(L4lower,L4),max(L4upper,L4)),
	type="n",xlab="distance",ylab="L4-hat")
points(s,L4,pch="*")
lines(s,L4)
lines(s,L4upper,lty=2)
lines(s,L4lower,lty=2)


#### Simulate a Matern(I) process
z = rMaternI(100, 0.02)
par(mfrow=c(1,1))
plot(z,pch=".")
z1 = cbind(z$x,z$y)

### Simulate an SSI process
z = rSSI(.01,500,giveup=10000)
par(mfrow=c(1,1))
plot(z,pch=".")
z1 = cbind(z$x,z$y)


### K-function and L-function for SSI, with confidence bounds:
par(mfrow=c(2,1)) ## makes a 2x1 grid of plots
s = seq(.001,0.03,length=50)  ## change the 0.3 here to 0.03, to zoom in
k4 = khat(z1,bdry,s)
L4 = sqrt(k4/pi)-s
k4conf = Kenv.csr(npts(z1), bdry, 100, s)
plot(c(0,max(s)),c(0,max(k4conf$upper,k4)),
	type="n",xlab="distance",ylab="K4-hat")
points(s,k4,pch="*")
lines(s,k4)
lines(s,k4conf$upper,lty=2)
lines(s,k4conf$lower,lty=2)
L4upper = sqrt(k4conf$upper/pi) - s
L4lower = sqrt(k4conf$lower/pi) - s
plot(c(0,max(s)),c(min(L4lower,L4),max(L4upper,L4)),
	type="n",xlab="distance",ylab="L4-hat")
points(s,L4,pch="*")
lines(s,L4)
lines(s,L4upper,lty=2)
lines(s,L4lower,lty=2)


### Back to Spatial Data 
data(lennon)
image(lennon,col=grey(c(20:64)/64))

## Data as a list
mydata1 = 
data.frame(x=rep(1:50,times=rep(50,50)),y=rep(1: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))

par(mfrow=c(1,1))
image(mydata2,zlim=range(mydata2),col=grey(c(64:20)/64))


## variogram estimates:
par(mfrow=c(1,1))
my.ls0 = surf.ls(0,data2)
sill1 = var(data2$z)
c1 = correlogram(my.ls0,nint=100,plot=F) 
m1 = max((1:length(c1$x))[c1$x<.5]) ## to only go up to distance 0.5 in what follows
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="variograms") 
lines(a1, 2*sill1 - 2*expcov(a1,0.08)*sill1,col=3)
lines(a1, 2*sill1 - 2*gaucov(a1,0.12)*sill1,col=4)
lines(a1, 2*sill1 - 2*sphercov(a1,0.18)*sill1,col=5)
legend(0.16,2000,c("nonparametric","exponential (d = .08)", 
	"Gaussian (d = .12)", "spherical (d = .18)"),
    col=c(1,3,4,5),lty=c(2,1,1,1))


par(mfrow=c(3,1))
## 2) Fitting surfaces
zz = mydata2
zmin = min(zz)
zmax = max(zz)
n1 = 50
n2 = 50
x = seq(0,1,length=n1)
y = seq(0,1,length=n2)
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
image(zz,zlim=c(.8*zmin,1.2*zmax),col=grey(c(64:20)/64))
image(my.matrix.ls,zlim=c(.8*zmin,1.2*zmax),col=grey(c(64:20)/64))
my.gls = surf.gls(2,gaucov,data2,d=.12)
## my.gls is the fit of a quadratic surface to the data by GLS,
## estimating the variogram exponentially with range d=.1. (change the .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(.8*zmin,1.2*zmax), col=grey(c(64:20)/64))