Azimuthal Equidistant map - gridline placement, labeling magnetic declination, adjusting position

johnm · June 9, 2025, 9:21pm

Hi All,

Below is a shell script I wrote to generate azimuthal equidistant maps given a particular latitude and longitude. Couple Questions:

Setting -Byg20 in the psbasemap seems to place the first line of latitude at 10° out from center, and subsequent lines in 20° increments, whereas setting -Byg10 places the first line at 10° and subsequent lines at 10°. How to place the first circle out from center at 20°? Here’s a screenshot using -Byg20:

How to include a magnetic declination diagram?
How to adjust position on page (e.g. left margin)?

Also, I’d like to come up with a more comprehensive list of city names and coordinates than what I put together manually in the cities.txt file pasting lon/lat values from Wikipedia. Is there a file somewhere or a recommended approach to accurately compile city names & coordinates?

Thanks in advance for any suggestions, corrections or pointers.

Cheers,

John

#!/bin/bash
#
# A script utilizing GMT ( https://www.generic-mapping-tools.org/ ) to generate an
# Azimuthal Equidistant projection map given latitude and longitude coordinates.
# With a text file listing City names and coordinate, locations of cities will be shown.
#
# I used this tutorial as a guide:
# http://www.ceri.memphis.edu/people/egdaub/datanotes/_build/html/gmt2.html#a-more-complicated-map
#
# See also:
# https://docs.generic-mapping-tools.org/6.5/reference/map-projections.html#azimuthal-projections
# 
# Last edit: 2025/06/06 Fri 10:26 -0700 
# Copyright © 2025 John Magolske | released under the MIT license

usage="Usage:
    `basename $0` \"latitude\" \"longitude\" [ letter, tabloid ]
will generate an Azimuthal Equidistant projection map centered on location specified by
\"latitude\" and \"longitude\". Values are entered in numerical degrees and quoted so negative
numbers aren't interpreted as options. e.g. Buenos Aires is located lat -34.603, lon 301.618
so the command would be: `basename $0` \"-34.603\" \"301.618\""

[ "$#" -eq 0 ] && echo "$usage" && exit 0

# Set default media to letter size with option to specify tabloid size.
# To add additional options, see PS_MEDIA parameter in gmt.conf documentation.
if [ -z "$3" ] || [ "$3" = "letter" ] 
then
    media="letter"
    width="6.5i"
elif [ "$3" = "tabloid" ]
then
    media="tabloid"
    width="9.5i"
else
    echo " \"$3\" is not a valid media option, consider editing the `basename $0` script to add more options."
    exit 1
fi

# Set size of media
gmt set PS_MEDIA $media

lat="$1"
lon="$2"
outfile="AziEqui_$1_$2.ps"

# create enclosing circle, 6.5 inch wide, label Azimuthal Equidistant with coordinates
gmt psbasemap -R0/360/-90/0 -JS0/-90/90/$width -Bxa30 -B+t"Azimuthal Equidistant"+s"$1@., $2@." -K -P > $outfile

# generate an Azimuthal Equidistant projection map, horizon at 175° (5° shy of antipode),
# 6.5 inches wide, permanent major rivers blue, shorelines black, shade landmass light gray
gmt pscoast -Rd -JE$lon/$lat/175/$width -Df -A50 -I1/.01p,blue -W0.02p,black -Glightgray -K -O -P  >> $outfile

# add gridlines at every 10° longitude & every 20° latitude.
gmt psbasemap -Rg -JE0/-90/175/$width -Bxg10 -Byg20 -K -O -P >> $outfile

# place triangular geometric symbol at city locations listed in cities.txt file
gmt psxy -Rd -JE$lon/$lat/175/$width -St.03i -H0.5 cities.txt -K -O -P >> $outfile

# place names for cities listed in cities.txt file
gmt pstext -Rd -JE$lon/$lat/175/$width -F+f4p,Helvetica,black+j,right cities.txt -O -P >> $outfile

# convert postscript file to pdf
ps2pdf "$outfile"

PlanetGus · June 10, 2025, 7:45am

Hi there,

You can use a “phase” shift : -byg20+10 : this will shift by 10°. So 20, 40, 60, … Or you can specify a file with the exact value you want.
You mean add another layer or just a rose ? First case : you need the file and use plot. Second case there’s an option implemented in coast (and several other modules) : -T if my memory is correct.
You can shift the whole plot using -X2c to move it to the right by 2cm. Though, using modern mode, the figure will be adjusted anyway giving a 0 shift apparent result. The workaround is given in another thread using psconvert.

Best,

johnm · June 13, 2025, 7:28pm

Thanks – your reply answered my questions.

Using -byg20+10 bumped that first line of latitude out from 10° to 20°. I was able to create a magnetic declination rose using the -T option with psbasemap and position the plot using the -X and -Y options. See below for the updated script.

A partial screenshot of the output:

One thing I’d still like to do is create a magnetic declination diagram something a bit more like those found on topo maps:

…will have to look into that when I get more time.

Regards,

John

The cities.txt file used by the script below: cities.txt (772 Bytes)

