Extending geofacet to algorithmically generate grids from lat/long data

This morning I learned about the very cool geofacet R package for making geographic facets with ggplot2. Unfortunately, my cities aren’t one of the existing grids. There is a neat web-based tool for arranging cities, but after 5 minutes, I was feeling frustrated by my lack of relative geographical knowledge (“Is Chicago really south of Boston? Seems wrong.”).

Anyway, I forked it and wrote a function for generating new grids based on lat/long. The resultant grids are intended to be (a) small and (b) not violate any cardinal ordering (e.g., a city that is south of one city will never shown as strictly north, though they can be on the same row).

To see how it works, suppose we have a collection of 50 cities w/ lat and long data (here plotted at their actual position):

If we have $n$ cities, we can always put them in an $n \times n$ grid such that all cities are in the correct relative position to each other (by sorting east-west and north-south):

Now we want to start squishing this matrix to make it smaller. What my algorithm does is pick the long side and then find columns/rows of that side it can combine without causing two cities to occupy the same cell. Every time it finds a pair, it combines them in removes a row or column, as the case may be. The algorithm keeps doing this until it cannot squish the matrix down any more in either direction. Here’s an example squished version of the plot from above:

As it runs quickly, I just brute force this some number of times and take the most compact one. There’s definitely more that could be done, say relaxing some constraints for greater compactness. Even this one, I’d hand-tune it a bit, but it was already much better than what I could generate by hand, starting from scratch.

Here’s my code for this example: