#=========================================================================
#
# File:		Ripley.R 	
# Description:	Imports data and conducts Ripley's Univariate K analyses
# on point pattern data
# 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 file containing data, defines the area of interest (via 
# its corners), subsets the data appropriately, and creates spatial point 
# pattern (spp) data sets
#
#=========================================================================
	
	# Load the SPLANCS library
	library(splancs)

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

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


#=========================================================================
#
# Univariate Analysis
# 
# Performs Ripley's K univariate analysis, diplaying a figure and writing 
# the results to a .csv file of your choosing.
#
#=========================================================================

	# L function to ease interpretation by stabilizing the variance
	l<-function(k, s)
	{
		sqrt(k/pi)-s
	}
	
	# Univariate Ripley's K envelope function (Same as above)
	lhat.env<-function(data.pp, data.s, poly, nsim)
	{
		nloc<-length(data.pp[,1])
		lall<-matrix(0,nrow=length(data.s),ncol=nsim)
		for (i in 1:nsim)
		{
			data.sim<-csr(poly,nloc)
			data.k<-khat(data.sim,poly,data.s)
			lall[,i]<-l(data.k,data.s)
		}
	return(lall)
	}

	# Univariate Ripley's K function
	univariate.khat<-function(data.pp, spoly, title, nsim, ylim)
	{
		# Set h equal to half the shortest side
		h<-0.5*ceiling(min(max(spoly[,1])-min(spoly[,1]),
		max(spoly[,2])-min(spoly[,2])))
		# Sequence from 0 to h by 2
		data.s<-seq(0, h, by=2)
		# Create the data frame to house the output
		results <- data.frame(data.s)
		# Calculate K  for the data provided
		data.k<-khat(data.pp,spoly,data.s)

		#Transfor the K data to L data
		data.l<-l(data.k,data.s)
		results <- cbind(results,data.l)

		#Simulate envelopes for K
		lhat.all<-lhat.env(data.pp,data.s,spoly,nsim)
		Up<-apply(lhat.all,1,max)
		Down<-apply(lhat.all,1,min)
		mn<-apply(lhat.all,1,mean)
		# Determine the 95% Conf. Envelopes
		l.q<-apply(lhat.all,1,quantile,c(0.025,0.975))
		results <- cbind(results,l.q[1,])
		results <- cbind(results,l.q[2,])

		# Plot the resulting analyses
		plot(data.s,data.l,type='l',lwd=2,xlab="Lag (m)",
		ylab="L(t)-t", xlim=c(0,h),ylim=c(-(ylim),ylim))
		points(data.s,l.q[1,],type='l',lty=2,lwd=1)
		points(data.s,l.q[2,],type='l',lty=2,lwd=1)
		abline(0,0)
		title(title)
		# Write the analysis results to a .csv file
		write.table(results,file=file.choose(new=T), append=FALSE,
		sep = ",", col.names = c("lag","lhat","lower95","upper95"))
	}

# Perform Univariate K analysis specifying data, polygon, title, number of 
# simulations and the limit for the positive y-axis
univariate.khat(contemp.spp, s4a.poly, "Contemporary", 99, 12)