Ellipsoidal geodesy

The package contains a port for Dart language of ellipsoidal geodesy tools (more spefically Vincenty solutions of geodesics on the ellipsoid), originally written in JavaScript by Chris Veness. See the online form on the Movable Type Scripts web site and source code on GitHub.

These functions are based on an ellipsoidal earth model using direct and inverse solutions of geodesics on the ellipsoid devised by Thaddeus Vincenty.

Functions available:

inverse Vincenty calculations
- input: origin and destination geographic positions
- output: initial and final bearings in degrees, and distance in meters between input positions along the shortest geodetic arc (geodesic) on the ellipsoid surface
direct Vincenty calculations
- input: an origin geographic position, an initial bearing and a distance in meters to a destination point
- output: a destination geographic position and a final bearing
shortcuts to calculate just distance or bearings using inverse functions
shortcuts to calculate midpoints and intermediate points (with bearings that varies along the geodesic) between origin and destination geographic positions

The accuracy of distances calculated is approximately within 0.5 mm and bearings within 0.000015″. However even if this sound quite accurate, normally source positions (measured by GPS handheld devices etc.) is not within that accuracy but rather within meter scale tolerances. So it’s important to evaluate also the accuracy of source positions, not only the quality of calculations.

Anyway using these ellipsoidal functions should yield better results than using similar methods on spherical geodesy functions, that are also available on this library.

The package also provides methods to transform between ellipsoidal position representations

geocentric cartesian coordinates represented by ECEF (earth-centric earth-fixed) positions
geographic positions (latitude and longitude as geodetic coordinates)

🌐 Vincenty ellipsoidal calculations

Distances & bearings between points, and destination points calculated on an ellipsoidal earth model, along geodesics on the surface of a reference ellipsoid selected:

  // sample geographic positions
  final greenwich = Geographic.parseDms(lat: '51°28′40″ N', lon: '0°00′05″ W');
  final sydney = Geographic.parseDms(lat: '33.8688° S', lon: '151.2093° E');

  // decimal degrees (DD) and degrees-minutes (DM) formats
  const dd = Dms(decimals: 2);
  const dm = Dms.narrowSpace(type: DmsType.degMin, decimals: 2);

  // the shortest arc along the geodesic on the ellipsoid surface between points
  final arc1 = greenwich.vincenty().inverse(sydney);

  // prints (distance of the geodesic): 16983.3 km
  final distanceKm = arc1.distance / 1000.0;
  print('${distanceKm.toStringAsFixed(1)} km');

  // prints (bearing varies along the geodesic): 60.59° -> 139.15°
  final initialBearing = arc1.bearing;
  final finalBearing = arc1.finalBearing;
  print('${dd.bearing(initialBearing)} -> ${dd.bearing(finalBearing)}');

  // the shortest arc along the geodesic from greenwich to the initial direction
  // defined by `bearing` and the length by `distance`
  // prints: 51° 31.28′ N, 0° 07.48′ E - bearing: 61.10°
  final arc2 = greenwich.vincenty().direct(distance: 10000, bearing: 61.0);
  final dest = arc2.destination;
  final destBrng = arc2.finalBearing;
  print('${dest.latLonDms(format: dm)} - bearing: ${dd.bearing(destBrng)}');

There are also methods to calculate mid points and intermediate points along the geodesic.

  // mid point, prints: 28° 52.77′ N, 104° 48.82′ E
  final midPoint = greenwich.vincenty().midPointTo(sydney);
  print(midPoint.latLonDms(format: dm));

  // intermediate points along the geodesic between Greenwich and Sydney
  // prints 10 points with bearings on intermediate geographic positions:
  // 0.0: 51° 28.67′ N, 0° 00.08′ W - bearing: 60.59°
  // 0.1: 56° 39.07′ N, 24° 34.88′ E - bearing: 80.62°
  // 0.2: 56° 03.20′ N, 52° 13.09′ E - bearing: 103.76°
  // 0.3: 49° 56.60′ N, 75° 34.27′ E - bearing: 122.53°
  // 0.4: 40° 19.81′ N, 92° 27.98′ E - bearing: 134.59°
  // 0.5: 28° 52.77′ N, 104° 48.82′ E - bearing: 141.66°
  // 0.6: 16° 30.55′ N, 114° 36.83′ E - bearing: 145.47°
  // 0.7: 3° 43.29′ N, 123° 12.60′ E - bearing: 146.99°
  // 0.8: 9° 09.07′ S, 131° 33.47′ E - bearing: 146.59°
  // 0.9: 21° 48.54′ S, 140° 31.88′ E - bearing: 144.18°
  // 1.0: 33° 52.13′ S, 151° 12.56′ E - bearing: 139.15°
  for (var fr = 0.0; fr < 1.0; fr += 0.1) {
    final ip = greenwich.vincenty().intermediatePointTo(sydney, fraction: fr);
    final point = ip.origin;
    final pointBrng = ip.bearing;
    print('${fr.toStringAsFixed(1)}: ${point.latLonDms(format: dm)}'
        ' - bearing: ${dd.bearing(pointBrng)}');
  }

