#!/bin/bash
#

# A script used to generate Azimuthal Equidistant projection maps centered on location
# given latitude and longitude coordinates. Used in conjunction with a text file listing
# city names and coordinate, the locations of cities will be shown. The cities.txt file
# must be copied / symlinked into the working directory where the command is issued and the
# PDF is generated. This utilizes GMT, which must be installed on your Linux/Unix system:
# https://www.generic-mapping-tools.org/
#
# Generate one of these maps with your location at center and use it to find the path to
# another location on Earth from your location -- both the "as the bird flies" surface
# distance and the "as the neutrino flies" straight-line distance through the Earth.
#
# I used this tutorial as a guide in creating this script:
# 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/25 Wed 10:57 -0700 
# Copyright © 2025 John Magolske | released under the MIT license

usage="USAGE
    `basename $0` [-m {letter|tabloid|a4|A4|a3|A3}] [-d <degrees>] [-c] <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, tabloid, a4, a3.
    -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.
    -c  Do not draw cities (proceed without cities.txt file)
    -h Print this usage command and exit.
    <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
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"

# Initialize variables
skip_cities=false

# 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 ;;
        -c) skip_cities=true ; shift ;;
        -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

# Check if cities.txt file exists
if [ "$skip_cities" = false ] && [ ! -e cities.txt ]; then
    echo -e "cities.txt file does not exist. Use -c option to generate map without cities labeled."
    exit 1
fi

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

# Set options to adjust size and position of plot elements given media size
if [ "$media" = "letter" ]; then
    width="6.5i"
    d_x="-X3.1i"
    d_y="-Y3.2i"
    p_x=""
    p_y="-Y2.25i"
    t_x="-X0.0i"
    t_y="-Y-3.8i"
elif [ "$media" = "tabloid" ]; then
    width="9.5i"
    d_x="-X3.6i"
    d_y="-Y5.3i"
    p_x="-X0.75i"
    p_y="-Y4.25i"
    t_x="-X1.1i"
    t_y="-Y-6.3i"
elif [ "$media" = "a4" -o "$media" = "A4" ] ; then
    width="6.5i"
    d_x="-X2.75i"
    d_y="-Y3.5i"
    p_x="-X0.89i"
    p_y="-Y2.5i"
    t_x="-X0.48i"
    t_y="-Y-4.3i"
elif [ "$media" = "a3" -o "$media" = "A3" ]; then
    width="9.8i"
    d_x="-X3.6i"
    d_y="-Y5.3i"
    p_x="-X.95i"
    p_y="-Y3.8i"
    t_x="-X1.4i"
    t_y="-Y-6.3i"
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}_$media.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"

if [ "$skip_cities" = false ]; then
    # 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"
fi

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

# Print explanatory text below map
{ echo "0 0 Use this map to determine the straight-line direction to other locations on Earth while situated at"
  echo "0 -0.2 $lat@., $lon@.. With the map oriented relative to North, observe the azimuthal heading to"
  echo "0 -0.4 another location. While maintaining this azimuthal heading, look down at the angle of declination"
  echo "0 -0.6 below horizontal to visualize the distant location; each concentric ring represents a 10@. increment"
  echo "0 -0.8 of declination, so a distant city on the 3rd circle out would be 30@. below the horizon."
} | gmt pstext -R0/1/0/1 -JX1i -F+f11p,+jCM $t_x $t_y -N -K -O -P >> "$outfile"

{ echo "0 -1.1 https://earthdirections.org/ideas/map"
} | gmt pstext -R0/1/0/1 -JX1i -F+f13p,Courier+jCM -N -O -P >> "$outfile"

# convert postscript file to pdf
ps2pdf "$outfile"

# Delete postscript file
rm "$outfile"
