Measure Ellipsoidal Routes with geod and invgeod

Work out the distance and bearings between two latitude and longitude pairs with geod, then run it backwards to find a destination. You will use the installed PROJ 9.4.0 command, feed it from standard input and change no files, so it takes about 10 minutes to try.

Two results to look for:

• Inverse calculation: two azimuths and a distance.
• Direct calculation: the destination latitude, longitude and back azimuth.

1. Check the command and choose an ellipsoid

geod performs direct calculations. invgeod is the inverse-mode alias this installation supplies. The command works on an ellipsoid, and the distance unit is normally the unit of that ellipsoid's major axis, which is metres for the usual definitions. The installed package is proj-bin 9.4.0-1build2.

$ geod --version
Rel. 9.4.0, March 1st, 2024
$ command -v invgeod
/usr/bin/invgeod

The man page examples use +ellps=clrk66, a historical Clarke 1866 ellipsoid. That makes their results reproducible, but it is no recommendation for new work. Pick the ellipsoid that belongs to your data and write that choice down. List the names available on this host before you guess one:

$ geod -le | head -5
    MERIT a=6378137.0      rf=298.257       MERIT 1983
    SGS85 a=6378136.0      rf=298.257       Soviet Geodetic System 85
    GRS80 a=6378137.0      rf=298.257222101 GRS 1980(IUGG, 1980)
    IAU76 a=6378140.0      rf=298.257       IAU 1976
    airy a=6377563.396    rf=299.3249646   Airy 1830

2. Solve the inverse problem

The inverse problem takes the latitude and longitude of the initial point, then those of the terminal point. The -I option switches geod to inverse mode. Input can be decimal degrees or degrees, minutes and seconds. This example follows the installed man page: Boston to Portland, in US statute miles.

$ printf "42d15'N 71d07'W 45d31'N 123d41'W\n" | \
    geod +ellps=clrk66 -I +units=us-mi
-66d31'50.141"    75d39'13.083"    2587.504

The fields are the forward azimuth from the initial point, the back azimuth from the terminal point, and the distance. The negative first azimuth is a signed DMS value. If your consumer expects azimuths from 0 to 360 degrees, add -p:

$ printf "42d15'N 71d07'W 45d31'N 123d41'W\n" | \
    geod +ellps=clrk66 -I -p +units=us-mi
293d28'9.859"    75d39'13.083"    2587.504

3. Run the result as a direct calculation

The direct problem takes an initial latitude and longitude, a forward azimuth and a distance. Feed the first inverse output back in as the azimuth and distance, with the same ellipsoid and units:

$ printf "42d15'N 71d07'W -66d31'50.141\" 2587.504\n" | \
    geod +ellps=clrk66 +units=us-mi
45d31'0.003"N    123d40'59.985"W    75d39'13.094"

The output is the destination latitude, destination longitude and back azimuth. It lands close to the original Portland coordinate but not character for character, because the distance was rounded to three decimal places. More input precision gives a more precise destination.

4. Set units and precision

Run geod -lu to see the unit names this installation accepts. The default is metres, so the +units=us-mi setting in the examples matters: leave it out and the distance field changes unit, not just its display.

$ geod -lu | head -10
          mm 0.001                millimetre
          cm 0.01                 centimetre
           m 1                    metre
          ft 0.3048               foot
     us-ft 0.304800609601219    US survey foot
        fath 1.8288               fathom
        kmi 1852                 nautical mile
      us-ch 20.1168402336805     US survey chain
      us-mi 1609.34721869444     US survey mile
          km 1000                 kilometre

Geographic output defaults to DMS, and distance defaults to three decimal places. Two switches adjust the coordinates:

• -w: changes the significant fractional seconds shown.
• -W: does the same with constant-width fields and leading zeroes.

Both affect presentation only. They cannot improve coordinates that were already rounded before you supplied them.

5. Read files or standard input

With no input file, geod reads standard input. You can name a file instead, and several files are processed from left to right. A lone - names standard input explicitly. That makes a saved input handy for a repeatable check:

$ printf "51d30'N 0d08'W 48d51'N 2d21'E\n" > /tmp/geod-route.txt
$ geod +ellps=WGS84 -I +units=km /tmp/geod-route.txt
147d59'57.761"    -30d5'34"    343.982

That creates a temporary input file and nothing else. Remove it when you have finished, or keep it if the coordinates are deliberately non-sensitive test data:

$ rm -- /tmp/geod-route.txt

6. Generate points along a geodesic

geod can also emit intermediate points. Give it the initial point with +lat_1 and +lon_1, then either a terminal point with +lat_2 and +lon_2, or a distance and azimuth with +S and +A. Use +n_S for a number of intermediate points or +del_S for a distance increment.

$ geod +ellps=WGS84 +lat_1=51.5 +lon_1=-0.133333 \
    +lat_2=48.85 +lon_2=2.35 +n_S=3
51d30'N    0d7'59.999"W
50d37'21.947"N    0d43'31.058"E
49d44'21.469"N    1d33'8.803"E
48d51'N    2d21'E

Here +n_S=3 gives the two endpoints plus two points between them: four lines in total. Check this specialised mode with a small, known input before you build a data pipeline on it. If all you need is one route measurement, stop after the direct or inverse calculation.

7. Diagnose the common errors

• Odd distance or non-zero result. Usually the input order, hemisphere, ellipsoid or unit is wrong. Re-run with one line of known test data.
• Latitude outside -90 to 90 degrees. That is invalid. Check whether longitude and latitude got swapped.
• No invgeod on another machine. Use geod -I, the documented fallback.
• A format change breaks a parser. Go back to default DMS output and change one option at a time. Display precision is not calculation precision.
• Reading it as a road route. A geodesic distance is not a road route, a travel time or a local survey result. It follows the selected ellipsoid, not a transport network or every local datum correction.

Apart from the temporary file, these commands are read-only, so there is no configuration change to undo. If you created that file, the recovery is the explicit rm -- /tmp/geod-route.txt command above.