# How can I use Mathematica to sort the US Congressional Districts by Perimeter÷Area

The United States is divided into 435 congressional districts. Some of their shapes can be really strange due to gerrymandering. I think that I can find very gerrymandered congressional districts if I sort them by the ratio of the perimeter/area.

The first step would be to import the perimeters and areas of all the congressional districts. How can I obtain this data directly from Mathematica?

• Try Entity["USCongressionalDistrict", "Georgia:District9"]["Polygon"] – Michael E2 Jul 3 '16 at 20:44

The US congressional districts are available as geographic entities in the form

Entity["USCongressionalDistrict", "Illinois:District7"]


As such, in principle you can enumerate them using

GeoEntities[Entity["Country", "UnitedStates"], "USCongressionalDistrict"]


but this is liable to timing out because you're asking for a lot of data, so instead you can break down the query into per-state lookups, giving you something like

Length[
districts = Flatten[GeoEntities[#, "USCongressionalDistrict"] & /@
GeoEntities[Entity["Country", "UnitedStates"], "USState"]]
]
(*449*)


(Or, as Michael pointed out, using the built-in EntityList["USCongressionalDistrict"] is way better.)

To get the data you want, you can use GeoArea. To get the perimeter you can apply GeoDistance to the data but you need to fish it out first. You can get the shape of the district out using EntityValue[#, "Polygon"] and this will return things of the form

Polygon[ GeoPosition[ {{lat1, lon1}, {lat2, lon2}, …, } ] ]

and you can't run GeoDistance on the Polygon - only on the GeoPosition object inside it or the lists of latitudes and longitudes inside that. Thus you can use

perimeter[polygon_] := GeoDistance[polygon[[1]]]


and to match it up you can define

area[district_] := EntityValue[district, "Area"]


To do the calculation, then, you can run

polygons = EntityValue[districts, "Polygon"]


which is relatively fast.

There is a slight complication with the data, though, and it's that not all districts appear as a single polygon. In fact, if you define districts = EntityList["USCongressionalDistrict"] and polygons as above, then there are

Select[polygons, Length[#] > 1 &] // Length


= 15 districts with multiple polygons. This is a bit surprising, since US congressional districts are meant to be territorially contiguous (and hence the gerrymandering stuff), but 11 of these have islands, which is fairly reasonable. Nevertheless, there are districts like the Tennessee 2nd district which are not, in fact, contiguous, even though their exclaves are separated only by land. To visualize these 15 districts, use

multiPolygonDistricts = Pick[districts, Length[#] > 1 & /@ polygons]
GeoGraphics[EntityValue[#, "Polygon"], PlotLabel -> #,
ImageSize -> {{250}, {250}}] & /@ multiPolygonDistricts


In addition to this, there are some districts which don't have polygon data available, such as

EntityValue[Entity["USCongressionalDistrict","Pennsylvania:District19"],"Polygon"]
(*Missing["NotAvailable"]*)


There are 12 such districts and you can find them using

Pick[districts, MissingQ /@ polygons]


These turn out to be historical districts (like the aforementioned PA 19th) that are no longer in use, and as of November 2016 and MM v11 they are no longer present in the data.

There are also some irregularities in the formatting of the polygon data, so for example

EntityValue[Entity["USCongressionalDistrict", "Tennessee:District9"], "Polygon"]


has an output of the form {Polygon[_]}, but

EntityValue[Entity["USCongressionalDistrict", "Texas:District1"], "Polygon"]


has an output of the form {{Polygon[_]}} with a deeper list, whereas

EntityValue[Entity["USCongressionalDistrict", "Alaska:District0"], "Polygon"]


returns a straight Polygon[_], though its insides are all messed up (i.e. it contains multiple polygons inside it, instead of a list of polygons, which will definitely mess with the perimeter calculator above); either way, if you're studying gerrymandering you presumably want to leave out this district as it's the entirety of Alaska. Again, this is buggy data which I'll leave to you to report.

EDIT: some issues with the data seem to have been resolved, but it's not necessarily all good yet. It's also not clear to me whether the fix is version-dependent (or fixed upstream at the server), either. So: watch out for this if you want to implement this.

To get clean data, then, you can do

cleanPolygons = Flatten /@ DeleteCases[DeleteMissing[polygons], Polygon[_]];
cleanDistricts = DeleteCases[Pick[districts, Not@*MissingQ /@ polygons],


You can now calculate the perimeters using

perimeters = Total /@ Map[perimeter, cleanPolygons, {2}]


which is now quite fast as well (the previous implementation was slowing the whole thing down). Calculating the areas as

areas = EntityValue[cleanDistricts, "Area"]


you finally have enough to get your rankings: doing

TableForm[
(data = SortBy[{EntityValue[cleanDistricts, "Name"], perimeters^2/areas
}\[Transpose], Minus@*Last])[[1 ;; 5]]
]


gets you for the top five

District                                        p^2/A
"North Carolina 12th Congressional District"    365.611
"Florida 5th Congressional District"            244.518
"North Carolina 1st Congressional District"     234.952
"North Carolina 4th Congressional District"     228.47
"Maryland 3rd Congressional District"           226.788


And, indeed, you get some pretty heavy gerrymandering going on:

As an interesting trivia, note that the fourth place on this list, Maryland's 3rd district, really stretches the definition for "contiguous" - if you zoom in on the map into the Baltimore harbour, there's not a lot of district land connecting Federal Hill and Little Italy.

It's also important to note that it matters what year the data is from: Mathematica seems to be using the 2014 electoral districts, which have since changed in some states (North Carolina in particular, but also Florida and Virginia) precisely because they were deemed to be too gerrymandered. For a good illustration of the changes see this NPR animation for the NC redistricting.

• See @Michael E2 answer for the suggestion of using perimeter^2 / area (so that the answer is dependent only on the shape and not the size) – mikado Jul 3 '16 at 21:52
• @mikado That's a different quantity that will yield a different ranking, and it's for the OP to decide. Having a dimensional answer is not particularly problematic, as the ranking will be the same regardless of the unit system in use. – Emilio Pisanty Jul 3 '16 at 22:17
• @EmilioPisanty The ranking could be different as they are different calculations, but I agree about leaving it up to the OP to ask. Just a side note, you shouldn't have to write Apply[x] /@ list when x@@@list does the same thing. – Jason B. Jul 4 '16 at 1:10
• I think @MichaelE2 's suggestion is not concerning unit system, but rather since area scales like $[\text{length}]^2$ and perimeter like $[\text{length}]^1$, my original ratio perimeter/area will simply yield smaller numbers for large districts with exactly the same shape as small districts. This may cause me to inadvertently miss heavily gerrymandered big districts. The dimensionless ratio perimeter^2/area is definitely better in that regard. – QuantumDot Jul 4 '16 at 11:04
• Too bad I can't give another +1 for the NPR link… – J. M.'s technical difficulties Jul 4 '16 at 18:13

First (update: I guess the missing ones should be deleted),

alldistricts = DeleteMissing@ EntityList["USCongressionalDistrict"];
allpolygons = #["Polygon"] & /@ alldistricts (* // DeleteMissing *);


Then find the GeoArea[poly] and perimeter* of them. Note that you might want to compute a dimensionless quantity such as perimeter^2/area.

Update: The perimeter was a little tricker than my original blithe remark led one to believe. Here is an improvement over GeoDistance, at least in terms of reliability:

perimeters = Total[allpolygons /. Polygon :> GeoLength@*GeoPath, {2}]
areas = Total[GeoArea /@ allpolygons, {2}]  (* why not put in the areas, too *)

• I'm not sure GeoDistance works here - trying GeoDistance[Entity[ "USCongressionalDistrict", "Alabama:District1"][ "Polygon"]] returns zero. – Emilio Pisanty Jul 4 '16 at 10:13
• @EmilioPisanty GeoDistance[] works on the points of the polygon, though. I guess I should clarify. – Michael E2 Jul 4 '16 at 17:26
• I agree. As is, GeoDistance[poly] produces either zeroes or an error. Do run it through to the end for the full list - the data is much less homogeneous than you'd assume at first brush. – Emilio Pisanty Jul 4 '16 at 17:26
• @EmilioPisanty Surely the OP can deal with the problematic data. SE isn't a consulting service, and I don't get paid for participating. :) – Michael E2 Jul 4 '16 at 17:42
• That's up to you. But as is, your answer doesn't get the OP anywhere closer to their stated goal. None of the polygons from your method yield useful answers with your code. – Emilio Pisanty Jul 4 '16 at 17:44