Visualizing Polygon Areas of One Variable

ArcMap

  1. Single symbol (0:00)
  2. Categorical choropleths (0:45)
    • Transparency (1:05)
  3. Quantity choropleths (1:55)
    • Normalization (2:30)
    • Add field / field calculator (3:00)
  4. Classification
    • Jenks natural breaks classification (4:30)
    • Equal interval classification (5:00)
    • Geometric interval classification (5:35)
    • Quantile classification (5:55)
    • Standard deviation classification (6:10)
    • Manual classification (6:50)
    • Lying with maps (7:30)
  5. Graduated symbol maps (7:50)
  6. Dot density maps (8:50)
  7. Pie chart maps (10:00)

Visualizing Polygon Areas in R

Spatial polygons can be loaded with the readOGR() function from the rgdal library. For a shapefile of polygons, the function returns a SpatialPolygonsDataFrame object that holds both the polygon shapes and attributes associated with those polygons. Zipped shapefiles of US data are available here

Loading Polygons

Polygons can be accessed and subsetted as if they were arrays. In the example below, the states of Alaska and Hawaii are removed to make the example maps just show the continental USA.

library(rgdal)
counties = readOGR(dsn=".", layer="2017-county-data", stringsAsFactors=F)

# Subset to remove Alaska, Hawaii and other territories so just the continential US is shown
counties = counties[!(counties$ST %in% c('AK', 'HI', NA)),]

A simple choropleth can be created with plot() or spplot() functions. The zcol parameter is the name of the field used to determine the colors for the polygons, in this case median household income:

spplot(counties, zcol="MEDHHINC")
An spplot() Choropleth

Plotting over an OpenStreetMap Base Map

While the spplot() visualization may be adequate for a quick exploration of the data, the appearance of the map will usually be improved if a base map is placed under the polygons to give them some context. The OpenStreetMap library has a number of functions for downloading and manipulating OpenStreetMaps that can be used as base maps. The openmap() function creates an object that can be plotted with plot():

library(OpenStreetMap)

upperLeft = c(55, -130)
lowerRight = c(20, -60)
basemap = openmap(upperLeft, lowerRight, type="osm")
plot(basemap)
An OpenStreetMap Base Map

The openmap() function must be passed the WGS 84 latitude/longitude of the upper left and bottom right corners of the desired map area. You can get the bounding box of a SpatialPolygonsDataFrame by examining the bbox slot. The actual upperLeft/lowerRight corners that you use should pad these values a bit to give a margin around the polygons.

> counties@bbox
         min       max
x -124.72584 -66.94989
y   24.49813  49.38436

To plot the polygons over the base map, you must transform the polygons to the projection used by the spherical mercator projection used by OpenStreetMap object. This is the same projection used by almost all web maps. The spTransform() function can be used to transform any Spatial* object as the first parameter, with the osm() passed as the second parameter specifying spherical mercator projection used by OpenStreetMap.

counties = spTransform(counties, osm())
plot(counties, col="gray", add=T)
Single Color Polygons

plot() can be used to plot a variety of different choropleths based on color parameter, which contains the colors used to draw each of the polygons.

Plotting a Choropleth with Dichotomous Data

To plot dichotomous data, the ifelse() function can be used to select colors based on whether a field does or does not have a specific value. The add=T parameter indicates that the polygons should be drawn on top of the existing plot. The border=NA parameter indicates not to draw any borders, which would clutter the map, especially in areas with alot of small counties.

# Dichotomous data
ramp = c("navy", "red2")
colors = ifelse(counties$WIN2012 == "Obama", ramp[1], ramp[2])

plot(basemap)
plot(counties, col=colors, border=NA, add=T)
Choropleth of Dichotomous Data

As with any choropleth, a legend() should be provided to indicate what the colrs mean. The legend parameter is the list of names for the legend. The fill parameter indicates the color of the boxes to be drawn beside the names. The bg="white" parameter indicates to fill the box with a white background so it stands above whatever is under it. The inset=0.05 parameter indicates that the legend should be drawn 5% inside of the edge.

The mtext() function can be used to draw a title.

legend(x = "bottomright", legend=c("Obama", "Romney"), fill = ramp, bg="white", inset=0.05)
mtext("2012 Election Winner By County", cex=1.5, font=2)
Choropleth With Legend and Title

Plotting a Choropleth With Categorical Data

To plot categorical data:

# Categorical data
ramp = rainbow(length(unique(counties$ST)))
colors = ramp[as.factor(counties$ST)]

plot(basemap)
plot(counties, col=colors, border=NA, add=T)
Choropleth of Categorical Data

Plotting A Choropleth From Continuous Data (Interval or Ratio)

Quantile Classification

For making a choropleth from a quantitative field, the values must be classified into ranges that are represented with different colors.

In the example below, four colors are interpolated between skyblue1 and navy (dark blue) to create a four-color color ramp of shades of blue:

# Quantile Classification

breaks = quantile(counties$MEDHHINC, na.rm=T)
categories = as.numeric(cut(counties$MEDHHINC, breaks))

palette = colorRampPalette(c("skyblue1", "navy"))
ramp = palette(4)
colors = ramp[categories]

plot(basemap)
plot(counties, col=colors, border=NA, add=T)

A legend for a choropleth where the colors represent ranges of values can be constructed with the sapply() function, which pastes together successive pairs of breakpoints to give labels reflecting those ranges:

