### 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") 
## 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)
## Also if you want you can do
library(help=spatstat)      ### briefly describes functions in spatstat.
library()		   ### lists available installed libraries.
################################################

#####  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)

##### Quadrat counts

x = matrix(0,ncol=10,nrow=10)
for(i in 1:10){
for(j in 1:10){
for(k in 1:n){
if((x1[k]<i/10) && (x1[k] >= (i-1)/10) &&
(y1[k]<j/10) && (y1[k] >= (j-1)/10)) x[i,j] = x[i,j] + 1
}
}}
sum(x) ## should = n, which here is 23

##### Plot the quadrat counts
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 counts of \n Celebrity homes")
image(x=c(0:9)/10+.05,y=c(0:9)/10+.05,z=x,col=grey(c(64:20)/64),add=T)
## contour(x=c(0:9)/10+.05,y=c(0:9)/10+.05,z=x,nlevels=3,add=T)
points(x1,y1,pch=".")
######### LEGEND:
zmin = min(x)
zmax = max(x)
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,"Values",srt=-90)

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

######### K-function:
### (Remember to LOAD SPLANCS first!!!)
b1 = as.points(x1,y1)
n = npts(b1)
bdry = matrix(c(0,0,1,0,1,1,0,1,0,0),ncol=2,byrow=T)
plot(c(0,1),c(0,1),type="n",xlab="x",ylab="y")
lines(bdry)
points(b1,pch="*")

par(mfrow=c(2,1)) ## makes a 2x1 grid of plots
s = seq(.001,.3,length=50)
k4 = khat(b1,bdry,s)
plot(s,k4,xlab="distance",ylab="K4-hat",pch="*")
lines(s,k4)
lines(s,pi*s^2,lty=2)
L4 = sqrt(k4/pi)-s
plot(c(0,.3),c(-.02,.1),type="n",xlab="lag, h",ylab="L4-hat(h)")
points(s,L4,pch="*")
lines(s,L4)
lines(s,rep(0,50),lty=2)


### CONFIDENCE BOUNDS FOR K-FUNCTION

k4conf = Kenv.csr(npts(b1), bdry, 1000, 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=3)
lines(s,k4conf$lower,lty=3)

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)


###  Kernel smoothing
bdw = sqrt(bw.nrd0(x1)^2+bw.nrd0(y1)^2)  ## default bandwidth
z = kernel2d(b1,bdry,bdw)
par(mfrow=c(1,2))
image(z,col=gray((64:20)/64),xlab="x",ylab="y")
points(b1)
x4 = (0:100)/100*(max(z$z)-min(z$z))+min(z$z)
plot(c(0,10),c(.8*min(x4),1.2*max(x4)),type="n",axes=F,xlab="",ylab="")
image(c(-1:1),x4,matrix(rep(x4,2),ncol=101,byrow=T),add=T,col=gray((64:20)/64))
text(2,min(x4),as.character(signif(min(x4),2)),cex=1)
text(2,(max(x4)+min(x4))/2,as.character(signif((max(x4)+min(x4))/2,2)),cex=1)
text(2,max(x4),as.character(signif(max(x4),2)),cex=1)
mtext(s=3,l=-3,at=1,"Rate (pts per unit area)")
## repeat the above, trying other values of bdw, for more or less smoothing









#### 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)
k4 = khat(z1,bdry,s)
L4 = sqrt(k4/pi)-s
k4conf = Kenv.csr(npts(z1), bdry, 1000, 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=".")


### 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)



## Brownian sheet

m = 101
m1 = m-1
y = matrix(rep(0,m*m),ncol=m)
for(i in c(2:m)) y[1,i] = y[1,(i-1)] + rnorm(1)/2
for(i in c(2:m)) y[i,1] = y[(i-1),1] + rnorm(1)/2
for(i in c(2:m)){
for(j in c(2:m)){
y[i,j] = (y[(i-1),j] + rnorm(1))/2 + (y[i,(j-1)] + rnorm(1))/2
}
cat(i,"\n")
}
x1 = c(0:m1)
x2 = c(0:m1)
x3 = y[(2:m),(2:m)]-y[1,1]
x4 = (0:100)/100*(max(x3)-min(x3))+min(x3)
par(mfrow=c(1,2))
image(x1,x2,x3,col=gray((64:0)/64),xlab="column",ylab="row")
plot(c(0,10),c(1.2*min(x4),1.2*max(x4)),type="n",axes=F,xlab="",ylab="")
image(c(-1:1),x4,matrix(rep(x4,2),ncol=101,byrow=T),add=T,col=gray((64:0)/64))
text(2,min(x4),as.character(signif(min(x4),2)),cex=1)
text(2,(max(x4)+min(x4))/2,as.character(signif((max(x4)+min(x4))/2,2)),cex=1)
text(2,max(x4),as.character(signif(max(x4),2)),cex=1)
## I use the function signif() here so that it only prints 2 sign. digits.