# Mapping a few Alabama block-groups from 2000 vs 2010
# and showing how they overlap each other
#
# Jerzy Wieczorek
# June 2012
# civilstat.com


# I downloaded and unzipped the shapefiles from
# http://www.census.gov/cgi-bin/geo/shapefiles2010/main
#
# All are in Alabama, in Autauga and Mobile counties:
# Block-groups 010010201001-2: county 001, AUTAUGA
# Block-groups 010970011002-4: county 097, MOBILE
#
# I also use the Autauga county 2010 blocks (not just block-groups).


# Set working directory; load packages
mydir <- "C:/BlockGroupMaps"
setwd(mydir)
library(sp)
library(rgdal)
library(OpenStreetMap)



############################################################################
# Autauga County, Alabama
# Block groups 010010201001-2

# Read in the 2000 and 2010 Block-Groups, and 2010 Blocks
Autauga2000.BG <- readOGR(mydir, "tl_2010_01001_bg00")
plot(Autauga2000.BG)
Autauga2010.BG <- readOGR(mydir, "tl_2010_01001_bg10")
plot(Autauga2010.BG)
Autauga2010.B <- readOGR(mydir, "tl_2010_01001_tabblock10")
plot(Autauga2010.B)

# Subset just the Block-Groups and Block we need
Autauga2000.BG <- subset(Autauga2000.BG, Autauga2000.BG$BKGPIDFP00 == "010010201001")
Autauga2010.BG <- subset(Autauga2010.BG, Autauga2010.BG$GEOID10 == "010010201002")
Autauga2010.B <- subset(Autauga2010.B, Autauga2010.B$GEOID10 == "010010201001000")

# Plot quickly to see how it looks
plot(Autauga2000.BG, border='blue')
plot(Autauga2010.BG, border='red', add=T, lty=2)
plot(Autauga2010.B, border='black', lty=3, add=T)

# Grab the OpenStreetMap to put underneath
bb <- array(dim=c(2,2,3))
bb[,,1] = bbox(spTransform(Autauga2000.BG, CRS("+init=epsg:4326")))
bb[,,2] = bbox(spTransform(Autauga2010.BG, CRS("+init=epsg:4326")))
bb[,,3] = bbox(spTransform(Autauga2010.B, CRS("+init=epsg:4326")))
bb.Autauga <- matrix(c( min(bb[1,1,]),min(bb[2,1,]),max(bb[1,2,]),max(bb[2,2,]) ), 2,2)
map.Autauga = openmap(c(bb.Autauga[2, 2], bb.Autauga[1, 1]), c(bb.Autauga[2, 1], bb.Autauga[1, 2]),
    type = "osm")

# Plot nicely to a PDF or PNG image
#png("Autauga.png",480,720) # width&height in pixels
pdf("Autauga.pdf",4.7,7) # width&height in inches
plot(map.Autauga,removeMargin=TRUE)
plot(spTransform(Autauga2010.BG, osm()), add=T, lwd=2, border="red")
plot(spTransform(Autauga2000.BG, osm()), add=T, lwd=4, lty=2, border="blue")
plot(spTransform(Autauga2010.B, osm()), add=T, lwd=6, lty=3, border='black')
legend("bottomleft", legend=c("Census 2000 Block-Group 1","Census 2010 Block-Group 2","Census 2010 Block 1000"),
	lty=c(2,1,3), col=c("blue","red",'black'), lwd=c(4,2,6), bg="white")
dev.off()



############################################################################
# Mobile County, Alabama
# Block groups 010970011002-4

# Read in the 2000 and 2010 Block-Groups
Mobile2000.BG <- readOGR(mydir, "tl_2010_01097_bg00")
plot(Mobile2000.BG)
Mobile2010.BG <- readOGR(mydir, "tl_2010_01097_bg10")
plot(Mobile2010.BG)

# Subset just the Block-Groups and Block we need
Mobile2000.BG4 <- subset(Mobile2000.BG, Mobile2000.BG$BKGPIDFP00 == "010970011004")
Mobile2000.BG3 <- subset(Mobile2000.BG, Mobile2000.BG$BKGPIDFP00 == "010970011003")
Mobile2000.BG2 <- subset(Mobile2000.BG, Mobile2000.BG$BKGPIDFP00 == "010970011002")
Mobile2010.BG <- subset(Mobile2010.BG, Mobile2010.BG$GEOID10 == "010970011003")

# Plot quickly to see how it looks
plot(Mobile2000.BG2, border='green')
plot(Mobile2000.BG3, add=T, border='orange', lty=2)
plot(Mobile2000.BG4, add=T, border='hotpink', lty=2)
plot(Mobile2010.BG, add=T, border='blue', lty=2)

# Grab the OpenStreetMap to put underneath
bb <- array(dim=c(2,2,4))
bb[,,1] = bbox(spTransform(Mobile2000.BG2, CRS("+init=epsg:4326")))
bb[,,2] = bbox(spTransform(Mobile2000.BG3, CRS("+init=epsg:4326")))
bb[,,3] = bbox(spTransform(Mobile2000.BG4, CRS("+init=epsg:4326")))
bb[,,4] = bbox(spTransform(Mobile2010.BG, CRS("+init=epsg:4326")))
bb.Mobile <- matrix(c( min(bb[1,1,]),min(bb[2,1,]),max(bb[1,2,]),max(bb[2,2,]) ), 2,2)
map.Mobile = openmap(c(bb.Mobile[2, 2], bb.Mobile[1, 1]), c(bb.Mobile[2, 1], bb.Mobile[1, 2]),
    type = "osm")

# Using diverging colors chosen from ColorBrewer
# http://colorbrewer2.org/
MCols <- c("#7570B3","#1B9E77","#D95F02","#E7298A")
# You could also get them from RColorBrewer package:
#library(RColorBrewer)
#brewer.pal(4,"Dark2")[c(3,1,2,4)]

# Plot nicely to a PDF or PNG image
#png("Mobile.png",480,480) # width&height in pixels
pdf("Mobile.pdf",7,7) # width/height in inches
plot(map.Mobile, removeMargin=TRUE)
plot(spTransform(Mobile2010.BG, osm()),  add=T, lwd=3, lty=1, border=MCols[1])
plot(spTransform(Mobile2000.BG2, osm()), add=T, lwd=1, lty=1, border=MCols[2])
plot(spTransform(Mobile2000.BG3, osm()), add=T, lwd=5, lty=2, border=MCols[3])
plot(spTransform(Mobile2000.BG4, osm()), add=T, lwd=6, lty=3, border=MCols[4])
legend("topleft", legend=c("Census 2010 Block-Group 3","Census 2000 Block-Group 2",
	"Census 2000 Block-Group 3","Census 2000 Block-Group 4"),
	lwd=c(3,1,5,6), lty=c(1,1,2,3), col=MCols, bg="white")
dev.off()