hepplerj/choro.R

## choro.R
library(ggmap)
library(rgdal)
library(tidyr)
library(dplyr)
library(ggplot2)
library(RColorBrewer)

# Read data
# The shapefile comes from NHGIS, modified to include only California census tracts.
tract_1960_shp <- readOGR("/Users/jheppler/data/census-transformed/1960-census/", "1960-ca-tracts")

names(tract_1960_shp)
# [1] "NHGISST"    "NHGISCTY"   "GISJOIN"    "GISJOIN2"   "SHAPE_AREA" "SHAPE_LEN"  "nhgis0007_" "nhgis000_1"
# [9] "nhgis000_2" "nhgis000_3" "nhgis000_4" "nhgis000_5" "nhgis000_6" "nhgis000_7" "nhgis000_8" "nhgis000_9"
#[17] "nhgis00010" "nhgis00011" "nhgis00012" "nhgis00013" "nhgis00014" "nhgis00015" "nhgis00016" "nhgis00017"
#[25] "nhgis00018" "nhgis00019" "nhgis00020" "nhgis00021" "nhgis00022" "nhgis00023" "nhgis00024" "nhgis00025"
#[33] "nhgis00026" "nhgis00027" "nhgis00028" "nhgis00029" "nhgis00030" "nhgis00031" "nhgis00032" "nhgis00033"
#[41] "nhgis00034" "nhgis00035" "nhgis00036" "nhgis00037" "nhgis00038" "nhgis00039" "nhgis00040" "nhgis00041"
#[49] "nhgis00042" "nhgis00043" "nhgis00044" "nhgis00045" "nhgis00046" "nhgis00047" "nhgis00048" "nhgis00049"
#[57] "nhgis00050" "nhgis00051" "nhgis00052" "nhgis00053" "nhgis00054" "nhgis00055" "nhgis00056" "nhgis00057"
#[65] "nhgis00058" "nhgis00059" "nhgis00060" "nhgis00061" "nhgis00062" "nhgis00063" "nhgis00064" "nhgis00065"
#[73] "nhgis00066" "nhgis00067" "nhgis00068" "nhgis00069" "nhgis00070" "nhgis00071" "nhgis00072" "nhgis00073"
#[81] "nhgis00074" "nhgis00075" "nhgis00076" "nhgis00077" "nhgis00078" "nhgis00079" "nhgis00080" "nhgis00081"
#[89] "nhgis00082" "nhgis00083" "nhgis00084" "nhgis00085" "nhgis00086"