#!/bin/bash
#
# A script utilizing GMT ( https://www.generic-mapping-tools.org/ ) to generate an
# Azimuthal Equidistant projection map given latitude and longitude coordinates.
# With a text file listing City names and coordinate, locations of cities will be shown.
#
# Used this tutorial as a guide:
# http://www.ceri.memphis.edu/people/egdaub/datanotes/_build/html/gmt2.html#a-more-complicated-map
#
# See also:
# https://docs.generic-mapping-tools.org/6.5/reference/map-projections.html#azimuthal-projections
# 
# Last edit: 2025/06/11 Wed 07:23 -0700 
# Copyright © 2025 John Magolske | released under the MIT license

usage="USAGE
    `basename $0` [-m {letter|tabloid}] [-d <degrees>] <latitude> <longitude>
Generates an Azimuthal Equidistant projection map centered on location specified by <latitude> and <longitude>.
DESCRIPTION
    -m <media>
        Enter media size. Default is letter, accepted values are letter and tabloid.
    -d <degrees>
        Magnetic declination in decimal degrees used to create magnetic rose. Note that
        declination East values are positive numbers and declination West values are negative.
    <latitude> <longitude>
        Coordinates in decimal degrees of map center. Note that latitude North values are
        positive, latitude South values are negative, longitude East values are positive,
        longitude West values are negative
    -h Print this usage command and exit.
EXAMPLE
    San Francisco is at 37.778° North latitude, 122.416° West longitude with a magnetic
    declination of 12.95° East. To Print this to tabloid size the command would be:
    `basename $0` -m tabloid -d 12.95 37.778 -122.416"

# Process arguments manually to handle negative coordinate values
temp_args=()
while [[ $# -gt 0 ]]; do
    case $1 in
        -m) media="$2" ; shift 2 ;;
        -d) decl="$2" ; shift 2 ;;
        -h) echo "$usage" ; exit 0 ;;
        -*) # Check if this looks like a negative number
            if [[ $1 =~ ^-[0-9]+\.?[0-9]*$ ]]; then
                temp_args+=("$1") ; shift
            else
                echo -e "Error: unknown option $1 \n $usage"
                exit 1
            fi ;;
        *) temp_args+=("$1") ; shift ;;
    esac
done

# Restore positional parameters
set -- "${temp_args[@]}"

# Check for required arguments
[ "$#" -eq 0 ] && echo "$usage" && exit 0
[ "$#" -ne 2 ] && echo -e "Wrong number of arguments \n $usage" && exit 1

# Set default media to letter size with option to specify tabloid size.
# To add additional options see PS_MEDIA parameter in gmt.conf documentation.
[ "$media" = "" ] && media="letter"

# Set options to adjust size and position of plot given media size
if [ "$media" = "letter" ]; then
    width="6.5i"
    d_x="-X2.7i"
    d_y="-Y3.8i"
    p_x=""
    p_y="-Y1.25i"
elif [ "$media" = "tabloid" ]; then
    width="9.5i"
    d_x="-X3.6i"
    d_y="-Y5.3i"
    p_x="-X0.75i"
    p_y="-Y4.25i"
else
    echo "\"$media\" is not a valid media option, consider editing the `basename $0` script to add more options."
    exit 1
fi

# Set size of media
gmt set PS_MEDIA $media

# variables for latitude and longitude
lat="$1"
lon="$2"

# variable for output file
outfile="AziEqui_${lat}_${lon}.ps"

# create enclosing circle, 6.5 inch wide, label Azimuthal Equidistant with coordinates
gmt psbasemap -R0/360/-90/0 -JS0/-90/90/$width -Bxa30 -B+t"Azimuthal Equidistant"+s"$lat@., $lon@." $p_x $p_y -K -P > "$outfile"

# generate an Azimuthal Equidistant projection map, horizon at 175° (5° shy of antipode),
# permanent major rivers blue, shorelines black, shade landmass light gray
gmt pscoast -Rd -JE$lon/$lat/175/$width -Df -A50 -I1/.01p,blue -W0.02p,black -Glightgray -K -O -P >> "$outfile"

# add gridlines at every 10° longitude & every 20° latitude.
gmt psbasemap -Rg -JE0/-90/175/$width -Bxg10 -Byg20+10 -K -O -P >> "$outfile"

# draw triangular geometric symbol at city locations listed in cities.txt file
gmt psxy -Rd -JE$lon/$lat/175/$width -St.03i -H0.5 -G0 cities.txt -K -O -P >> "$outfile"

# draw names for cities listed in cities.txt file
gmt pstext -Rd -JE$lon/$lat/175/$width -F+f4p,Helvetica,black+j,right cities.txt -K -O -P >> "$outfile"

# draw magnetic declination rose
[ "$decl" != "" ] && gmt psbasemap -Rg -JE$lon/$lat/175/$width \
 -Tmg$lon/$lat+w2.4i+d$decl+t360/90/90/9450/90/90+i0.5p,blue+l,,,* $d_x $d_y -O -P >> "$outfile"

# convert postscript file to pdf
ps2pdf "$outfile"

Andreas · June 14, 2025, 9:30pm

Don’t have anything to contribute with, but I really like your ideas/website

PlanetGus · June 15, 2025, 7:37am

Well for that, I can only suggest to use plot (eventually inside an inset environment).

Joaquim · June 15, 2025, 6:08pm

I guess that a custom symbol could be made as well.

johnm · September 10, 2025, 4:38pm

Thank you all for the helpful replies. Got busy and forgot to mention that I put the script along with example PDFs up on this page: https://earthdirections.org/ideas/map

– John

PlanetGus · September 11, 2025, 7:46am

Fun website, I liked it