4

I have created state and county choropleths in R using ggplot2. The method I have used for this came from Winston Chang's "R Graphics Cookbook". For those choropleths the map was already represented as a data.frame in R (as opposed to a shapefile).

I now need to create choropleths of the continental US at the ZIP code level. The boundaries that I would like to use are the US Census ZCTA boundaries which you can download as a shapefile here (select "ZIP Code Tabulation Areas" from the dropdown). The shapefile is large, at over 500MB zipped.

I do not have a background in GIS, and attempted to create a choropleth from the shapefile in R using ggplot2 by following the advice on the ggplot2 wiki: Plotting Polygon Shapefiles.

However, when I attempt to follow the advice listed there the program just hangs:

require("rgdal") # requires sp, will use proj.4 if installed
require("maptools")
require("ggplot2")
require("plyr")

setwd("~/Desktop/tl_2013_us_zcta510")
zcta = readOGR(dsn=".", layer="tl_2013_us_zcta510")
zcta@data$id = rownames(zcta@data)

zcta.points = fortify(zcta, region="id")

The call to fortify hangs - even running it overnight it did not complete. I have dug into the code for fortify a bit and can see the line that it causing it to hang (coords <- ldply(model@polygons,fortify)), but am not sure how to fix it. I have also tried to follow the advice I have seen on numerous other websites about reading in shapefiles into R and ggplot2, but they don't seem to work for me either. I can't tell if it is because of the size of this shapefile or for some other reason which I do not understand because of my lack of knowledge about GIS.

It would be great if someone could help me with this.

Improve this question

2 Answers 2

Reset to default
2

Yes. My guess too. Try the subset suggested. If that works you can either simplify the shape file and try again or loop through subsets and rbind them together.

To simplify I'd look at the thinnedSpatialPoly() function in the maptools library. I don't have much guidance regarding the right threshold setting but something like

require("rgdal") # requires sp, will use proj.4 if installed
require("maptools")
require("ggplot2")
require("plyr")

Map <- readShapeSpatial('~/Desktop/tl_2013_us_zcta510/tl_2013_us_zcta510.shp')
#simplify the spatial polygon:  (not sure about ideal tolerance param)
zcta <- thinnedSpatialPoly(Map, tolerance=0.05)

zcta@data$id = rownames(zcta@data)

zcta.points = fortify(zcta, region="id")

Even if this works you will have to wait a long time for ggplot to generate a nationwide choropleth at zip level. Your code looks right to me

Improve this answer
3
  • I did the subset suggested and then rbind'd them together as suggested. The result is a large data.frame which takes ages to do any operations with (such as merging or rendering). I like your suggestion of simplifying the shapefile. I don't have prior GIS experience - can you explain how to do that? My use case is for exploratory data analysis of business data (which will eventually be mashed up with census data - such as sales per capita). So in my case precise boundaries are much less important being able to rapidly create multiple maps in a short amount of time.
    – Ari
    Commented Feb 21, 2014 at 19:10
  • Additionally, is there a way to remove ZCTAs from the shapefile that are not part of the continental US? That would have the benefit of a) likely reducing the size of the shapefile quite a bit and b) making the map more analogous to the existing state and county maps already present in R's maps package
    – Ari
    Commented Feb 21, 2014 at 19:55
  • I'd inspect the dataframe embedded in the spatialpolygon 'zcta'. Use head(zcta@data) and str(zcta@data) to get a sense of what is available within the zcta shapefile. If you don't have anything obvious to subset with you might have to use another shapefile to mask the zcta one... I don't know how this is done in R (though I imagine it is possible).
    – Michael Williams
    Commented Feb 22, 2014 at 4:36
0

Are you sure you're not running out of memory? Did you try subsetting the data first? 

zcta_subset <- zcta[as.integer(zcta@data$id) < 10,]
Improve this answer
1
  • Thank you - subsetting seemed to work. I then ran the loop over 1% of the shapefil at a time and rbind'd the final data.frame. But the resulting data.frame is so large that it takes ages to work with.
    – Ari
    Commented Feb 21, 2014 at 19:16

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.