Mapping Income by Neighborhood with ggplot2

Question

I have the mean income data for different section of Oklahoma City and I'm trying to create a choropleth of the city. Here is the data.

library(acs)
library(ggplot2)
library(ggmap)
library(UScensus2010)
library(spdep)
library(RColorBrewer)

# Use your census API key (You'll need to get your own API key)
#http://api.census.gov/data/key_signup.html
api.key.install(key="c369cd6ed053a84332caa62301eb8afe98bed825")

# Load in Shape File (You'll need to download this file from the census)
#ftp://ftp2.census.gov/geo/tiger/TIGER2013/TRACT/tl_2013_40_tract.zip
geodat<-readShapePoly((""), proj4string=CRS('+proj=longlat +datum=NAD83'))

geodat<-geodat[geodat$COUNTYFP==109,]

# Tract Level Data
mdat<-geodat@data

# American Community Survey Data: Median HH Income for OK Census Tracts
ok.counties=geo.make(state="OK", county="Oklahoma", tract="*")
ok.income<-acs.fetch(geography=ok.counties, table.number="B19013", endyear=2013)

# Merge Data Sets 
geo_dat<-geography(ok.income)
var_dat<-as.data.frame(estimate(ok.income))
acs_data<-cbind(geo_dat,var_dat)
mdat2<-merge(mdat,acs_data,by.x="TRACTCE", by.y= "tract")

#cut data
mdat2$cut <- cut(mdat2$B19013_001, breaks = c(5000, 200000, by = 10000))
map <- mdat2[, 11:12]

#create map
 okc <- ggplot(mdat2, aes(x = INTPTLON, y = INTPTLAT, fill = B19013_001)) +
                  geom_map(aes(fill=cut), map=map)

I'm new to ggplot2 and can't seem to get it to work properly. I was able to create a choropleth using spplot but I would like to make one with ggplot so that I can layer it with other maps.

edit: Forgot to include that I am attempting to create a choropleth with median income as the fill. Median income is B19013_001 in the mdat2 data frame. The INTPTLON vector is the longitude coordinates for the different sections of the city and INTPTLAT is the latitude coordinates. Here is an example of a plot I created in sppolot. I am trying to make something similar in ggplot

Here are the first 10 observations from mdat2:

  TRACTCE STATEFP COUNTYFP       GEOID NAME.x          NAMELSAD MTFCC FUNCSTAT   ALAND AWATER    INTPTLAT
1  100100      40      109 40109100100   1001 Census Tract 1001 G5020        S 2141969      0 +35.5207581
2  100200      40      109 40109100200   1002 Census Tract 1002 G5020        S 3349288      0 +35.5039838
3  100300      40      109 40109100300   1003 Census Tract 1003 G5020        S 2098254      0 +35.5123370
4  100400      40      109 40109100400   1004 Census Tract 1004 G5020        S 2593440      0 +35.5004315
5  100500      40      109 40109100500   1005 Census Tract 1005 G5020        S 2734867      0 +35.5006179
6  100600      40      109 40109100600   1006 Census Tract 1006 G5020        S  551754      0 +35.5040878
      INTPTLON                                       NAME.y state county B19013_001
1 -097.5315537 Census Tract 1001, Oklahoma County, Oklahoma    40    109      33440
2 -097.5546474 Census Tract 1002, Oklahoma County, Oklahoma    40    109      29338
3 -097.5259348 Census Tract 1003, Oklahoma County, Oklahoma    40    109      85592
4 -097.4855750 Census Tract 1004, Oklahoma County, Oklahoma    40    109      21084
5 -097.5037840 Census Tract 1005, Oklahoma County, Oklahoma    40    109      23411
6 -097.5172200 Census Tract 1006, Oklahoma County, Oklahoma    40    109      62857

Answer 1

Here's a solution that works with ggplot:

library(acs)
library(ggplot2)
library(ggmap)
library(UScensus2010)
library(spdep)
library(RColorBrewer)


library(dplyr)
library(scales)

# Use your census API key (You'll need to get your own API key)
#http://api.census.gov/data/key_signup.html
api.key.install(key="c369cd6ed053a84332caa62301eb8afe98bed825")

# Load in Shape File (You'll need to download this file from the census)
#ftp://ftp2.census.gov/geo/tiger/TIGER2013/TRACT/tl_2013_40_tract.zip

## load, subset shapefile
geodat<-readShapePoly("/home/aksel/Downloads/ohio//tl_2013_40_tract.shp", proj4string=CRS('+proj=longlat +datum=NAD83'))
geodat<-geodat[geodat$COUNTYFP==109,]

## fortify for ggplot digestion
geodat.f<-fortify(geodat,region="GEOID")

# American Community Survey Data: Median HH Income for OK Census Tracts
ok.counties=geo.make(state="OK", county="Oklahoma", tract="*")
ok.income<-acs.fetch(geography=ok.counties, table.number="B19013", endyear=2013)


# Merge Data Sets 
geo_dat<-geography(ok.income)
var_dat<-as.data.frame(estimate(ok.income))
acs_data<-cbind(geo_dat,var_dat)
acs_data$id<- paste("40109", acs_data$tract, sep = "")

## from dplyr
mapdata<-left_join(geodat.f,acs_data)

ggplot() +
  geom_polygon(data = mapdata, aes(x = long, y = lat, group = group,
                                    fill = B19013_001), color = "black", size = 0.5)+
  scale_fill_distiller(palette = "Reds", labels = comma,
                       breaks = pretty_breaks(n = 10), values = c(1,0)) +
  guides(fill = guide_legend(reverse = TRUE)) +
  theme_nothing(legend = TRUE) +
ggtitle('Map of 40109')  

Which yields: