Bivariate Choropleth Maps with Arcpy

In my previous post, I showed how to prepare the data for a bivariate choropleth map using PostGIS and QGIS. I also indicated that there is a website that shows an ArcGIS tool to do it. But, this actually turns into a good opportunity to illustrate some Python, and how to create the bivariate data using Arcpy.

Arcpy is certainly not as terse as SQL, but it does get the job done, and rather easily. We just have to think about the project a little differently. The code below is a Script tool that I created.

import arcpy, math, numpy
fc = arcpy.GetParameterAsText(0)

numrecs = int(arcpy.GetCount_management(fc).getOutput(0))


fields = arcpy.ListFields(fc, "bimode")
if len(fields) != 1:
    arcpy.AddField_management(fc, "bimode", "text", 3)

f1 = arcpy.GetParameterAsText(1)
f2 = arcpy.GetParameterAsText(2)
fields = ['bimode',f1,f2]

var1 = arcpy.UpdateCursor(fc, sort_fields=f1)

i=1
for row in var1:
    row.setValue("bimode",str(int(math.ceil((float(i) / float(numrecs)) * 3.0))))
    var1.updateRow(row)
    i=i+1

var2 = arcpy.UpdateCursor(fc, sort_fields=f2)

i=1
for row in var2:
    row.setValue("bimode",row.getValue("bimode") + "." + str(int(math.ceil((float(i) / float(numrecs)) * 3.0))))
    var2.updateRow(row)
    i=i+1
 

So let’s think about what we are doing. If I want to break the data up into 3 sections, I have to know how many records we have, sort the records, and then determine which group (1, 2, or 3) each record falls in. Let’s assume we have 10 records, all sorted. What group is the nth record in. The formula is actually easy:

Group = ceil(N / i * 3)
when i = 2: (2/10) * 3 = .6 = ceil(.6) = 1
when i = 7: (7/20) * # = 2.1 = ceil(2.1) = 3
 

So, you’ll see from the code that the first thing we do is get N,

numrecs = int(arcpy.GetCount_management(fc).getOutput(0))
 

Then, we create a recordset that allows updating the data, but we do it for one of the fields, and sort the cursor by that field:

var1 = arcpy.UpdateCursor(fc, sort_fields=f1)
 

Now, with the field sorted, we can loop through the cursor and apply our formula:

i=1
for row in var1:
    row.setValue("bimode",str(int(math.ceil((float(i) / float(numrecs)) * 3.0))))
    var1.updateRow(row)
    i=i+1
 

I do a little sleight of hand with converting the data to an Integer, and then to String, but that is mostly to get the data into the form I want. From there, you can symbolize your layer based on the bimode field.