x = c("ancheta", "anderson", "brett-esborn", "cassel", "cetinkaya",
"chen","diez", "enos", "firestine", "jiang",
"li", "mejgaard", "nesbitt", "nguyen", "parfenov",
"pfeiffe", "rundel", "shanata", "snyder", "sun",
"tsai", "wang", "wong", "wu", "yana","zes")
n = length(x) ## 26
y = matrix(x[sample(26)],ncol=2,nrow=13)
y


      [,1]        [,2]          
 [1,] "anderson"  "mejgaard"    
 [2,] "wang"      "li"          
 [3,] "snyder"    "brett-esborn"
 [4,] "enos"      "jiang"       
 [5,] "wong"      "nesbitt"     
 [6,] "shanata"   "rundel"      
 [7,] "zes"       "sun"        
 [8,] "pfeiffe"   "chen"        
 [9,] "parfenov"  "ancheta"     
[10,] "tsai"       "cetinkaya"   
[11,] "yana"      "nguyen"      
[12,] "wu"        "diez"        
[13,] "firestine" "cassel"  



### 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 = 103
plot(c(0,1),c(0,1),type="n")
x1 = rep(0,n)
y1 = rep(0,n)
for(i in 1:n){
z1 = locator(1)
x1[i] = z1$x
y1[i] = z1$y
}


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

##### 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)
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(2,1))
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.


######## Pocket plot
## e = unit vector pointing north
y = matrix(rep(0,(m1*m1)),ncol=m1)
for(i in c(1:99)){
cat(i,"\n")
for(j in c((i+1):100)){
y[i,j] = mean(sqrt(abs(x3[i,] - x3[j,])))
y[j,i] = y[i,j]
}
}


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


p = matrix(rep(0,(m1*m1)),ncol=m1)
for(i in c(1:99)){
cat(i,"\n")
for(j in c((i+1):100)){
k = j-i
p[i,j] = y[i,j] - z[k]
p[j,i] = p[i,j]
}
}

p1 = p[5,c(1:4,6:m1)]
for(i in seq(15,95,by=10)){
p1 = cbind(p1,p[i,c(1:(i-1),(i+1):m1)])
}



plot(c(1,10),c(-.5,2.5),xlab="Row number",ylab="P_jk",type="n",xaxt="n")
boxplot(p1[,1],p1[,2],p1[,3],p1[,4],p1[,5],p1[,6],p1[,7],p1[,8],p1[,9],p1[,10],
names=as.character(seq(5,95,by=10)),add=T)





## First, suppose regular spatial data on [0,50]x[0,100]:

x = seq(0,50,by=2.5)
y = seq(0,100,by=5)
z = matrix(rnorm(21*21)+.05*c(1:(21*21)),ncol=21)
z[1,2] = 6.5
data1 = list(x=x,y=y,z=z)
par(mfrow=c(3,3))
image(data1)
#persp(data1)
#contour(data1)
data2 = data.frame(x=rep(x,times=rep(21,21)),y=rep(y,21),z=c(t(z)))
## data2 has 3 long vectors: the x-coordinates, the y's, and the z's.
my.ls = surf.ls(2,data2)
## my.ls is the least squares fit of a quadratic surface to the data
my.matrix.ls = trmat(my.ls,0,50,0,100,21)
image(my.matrix.ls)
my.gls = surf.gls(2,expcov,data2,d=38)
## my.gls is the fit of a quadratic surface to the data by GLS,
## estimating the variogram exponentially with range d=38
my.matrix.gls = trmat(my.gls,0,50,0,100,21)
image(my.matrix.gls)
my.krig = prmat(my.gls,0,50,0,100,21)
image(my.krig)
my.ses = semat(my.gls,0,50,0,100,21)
image(my.ses)

correlogram(my.ls,38)
c1 = correlogram(my.gls,38)
lines(c1$x,sphercov(c1$x,38))
variogram(my.ls,25)  ## fits variogram to LS residuals
variogram(my.gls,25) ## fits variogram to GLS residuals

## or, if you want to set the parameters of the plot yourself:
par(mfrow=c(1,1))
b5 = variogram(my.gls,25,plot=F) ## fits variogram to GLS residuals
plot(c(0,110),c(-.5,3),type="n",xlab="h",ylab="gamma(h)")
lines(b5)
title(main="variogram of GLS residuals")