Maidenhead Grid Square Maps using Python
About a week ago, I started working FT8 on 40 meters. It’s been fun making contacts all over the country, and I’ve even managed a few DX contacts in Europe, South America, and the Pacific. As the grid squares started piling up, I started tracking them by coloring in the squares on a sheet of paper.
Just like I would have done 40 years ago. What the heck was I thinking? There must be an easier way. I’d previously used the matplotlib toolkit Basemap, so I thought I might be able to scrape something together.
The first problem to solve was how to pull in my QSO data. After installing the adif_io package I was able to load the WSJT-X ADIF log file into a Pandas DataFrame on the first attempt!
qsos_raw, adif_header = adif_io.read_from_file("/Full/Path/wsjtx_log.adi")
qso_df = pd.DataFrame(qsos_raw)

With the log data in hand, the next step was to draw a basic map. Here is the code for doing a barebones map:
# Define the map boundaries (lat, lon)
swcorner = (24, -126)
necorner = (50, -64)
# Create the figure
fig = plt.figure(figsize=(12,6), dpi=200)
ax = fig.add_subplot(1, 1, 1, aspect=1.0)
map = Basemap(
Next, I needed a way to shade in the grid squares that appear in the log. Iterating over the map area and drawing 2º × 1º rectangles was very easy.
for gs_lon in np.arange(swcorner[1], necorner[1], 2):
    for gs_lat in np.arange(swcorner[0], necorner[0], 1):
        cx_lat, cx_lon = gs_lat+.5, gs_lon+1. # compute grid square center
        gridid = mh.to_maiden(cx_lat, cx_lon, precision=2) # get grid id
        if np.count_nonzero(qso_df['GRIDSQUARE']==gridid) > 0:
            x1, y1 = map(gs_lon, gs_lat)
            x2, y2 = map(gs_lon+2., gs_lat+1.)
            rect = Polygon([(x1,y1), (x1,y2), (x2,y2), (x2, y1)])
            ax.add_patch(rect) # add to the plot
The result was a modest success.

Fig. 1. My first attempt using the Mercator projection

Not bad for a first attempt, but not without some problems. There are no labels or boundary lines, so it's hard to tell which grid squares have been worked and which are still needed. There's also the matter of the Mercator projection. Although it is (or once was) useful for maritime navigation, it grotesquely distorts the shape and—more importantly—size of land masses. In my judgement, it’s just not an attractive projection to use.
I updated the code to use the Cassini projection, in which both meridians and parallels are curved rather than straight lines. Here is the code to add meridians and parallels:
map.drawparallels(np.arange(-90,90,1.), color='#CC0000', linewidth=0.1)
map.drawmeridians(np.arange(-180,180,2.), color='#CC0000', linewidth=0.1)
I was lazy and added them for the entire globe. That way, I am sure to get the grid lines edge to edge. I don't have the problem of lines appearing outside the plot area because Basemap takes care of that for me. It adds a few seconds to the time it takes to generate the map, but I'm okay with that.
I supposed these changes would make my map a little more pleasing to the eye. And so it did.

Fig. 2. Grid squares are misshapen after switching to the Cassini projection.

But I introduced a couple of new problems. First and most obvious is the misalignment and distortion of the shaded grid squares. My code correctly translates the coordinates of the lower-left and upper-right corners, but still plots the grid squares as rectangles in the plot space. This problem is more clearly visible when applied to a larger area, such as a Maidenhead Grid Field.

Fig. 3. The grid problem: attempting to plot rectangular spaces in the map space using coordinates in the plot space gives incorrect results.

To solve this problem, I wrote a function that interpolates points all along the edges of the rectangles in the map space (latitude and longitude) and maps them into the plot space (X and Y coordinates). The function returns a complex polygon in the plot space that can be inserted into the plot.
def draw_map_rect(lllat, lllon, urlat, urlon, m, fc='#CC0000', ec=None, a=0.25, lw=0):
    L = 10
    lats1 = np.linspace(lllat, urlat, L)
    lons1 = np.ones(L) * lllon
    lats2 = np.ones(L) * urlat
    lons2 = np.linspace(lllon, urlon, L)
    lats3 = np.flip(np.linspace(lllat, urlat, L))
    lons3 = np.ones(L) * urlon
    lats4 = np.ones(L) * lllat
    lons4 = np.flip(np.linspace(lllon, urlon, L))
    lats = np.concatenate((lats1,lats2,lats3,lats4))
    lons = np.concatenate((lons1,lons2,lons3,lons4))
    XX = np.stack(m(lons,lats), axis=-1)
    return Polygon(XX, edgecolor=ec, facecolor=fc, alpha=a, linewidth=lw)
Notice also that I included a few additional parameters that can be used to tweak the colors and opacity.

Fig. 4. The grid problem corrected

I replaced the Polygon(…) and add_patch(…) lines in the code above with the following:
poly = draw_map_rect(gs_lat, gs_lon, gs_lat+1, gs_lon+2, map)
The result was more like what I was shooting for.

Fig. 5. Grid squares are now correctly shaded.

Two final flourishes remain. First, to add the grid field borders and labels:
map_ctr_lon = (necorner[1]+swcorner[1])/2
for gs_lon in np.arange(np.floor((map_ctr_lon-45)/20)*20, np.ceil((map_ctr_lon+45)/20)*20, 20):
    for gs_lat in np.arange(-80, 80, 10):
        cx_lat, cx_lon = gs_lat+5., gs_lon+10.
        cx, cy = map(cx_lon, cx_lat)
        gridid = mh.to_maiden(cx_lat, cx_lon, precision=1)
        plt.text(cx, cy, gridid, ha='center', va='center', color='k', alpha=0.25, fontsize=40, fontweight='bold').set_clip_on(True)
        poly = draw_map_rect(gs_lat, gs_lon, gs_lat+10, gs_lon+20, map, fc='None', ec='#FFCC00', lw=2, a=1.)

Of particular interest here is the use of the set_clip_on(True) call. Without it, those labels extending beyond the borders of the map area will still be included in the plot. Trust me, you don't want that. The result is very satisfying.

Fig. 6. Grid Fields are identified.

All that's left are the individual grid square labels. This is also pretty straightforward. I replaced my grid square loop as follows:
for swlat1 in np.arange(-90, 90, 10):
  for swlon1 in np.arange(-180, 180, 20):
    cx_lat1, cx_lon1 = swlat1 + 5, swlon1 + 10
    grid1 = mh.to_maiden(cx_lat1, cx_lon1, precision=1)
    for gx in range(10):
      swlon2 = swlon1 + (gx*2.)
      for gy in range(10):
        swlat2 = swlat1 + gy
        grid2 = '{}{}'.format(gx, gy)
        gridid = '{}{}'.format(grid1, grid2)
        if np.count_nonzero(qso_df['GRIDSQUARE']==gridid) > 0:
          poly = draw_map_rect(swlat2, swlon2, swlat2+1., swlon2+2., map)
        cx, cy = map(swlon2+1., swlat2+0.5)
        plt.text(cx, cy, grid2, ha='center', va='center', color='black', fontsize=4).set_clip_on(True)

Fig. 7. Mostly complete map. The colors and map boundaries need a little tweaking.

All the includes needed for the code above are here:
import adif_io
import maidenhead as mh
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from matplotlib.patches import Polygon
from mpl_toolkits.basemap import Basemap

You may also like