glimpse(tract_1960_shp)
#Observations: 582
#Variables:
#$ NHGISST    (fctr) 060, 060, 060, 060, 060, 060, 060, 060, 060, 060, 060, 060, 060, 060, 060, 060, 060, 0...
#$ NHGISCTY   (fctr) 0010, 0010, 0010, 0010, 0010, 0010, 0010, 0010, 0010, 0010, 0010, 0010, 0010, 0010, 00...
#$ GISJOIN    (fctr) G0600010CV0022, G0600010HA0049, G0600010HA0051, G0600010HA0048, G0600010HA0053, G06000...
#$ GISJOIN2   (fctr) 0600010CV0022, 0600010HA0049, 0600010HA0051, 0600010HA0048, 0600010HA0053, 0600010OK00...
#$ SHAPE_AREA (dbl) 4.827141e+06, 2.035771e+06, 1.785122e+06, 1.307165e+06, 1.810945e+06, 4.195454e+05, 6.9...
#$ SHAPE_LEN  (dbl) 14137.0237, 7722.5626, 5902.0262, 7144.6136, 6301.1831, 3107.9843, 4515.7960, 3519.7785...
#$ nhgis0007_ (fctr) 1960, 1960, 1960, 1960, 1960, 1960, 1960, 1960, 1960, 1960, 1960, 1960, 1960, 1960, 19...
#$ nhgis000_1 (fctr) California, California, California, California, California, California, California, Ca...
#$ nhgis000_2 (fctr) 06, 06, 06, 06, 06, 06, 06, 06, 06, 06, 06, 06, 06, 06, 06, 06, 06, 06, 06, 06, 06, 06...
#$ nhgis000_3 (fctr) Alameda, Alameda, Alameda, Alameda, Alameda, Alameda, Alameda, Alameda, Alameda, Alame...
#$ nhgis000_4 (fctr) 001, 001, 001, 001, 001, 001, 001, 001, 001, 001, 001, 001, 001, 001, 001, 001, 001, 0...
#$ nhgis000_5 (fctr) CV, HA, HA, HA, HA, OK, OK, OK, OK, OK, OK, OK, OK, OK, OK, OK, OK, OK, OK, OK, OK, OK...
#$ nhgis000_6 (fctr) 0022, 0049, 0051, 0048, 0053, 0033, 0045, 0047, 0046, 0008, 0010, 0013, 0012, 0014, 00...
#$ nhgis000_7 (fctr) NA, NA, NA, NA, NA, B, A, A, A, NA, NA, NA, NA, A, B, A, NA, NA, NA, NA, NA, NA, NA, B...
#$ nhgis000_8 (fctr) TRACT CV-22 PT CASTRO VALLEY DIV, TRACT HA-49 PT HAYWARD DIV, TRACT HA-51 PT HAYWARD C...
#$ nhgis000_9 (fctr) 253, 399, 341, 204, 442, 114, 100, 89, 148, 110, 69, 43, 25, 16, 11, 60, 91, 132, 120,...
#$ nhgis00010 (fctr) 229, 381, 333, 212, 403, 103, 89, 94, 158, 108, 66, 39, 26, 10, 12, 57, 79, 132, 98, 2...
#$ nhgis00011 (fctr) 319, 206, 212, 303, 375, 61, 109, 82, 77, 69, 33, 31, 21, 7, 3, 40, 77, 96, 95, 215, 1...
#$ nhgis00012 (fctr) 344, 234, 193, 269, 388, 62, 88, 60, 82, 71, 51, 42, 16, 11, 3, 52, 94, 78, 81, 213, 1...
#$ nhgis00013 (fctr) 359, 236, 185, 270, 338, 63, 129, 88, 85, 80, 39, 36, 16, 7, 11, 39, 88, 65, 66, 194, ...
#$ nhgis00014 (fctr) 324, 254, 217, 284, 340, 50, 138, 69, 85, 84, 27, 34, 27, 2, 7, 38, 101, 60, 89, 219, ...
#$ nhgis00015 (fctr) 188, 195, 137, 161, 159, 45, 93, 70, 64, 59, 38, 26, 28, 4, 6, 20, 59, 53, 79, 182, 11...
#$ nhgis00016 (fctr) 199, 263, 178, 176, 195, 66, 118, 100, 177, 67, 44, 31, 67, 1, 2, 35, 92, 133, 101, 17...
#$ nhgis00017 (fctr) 42, 293, 204, 74, 151, 87, 57, 52, 125, 70, 49, 46, 35, 7, 9, 23, 65, 131, 125, 203, 1...
#$ nhgis00018 (fctr) 55, 320, 222, 66, 199, 118, 58, 82, 229, 95, 74, 42, 53, 9, 5, 22, 90, 264, 151, 197, ...
#$ nhgis00019 (fctr) 69, 239, 180, 73, 210, 129, 57, 69, 199, 78, 66, 53, 29, 12, 11, 28, 80, 185, 127, 174...
#$ nhgis00020 (fctr) 92, 191, 161, 91, 226, 97, 64, 49, 153, 60, 64, 44, 25, 4, 1, 27, 93, 158, 109, 172, 1...
#$ nhgis00021 (fctr) 146, 191, 129, 123, 242, 85, 58, 58, 109, 71, 70, 52, 36, 7, 6, 30, 93, 177, 139, 170,...
#$ nhgis00022 (fctr) 222, 150, 123, 159, 255, 77, 72, 63, 117, 63, 48, 43, 26, 5, 12, 50, 97, 160, 110, 179...
#$ nhgis00023 (fctr) 273, 180, 132, 222, 285, 87, 90, 60, 95, 85, 48, 57, 35, 8, 10, 45, 97, 135, 91, 158, ...
#$ nhgis00024 (fctr) 297, 217, 148, 255, 291, 81, 114, 70, 100, 89, 53, 42, 23, 4, 10, 53, 118, 125, 115, 1...
#$ nhgis00025 (fctr) 262, 198, 131, 231, 214, 67, 92, 56, 90, 63, 65, 61, 28, 1, 12, 38, 80, 110, 101, 158,...
#$ nhgis00026 (fctr) 253, 209, 136, 214, 187, 82, 122, 91, 97, 81, 79, 50, 35, 6, 7, 35, 147, 135, 109, 184...
#$ nhgis00027 (fctr) 178, 181, 97, 166, 145, 64, 126, 73, 89, 108, 70, 66, 45, 6, 17, 50, 117, 116, 113, 19...
#$ nhgis00028 (fctr) 157, 197, 129, 133, 150, 100, 148, 88, 101, 125, 104, 53, 54, 7, 15, 47, 149, 158, 168...
#$ nhgis00029 (fctr) 109, 182, 95, 103, 120, 68, 110, 79, 90, 111, 83, 97, 61, 14, 17, 71, 113, 120, 111, 1...
#$ nhgis00030 (fctr) 89, 205, 102, 83, 102, 103, 155, 100, 108, 139, 92, 83, 79, 13, 15, 48, 148, 144, 161,...

