Easy bivariate map with Postgres

A friend recently asked me about the cool looking bivariate maps produced with ArcOnline, lamenting that the capability seemed lacking in ArcGIS. Well, it turns out that ESRI has a .dll you can use, and there is a good article here. So, if you want to create these great looking maps in ArcGIS, it shouldn’t be a problem. The website will allow you to download the .dll, and you can also watch the video on how to use it. Well worth the time.

So of course that got me thinking: could we make the same map with spatial SQL. Well, sure, and it is super easy. If you want to get spun up on what these 9 color bivariate choropleth maps are, and the theory behind it, have a look at this great site. Josh does a great job explaining how this works, but it is a little cumbersome if you want to pump out map after map. But, with a very little bit of SQL, you can easily pull it off.

Let’s start with our data: I have a Postgres table of United States County boundaries, with the attributes percentobese and percentdiabetes, along with a FIPS code and a geometry column.

To prepare the data for the map, we simply issue this SQL query:

SELECT 
       fips,geom,
       ntile(3) over (order by percentobese) || '.' ||
       ntile(3) over (order by percentdiabetes) AS bimode
  INTO qlayer
  FROM ushealthrisk 

yes, that’s it. Really.

I’m not kidding. We’re done, folks. Go home, nothing left to see.

Well, if you insist on reading, I’ll tell you what the SQL does, and how to actually visualize the data.

The real magic here is in the ntile function:

      ntile(3) over (order by percentobese) || '.' || 

the ntile function ntile(num_buckets integer) returns an integer ranging from 1 to the argument value (in our case the number breaks, dividing the partition as equally as possible. Therefore, when we use 3, we are breaking our data into 3 sections, and the lowest third is labeled 1, the middle third is labeled 2, and the highest third is labeled 3. Ordering by percent obese ensures that our lowest obese records are in group 1, and the highest are in group 3 like you can see here:

Similarly, we do the same thing for percentdiabetes, and combine both results into a single field where we separate the group numbers by a decimal point:

So now we know which quadrat each record falls in as Josh demonstrated with his 3×3 matrix :

Josh was also great to provide a number of palettes with the color codes associated with them:

and I used one of them to create a symbology in QGIS which I saved as a .SLD:

and that yielded our final map:

I’m going to get one of my students to create the rest of.SLD files, and if anyone (maybe Josh) wanted to create a few more example matrices, I can get the student to create those.

So there you have it:

1. with a very short SQL statement, you’ll have the 9 categories for the two-variable map.

2. reading the layer into QGIS and assigning a .SLD symbology to the file will give you the nice display.

3. with a little more effort you can convert my SQL query to a Function, and simply pass it information to continue pumping out more maps.