Measure Ellipsoidal Routes with geod and invgeod

You will use the installed PROJ 9.4.0 command to calculate a route between two latitude and longitude pairs, then use the result to find a destination from a starting point, azimuth and distance. The examples use standard input and change no files, so they take about 10 minutes to try.

1. Check the installed command and choose an ellipsoid

geod performs direct calculations. invgeod is the inverse-mode alias supplied by this installation. The command uses 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 manpage examples use +ellps=clrk66, a historical Clarke 1866 ellipsoid. That makes their result reproducible, but it is not a general recommendation for new work. Select the ellipsoid that belongs to your data and document that choice. List the names available on this host before guessing 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

Do not silently swap ellipsoids when comparing results. A plausible-looking number can still be wrong for the reference surface your coordinates require.

2. Solve the inverse problem between two points

The inverse problem starts with the latitude and longitude of the initial point, followed by the latitude and longitude of the terminal point. The -I option tells geod to use inverse mode. The input accepts decimal degrees or degrees, minutes and seconds. This example follows the installed manpage, using Boston and Portland and 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. Use the inverse result for a direct calculation

The direct problem takes an initial latitude and longitude, a forward azimuth and a distance. Feed the first inverse output back as the azimuth and distance, keeping 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 is close to the original Portland coordinate, but not character-for-character identical because the distance was rounded to three decimal places. More input precision gives a more precise destination.

For scripts, prefer a stable numeric format instead of parsing the default DMS text. The -f option controls geographic coordinate formatting and -F controls the distance format. First test the exact format your downstream program accepts; these are printf-style strings, not arbitrary labels.

4. Make units and precision explicit

Use geod -lu to see unit names accepted by this installation. The default is metres. The +units=us-mi setting in the examples is therefore significant, and omitting it changes the distance field's unit rather than merely changing 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, while distance defaults to three decimal places. -w changes the significant fractional seconds shown in coordinate output. -W does the same with constant-width fields and leading zeroes. These switches affect presentation; they do not improve the information in coordinates that were already rounded before you supplied them.

5. Process files or standard input without surprises

With no input file, geod reads standard input. A named file can be supplied instead, and multiple files are processed from left to right. A lone - explicitly names standard input. This makes a saved input suitable 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

This example creates a temporary input file only. Remove it when you have finished, or keep it if the coordinates are intentionally non-sensitive test data:

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

Do not put confidential coordinates in shell history, shared temporary files or copied transcripts. The command itself does not require elevated privileges. No sudo step is needed for calculations.

6. Generate points along a geodesic when needed

geod can also emit intermediate points. Supply 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 produces three intermediate points plus the two endpoints. Check this specialised mode with a small, known input before building a data pipeline around it. If your task is simply one route measurement, stop after the direct or inverse calculation.

7. Diagnose the common errors

• A non-zero result or an implausible distance usually means the input order, hemisphere, ellipsoid or unit is wrong. Re-run with one line of known test data.
• A latitude outside -90 to 90 degrees is invalid. Check whether longitude and latitude were accidentally swapped.
• If invgeod is absent on another machine, use geod -I. This is the documented fallback.
• If a format change breaks a parser, return to default DMS output and alter one option at a time. Do not mistake display precision for calculation precision.
• Do not treat a geodesic distance as a road route, travel time or local survey result. It follows the selected ellipsoid, not a transport network or every local datum correction.

These commands are read-only apart from the temporary file example. There is no configuration change to undo. If you created that file, the recovery action is the explicit rm -- /tmp/geod-route.txt command above.