Here is the python function that I came up with. Far from perfect, but seems to get the job done:
def plot_usgs_scale_bar(fig, region, proj, km_tick, mi_tick, xshift='w*0.6i', yshift='-0.1i'):
"""
Plots a USGS-style scale bar on a map figure, displaying distances in kilometers and miles. The function
utilizes map geometry and projection information to appropriately size and position the scale bar
relative to the map region. It also calculates appropriate scale tick lengths based on map distance
and projection. The scale bar is rendered with tick marks and labels both above and below the baseline.
:param fig: The PyGMT figure on which the scale bar will be plotted.
:type fig: pygmt.Figure
:param region: A list containing the geographic region covered by the map in the format
[west, east, south, north].
:type region: list[float]
:param proj: The map projection string used for the map, e.g., "M10i" or "X8i", where the last character
of the string indicates the projection unit.
:type proj: str
:param km_tick: A list of numerical tick distances in kilometers for the scale bar.
These will define the position of ticks above the baseline on the scale bar.
:type km_tick: list[float]
:param mi_tick: A list of numerical tick distances in miles for the scale bar.
These will define the position of ticks below the baseline on the scale bar.
:type mi_tick: list[float]
:param xshift: The horizontal shift of the scale bar relative to the map, specified in map units, defaulting
to 'w*0.6i' for positioning relative to the map’s width.
:type xshift: str, optional
:param yshift: The vertical shift of the scale bar relative to the map, specified in map units, defaulting
to '-0.1i' for positioning directly beneath the map.
:type yshift: str, optional
:return: The modified PyGMT figure with the scale bar added.
:rtype: pygmt.Figure
"""
units = proj[-1]
# This is a wrapper around a "mapproject" function to get the height and width of the map plot.
geometry = pygmtPlots.getTriPanelGeometry(proj, region, crossSectionWidth=0, buffer=0.0, units=units)
scale_width = geometry["width"] * 0.375
# Get length of map at center latitude
dist, _, _ = gps2dist_azimuth(region[2], region[0], region[2], region[1])
# Get inch/distance
inch_km = geometry["width"] / (dist / 1000)
inch_mi = geometry["width"] / ((dist / 1000) * 0.621371)
# Get appropriate tick lengths
km_tick_length = []
mi_tick_length = []
for tick in km_tick:
km_tick_length.append(tick * inch_km)
for tick in mi_tick:
mi_tick_length.append(tick * inch_mi)
# Shift proportionally below plot
with fig.shift_origin(xshift=xshift, yshift=yshift):
max_length = max(km_tick_length + mi_tick_length)
scale_proj = f'X{scale_width}i/{0.3 * scale_width}i'
region = [0, scale_width, -1, 1]
tick_top = 0.2
# Baseline
fig.plot(x=[0, max_length], y=[0, 0], pen='1p,black', projection=scale_proj, region=region)
# Zero ticks (top and bottom)
fig.plot(x=[0, 0], y=[0, tick_top], pen='1p,black', projection=scale_proj, region=region)
fig.text(text='0', x=0, y=tick_top, font='12p,ArialNarrow,black', justify='BC', offset='0p/5p', no_clip=True)
fig.plot(x=[0, 0], y=[0, -tick_top], pen='1p,black', projection=scale_proj, region=region)
fig.text(text='0', x=0, y=-tick_top, font='12p,ArialNarrow,black', justify='TC', offset='0p/-5p', no_clip=True)
# Loop over kilometers
for t, t_len in zip(km_tick, km_tick_length):
fig.plot(x=[t_len, t_len], y=[0, tick_top], pen='1p,black', projection=scale_proj, region=region)
if t == km_tick[-1]:
string = f'{t} KILOMETERS'
justify = 'BL'
else:
string = str(t)
justify = 'BC'
fig.text(text=string, x=t_len, y=tick_top, font='12p,ArialNarrow,black', justify=justify, offset='0p/5p',
no_clip=True)
# Loop over miles
for t, t_len in zip(mi_tick, mi_tick_length):
fig.plot(x=[t_len, t_len], y=[0, -tick_top], pen='1p,black', projection=scale_proj, region=region)
if t == mi_tick[-1]:
string = f'{t} MILES'
justify = 'TL'
else:
string = str(t)
justify = 'TC'
fig.text(text=string, x=t_len, y=-tick_top, font='12p,ArialNarrow,black', justify=justify, offset='0p/-5p',
no_clip=True)
return fig