Plotting Maps with R

Install and load the following R packages.

library(maptools)
library(RColorBrewer)
library(classInt)
library(rgdal)

First download the shapefile and save it on your computer. Copy the path to the directory containing the shapefile.

Set the working directory to the directory where you have save the shapefile. Session > Set Working Directory > Choose Directory. Read in the shape file for the Bay Area.

## set the working directory.  Windows
# setwd("C:\\Users\\Eric A. Suess\\Downloads\\bayarea_zipcodes\\")
## set the working directory.  Mac or Linux
setwd("~/classes/2016-2017/02 Winter 2017/Stat 6610/Activity/BayArea")
## load the shapefile
zip <- readOGR("bayarea_zipcodes.shp")

OGR data source with driver: ESRI Shapefile 
Source: "bayarea_zipcodes.shp", layer: "bayarea_zipcodes"
with 187 features
It has 5 fields

Add the median household incomes for each zipcode to the .shp file.

med.income <- c(55995, 54448, 50520, 43846, 50537, 67705, 45142, 36518, NA, 65938, 50500, 61022, 65959, 17188, 53881, 66970, 35699, 76194, 63777, 63838, 66010, 48523, 70758, 41002, 60082, 48672, 61429, 60402, 75707, 58333, 60769, 60375, 76627, 68515, 64485, 60971, 60833, 54732, 60804, 33962, 57601, 43649, 59889, NA, 61494, 67824, 64429, 82528, 101555, 96658, 50300, 41573, 75747, 53750, 56905, 85479, 77455, 64389, 91283, 119832, 103791, 100590, 85109, 88184, 57153, 34951, 51418, 38613, 43640, 19750, 139997, 39290, NA, 76808, 98525, 106492, 68112, 68853, 77952, 34398, 75026, 33556, NA, 109771, 142459, 80959, 21124, 49066, 95588, 55321, 20034, 57976, 40990, 73571, 84710, 43444, 32273, 56569, 54174, 105393, 51896, 31542, 33152, 54879, 88976, 14609, 61609, 61776, 22351, 31131, 47288, NA, 63983, 60733, 29181, 75727, 61362, 53795, 76044, 34755, 66627, 37146, 92644, 87855, 95313, 50888, 55000, 57629, 54342, 77122, 44723, 64534, 65658, 60711, 57214, 54594, 48523, 90107, 69014, 49452, 72288, 56973, 81923, 61289, 71863, 61939, 0, 82735, 68067, 82188, 70026, 101977, 55112, 84442, 82777, 82796, 92989, 67152, 68121, 69350, 104958, 49279, 80973, 89016, 96677, 89572, 64256, 84565, 16250, 64839, 200001, 82072, 58304, 66807, 97758, 68721, 77539, 41313, NA, 82314, 164479, 69087, 145425, NA, 71056, 128853, 84856)
zip$INCOME <- med.income

Select color palette and the number colors (levels of income) to represent on the map.

colors <- brewer.pal(9, "YlOrRd")

Set breaks for the 9 colors.

brks <- classIntervals(zip$INCOME, n=9, style="quantile")

var has missing values, omitted in finding classes

brks <- brks$brks

Plot the map. Add a title. Add a legend.

plot(zip, col=colors[findInterval(zip$INCOME, brks,all.inside=TRUE)])
title(paste ("SF Bay Area Median Household Income"))

legend(x=6298809, y=2350000, legend=leglabs(round(brks)), fill=colors, bty="n", x.intersp = .8, y.intersp = .8)