The package provides Ellipsoid.WGS84 and Ellipsoid.GRS80 ellipsoids as constant value objects. These ellipsoids should be enough for most modern geodetic calculations:

  // to use alternative ellipsoids set an optional argument on `vincenty` method
  final distanceGRS80 =
      greenwich.vincenty(ellipsoid: Ellipsoid.GRS80).distanceTo(sydney);

However if an older ellipsoid is needed it’s possible to define such manually using known ellipsoidal parameters. The following sample defines the Airy 1830 ellipoid and compares distance calculated along geodesics of each ellipsoid surfaces:

  // custom ellipsoids can be used also
  final airy = Ellipsoid(
    id: 'airy',
    name: 'Airy 1830',
    a: 6377563.396,
    b: 6356256.909,
    f: 1.0 / 299.3249646,
  );
  final distanceAiry = greenwich.vincenty(ellipsoid: airy).distanceTo(sydney);

  // Distances printed: 16983280.66025 m, 16983280.66013 m, 16981837.55212 m
  // (Note only very small difference between WGS84 and GRS80 ellipsoids,
  //  however this level of "accuracy" is out of nominal accuracy of measured
  //  points and Vincenty calculations with 0.5 mm expected accuracy)
  print('Distance WGS84: ${(distanceKm * 1000.0).toStringAsFixed(5)} m');
  print('Distance GRS80: ${distanceGRS80.toStringAsFixed(5)} m');
  print('Distance Airy1830: ${distanceAiry.toStringAsFixed(5)} m');

🗺️ Transform between geographic and geocentric coordinates

This sample shows how to transform coordinate values between geographic positions (latitude and longitude as geodetic coordinates) and geocentric cartesian coordinates represented by ECEF (earth-centric earth-fixed) positions:

  // a sample geographic position with geodetic latitude and longitude
  const geographic1 = Geographic(lat: 51.4778, lon: -0.0014, elev: 45.0);

  // same as ECEF (earth-centric earth-fixed) geocentric cartesian coordinates
  final geocentric1 = geographic1.toGeocentricCartesian();

  // returned object is of type `Position` with x, y and z cartesian coordinates
  // prints (X, Y, Z): 3980609.2373, -97.2646, 4966859.7285
  print(geocentric1.toText(decimals: 4, delimiter: ', '));

  // let's try inverse, first create a geocentric cartesian position (ECEF)
  final geocentric2 =
      Position.create(x: 3980609.2373, y: -97.2646, z: 4966859.7285);

  // convert this to a geographic position with geodetic latitude and longitude
  final geographic2 = EllipsoidalExtension.fromGeocentricCartesian(geocentric1);

  // returned object is of the type `Geographic`
  // prints (longitude, latitude, elevation): -0.0014, 51.4778, 45.0000
  print(geographic2.toText(decimals: 4, delimiter: ', ', compactNums: false));

Samples above used the WGS84 reference ellipsoid for ellipsoidal calculations. You can also use other ellipsoids on toGeocentricCartesian and fromGeocentricCartesian methods:

  // create a geocentric cartesian coordinates based on the GRS80 ellipsoid from
  // the same geodetic coordinate values as specified on `geographic1`
  final geocentricGRS80 =
      geographic1.toGeocentricCartesian(ellipsoid: Ellipsoid.GRS80);

  // prints (X, Y, Z): 3980609.2373, -97.2646, 4966859.7284
  print(geocentricGRS80.toText(decimals: 4, delimiter: ', '));

  // As WGS84 and GRS80 ellipsoids are very close to each other you may notice
  // only a small difference compared to values printed from `geocentric1`.