library(splancs)
try(data(package = "splancs") )

nests<-read.table("nests.txt",header=T)

dim(nests)
names(nests)
nests[1:10,]

#nest.pts<-as.points(nests)
nest.pts<-as.points(nests[,1],nests[,2])

npts(nest.pts)

summary(nest.pts)

#Plotting polygons and points, other utilities:
#create a dot-plot of the nests data set
pointmap(nest.pts,pch="+")

#plot the data within a bounding polygon with corners
#(0,0),(30,0), (30,20) and (0,20)
nest.poly<-as.points(c(0,30,30,0),c(0,0,20,20))
#or, you can do this interactively. Type
nest.poly<-getpoly(quiet=TRUE)
#then click on the graph to define the "corners" of the polygon

#then draw
polymap(nest.poly)
pointmap(nest.pts,add=T,pch=20)


#estimating and mapping intensity
#splancs provides: kernel2d(),mse2d()
#use mse2d to decide on the degree of smoothing
#choose h so that mse is smallest
nest.mse<-mse2d(nest.pts,nest.poly,nsmse=20,range=10)
plot(nest.mse$h,nest.mse$mse,pch=20,ylab="MSE",xlab="h")

nest.lambda<-kernel2d(nest.pts,nest.poly,h0=3.5,nx=60,ny=40)
polymap(nest.poly)
image(nest.lambda,add=T,col=topo.colors(12))
pointmap(nest.pts,add=T,pch=20)
contour(nest.lambda,add=T)

#Second order properties
#I need to define the set of distances
nest.all.dists<-nndistG(nest.pts)$dist
#could plot a boxplot of nearest neighbor distances (to see if there are any "excess" small/large distances)
boxplot(nndistG(nest.pts)$dist)


#Estimate Ghat
#First we need to extract the set of distances
nest.all.dists<-nndistG(nest.pts)$dist
#and get rid of duplicates
nest.su<-sort(unique(nest.all.dists))
#Finally, compute Ghat 
nest.ghat<-Ghat(nest.pts,nest.su)

#plot Ghat vs. Distance
plot(nest.su,nest.ghat,xlab="Distance",ylab="Ghat",pch=20,cex=.75)
title('Ghat for Duck Nests')


#Estimate K function
#first, define the array of distances for which 
#we will compute khat
nest.k.dists<-seq(from=0.25,to=10,by=0.25)
nest.k<-khat(nest.pts,nest.poly,nest.k.dists)
plot(nest.k.dists,nest.k,type='l')
title("Ripley's K for Duck Nests")

#plot the theoretical K function for Homogeneour Poisson
lines(nest.k.dists,pi*(nest.k.dists)^2,col=2)

#compute simulated envelopes using the Kenv.csr function
nests.K.env<-Kenv.csr(n.nests,nest.poly,100,nest.k.dists)

#plot them on the same graph as Khat
plot(nest.k.dists,nest.k,type='l', ylab="Estimated L")
lines(nest.k.dists,nests.K.env$up,col=2,lty=2)
lines(nest.k.dists,nests.K.env$low,col=3,lty=2)


#estimate and plot L function
nest.l<-sqrt(nest.k/pi)-nest.k.dists
#plot estimated L function
plot(nest.k.dists,nest.l,type="l",ylim=c(-0.3,0.7))
#plot simulated envelopes for L
lines(nest.k.dists,sqrt(nests.K.env$up/pi)-nest.k.dists,col=2,lty=2)
lines(nest.k.dists,sqrt(nests.K.env$low/pi)-nest.k.dists,col=3,lty=2)


###########################################################
#in case you would like to get a bit more adventurous and
#simulate some spatial point patterns, here's how you'd do it:
nests.rand<-csr(nest.poly,n.nests)
pointmap(nests.rand)

#constructing simulated envelopes

#construct a function that simulates a CSR process
#for the given region, consisting of the SAME number
#of points as the original data set

 ghat.env<-function(data.p,data.s,polyg,nsim)
 {
  nloc<-length(data.p[,1])
  gall<-matrix(0,nrow=length(data.s),ncol=nsim)
  for (i in 1:nsim){
   data.sim<-csr(polyg,nloc)
   gall[,i]<-Ghat(data.sim,data.s)
  }
  return(gall)
 }

#Call this function for your data set, for a number of sim.

ghat.all<-ghat.env(nest.pts,nest.su,nest.poly,99)

#calculate max, min and mean for each distance
Up<-apply(ghat.all,1,max)
Down<-apply(ghat.all,1,min)
mn<-apply(ghat.all,1,mean)

#plot Ghat vs. dist
plot(nest.su,nest.ghat,type='l',lwd=2,xlab="Distance",ylab="Ghat")
points(nest.su,Up,type='l',col=2,lwd=2,lty=2)
points(nest.su,Down,type='l',col=3,lwd=2,lty=2)
title("Simulated envelopes for Ghat")


#plot Ghat vs. mean, Up vs. mean, and Down vs. mean
plot(mn,nest.ghat,type='l',lwd=2,xlab="Mean",ylab="Ghat")
points(mn,Up,type='l',col=2,lwd=2,lty=2)
points(mn,Down,type='l',col=3,lwd=2,lty=2)
title("Simulated envelopes for Ghat")