#=========================================================================
#
# File:		CENN.R 	
# Description:	Imports data, creates a stem map, and conducts first
#		order Nearest Neighbor Distance Analyses, including nnd 
#		summary statistics and Clark and Evans' R
# Written by: 	Andrew J. Sánchez Meador
# Date: 	November 25, 2006
# Files needed: s4acorners.csv and s4ac.csv
#
# Modified by JD Bakker on 080304.
#
#=========================================================================

#=========================================================================
#
# Data Set-Up
# 
# Reads the csv files containing data, defines the area of interest,
# subsets the data appropriately, and creates spatial point pattern (spp) 
# data sets
#
#=========================================================================

	# Read a csv data file in table format and create a data frame from
	# it for analysis purposes
	corners <- read.csv(file.choose(), header = TRUE, sep = ",")
	contemp.data <- read.csv(file.choose(), header = TRUE, sep = ",")

	# Create a polygon defining the boundary of the point pattern, 
	# arranged as two vectors (x's and y's)
	s4a.poly <- as.points(corners$X, corners$Y)

	# Subset the data frame as needed - two populations (presettlement
	# (os) and postsettlement (us)) based on population, which reflects
	# tree age and diameter
	contemp.us<-subset(contemp.data, Pop==0)
	contemp.os<- subset(contemp.data, Pop==1)

	# create datasets defining the points of interest, arranged as two 
	# vectors (x's and y's)
	contemp.us.spp <- as.points(contemp.us$X, contemp.us$Y)
	contemp.os.spp <- as.points(contemp.os$X, contemp.os$Y)
	contemp.spp <- as.points(contemp.data$X, contemp.data$Y)
	colnames(contemp.spp) <- c("X", "Y")
	contemp.spp <- as.data.frame(contemp.spp)

#=========================================================================
#
# Create Stem Map 
#
# Display a stem map for the contemporary live trees
#
#=========================================================================

	# Creates a figure, showing the plot boundary, the presettlement and 
	# the postsettlement trees (scaled to size based on diameter)
	library(splancs)
	polymap(s4a.poly, lty=2, xlab="X Coordinate (m)", ylab="Y Coordinate (m)")
	symbols(contemp.os.spp, circles=contemp.os$Dia/25.4, inches=FALSE, fg=2, bg=2, add=TRUE)
	symbols(contemp.us.spp, squares=contemp.us$Dia/25.4, inches=FALSE, fg=1, bg=1, add=TRUE)
	# Detach the SPLANCS library
	detach("package:splancs")

#=========================================================================
#
# Calculate NND Statistic and C&E Index
#
# Calculate statistics on nearest neighbor distances (nnd), load the
# Clark and Evans' R (CENN; Clark and Evans 1954) function, and use it to
# analyze the point pattern data.
#
#=========================================================================

	# Load the SPATSTAT library
	library(spatstat)

	# Calculate the individual first nearest neighbor distance values
	contemp.nnd<-nndist(contemp.spp)

	# Load the generalize CENN function from Reich and Davis's online 
	# spatial package (Reich & Davis 1998) 
	cenn<-function(data, poly, n, graph = T, sim = F)
	{
		# Load the SPLANCS library
		library(splancs)

		# calculating estimate of density (D) and polygon area (area)
		num.ind <- length(data$X)
		area <- areapl(poly)
		D <- num.ind/area

		#calculating the eXpected mean distance under CSR
		exp <- 1/(sqrt(D) * 2)

		#locating n random objects
		index <- sample(num.ind, n, replace = F)
		point.x <- data$X[index]
		point.y <- data$Y[index]
		coordin <- list(a = point.x, b = point.y)

		# making border adjustments to correct for edge effects
		plusx <- data$X + poly[3,1] - poly[1,1]
		plusy <- data$Y + poly[3,2] - poly[1,2]
		minusx <- data$X - (poly[3,1] - poly[1,1])
		minusy <- data$Y - (poly[3,2] - poly[1,2])
		new <- list()
		new$X <- rep(c(minusx, data$X, plusx), 3)
		new$Y <- c(rep(plusy, 3), rep(data$Y, 3), rep(minusy, 3))
		index <- length(new$X)

		#finding the minimum nearest neighbor distance
		dist <- list()
		neighbor <- list()
		place <- numeric()
		w <- numeric()
		n.x <- numeric()
		n.y <- numeric()
		for(i in 1:n) {
			dist[[i]] <- (data$X - coordin$a[[i]])^2 + (data$Y - coordin$b[[i]])^2
			test <- sort(dist[[i]])
			w[i] <- test[2]

			#Locating nearest neighbor
			place[i] <- (1:num.ind)[dist[[i]] == w[i]]
			n.x[i] <- data$X[place[i]]
			n.y[i] <- data$Y[place[i]]
		}
		neighbor <- list(x = n.x, y = n.y)

		#finding the index of non-randomness: R 
		w <- sqrt(w)
		avg.w <- mean(w)
		R <- avg.w/exp
		z <- abs(R - 1)/sqrt((4 - pi)/(pi * n))
		pvalue <- 2 * (1 - pnorm(z, 0, 1))
		if(!sim) {
			cat("\n Number of trees = ", num.ind)
			cat("\n Density (Number of trees / plot area) = ", round(D, 6))
			cat("\n Expected mean distance under CSR = ", round(exp, 6))
			cat("\n Average distance to nearest neighbor = ", round(avg.w, 6))
			cat("\n Clark & Evans Nearest Neighbor Index = ", round(R, 6))
			cat("\n P-value for testing null hypothesis of CSR (index=1) = ",
				round(pvalue, 6), "\n")
			#plotting the scene
			if(graph) {
				plot(data$X, data$Y, pch = 0)
				title("Clark & Evans Nearest Neighbor Index")
				points(coordin$a, coordin$b, type = "p", pch = 15)
				segments(coordin$a, coordin$b, neighbor$x, neighbor$y)
			}
		}
		invisible(list(index = R, pvalue = pvalue))

		# Detach the SPLANCS library
		detach("package:splancs")
	}

	# Use the CENN function to calculate R and test against CSR
	cenn(contemp.spp, s4a.poly, nrow(contemp.spp), graph = F, sim = F)
	# Turn graph argument on or off as desired

	# Compute the mean, sd, median, min, and max nnd
	# Plot the contemporary data
		mean(contemp.nnd)
		sd(contemp.nnd)
		median(contemp.nnd)
		min(contemp.nnd)
		max(contemp.nnd)

	# Detach the SPATSTAT library
	detach("package:spatstat")