Package org.apache.sis.referencing

# Class GeodeticCalculator

Object
GeodeticCalculator

public class GeodeticCalculator extends Object
Performs geodetic calculations on a sphere or an ellipsoid. This class computes the distance between two points, or conversely the point located at a given distance from another point when navigating in a given direction. The distance depends on the path (or track) on Earth surface connecting the two points. The track can be great circles (shortest path between two points) or rhumb lines (path with constant heading).

This class uses the following information:

- The start point, which is always considered valid after the first call to `setStartPoint(…)`. Its value can only be changed by another call to `setStartPoint(…)`.
- One of the followings (the latest specified properties override other properties and determines what will be calculated):
  - the end point, or
  - the azimuth at start point together with the geodesic distance from that point.

## Algorithms

`GeodeticCalculator` uses two set of formulas, depending if the figure of the Earth is a sphere or an ellipsoid. Publications relevant to this class are:

- Wikipedia, Great-circle navigation for spherical formulas.
- Wikipedia, Rhumb line for spherical formulas.
- Charles F. F. Karney (2013), Algorithms for geodesics for ellipsoidal formulas.
- G.G. Bennett, 1996. Practical Rhumb Line Calculations on the Spheroid for ellipsoidal formulas.
- Charles F. F. Karney (2010), Test set for geodesics for `GeodeticCalculator` tests.
- Charles F. F. Karney, GeographicLib for the reference implementation.

## Accuracy

`GeodeticCalculator` aims for a positional accuracy of one centimetre. The accuracy is often better (about one millimetre), but not everywhere. Azimuthal accuracy corresponds to an error of one centimetre at a distance of one kilometer, except for nearly antipodal points (less than 1° of longitude and latitude from antipode) and points close to the poles where the azimuthal errors are larger. Karney's GeographicLib should be used if better accuracy is desired. Apache SIS accuracy does not go as far as GeographicLib because the rest of Apache SIS library (map projections, etc.) aims for an one centimetre accuracy anyway.

## Limitations

Current implementation can not compute the geodesics in some cases. In particular, calculation may fail for antipodal points on an ellipsoid. Karney's algorithm should cover those cases, but this `GeodeticCalculator` implementation may not be sufficiently tuned. See SIS-467 for more information.

## Thread safety

This class is not thread-safe. If geodetic calculations are needed in a multi-threads environment, then a distinct instance of `GeodeticCalculator` needs to be created for each thread.

Since:
  1.0

Defined in the `sis-referencing` module

- ## Method Summary

  Modifier and Type
  Method
  Description
  `static GeodeticCalculator`
  `create(CoordinateReferenceSystem crs)`
  Constructs a new geodetic calculator expecting coordinates in the supplied CRS.
  `Shape`
  `createGeodesicCircle2D(double tolerance)`
  Creates an approximation of the curve at a constant geodesic distance around the start point.
  `Shape`
  `createGeodesicPath2D(double tolerance)`
  Creates an approximation of the geodesic track from start point to end point as a Java2D object.
  `MathTransform`
  `createProjectionAroundStart()`
  Creates an Azimuthal Equidistant projection centered on current starting point.
  `double`
  `getConstantAzimuth()`
  Computes the angular heading of a rhumb line path.
  `Unit<Length>`
  `getDistanceUnit()`
  Returns the unit of measurement of all distance measurements.
  `double`
  `getEndingAzimuth()`
  Computes the angular heading at the ending point of a geodesic path.
  `DirectPosition`
  `getEndPoint()`
  Returns or computes the destination in the CRS specified at construction time.
  `double`
  `getGeodesicDistance()`
  Returns or computes the shortest distance from start point to end point.
  `GeographicCRS`
  `getGeographicCRS()`
  Returns the coordinate reference system for all methods expecting (φ,λ) as `double` values.
  `CoordinateReferenceSystem`
  `getPositionCRS()`
  Returns the Coordinate Reference System (CRS) in which `Position`s are represented, unless otherwise specified.
  `double`
  `getRhumblineLength()`
  Returns or computes the length of rhumb line (part of constant heading) from start point to end point.
  `double`
  `getStartingAzimuth()`
  Returns or computes the angular heading at the starting point of a geodesic path.
  `DirectPosition`
  `getStartPoint()`
  Returns the starting point in the CRS specified at construction time.
  `void`
  `moveToEndPoint()`
  Sets the start point and starting azimuth to the current end point and ending azimuth values.
  `void`
  `setEndGeographicPoint(double latitude, double longitude)`
  Sets the destination as geographic (latitude, longitude) coordinates.
  `void`
  `setEndPoint(Position position)`
  Sets the destination as coordinates in arbitrary reference system.
  `void`
  `setGeodesicDistance(double distance)`
  Sets the geodesic distance from the start point to the end point.
  `void`
  `setStartGeographicPoint(double latitude, double longitude)`
  Sets the starting point as geographic (latitude, longitude) coordinates.
  `void`
  `setStartingAzimuth(double azimuth)`
  Sets the angular heading at the starting point of a geodesic path.
  `void`
  `setStartPoint(Position point)`
  Sets the starting point as coordinates in arbitrary reference system.
  `String`
  `toString()`
  Returns a string representation of start point, end point, azimuths and distance.

  ### Methods inherited from class Object

  `clone, equals, finalize, getClass, hashCode, notify, notifyAll, wait, wait, wait`

