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)

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)

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)

fig = plt.figure(figsize=(12,6), dpi=200)

ax = fig.add_subplot(1, 1, 1, aspect=1.0)

map = Basemap(

projection='merc',

llcrnrlat=swcorner[0],

llcrnrlon=swcorner[1],

urcrnrlat=necorner[0],

urcrnrlon=necorner[1],

resolution='l',

ax=ax

)

projection='merc',

llcrnrlat=swcorner[0],

llcrnrlon=swcorner[1],

urcrnrlat=necorner[0],

urcrnrlon=necorner[1],

resolution='l',

ax=ax

)

map.drawmapboundary()

map.fillcontinents()

map.drawcoastlines()

map.drawstates()

map.drawcountries()

map.fillcontinents()

map.drawcoastlines()

map.drawstates()

map.drawcountries()

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

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.

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)

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.

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.

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)

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

lons3 = np.ones(L) * urlon

lats4 = np.ones(L) * lllat

lons4 = np.flip(np.linspace(lllon, urlon, L))

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)

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.

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)

ax.add_patch(poly)

ax.add_patch(poly)

The result was more like what I was shooting for.

Two final flourishes remain. First, to add the grid field borders and labels:

# DRAW 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.)

ax.add_patch(poly)

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.)

ax.add_patch(poly)

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.

All that's left are the individual grid square labels. This is also pretty straightforward. I replaced my grid square loop as follows:

# PLOT GRID SQUARES

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)

ax.add_patch(poly)

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)

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)

ax.add_patch(poly)

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)

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

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

from mpl_toolkits.basemap import Basemap