###########################################################
#
# Continuous Education Seminar :
#
# Introduction of geospatial data visualization 
#                 and geographically weighted regression (GWR)
#
#                                               8/16/2012
#   Frank Fan 
#
###########################################################

rm(list=ls())
library(maps)
library(maptools)

###########################################################
#
# Setting Color Scales
#
###########################################################

color_scales <- colorRampPalette(c("orange","red","green"))
col.lists <- color_scales(21)
plot(1:21,rep(1,21),col=col.lists,pch=15,cex=4,xaxt="n",
      yaxt="n",xlab="",ylab="",axes=FALSE)
for(i in 1:21){
text(i,0.8,i)
}

######################################
### Showing Examples from here !!! ###
######################################

#############################################################
#
# Maps of United States including Alaska (AK) & Hawaii (HI)
#
#                                          Frank 2012/03/09
# 
#############################################################

# Remember to set working directory !!

load("coords of AK and HI.RData")

# str(ak)
# str(hi)

# pdf("US Map.pdf",height=8,width=10)
# jpeg("US Map.jpeg",height=800,width=1000,quality=100)

layout(rbind(c(1,1),c(2,3)),heights=c(5,2))

# layout 1
par(mar=c(0,0,3,0))
map("state")
mtext("United States",line=1)

# layout 2
par(mar=c(0,2,0,8))
plot(NULL,xlim=c(-179,-131.5),ylim=ak@bbox[2,],
     xaxt="n",yaxt="n",bty="n",xlab="",ylab="")
for(i in 1:length(ak@polygons[[1]]@Polygons)){
  coords <- ak@polygons[[1]]@Polygons[[i]]@coords
  polygon(coords[,1],coords[,2])
}
mtext("Alaska")

# layout 3
par(mar=c(5,10,5,12))
plot(NULL,xlim=c(-161.5,-154.7579),ylim=c(18.86546,22.5),
     xaxt="n",yaxt="n",bty="n",xlab="",ylab="")
for(i in 1:length(hi@polygons[[1]]@Polygons)){
  coords <- hi@polygons[[1]]@Polygons[[i]]@coords
  polygon(coords[,1],coords[,2])
}
mtext("Hawaii",line=2)

# dev.off()

###########################################################
#
# Plot 
#
###########################################################

load("map.RData")

######################################################
# Zip code vs. Latitude & Longitude Dictionary 
######################################################

TN.zip.map@polygons[964]

# Get "center" coordinates of the zip coed
labelpos <- data.frame(do.call(rbind, lapply(TN.zip.map@polygons, function(x) x@labpt)))

# Get zip code "in order"
zip.name <- TN.zip.map@data

# Merge as a dataframe 
zip.dic <- data.frame(labelpos,zip.name$NAME,block=1:964)
names(zip.dic) <- c("longitude.x","latitude.y","Zip.Code","Block")

# str(zip.dic)
# head(zip.dic)

# Select the zip code area of interests 
data <- merge(pi.data,zip.dic,all.x=TRUE)

# str(data)
# head(data)

#######################################################
# Data Correction
#######################################################

data$Zip.Code[which(is.na(data$Block))]

# http://maps.huge.info/zip.htm
# 37132 belongs to 37130 ; 37240 belongs to 37203

data <- data[-which(is.na(data$Block)),]

#######################################################
# Color Scale Function
#######################################################

cs <- function(x){
	cutoff <- unname(quantile(x,probs=seq(0,1,0.2)))
	y <-  ifelse( x < cutoff[2], 1,
		ifelse( x < cutoff[3], 2,
		ifelse( x < cutoff[4], 3,
		ifelse( x < cutoff[5], 4, 5 ))))
	return(y)
}

#######################################################
# Assign Color Scale 
#######################################################

mat <- data.frame(apply(data[,4:5],2,cs))
data <- cbind(data,mat)

# head(data)

#######################################################
# Color Palette
#######################################################

##################### light to dark !!!! ###################