labels = sapply(1:4, function(z) paste(round(breaks[z]), '-', round(breaks[z + 1])))
legend(x = "bottomright", legend=labels, fill=ramp, bg="white", inset=0.05, cex=0.8)
mtext("Median Household Income (ACS 2016)", cex=1.5, font=2)
Choropleth of Continuous Data With a Quantile Classification

Equal Interval Classification

Since quantiles can exaggerate differences, sometimes an equal interval classification is more appropriate. For example, the map below reinforces a narrative of the US with high concentrations of wealth in a very limited number of wealthy urban counties, surrounded by vast, considerably less-well-off rural areas.

# Equal Interval Classification

breaks = qunif(seq(0, 1, 0.25), min = min(counties$MEDHHINC, na.rm=T), max = max(counties$MEDHHINC, na.rm=T))

categories = as.numeric(cut(counties$MEDHHINC, breaks))
categories = cut(counties$MEDHHINC, breaks)

palette = colorRampPalette(c("skyblue1", "navy"))
ramp = palette(4)
colors = ramp[categories]

plot(basemap)
plot(counties, col=colors, border=NA, add=T)

labels = sapply(1:4, function(z) paste(round(breaks[z]), '-', round(breaks[z + 1])))
legend(x = "bottomright", legend=labels, fill=ramp, bg="white", inset=0.05, cex=0.8)
mtext("Median Household Income (ACS 2016)", cex=1.5, font=2)
Choropleth of Continuous Data With a Equal-Interval Classification

Zero-Based Equal Interval Classification

Similarly, stretching the colors between the minimum and maximum can also exaggerate differences. Starting with a minimum of zero creates the map with the least contrast, which accentuates a narrative of the US as a wealthy country with a broadly dispersed middle-class, especially in comparison to many other less-wealthy countries. Whether that story is appropriate is a value judgement.

# Equal Interval Classification

breaks = qunif(seq(0, 1, 0.25), min = 0, max = max(counties$MEDHHINC, na.rm=T))

categories = as.numeric(cut(counties$MEDHHINC, breaks))
categories = cut(counties$MEDHHINC, breaks)

palette = colorRampPalette(c("skyblue1", "navy"))
ramp = palette(4)
colors = ramp[categories]

plot(basemap)
plot(counties, col=colors, border=NA, add=T)

labels = sapply(1:4, function(z) paste(round(breaks[z]), '-', round(breaks[z + 1])))
legend(x = "bottomright", legend=labels, fill=ramp, bg="white", inset=0.05, cex=0.8)
mtext("Median Household Income (ACS 2016)", cex=1.5, font=2)
Choropleth of Continuous Data With a Equal-Interval Zero-Based Classification

Transparency

Colors can be specified in the same type of RGB hexadecimal codes used for specifying colors on web pages. A fourth byte can be added for alpha (#RRGGBBAA), with 00 meaning transparent and FF meaning opaque. Specifying colors as semi-transparent will allow the base map features to be visible through the choropleth.

ramp = c("#0000FF80", "#FF000080")
colors = ifelse(counties$WIN2012 == "Obama", ramp[1], ramp[2])

plot(basemap)
plot(counties, col=colors, border=NA, add=T)
Semi-Transparent Categorical Choropleth

With color ramps, the alpha=T parameter needs to be passed to colorRampPalette() so transparency in included in interpolated colors.

breaks = quantile(counties$MEDHHINC, na.rm=T)
categories = as.numeric(cut(counties$MEDHHINC, breaks))

palette = colorRampPalette(c("#87CEFF80", "#00008080"), alpha=T)
ramp = palette(4)
colors = ramp[categories]

plot(basemap)
plot(counties, col=colors, border=NA, add=T)

labels = sapply(1:4, function(z) paste(round(breaks[z]), '-', round(breaks[z + 1])))
legend(x = "bottomright", legend=labels, fill=ramp, bg="white", inset=0.05, cex=0.8)
mtext("Median Household Income (ACS 2016)", cex=1.5, font=2)
Semi-Transparent Quantile Choropleth

Other Base Maps

The openmap() function can load base maps from other sources besides OpenStreetMap. See ?openmap documentation page for a list of values for the type parameter. One attractive map is the stamen-watercolor map from Stamen Design:

basemap = openmap(upperLeft, lowerRight, type="stamen-watercolor")
plot(basemap)
Stamen Watercolor Base Map

Other Projections

While it is usually safest when using OpenStreetMap to just use the default spherical mercator projection, the openproj() function can be used to transform the map and plot to another projection. Note that you must also spTransform() your polygons to the same projection to plot.

With the conic projection below, the edges of the base map are no longer at right angles, which may be undesirable and would require loading and expanded area and cropping to restore.

# EPSG: 102008 North America Albers Equal Area Conic

usa_albers = CRS("+proj=aea +lat_1=20 +lat_2=60 +lat_0=40 +lon_0=-96 +x_0=0 +y_0=0 
	+datum=NAD83 +units=m +no_defs")

basemap = openproj(basemap, usa_albers)
counties = spTransform(counties, usa_albers)

ramp = c("#2693a0", "#a03626")
colors = ifelse(counties$WIN2012 == "Obama", ramp[1], ramp[2])

plot(basemap)
plot(counties, col=colors, border=NA, add=T)
legend(x = "bottomright", legend=c("Obama", "Romney"), fill=ramp, bg="white", inset=0.05)
Albers Equal Area Conic Projection