Harlinn.GeographicLib.Net - for solving geodesic problems.

The library can be used to solve geodesic and rhumb line calculations, convert between geographic, UTM, UPS, MGRS, geocentric, and local cartesian coordinates, and solve gravity and geomagnetic field calculations.

Implemented in C++/CLI to provide access to the functionality of the GeographicLib library.

Harlinn.GeographicLib.Net works with the release build of Harlinn.GeographicLib.

The information below was copied from the original GeographicLib documentation and altered slightly, as Harlinn.GeographicLib.Net will mostly be used from apps written in C#.

Geodesic calculations

The shortest path between two points on an ellipsoid at $(l a t_{1}, l o n_{1})$ and $(l a t_{2}, l o n_{2})$ is called the geodesic. Its length is $s_{12}$ and the geodesic from point 1 to point 2 has azimuths $a z i_{1}$ and $a z i_{2}$ at the two end points. (The azimuth is the heading measured clockwise from north. $a z i_{2}$ is the “forward” azimuth, i.e., the heading that takes you beyond point 2 not back to point 1.) In the figure below, latitude is labeled $φ$ , longitude $λ$ (with $λ_{12} = λ_{2} - λ_{1}$ ), and azimuth $α$ .

Given $l a t_{1}$ , $l o n_{1}$ , $a z i_{1}$ , and $s_{12}$ , we can determine $l a t_{2}$ , $l o n_{2}$ , and $a z i_{2}$ . This is the direct geodesic problem and its solution is given by the function Geodesic.Direct. (If $s_{12}$ is sufficiently large that the geodesic wraps more than halfway around the earth, there will be another geodesic between the points with a smaller $s_{12}$ .)

Given $l a t_{1}$ , $l o n_{1}$ , $l a t_{2}$ , and $l o n_{2}$ , we can determine $a z i_{1}$ , $a z i_{2}$ , and $s_{12}$ . This is the inverse geodesic problem, whose solution is given by Geodesic.Inverse. Usually, the solution to the inverse problem is unique. In cases where there are multiple solutions (all with the same $s_{12}$ , of course), all the solutions can be easily generated once a particular solution is provided.

The standard way of specifying the direct problem is the specify the distance $s_{12}$ to the second point. However it is sometimes useful instead to specify the arc length $a_{12}$ (in degrees) on the auxiliary sphere. This is a mathematical construct used in solving the geodesic problems. The solution of the direct problem in this form is provided by Geodesic.ArcDirect. An arc length in excess of $180 °$ indicates that the geodesic is not a shortest path. In addition, the arc length between an equatorial crossing and the next extremum of latitude for a geodesic is $90 °$ .

This class can also calculate several other quantities related to geodesics. These are:

reduced length. If we fix the first point and increase $a z i_{1}$ by $d a z i_{1}$ (radians), the second point is displaced $m_{12} d a z i_{1}$ in the direction $a z i_{2} + 90 °$ . The quantity $m_{12}$ is called the “reduced length” and is symmetric under interchange of the two points. On a curved surface the reduced length obeys a symmetry relation, $m_{12} + m_{21} = 0$ . On a flat surface, we have $m_{12} = s_{12}$ . The ratio $s_{12} / m_{12}$ gives the azimuthal scale for an azimuthal equidistant projection.
geodesic scale. Consider a reference geodesic and a second geodesic parallel to this one at point 1 and separated by a small distance $d t$ . The separation of the two geodesics at point 2 is $M_{12} d t$ where $M_{12}$ is called the “geodesic scale”. $M_{21}$ is defined similarly (with the geodesics being parallel at point 2). On a flat surface, we have $M_{12} = M_{21} = 1$ . The quantity $1 / M_{12}$ gives the scale of the Cassini-Soldner projection.
area. The area between the geodesic from point 1 to point 2 and the equation is represented by $S_{12}$ ; it is the area, measured counter-clockwise, of the geodesic quadrilateral with corners $(l a t_{1}, l o n_{1})$ , $(0, l o n_{1})$ , $(0, l o n_{2})$ , and $(l a t_{2}, l o n_{2})$ . It can be used to compute the area of any geodesic polygon.

Overloaded versions of Geodesic.Direct, Geodesic.ArcDirect, and Geodesic.Inverse allow these quantities to be returned. In addition there are general functions Geodesic.GenDirect, and Geodesic.GenInverse which allow an arbitrary set of results to be computed. The quantities $m_{12}$ , $M_{12}$ , $M_{21}$ which all specify the behavior of nearby geodesics obey addition rules. If points 1, 2, and 3 all lie on a single geodesic, then the following rules hold:

\begin{array}{r} s_{13} = s_{12} + s_{23} \\ a_{13} = a_{12} + a_{23} \\ S_{13} = S_{12} + S_{23} \\ m_{13} = m_{12} M_{23} + m_{23} M_{21} \\ M_{13} = M_{12} M_{23} - (1 - M_{12} M_{21}) m_{23} / m_{12} \\ M_{31} = M_{32} M_{21} - (1 - M_{23} M_{32}) m_{12} / m_{23} \end{array}

Additional functionality is provided by the GeodesicLine class, which allows a sequence of points along a geodesic to be computed.

Geodesic example

using Harlinn.GeographicLib.Net;