#############
# mortality #
#############

mort.scl <- colorRampPalette(c("#9999FF","#3333FF"))
mort.col.lists <- mort.scl(5)
data$mort.col.sets <- sapply(mat$Female.Mortality.Rate.per.100K,function(x) return(mort.col.lists[x]))

##########
# Income #
##########

income.scl <- colorRampPalette(c("#99FF99","#009900"))
income.col.lists <- income.scl(5)
data$income.col.sets <- sapply(mat$Median.Houshold.Income,function(x) return(income.col.lists[x]))

#########################################################
# Color of the Zip Code label
#########################################################

zip.col <- c(rgb(254,196,79,max=255),rgb(255,20,147,max=255))

#########################################################
# Plot !!!
#########################################################

layout(1)

#############
# mortality #
#############

par(mar=c(2,1,2,1))
plot(NULL,xlim=c(-87.8,-85.8),ylim=c(35.4,36.7),xaxt="n",yaxt="n",bty="n",xlab="",ylab="")

for(i in data$Block){
ind <- which(data$Block==i)
coords <- TN.zip.map@polygons[[i]]@Polygons[[1]]@coords
polygon(coords[,1],coords[,2],col=data$mort.col.sets[ind])
}

for( j in 1:nrow(data)){
text(data$longitude.x[j],data$latitude.y[j],data$Zip.Code[j],cex=.5,col=zip.col[1])
}

mtext("Greater Nashville Affiliate",1,0,font=2,cex=1.2)
mtext("Female Mortality Rate Per 100K Female Population",3,0,cex=1.5,font=3)

legend(-86.2,35.7,legend=c("  6.50 - 20.05","20.06 - 23.12","23.13 - 26.63","26.64 - 29.56","29.57 - 47.76")
		,col=mort.col.lists, pch=15,pt.cex=1.5, title="Mortality Rate",bty="n")

##########
# Income #
##########

par(mar=c(2,1,2,1))
plot(NULL,xlim=c(-87.8,-85.8),ylim=c(35.4,36.7),xaxt="n",yaxt="n",bty="n",xlab="",ylab="")

for(i in data$Block){
ind <- which(data$Block==i)
coords <- TN.zip.map@polygons[[i]]@Polygons[[1]]@coords
polygon(coords[,1],coords[,2],col=data$income.col.sets[ind])
}

for( j in 1:nrow(data)){
text(data$longitude.x[j],data$latitude.y[j],data$Zip.Code[j],cex=.5,col=zip.col[2])
}

mtext("Greater Nashville Affiliate",1,0,font=2,cex=1.2)
mtext("Median Household Income",3,0,cex=1.5,font=3)

legend(-86.2,35.7,legend=c("13,523 -   43,102","43,103 -   48,770","48,771 -   53,982","53,983 -   61,586","61,587 - 108,498")
		,col=income.col.lists, pch=15,pt.cex=1.5, title="Median Income",bty="n")

#########################################################
# geographically weighted regression
#########################################################

library(spgwr)
data(columbus)

# Simple linear model 
col.lm <- lm(crime ~ income + housing, data=columbus)
summary(col.lm)

# GWR with Gauss
col.bw <- gwr.sel(crime ~ income + housing, data=columbus,
                  coords=cbind(columbus$x, columbus$y))

(col.gauss <- gwr(crime ~ income + housing, data=columbus,
                 coords=cbind(columbus$x, columbus$y), bandwidth=col.bw))


# Distribution of betas

d <- cbind(col.gauss$SDF$income,col.gauss$SDF$housing)

par(mar=c(3,4,2,2))
boxplot(d,xaxt="n",yaxt="n",pars=list(boxwex=0.3))
axis(1,at=1:2,label=c("Income","Housing"))
axis(2,at=seq(-4,2,.2),las=1)
abline(h=0,lty="4343",col="#7E7E7E")
mtext("Beta i",2,line=3)

# Information of original data and betas 

col.gauss$SDF