# `fortify` fails here. RStudio gives me a fatal error. Running on the command line tells me there's a segmentation fault.
tract_1960_df <- fortify(tract_1960_shp, region = "GISJOIN")

data <- read.csv("/Users/jheppler/data/census-data/nhgis0007_ds92_1960_tract.csv", sep=",", header=T)

# Join data
full.data <- left_join(tract_1960_df@data, data, by="GISJOIN")

# Change column names
names(full.data)[names(full.data) == 'B7B002'] <- 'black'
names(full.data)[names(full.data) == "B7B001"] <- 'white'
names(full.data)[names(full.data) == "B7B003"] <- 'other'
glimpse(full.data) #check

# Set color scheme for choropleth
pal <- brewer.pal(7, "OrRd")

# Map!
	library(ggmap)
	library(rgdal)
	library(tidyr)
	library(dplyr)
	library(ggplot2)
	library(RColorBrewer)

	# Read data
	# The shapefile comes from NHGIS, modified to include only California census tracts.
	tract_1960_shp <- readOGR("/Users/jheppler/data/census-transformed/1960-census/", "1960-ca-tracts")

	names(tract_1960_shp)
	# [1] "NHGISST" "NHGISCTY" "GISJOIN" "GISJOIN2" "SHAPE_AREA" "SHAPE_LEN" "nhgis0007_" "nhgis000_1"
	# [9] "nhgis000_2" "nhgis000_3" "nhgis000_4" "nhgis000_5" "nhgis000_6" "nhgis000_7" "nhgis000_8" "nhgis000_9"
	#[17] "nhgis00010" "nhgis00011" "nhgis00012" "nhgis00013" "nhgis00014" "nhgis00015" "nhgis00016" "nhgis00017"
	#[25] "nhgis00018" "nhgis00019" "nhgis00020" "nhgis00021" "nhgis00022" "nhgis00023" "nhgis00024" "nhgis00025"
	#[33] "nhgis00026" "nhgis00027" "nhgis00028" "nhgis00029" "nhgis00030" "nhgis00031" "nhgis00032" "nhgis00033"
	#[41] "nhgis00034" "nhgis00035" "nhgis00036" "nhgis00037" "nhgis00038" "nhgis00039" "nhgis00040" "nhgis00041"
	#[49] "nhgis00042" "nhgis00043" "nhgis00044" "nhgis00045" "nhgis00046" "nhgis00047" "nhgis00048" "nhgis00049"
	#[57] "nhgis00050" "nhgis00051" "nhgis00052" "nhgis00053" "nhgis00054" "nhgis00055" "nhgis00056" "nhgis00057"
	#[65] "nhgis00058" "nhgis00059" "nhgis00060" "nhgis00061" "nhgis00062" "nhgis00063" "nhgis00064" "nhgis00065"
	#[73] "nhgis00066" "nhgis00067" "nhgis00068" "nhgis00069" "nhgis00070" "nhgis00071" "nhgis00072" "nhgis00073"
	#[81] "nhgis00074" "nhgis00075" "nhgis00076" "nhgis00077" "nhgis00078" "nhgis00079" "nhgis00080" "nhgis00081"
	#[89] "nhgis00082" "nhgis00083" "nhgis00084" "nhgis00085" "nhgis00086"

	glimpse(tract_1960_shp)
	#Observations: 582
	#Variables:
	#$ NHGISST (fctr) 060, 060, 060, 060, 060, 060, 060, 060, 060, 060, 060, 060, 060, 060, 060, 060, 060, 0...
	#$ NHGISCTY (fctr) 0010, 0010, 0010, 0010, 0010, 0010, 0010, 0010, 0010, 0010, 0010, 0010, 0010, 0010, 00...
	#$ GISJOIN (fctr) G0600010CV0022, G0600010HA0049, G0600010HA0051, G0600010HA0048, G0600010HA0053, G06000...
	#$ GISJOIN2 (fctr) 0600010CV0022, 0600010HA0049, 0600010HA0051, 0600010HA0048, 0600010HA0053, 0600010OK00...
	#$ SHAPE_AREA (dbl) 4.827141e+06, 2.035771e+06, 1.785122e+06, 1.307165e+06, 1.810945e+06, 4.195454e+05, 6.9...
	#$ SHAPE_LEN (dbl) 14137.0237, 7722.5626, 5902.0262, 7144.6136, 6301.1831, 3107.9843, 4515.7960, 3519.7785...
	#$ nhgis0007_ (fctr) 1960, 1960, 1960, 1960, 1960, 1960, 1960, 1960, 1960, 1960, 1960, 1960, 1960, 1960, 19...
	#$ nhgis000_1 (fctr) California, California, California, California, California, California, California, Ca...
	#$ nhgis000_2 (fctr) 06, 06, 06, 06, 06, 06, 06, 06, 06, 06, 06, 06, 06, 06, 06, 06, 06, 06, 06, 06, 06, 06...
	#$ nhgis000_3 (fctr) Alameda, Alameda, Alameda, Alameda, Alameda, Alameda, Alameda, Alameda, Alameda, Alame...
	#$ nhgis000_4 (fctr) 001, 001, 001, 001, 001, 001, 001, 001, 001, 001, 001, 001, 001, 001, 001, 001, 001, 0...
	#$ nhgis000_5 (fctr) CV, HA, HA, HA, HA, OK, OK, OK, OK, OK, OK, OK, OK, OK, OK, OK, OK, OK, OK, OK, OK, OK...
	#$ nhgis000_6 (fctr) 0022, 0049, 0051, 0048, 0053, 0033, 0045, 0047, 0046, 0008, 0010, 0013, 0012, 0014, 00...
	#$ nhgis000_7 (fctr) NA, NA, NA, NA, NA, B, A, A, A, NA, NA, NA, NA, A, B, A, NA, NA, NA, NA, NA, NA, NA, B...
	#$ nhgis000_8 (fctr) TRACT CV-22 PT CASTRO VALLEY DIV, TRACT HA-49 PT HAYWARD DIV, TRACT HA-51 PT HAYWARD C...
	#$ nhgis000_9 (fctr) 253, 399, 341, 204, 442, 114, 100, 89, 148, 110, 69, 43, 25, 16, 11, 60, 91, 132, 120,...
	#$ nhgis00010 (fctr) 229, 381, 333, 212, 403, 103, 89, 94, 158, 108, 66, 39, 26, 10, 12, 57, 79, 132, 98, 2...
	#$ nhgis00011 (fctr) 319, 206, 212, 303, 375, 61, 109, 82, 77, 69, 33, 31, 21, 7, 3, 40, 77, 96, 95, 215, 1...
	#$ nhgis00012 (fctr) 344, 234, 193, 269, 388, 62, 88, 60, 82, 71, 51, 42, 16, 11, 3, 52, 94, 78, 81, 213, 1...
	#$ nhgis00013 (fctr) 359, 236, 185, 270, 338, 63, 129, 88, 85, 80, 39, 36, 16, 7, 11, 39, 88, 65, 66, 194, ...
	#$ nhgis00014 (fctr) 324, 254, 217, 284, 340, 50, 138, 69, 85, 84, 27, 34, 27, 2, 7, 38, 101, 60, 89, 219, ...
	#$ nhgis00015 (fctr) 188, 195, 137, 161, 159, 45, 93, 70, 64, 59, 38, 26, 28, 4, 6, 20, 59, 53, 79, 182, 11...
	#$ nhgis00016 (fctr) 199, 263, 178, 176, 195, 66, 118, 100, 177, 67, 44, 31, 67, 1, 2, 35, 92, 133, 101, 17...
	#$ nhgis00017 (fctr) 42, 293, 204, 74, 151, 87, 57, 52, 125, 70, 49, 46, 35, 7, 9, 23, 65, 131, 125, 203, 1...
	#$ nhgis00018 (fctr) 55, 320, 222, 66, 199, 118, 58, 82, 229, 95, 74, 42, 53, 9, 5, 22, 90, 264, 151, 197, ...
	#$ nhgis00019 (fctr) 69, 239, 180, 73, 210, 129, 57, 69, 199, 78, 66, 53, 29, 12, 11, 28, 80, 185, 127, 174...
	#$ nhgis00020 (fctr) 92, 191, 161, 91, 226, 97, 64, 49, 153, 60, 64, 44, 25, 4, 1, 27, 93, 158, 109, 172, 1...
	#$ nhgis00021 (fctr) 146, 191, 129, 123, 242, 85, 58, 58, 109, 71, 70, 52, 36, 7, 6, 30, 93, 177, 139, 170,...
	#$ nhgis00022 (fctr) 222, 150, 123, 159, 255, 77, 72, 63, 117, 63, 48, 43, 26, 5, 12, 50, 97, 160, 110, 179...
	#$ nhgis00023 (fctr) 273, 180, 132, 222, 285, 87, 90, 60, 95, 85, 48, 57, 35, 8, 10, 45, 97, 135, 91, 158, ...
	#$ nhgis00024 (fctr) 297, 217, 148, 255, 291, 81, 114, 70, 100, 89, 53, 42, 23, 4, 10, 53, 118, 125, 115, 1...
	#$ nhgis00025 (fctr) 262, 198, 131, 231, 214, 67, 92, 56, 90, 63, 65, 61, 28, 1, 12, 38, 80, 110, 101, 158,...
	#$ nhgis00026 (fctr) 253, 209, 136, 214, 187, 82, 122, 91, 97, 81, 79, 50, 35, 6, 7, 35, 147, 135, 109, 184...
	#$ nhgis00027 (fctr) 178, 181, 97, 166, 145, 64, 126, 73, 89, 108, 70, 66, 45, 6, 17, 50, 117, 116, 113, 19...
	#$ nhgis00028 (fctr) 157, 197, 129, 133, 150, 100, 148, 88, 101, 125, 104, 53, 54, 7, 15, 47, 149, 158, 168...
	#$ nhgis00029 (fctr) 109, 182, 95, 103, 120, 68, 110, 79, 90, 111, 83, 97, 61, 14, 17, 71, 113, 120, 111, 1...
	#$ nhgis00030 (fctr) 89, 205, 102, 83, 102, 103, 155, 100, 108, 139, 92, 83, 79, 13, 15, 48, 148, 144, 161,...

	# `fortify` fails here. RStudio gives me a fatal error. Running on the command line tells me there's a segmentation fault.
	tract_1960_df <- fortify(tract_1960_shp, region = "GISJOIN")

	data <- read.csv("/Users/jheppler/data/census-data/nhgis0007_ds92_1960_tract.csv", sep=",", header=T)

	# Join data
	full.data <- left_join(tract_1960_df@data, data, by="GISJOIN")

	# Change column names
	names(full.data)[names(full.data) == 'B7B002'] <- 'black'
	names(full.data)[names(full.data) == "B7B001"] <- 'white'
	names(full.data)[names(full.data) == "B7B003"] <- 'other'
	glimpse(full.data) #check

	# Set color scheme for choropleth
	pal <- brewer.pal(7, "OrRd")

	# Map!