#########################################################
########################################################
##  G and F Functions.R.txt
##########################################################
#########################################################

library(splancs)
library(spatstat)

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

####Load K and L and G

source("http://reuningscherer.net/fes781/RScripts/SpatialPPFuncs.R.txt")


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

#plot two species
par(mfrow=c(3,2))
par(pty="s")
plot(paperbirch.ppp$x,paperbirch.ppp$y,pch=19,col='red',xlab="X",ylab="Y",main="Paper Birch")
plot(whitepine.ppp$x,whitepine.ppp$y,pch=19,col='red',xlab="X",ylab="Y",main="White Pine")

#plot L function for these species
L(paperbirch.env,paperbirch.K,275,30,"Paper Birch")
L(whitepine.env,whitepine.K,275,30,"White Pine")

#plot kernel smoothed plots
par(pty="s")
dens3paper<-density.ppp(paperbirch.ppp,sigma=70,dimyx=c(20,20))
contour(dens3paper,xlab="u",ylab="v",main="Paper Birch Sigma=70",lwd=1.5,col='blue')
dens3white<-density.ppp(whitepine.ppp,sigma=70,dimyx=c(20,20))
contour(dens3white,xlab="u",ylab="v",main="White Pine Sigma=70",lwd=1.5,col='blue')



############plots for made up data

par(mfrow=c(3,2),pty="s")

#100 random points on square (CSR)
x<-runif(100,1,100)
y<-runif(100,1,100)
plot(x,y,pch=19,col='red',xlab="",ylab="",xlim=c(0,100),ylim=c(0,100),pty="s")
CSR.ppp<-ppp(x,y,xrange=range(x),yrange=range(y))
CSR.env<-envelope(CSR.ppp,fun=Gest)
CSR.G<-Gest(CSR.ppp)
G(CSR.env,CSR.G,"Randomly Generated CSR")

#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))
grid.env<-envelope(grid.ppp,fun=Gest)
grid.G<-Gest(grid.ppp)
G(grid.env,grid.G,"Regular Grid")

#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
par(pty="s")
plot(movex,movey,pch=19,col='red',xlab="",ylab="")

Clust.ppp<-ppp(movex, movey,xrange=range(movex),yrange=range(movey))
Clust.env<-envelope(Clust.ppp,fun=Gest)
Clust.G<-Gest(Clust.ppp)
G(Clust.env,Clust.G,"Clustered Process")




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

par(mfrow=c(2,2))

plot(ug.pts[,1],ug.pts[,2],pch=20,cex=1.2,col="red",,main="Uganda Volcanos",xlab="", ylab=""
, xlim=c(0,3500))
polymap(ug.poly,add=T)

ug.env<-envelope(ugppp,fun=Gest)
ug.G<-Gest(ugppp)
G(ug.env,ug.G,"Uganda Volcano Data")


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


card.env<-envelope(cardppp,fun=Gest)
card.G<-Gest(cardppp)
G(card.env,card.G,"Cardiff Delinquents")

#####################################################
####  Tim's Maine Tree Data
#####################################################



# Now actually compute the G function

par(mfrow=c(3,2))
par(pty="s")
plot(maine$x,maine$y,pch=19,col='red',cex=0.2,xlab="",ylab="",main="All Trees")
maineppp=ppp(maine$x,maine$y,xrange=range(maine$x),yrange=range(maine$y))
maine.env<-envelope(maineppp,fun=Gest)
maine.G<-Gest(maineppp)
G(maine.env,maine.G,"All Trees")

par(pty="s")
plot(hemlock$x,hemlock$y,pch=19,col='red',cex=0.2,xlab="",ylab="",main="Hemlock")
hemlockppp=ppp(hemlock$x,hemlock$y,xrange=range(hemlock$x),yrange=range(hemlock$y))
hemlock.env<-envelope(hemlockppp,fun=Gest)
hemlock.G<-Gest(hemlockppp)
G(hemlock.env,hemlock.G,"Hemlock")

par(pty="s")
plot(redmaple$x,redmaple$y,pch=19,col='red',cex=0.2,xlab="",ylab="",main="Red Maple")
redmapleppp=ppp(redmaple$x,redmaple$y,xrange=range(redmaple$x),yrange=range(redmaple$y))
redmaple.env<-envelope(redmapleppp,fun=Gest)
redmaple.G<-Gest(redmapleppp)
G(redmaple.env,redmaple.G,"Red Maple")




par(mfrow=c(3,2))

#100 random points on square (CSR)
x<-runif(100,1,100)
y<-runif(100,1,100)
plot(x,y,pch=19,col='red',xlab="",ylab="",xlim=c(0,100),ylim=c(0,100))
CSR.ppp<-ppp(x,y,xrange=range(x),yrange=range(y))
CSR.env<-envelope(CSR.ppp,fun=Fest)
CSR.F<-Fest(CSR.ppp)
Ffunc(CSR.env,CSR.F,"Randomly Generated CSR")

#regular grid
x10<-rep(c(1:10),10)
y10<-sort(x10)
plot(x10,y10,pch=19,col='red',xlab="",ylab="")
grid.ppp<-ppp(x10,y10,xrange=range(x10),yrange=range(y10))
grid.env<-envelope(grid.ppp,fun=Fest)
grid.F<-Fest(grid.ppp)
Ffunc(grid.env,grid.F,"Regular Grid")

#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
plot(movex,movey,pch=19,col='red',xlab="",ylab="")

Clust.ppp<-ppp(movex, movey,xrange=range(movex),yrange=range(movey))
Clust.env<-envelope(Clust.ppp,fun=Fest)
Clust.F<-Fest(Clust.ppp)
Ffunc(Clust.env,Clust.F,"Clustered Process")










#made up examples

par(mfrow=c(3,2))

#100 random points on square (CSR)
x<-runif(100,1,100)
y<-runif(100,1,100)
par(pty="s")
plot(x,y,pch=19,col='red',xlab="",ylab="",xlim=c(0,100),ylim=c(0,100))
CSR.ppp<-ppp(x,y,xrange=range(x),yrange=range(y))
CSR.env<-envelope(CSR.ppp,fun=Jest)
CSR.J<-Jest(CSR.ppp)
Jfunc(CSR.env,CSR.J,"Randomly Generated CSR",xlims=c(0,8),ylims=c(0,3),legend=c(0,3))



#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))
grid.env<-envelope(grid.ppp,fun=Jest)
grid.J<-Jest(grid.ppp)
Jfunc(grid.env,grid.J,"Regular Grid",xlims=c(0,0.6),ylims=c(0,15),legend=c(0,15))



#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
par(pty="s")
plot(movex,movey,pch=19,col='red',xlab="",ylab="")

Clust.ppp<-ppp(movex, movey,xrange=range(movex),yrange=range(movey))
Clust.env<-envelope(Clust.ppp,fun=Jest)
Clust.J<-Jest(Clust.ppp)
Jfunc(Clust.env,Clust.J,"Clustered Process",xlims=c(0,0.6),ylims=c(0,2),legend=c(0,2))