public
Last active

Calculating and mapping district compactness

  • Download Gist
district_compactness.R
R
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73
doInstall <- TRUE # Change to FALSE if you don't want packages installed.
toInstall <- c("maptools", "rgdal", "ggplot2", "spatstat", "RColorBrewer")
if(doInstall){install.packages(toInstall, repos = "http://cran.r-project.org")}
lapply(toInstall, library, character.only = TRUE)
 
# Taking an online compressed shapefile, and opening it in R
# From http://stackoverflow.com/a/3053883
temp <- tempfile() # 110th & 111th Congressional District Shapefiles
download.file("http://www.census.gov/geo/cob/bdy/cd/cd110shp/cd99_110_shp.zip",
temp) # See http://www.census.gov/geo/www/cob/cd110.html#shp
unzip(temp) # Will put the unzipped files in the working directory
shapeFile <- readShapeSpatial("cd99_110.shp") # Load shapefile
 
mapObject <- fortify(shapeFile) # Convert to a data.frame
mapObject <- data.frame(mapObject, shapeFile@data[mapObject$id, ])
head(mapObject) # ^ Add some labels
 
# Reduce to just the Continental U.S.
mapObject <- mapObject[mapObject$long < -50, ]
mapObject <- mapObject[mapObject$long > -135, ]
mapObject <- mapObject[mapObject$lat < 50, ]
mapObject <- mapObject[mapObject$lat > 22, ]
 
# Basic map of congressional districts in the Continental U.S.
with(mapObject, plot(long, lat, pch = ".", asp = 1))
 
### Calculate district area and perimeter length ###
 
mapObject$piece <- as.character(mapObject$piece)
mapObject$stateCD <- with(mapObject, paste(STATE, CD))
mapObject$Area <- mapObject$Perimeter <- NA
uniqueCDs <- sort(unique(mapObject$stateCD))
 
for(cd in uniqueCDs){ # This loop will take some time
print(cd)
cdShape <- mapObject[mapObject$stateCD == cd, ]
cdPoly <- SpatialPolygons(list(Polygons(lapply(split(cdShape[, c("long", "lat")],
cdShape$piece), Polygon), ID = "b")))
owinObject <- try(as(cdPoly, "owin")) # owin is very finicky, and I won't try
if(class(owinObject) == "try-error"){next()} # to troubleshoot polygons here.
plot(owinObject)
Sys.sleep(1e-10)
mapObject[mapObject$stateCD == cd, "Area"] <- area.owin(owinObject)
mapObject[mapObject$stateCD == cd, "Perimeter"] <- perimeter(owinObject)
}
 
# A simple ratio measure of district compactness:
mapObject$Compactness <- with(mapObject, Area / Perimeter)
plot(density(unique(mapObject$Compactness), na.rm = T))
 
### Plot ###
new_theme_empty <- theme_bw() # Create our own, mostly blank, theme
new_theme_empty$line <- element_blank()
new_theme_empty$rect <- element_blank()
new_theme_empty$strip.text <- element_blank()
new_theme_empty$axis.text <- element_blank()
new_theme_empty$axis.title <- element_blank()
new_theme_empty$plot.margin <- structure(c(0, 0, -1, -1), unit = "lines",
valid.unit = 3L, class = "unit")
myPalette <- colorRampPalette(brewer.pal(11, "Spectral"))
 
mapPlot <- ggplot(mapObject)
mapPlot <- mapPlot + geom_polygon(aes(x = long, y = lat, group = group,
fill = log(Compactness)),
colour = "white", lwd = 1/9)
mapPlot <- mapPlot + coord_map(project="conic", lat0 = 30)
mapPlot <- mapPlot + new_theme_empty
mapPlot <- mapPlot + ggtitle("A measure of district compactness")
mapPlot <- mapPlot + scale_fill_gradientn(colours = myPalette(100),
na.value = "gray80") # I just can't quit the Spectral palette!
#print(mapPlot)
ggsave("District compactness map.png", mapPlot, h = 9, w = 16, type = "cairo-png")

Please sign in to comment on this gist.

Something went wrong with that request. Please try again.