namespace GeodesicExample
{
    internal class Program
    {
        static void Main(string[] args)
        {
            Geodesic geod = Geodesic.WGS84;

            // Sample direct calculation, travelling about NE from JFK
            double lat1 = 40.6, lon1 = -73.8, s12 = 5.5e6, azi1 = 51;
            double lat2, lon2;
            geod.Direct(lat1, lon1, azi1, out s12, out lat2, out lon2);

            Console.WriteLine($"Latitude: {lat2}, Longitude:{lon2}");

            // Sample inverse calculation, JFK to LHR

            // JFK Airport
            lat1 = 40.6;
            lon1 = -73.8;

            // LHR Airport
            lat2 = 51.6; 
            lon2 = -0.5;  
            
            geod.Inverse(lat1, lon1, lat2, lon2, out s12);

            Console.WriteLine($"Distance: {s12}");
        }
    }
}

Latitude: 51,88456450560677, Longitude:-1,1411728612008574
Distance: 5551759,400318679

GeodesicLine

GeodesicLine facilitates the determination of a series of points on a single geodesic. The starting point $(l a t_{1}, l o n_{1})$ and the azimuth $a z i_{1}$ are specified in the constructor; alternatively, the Geodesic::Line method can be used to create a GeodesicLine. GeodesicLine.Position returns the location of point 2 a distance $s_{12}$ along the geodesic. In addition, GeodesicLine.ArcPosition gives the position of point 2 an arc length $a_{12}$ along the geodesic.

You can register the position of a reference point 3 a distance (arc length), $s_{13}$ ( $a_{13}$ ) along the geodesic with the GeodesicLine.SetDistance (GeodesicLine.SetArc) functions. Points a fractional distance along the line can be found by providing, for example, 0.5 * Distance as an argument to GeodesicLine.Position. The Geodesic.InverseLine or Geodesic.DirectLine methods return GeodesicLine objects with point 3 set to the point 2 of the corresponding geodesic problem. GeodesicLine objects created with the public constructor or with Geodesic.Line have $s_{13}$ and $a_{13}$ set to NaNs.

The calculations are accurate to better than 15 nm (15 nanometers). See Sec. 9 of arXiv:1102.1215v1 for details. With exact = false (the default) in the constructor for the Geodesic object, the algorithms used by this class are based on series expansions using the flattening $f$ as a small parameter. These are only accurate for $| f | < 0.02$ ; however reasonably accurate results will be obtained for $| f | < 0.2$ . For very eccentric ellipsoids, set exact = true in the constructor for the Geodesic object; this will delegate the calculations to GeodesicLineExact.

GeodesicLine example

using Harlinn.GeographicLib.Net;

namespace GeodesicLineExample
{
    internal class Program
    {
        static void Main(string[] args)
        {
            Geodesic geod = Geodesic.WGS84;

            double
              // JFK airport
              lat1 = 40.640, lon1 = -73.779,
              // SIN airport
              lat2 = 1.359, lon2 = 103.989; 

            GeodesicLine line = geod.InverseLine(lat1, lon1, lat2, lon2);

            // Nominal distance between points = 500 km
            double ds0 = 500e3;

            // The number of intervals
            int num = (int)Math.Ceiling(line.Distance / ds0);

            double ds = line.Distance / num;
            for (int i = 0; i <= num; ++i)
            {
                double lat, lon;
                line.Position(i * ds, out lat, out lon);
                Console.WriteLine($"{i} -> Latitude: {lat}, Longitude: {lon} ");
            }
        }
    }
}

-> Latitude: 40,64, Longitude: -73,779
-> Latitude: 45,088624522948315, Longitude: -73,41639124870959
-> Latitude: 49,53243435521155, Longitude: -72,99265573206968
-> Latitude: 53,97096253429303, Longitude: -72,48430287040067
-> Latitude: 58,40353976434537, Longitude: -71,85489098898533
-> Latitude: 62,829059340035606, Longitude: -71,04452466778625
-> Latitude: 67,24547195373592, Longitude: -69,94724090357691
-> Latitude: 71,64856286681426, Longitude: -68,35623186483937
-> Latitude: 76,02845249253366, Longitude: -65,80772484891887
-> Latitude: 80,3567290132899, Longitude: -61,014890489628954
-> Latitude: 84,51420549018647, Longitude: -48,83492431665228
-> Latitude: 87,45474332671364, Longitude: 3,894014378845668
-> Latitude: 85,295259509238, Longitude: 71,94358213809562
-> Latitude: 81,22103811744249, Longitude: 87,60160059653242
-> Latitude: 76,91141529030632, Longitude: 93,21228048415551
-> Latitude: 72,53859298991506, Longitude: 96,06333035970746
-> Latitude: 68,13910398636213, Longitude: 97,79733028285946
-> Latitude: 63,72494286554836, Longitude: 98,97301258675924
-> Latitude: 59,301058336202594, Longitude: 99,83082330184723
-> Latitude: 54,86980145924297, Longitude: 100,49102303937627
-> Latitude: 50,432415581007945, Longitude: 101,02040715447234
-> Latitude: 45,98963410041732, Longitude: 101,45907366928466
-> Latitude: 41,54195282757638, Longitude: 101,83259447654174
-> Latitude: 37,089765877640964, Longitude: 102,15809755054522
-> Latitude: 32,633437629751434, Longitude: 102,44754031410375
-> Latitude: 28,17334176771984, Longitude: 102,70958072995484
-> Latitude: 23,709881876563497, Longitude: 102,95070239501462
-> Latitude: 19,243500881260854, Longitude: 103,17592172947425
-> Latitude: 14,774683272647815, Longitude: 103,38925110492858
-> Latitude: 10,303952429663326, Longitude: 103,59401483935383
-> Latitude: 5,831864516291208, Longitude: 103,79307473649479
-> Latitude: 1,3589999999999611, Longitude: 103,989