library(nlme)
library(spatstat)
library(spatial)
library(HistData)
library(rgeos)
library(maptools)
library(sp)
data(Snow.deaths)
data(Snow.pumps)
data(Snow.streets)
data(Snow.polygons)

# draw a rough approximation to Snow's map and data

# define some funtions to make the pieces re-usable
Sdeaths <- function(col="red", pch=15, cex=0.6) {
  # make sure that the plot limits include all the other stuff
  plot(Snow.deaths[,c("x","y")], col=col, pch=pch, cex=cex, 
       xlab="", ylab="", xlim=c(3,20), ylim=c(3,20),
       main="Snow's Cholera Map of London")
}
# function to plot and label the pump locations
Spumps <- function(col="blue", pch=17, cex=1.5)  {
  points(Snow.pumps[,c("x","y")], col=col, pch=pch, cex=cex)
  text(Snow.pumps[,c("x","y")], labels=Snow.pumps$label, pos=1, cex=0.8)
}

# function to draw the streets 
Sstreets <- function(col="gray") {
  slist <- split(Snow.streets[,c("x","y")],as.factor(Snow.streets[,"street"]))
  invisible(lapply(slist, lines, col=col))
}

# draw a scale showing distance in meters in upper left
mapscale <- function(xs=3.5, ys=19.7) {
  scale <- matrix(c(0,0, 4,0, NA, NA), nrow=3, ncol=2, byrow=TRUE)
  colnames(scale)<- c("x","y")
  # tick marks
  scale <- rbind(scale, expand.grid(y=c(-.1, .1, NA), x=0:4)[,2:1])
  lines(xs+scale[,1], ys+scale[,2])
  # value and axis labels
  stext <- matrix(c(0,0, 2,0, 4,0, 4, 0.1), nrow=4, ncol=2, byrow=TRUE)
  text(xs+stext[,1], ys+stext[,2], labels=c("0", "2", "4", "100 m."), pos=c(1,1,1,4), cex=0.8)
}

# draw the map with the pieces
Sdeaths()
Spumps()
Sstreets()
mapscale()


# draw the Thiessen polygon boundaries
starts <- which(Snow.polygons$start==0)
for(i in 1:length(starts)) {
  this <- starts[i]:(starts[i]+1)
  lines(Snow.polygons[this,2:3], col="blue", lwd=2, lty=2)
}


## overlay bivariate kernel density contours of deaths
Sdeaths()
Spumps()
Sstreets()
mapscale()
require(KernSmooth)
kde2d <- bkde2D(Snow.deaths[,2:3], bandwidth=c(0.5,0.5))
contour(x=kde2d$x1, y=kde2d$x2,z=kde2d$fhat, add=TRUE)


#######################################
#CSR Test
#########################################

# streets
slist <- split(Snow.streets[,c("x","y")],as.factor(Snow.streets[,"street"]))
Ll1 <- lapply(slist,Line)
Lsl1 <- Lines(Ll1,"Street")
Snow.streets.sp <- SpatialLines(list(Lsl1))
plot(Snow.streets.sp, col="gray")
title(main="Snow's Cholera Map of London (sp)")

# deaths
Snow.deaths.sp = SpatialPoints(Snow.deaths[,c("x","y")])
p<-as(Snow.deaths.sp,"ppp")
qt<-quadrat.test(p,nx=5,ny=5)
plot(Snow.deaths.sp, add=TRUE, col ='red', pch=15, cex=0.6)
# pumps
spp <- SpatialPoints(Snow.pumps[,c("x","y")])
Snow.pumps.sp <- SpatialPointsDataFrame(spp,Snow.pumps[,c("x","y")])
plot(Snow.pumps.sp, add=TRUE, col='blue', pch=17, cex=1.5)
text(Snow.pumps[,c("x","y")], labels=Snow.pumps$label, pos=1, cex=0.8)
plot(qt,add=T,cex=.3)
#test for significance. Null Hypothesis is complete spatial randomness
qt