Geodesics and planar approximations using PROJ and geod

The PROJ system supports computations with Geodesic calculations, through API calls or through the geod command line interface. This exercise will proceed mostly by having you copy sample commands from the text onto the command line, and noting the results.

The tutorial is written with the Windows command line in mind but it should be straight forward to use it on Unix-based systems as well.

Hint

Users of Unix-based systems can replace set with export and %VAR% with $VAR to follow the tutorial.

1. Objective

The objective of this exercise is to investigate how well a straight line in a projected plane approximates the geodesic between the same end points.

2. The geodesic

The geodesic we will consider runs from Helsinki to Tallinn, using the following coordinates:

Helsinki  60.171 N 24.938 E
Tallinn   59.437 N 24.745 E

So to save some typing, let's define some environment variables. On Windows:

set HEL="60.171  24.938"
set TAL="59.437  24.745"

export HEL="60.171  24.938"
export TAL="59.437  24.745"

(i.e. mark, copy, and paste these lines into the command prompt)

First, we want to solve the "inverse geodetic problem" for the geodesic, i.e. finding the azimuth and distance between the two points:

echo %HEL% %TAL% | geod -I +ellps=GRS80

where the -I option indicates "the inverse geodetic problem".

geod replies with these 3 numbers:

-172d22'12.772"    7d27'46.693"    82488.500

The first is the forward azimuth from Helsinki to Tallinn, the second is the return azimuth, while the third is the distance in meters.

We would actually rather have the azimuths in fractional degrees, so we provide an explicit output format, using the -f option to geod:

echo %HEL% %TAL% | geod -I -f "%.12f" +ellps=GRS80

getting us:

-172.370214337896    7.462970382137    82488.500

We assert that the planar approximation is worst at the midpoint of the geodesic, i.e. (82488.5 m) / 2  =  41244.25 m, from Helsinki. So now, we want to find the coordinates of this point, by solving the "direct geodetic problem":

echo %HEL%  -172.370214337896    41244.25 | geod -f "%.12f" +ellps=GRS80

resulting in:

59.804045696236    24.840439590098    7.545306054713

i.e. the latitude, the longitude, and the return azimuth of the mid point of the geodesic.

Now, for convenience, define this environment variable:

set MID="59.804045696236  24.840439590098"

Note

You can save a bit of manual copying by capturing the geod output into the clipboard, using the "clip" command on Windows, i.e. by appending | clip

to the geod command above.

1. The planar approximation

We will be working on ETRS89/UTM35 coordinates, so first define an environment variable, to save some typing:

set utm35=proj -rs +proj=utm +zone=35 +ellps=GRS80

Note

The -rs option switches the proj coordinate i/o order to latitude/longitude and northing/easting, respectively. We do this to comply with the expected order for the geod program.

Now, find the UTM coordinates of the two end points:

echo %HEL% | %utm35%
echo %TAL% | %utm35%

Note your results and compute the mean of the northings, and the mean of the eastings, to obtain the planar approximation of the mid point

Hint

you can use python as a makeshift command line calculator by saying:

python -c print((4+8)/2)

Now, compute the geographical coordinates corresponding to the UTM mid point:

echo <your northing   your easting> | %utm35% -f "%.12f" -I

(note the -I option for doing the inverse projection)

For convenience define:

set MID_APPROX=<your result here>

Finally, compute the distance between the geodesic mid point and its planar approximation, by stating it as another case of the inverse geodetic problem:

echo %MID%  %MID_APPROX% | geod -I +ellps=GRS80

resulting in:

107d0'16.673"    -72d59'43.193"    2.527

i.e. a deviation of 2.5 m over a stretch of about 85 km.

4. Suggested meditations

Consider these aspects:

  1. Given that geod is available, fast, and reliable - is it really worth the effort doing approximate calculations in the projected plane?

  2. geod includes functionality for computing intermediate points along a geodesic. Check the manual, especially the description of the n_S option, and try to compute the geodesic mid point directly, by setting n_S=2.

  3. One should actually expect the result of meditation 2 above to be slightly superior to the result obtained in the exercise. Why?

Answers

Helsinki,  UTM:  6672241.54   385592.95
Tallinn,   UTM:  6590881.40   372106.37
Mid point, UTM:  6631561.47   378849.66
Mid point, GEO:  59.804039062602   24.840482644156

Meditations:

  1. Probably not

  2. geod +lat_1=60.171 +lon_1=24.938 +lat_2=59.437 +lon_2=24.745 +n_S=2 +ellps=GRS80 -f "%.12f"

  3. Avoids a bit of truncation of the azimuth, and may use a slightly superior algorithm (check the code!)