#########################################################
########################################################
##  Clark Evans Test and Chi-Square Test
##########################################################
#########################################################


library(splancs)
library(spatstat)


#########################################################
########################################################
##  Clark Evans Test and Chi-Square Test
##########################################################
#########################################################


#made up examples

par(mfrow=c(1,1))

#100 random points on square (CSR)
x<-runif(100,1,100)
y<-runif(100,1,100)
CSR.ppp<-ppp(x,y,xrange=range(x),yrange=range(y))
clarkevans.test(CSR.ppp)
#actual z statistic
(clarkevans.test(CSR.ppp)$stat-1)/sqrt((4-pi)/CSR.ppp$n*pi)

#quatrat count test
CHSQ=quadrat.test(CSR.ppp, nx = 4, ny = 4)
CHSQ
plot(CSR.ppp, cols="red", pch=19)
plot(CHSQ, add = TRUE, cex=1)



#regular grid
x10<-rep(c(1:10),10)
y10<-sort(x10)
par(pty="s")
plot(x10,y10,pch=19,col='red',xlab="",ylab="")
grid.ppp<-ppp(x10,y10,xrange=range(x10),yrange=range(y10))
clarkevans.test(grid.ppp)
(clarkevans.test(grid.ppp)$stat-1)/sqrt((4-pi)/grid.ppp$n*pi)

#quatrat count test
CHSQ=quadrat.test(grid.ppp, nx = 4, ny = 4)
CHSQ
plot(grid.ppp, cols="red", pch=19)
plot(CHSQ, add = TRUE, cex=1)




#10 clusters of 10 points
xa<-runif(100)
ya<-runif(100)
movex<-runif(10,1,10)
movex<-sort(rep(movex,10))+xa
movey<-runif(10,1,10)
movey<-rep(movey,each=10)+ya
Clust.ppp<-ppp(movex, movey,xrange=range(movex),yrange=range(movey))
clarkevans.test(Clust.ppp)
(clarkevans.test(Clust.ppp)$stat-1)/sqrt((4-pi)/Clust.ppp$n*pi)

#quatrat count test
CHSQ=quadrat.test(Clust.ppp, nx = 4, ny = 4)
CHSQ
plot(Clust.ppp, cols="red", pch=19)
plot(CHSQ, add = TRUE, cex=1)



######################################################################
#  Clark Evans for real data
######################################################################





##################################################################
####  NOW SOME REAL DATA
##################################################################

###################################
#UGANDA VOLCANO DATA
###################################



 ugpoly=read.csv("http://reuningscherer.net/fes781/data/ugpoly.csv")
 ugpts=read.csv("http://reuningscherer.net/fes781/data/ugdata.csv")
 ug.pts=as.points(ugpts[,2],ugpts[,3])
 ug.poly=as.points(ugpoly[,2],ugpoly[,3])
 ugppp=ppp(ug.pts[,1],ug.pts[,2],xrange=range(ug.poly[,1]),yrange=range(ug.poly[,2]))


clarkevans.test(ugppp)
(clarkevans.test(ugppp)$stat-1)/sqrt((4-pi)/ugppp$n*pi)

#quatrat count test
CHSQ=quadrat.test(ugppp, nx = 4, ny = 4)
CHSQ
plot(ugppp, cols="red", pch=19)
plot(CHSQ, add = TRUE, cex=1)

###################################
#CARDIFF JUVENLE DELINQUENCY DATA
###################################

 cardpoly=read.csv("http://reuningscherer.net/fes781/data/cardpoly.csv")
 cardpts=read.csv("http://reuningscherer.net/fes781/data/cardpts.csv")
 card.pts=as.points(cardpts[,2],cardpts[,3])
 card.poly=as.points(cardpoly[,2],cardpoly[,3])
 #a plot of the data with a polygon around it.

cardppp=ppp(card.pts[,1],card.pts[,2],xrange=range(card.poly[,1]),yrange=range(card.poly[,2]))


plot(card.pts[,1],card.pts[,2],pch=20,cex=1.2,col="red",,main="Cardiff Delinquents",xlab="", ylab=""
,)
polymap(card.poly,add=T)

clarkevans.test(cardppp)




###################################
#Tim's MAINE tree data
###################################

# 1. import population and store population size and total stem weight

maine <- read.fwf("http://reuningscherer.net/fes781/data/maine.txt", skip=44,
     col.names=c("treecode","spp","crown", "dbh", "ht", "stemwt", "x", "y"),
     widths=c(7,3,3,8,8,8,8,8), strip.white=T)
maine$ba <- pi*(maine$dbh/2/12)^2

  maine$dbh.cm <- maine$dbh * 2.54;  maine$height.m <- maine$ht * .3048
  maine$ba.m2 <- maine$ba * .3048**2

  maine$tree.label <- 1:nrow(maine); maine$tree.count <- rep(1, nrow(maine))
  maine$stemwt.kg <- maine$stemwt / 2.2046
attach(maine)


N <- nrow(maine)       # store population size
N                      # display population size

spruce <- maine[spp == 1,]  ; spruce$Nsp <- nrow(spruce)
hemlock <- maine[spp == 2,]  ; hemlock$Nhe <- nrow(hemlock)
redmaple <- maine[spp == 3,]  ; redmaple$Nrm <- nrow(redmaple)
cedar <- maine[spp == 4,]  ; cedar$Nce <- nrow(cedar)
balsamfir <- maine[spp == 5,]  ; balsamfir$Nbf <- nrow(balsamfir)
paperbirch <- maine[spp == 6,]  ; paperbirch$Npb <- nrow(paperbirch)
whitepine <- maine[spp == 7,]  ; whitepine$Nwp <- nrow(whitepine)

paperbirch.ppp<-ppp(paperbirch$x,paperbirch$y,xrange=range(paperbirch$x),yrange=range(paperbirch$y))
whitepine.ppp<-ppp(whitepine$x,whitepine$y,xrange=range(whitepine$x),yrange=range(whitepine$y))



clarkevans.test(paperbirch.ppp)
clarkevans.test(whitepine.ppp)





######################################################################
#  Quadrat Count Chi-Square Test
######################################################################

######################################################################
#Example of non-rectangual quadrat counts - Baddeley 2010, section 13.2
#    Quadrats created using a potentially related covariate
######################################################################

#attach dataset for tropical rain forest
data(bei)

#Z is the information on slope (gradient)
Z <- bei.extra$grad

#Get quantiles - make 8 breaks
b <- quantile(Z, probs = (0:8)/8)

#Divide slope into 8 groups
Zcut <- cut(Z, breaks = b, labels = 1:8)

#V is a tesselation (i.e. divide into regions) based on slope

V <- tess(image = Zcut)
plot(V)
plot(bei, add = TRUE, pch = "+")

#Get quadrat counts based on 8 slope regions

qb <- quadratcount(bei, tess = V)
qb
#add numbers to the plot
plot(qb)

#quatrat count test using 4 x 2 grid
CHSQ=quadrat.test(bei, nx = 4, ny = 2)
CHSQ
plot(bei, pch=".")
plot(CHSQ, add = TRUE, cex=1)


#Do Chi-square quadrat count test using tesselation
CHSQ2=quadrat.test(bei, tess=V)
CHSQ2
plot(V)
plot(CHSQ2, add = TRUE, cex=1)