- ## Method Details

  - ### create

    public static GeodeticCalculator create(CoordinateReferenceSystem crs)
    Constructs a new geodetic calculator expecting coordinates in the supplied CRS. All `GeodeticCalculator` methods having a `Position` argument or return value will use that specified CRS. That CRS is the value returned by `getPositionCRS()`.

    #### Limitations

    Current implementation uses only spherical formulas. Implementation using ellipsoidal formulas will be provided in a future Apache SIS release.
    Parameters:
      `crs` - the reference system for the `Position` objects.
    Returns:
      a new geodetic calculator using the specified CRS.

  - ### getPositionCRS

    public CoordinateReferenceSystem getPositionCRS()
    Returns the Coordinate Reference System (CRS) in which `Position`s are represented, unless otherwise specified. This is the CRS of all `Position` instances returned by methods in this class. This is also the default CRS assumed by methods receiving a `Position` argument when the given position does not specify its CRS. This default CRS is specified at construction time. It is not necessarily geographic; it may be projected or geocentric.
    Returns:
      the default CRS for `Position` instances.

  - ### getGeographicCRS

    public GeographicCRS getGeographicCRS()
    Returns the coordinate reference system for all methods expecting (φ,λ) as `double` values. This CRS always has (latitude, longitude) axes, in that order and in degrees. The CRS may contain an additional axis for ellipsoidal height.
    Returns:
      the coordinate reference system of (φ,λ) coordinates.

  - ### getStartPoint

    public DirectPosition getStartPoint()
    Returns the starting point in the CRS specified at construction time. This method returns the last point given to a `setStartPoint(…)` method, transformed to the position CRS.
    Returns:
      the starting point represented in the CRS specified at construction time.
    Throws:
      `IllegalStateException` - if the start point has not yet been specified.
      `GeodeticException` - if the coordinates can not be transformed to position CRS.
    See Also:
      - `getEndPoint()`

  - ### setStartPoint

    public void setStartPoint(Position point)
    Sets the starting point as coordinates in arbitrary reference system. This method transforms the given coordinates to geographic coordinates, then delegates to `setStartGeographicPoint(double, double)`. If the given point is not associated to a Coordinate Reference System (CRS), then this method assumes the CRS specified at construction time.
    Parameters:
      `point` - the starting point in any coordinate reference system.
    Throws:
      `IllegalArgumentException` - if the given coordinates can not be transformed.
    See Also:
      - `setEndPoint(Position)`

  - ### setStartGeographicPoint

    public void setStartGeographicPoint(double latitude, double longitude)
    Sets the starting point as geographic (latitude, longitude) coordinates. The starting and ending azimuths, the end point, the geodesic distance and the rhumb line length are discarded by this method call; some of them will need to be specified again.
    Parameters:
      `latitude` - the latitude in degrees between -90.0° and 90.0°.
      `longitude` - the longitude in degrees.
    See Also:
      - `setEndGeographicPoint(double, double)`
      - `moveToEndPoint()`

  - ### getEndPoint

    public DirectPosition getEndPoint()
    Returns or computes the destination in the CRS specified at construction time. This method returns the point specified in the last call to a `setEndPoint(…)` method, unless the starting azimuth and geodesic distance have been set more recently. In the later case, the end point will be computed from the start point and the current azimuth and distance.
    Returns:
      the destination (end point) represented in the CRS specified at construction time.
    Throws:
      `IllegalStateException` - if the destination point, azimuth or distance have not been set.
      `GeodeticException` - if the coordinates can not be computed.
    See Also:
      - `getStartPoint()`

  - ### setEndPoint

    public void setEndPoint(Position position)
    Sets the destination as coordinates in arbitrary reference system. This method transforms the given coordinates to geographic coordinates, then delegates to `setEndGeographicPoint(double, double)`. If the given point is not associated to a Coordinate Reference System (CRS), then this method assumes the CRS specified at construction time.
    Parameters:
      `position` - the destination (end point) in any coordinate reference system.
    Throws:
      `IllegalArgumentException` - if the given coordinates can not be transformed.
    See Also:
      - `setStartPoint(Position)`

  - ### setEndGeographicPoint

    public void setEndGeographicPoint(double latitude, double longitude)
    Sets the destination as geographic (latitude, longitude) coordinates. The starting azimuth, ending azimuth geodesic distance and rhumb line length will be updated as an effect of this call.
    Parameters:
      `latitude` - the latitude in degrees between -90.0° and 90.0°.
      `longitude` - the longitude in degrees.
    See Also:
      - `setStartGeographicPoint(double, double)`

  - ### getStartingAzimuth

    public double getStartingAzimuth()
    Returns or computes the angular heading at the starting point of a geodesic path. Azimuth is relative to geographic North with values increasing clockwise. This method returns the azimuth normalized to [-180 … +180]° range given in last call to `setStartingAzimuth(double)` method, unless the `setEndPoint(…)` method has been invoked more recently. In the later case, the azimuth will be computed from the start point and the current end point.
    Returns:
      the azimuth in degrees from -180° to +180°. 0° is toward North and values are increasing clockwise.
    Throws:
      `IllegalStateException` - if the end point, azimuth or distance have not been set.
      `GeodeticException` - if the azimuth can not be computed.

  - ### setStartingAzimuth

    public void setStartingAzimuth(double azimuth)
    Sets the angular heading at the starting point of a geodesic path. Azimuth is relative to geographic North with values increasing clockwise. The ending azimuth, end point and rhumb line length will be updated as an effect of this method call.
    Parameters:
      `azimuth` - the starting azimuth in degrees, with 0° toward north and values increasing clockwise.
    See Also:
      - `setGeodesicDistance(double)`

  - ### getEndingAzimuth

    public double getEndingAzimuth()
    Computes the angular heading at the ending point of a geodesic path. Azimuth is relative to geographic North with values increasing clockwise. This method computes the azimuth from the current start point and end point, or from start point and the current starting azimuth and geodesic distance.
    Returns:
      the azimuth in degrees from -180° to +180°. 0° is toward North and values are increasing clockwise.
    Throws:
      `IllegalStateException` - if the destination point, azimuth or distance have not been set.
      `GeodeticException` - if the azimuth can not be computed.

  - ### getConstantAzimuth

    public double getConstantAzimuth()
    Computes the angular heading of a rhumb line path. Azimuth is relative to geographic North with values increasing clockwise.
    Returns:
      the azimuth in degrees from -180° to +180°. 0° is toward North and values are increasing clockwise.
    Throws:
      `IllegalStateException` - if the start point or end point has not been set.
      `GeodeticException` - if the azimuth can not be computed.

  - ### getGeodesicDistance

    public double getGeodesicDistance()
    Returns or computes the shortest distance from start point to end point. This is sometime called "great circle" or "orthodromic" distance. This method returns the value given in last call to `setGeodesicDistance(double)`, unless the `setEndPoint(…)` method has been invoked more recently. In the later case, the distance will be computed from the start point and current end point.
    Returns:
      the shortest distance in the unit of measurement given by `getDistanceUnit()`.
    Throws:
      `IllegalStateException` - if the start point or end point has not been set.
      `GeodeticException` - if the distance can not be computed.
    See Also:
      - `getDistanceUnit()`

  - ### setGeodesicDistance

    public void setGeodesicDistance(double distance)
    Sets the geodesic distance from the start point to the end point. The end point, ending azimuth and rhumb line length will be updated as an effect of this method call.
    Parameters:
      `distance` - the geodesic distance in unit of measurement given by `getDistanceUnit()`.
    See Also:
      - `setStartingAzimuth(double)`

  - ### getRhumblineLength

    public double getRhumblineLength()
    Returns or computes the length of rhumb line (part of constant heading) from start point to end point. This is sometime called "loxodrome". This is not the shortest path between two points. The rhumb line distance may be up to 50% longer than the geodesic distance.
    Returns:
      length of rhumb line in the unit of measurement given by `getDistanceUnit()`.
    Throws:
      `IllegalStateException` - if a point has not been set.

  - ### getDistanceUnit

    public Unit<Length> getDistanceUnit()
    Returns the unit of measurement of all distance measurements. This is the ellipsoid axis unit.
    Returns:
      the unit of measurement of all distance measurements.
    See Also:
      - `getGeodesicDistance()`

  - ### moveToEndPoint

    public void moveToEndPoint()
    Sets the start point and starting azimuth to the current end point and ending azimuth values. The ending azimuths, the geodesic distance and the end point are discarded by this method call; some of them will need to be specified again.
    Throws:
      `GeodeticException` - if the end point or ending azimuth can not be computed.
    See Also:
      - `setStartGeographicPoint(double, double)`

  - ### createGeodesicPath2D

    public Shape createGeodesicPath2D(double tolerance)
    Creates an approximation of the geodesic track from start point to end point as a Java2D object. The coordinates are expressed in the coordinate reference system specified at creation time. The approximation uses linear, quadratic or cubic Bézier curves. The returned path has the following characteristics:
    1. The first point is `getStartPoint()`.
    2. The beginning of the curve (more specifically, the tangent at starting point) is oriented toward the direction given by getStartingAzimuth(), adjusted for the map projection (if any) deformation at that location.
    3. The point B(½) in the middle of the Bézier curve is a point of the geodesic path.
    4. The end of the curve (more specifically, the tangent at ending point) is oriented toward the direction given by getEndingAzimuth(), adjusted for the map projection (if any) deformation at that location.
    5. The last point is `getEndPoint()`, potentially with 360° added or subtracted to the longitude.
    This method tries to stay within the given tolerance threshold of the geodesic track. The `tolerance` parameter should not be too small for avoiding creation of unreasonably long chain of Bézier curves. For example a value of 1/10 of geodesic length may be sufficient.
    Note: this method depends on the presence of `java​.desktop` module. This constraint may be addressed in a future Apache SIS version (see SIS-453). The "2D" suffix in the method name represents this relationship with Java2D. The `createGeodesicPath(…)` method name (without suffix) is reserved for a future version using ISO curves instead.
    Parameters:
      `tolerance` - maximal error between the approximated curve and actual geodesic track in the units of measurement given by `getDistanceUnit()`. This is approximate; the actual errors may vary around that value.
    Returns:
      an approximation of geodesic track as Bézier curves in a Java2D object.
    Throws:
      `IllegalStateException` - if some required properties have not been specified.
      `GeodeticException` - if some coordinates can not be computed.

  - ### createGeodesicCircle2D

    public Shape createGeodesicCircle2D(double tolerance)
    Creates an approximation of the curve at a constant geodesic distance around the start point. The returned shape is circlelike with the start point in its center. The coordinates are expressed in the coordinate reference system specified at creation time. The approximation uses cubic Bézier curves.
    Note: some authors define geodesic circle as the curve which enclose the maximum area for a given perimeter. This method adopts a different definition, the locus of points at a fixed geodesic distance from center point.
    This method tries to stay within the given tolerance threshold of the geodesic track. The `tolerance` parameter should not be too small for avoiding creation of unreasonably long chain of Bézier curves. For example a value of 1/10 of geodesic length may be sufficient.
    Note: this method depends on the presence of `java​.desktop` module. This constraint may be addressed in a future Apache SIS version (see SIS-453). The "2D" suffix in the method name represents this relationship with Java2D. The `createGeodesicCircle(…)` method name (without suffix) is reserved for a future version using ISO curves instead.
    Parameters:
      `tolerance` - maximal error in the units of measurement given by `getDistanceUnit()`. This is approximate; the actual errors may vary around that value.
    Returns:
      an approximation of circular region as a Java2D object.
    Throws:
      `IllegalStateException` - if some required properties have not been specified.
      `GeodeticException` - if some coordinates can not be computed.

  - ### createProjectionAroundStart

    public MathTransform createProjectionAroundStart()
    Creates an Azimuthal Equidistant projection centered on current starting point. On input, the `MathTransform` expects coordinates expressed in the position CRS. On output, the `MathTransform` produces coordinates in a `ProjectedCRS` having the following characteristics:
    - Coordinate system is a two-dimensional `CartesianCS` with (Easting, Northing) axis order and directions.
    - Unit of measurement is the same as position CRS if those units are linear, or `Units​.METRE` otherwise.
    - Projection of the start point results in (0,0).
    - Distances relative to (0,0) are approximately exact for distances less than 800 km.
    - Azimuths from (0,0) to other points are approximately exact for points located at less than 800 km.
    Given above characteristics, the following calculations are satisfying approximations when using (x, y) coordinates in the output space for D < 800 km:

    > D = √(x² + y²) — distance from projection center.
    > θ = atan2(y, x) — arithmetic angle from projection center to (x, y).
    > x = D⋅cos θ
    > y = D⋅sin θ — end point for a distance and angle from start point.

    The following calculations are not exacts, because distances and azimuths are approximately exacts only when measured from (0,0) coordinates:

    > D = √[(x₂ − x₁)² + (y₂ − y₁)²] — distances between points other then projection center are not valid.
    > θ = atan2(y₂ − y₁, x₂ − x₁) — azimuths between points other then projection center are not valid.
    > etc.

    This method can be invoked repetitively for doing calculations around different points. All returned `MathTransform` instances are immutable; changing `GeodeticCalculator` state does not affect those transforms.
    Returns:
      transform from position CRS to Azimuthal Equidistant projected CRS centered on current start point.
    Throws:
      `IllegalStateException` - if the start point has not been set.
      `GeodeticException` - if the projection can not be computed.
    Since:
      1.1
    See Also:
      - `moveToEndPoint()`
      - `AzimuthalEquidistant`

  - ### toString

    public String toString()
    Returns a string representation of start point, end point, azimuths and distance. The text representation is implementation-specific and may change in any future version. Current implementation is like below:

    >     Coordinate reference system: Unspecified datum based upon the GRS 1980 Authalic Sphere
    >     ┌─────────────┬─────────────────┬──────────────────┬─────────────┐
    >     │             │    Latitude     │    Longitude     │   Azimuth   │
    >     │ Start point │  9°39′06.1120″N │ 132°37′37.1248″W │  -17°10′37″ │
    >     │ End point   │ 70°32′45.0206″N │ 109°50′05.0533″E │ -119°03′12″ │
    >     └─────────────┴─────────────────┴──────────────────┴─────────────┘
    >     Geodesic distance: 9,967,530.74 m

    Overrides:
      `toString` in class `Object`
    Returns:
      a string representation of this calculator state.