Search code examples
gisgeotoolsazimuth

geotools GeodeticCalculator

Having two points and a distance, I am trying to compute the azimuth and then to recompute back one of the points.

However, the distance between the computed point and the original point is more then 50 meters, which is quite a big error.

Here is the code:

 public static void main(String[] args) {

 double startLongitude = -5.1085;
 double startLatitude = 40.6682667;
 double endLongitude = -4.000597497067124;
 double endLatitude = 41.49682079962159;
 double distance = 130947.51;

 try {
 CoordinateReferenceSystem crs = CRS.decode("EPSG:4326");
 GeodeticCalculator calculator = new GeodeticCalculator(crs);

 calculator.setStartingGeographicPoint(startLongitude, startLatitude);
 calculator.setDestinationGeographicPoint(endLongitude, endLatitude);
 double azimuth = calculator.getAzimuth();
 System.out.println("Azimuth=" + azimuth);

 calculator = new GeodeticCalculator(crs);
 calculator.setStartingGeographicPoint(startLongitude, startLatitude);
 calculator.setDirection(azimuth, distance);
 Point2D computedEndPoint = calculator.getDestinationGeographicPoint();
 System.out.println("computedEndPoint=" + computedEndPoint);

 calculator = new GeodeticCalculator(crs);
 calculator.setStartingGeographicPoint(endLongitude, endLatitude);
 calculator.setDestinationGeographicPoint(computedEndPoint);
 distance = calculator.getOrthodromicDistance();
 System.out.println("Distance=" + distance);

 } catch (FactoryException e) {
 e.printStackTrace();
 }

}

The output is:

Azimuth=44.97189638988797

computedEndPoint=Point2D.Double[-4.00014170719737, 41.49715519864095]

Distance=53.17698966547863

I expect the computedEndPoint to be quite similar (if not exactly) to declared end point from the beginning. And the distance between these two points to be close to zero.

Now my question is: what am I doing wrong? Or is there some bug in the GeodedicCalculator?

Solution

  • 50m over 130km is 0.04% error - that is pretty good for a round trip of an iterative numerical method.

    You are using the GeodedicCalculator which is using an approximation of the shape of the Earth for it's calculations. GeoTools uses the GeographicLib implementation of C. F. F. Karney, Algorithms for geodesics, J. Geodesy 87, 43–55 (2013), which explains the approximations used for solving the problem.

    This answer on gis.stackexchange.com explains the level of accuracy you can expect from various numbers of decimal places when using latitude and longitude is also summarized in this XKCD cartoon:

    enter image description here

    Your lowest precision point is 4DP so you can't expect much better than 10's of metres from the rest of the calculation. You are unlikely to get much better than 5 DP of actual measurement even in the aviation domain and more likely to be working with 3DP and 10s-100s metres accuracy.

    Update

    Further investigation suggests that it is your original distance value that is wrong.

      calculator.setStartingGeographicPoint(startLongitude, startLatitude);
  calculator.setDestinationGeographicPoint(endLongitude, endLatitude);
  double azimuth = calculator.getAzimuth();
  System.out.println("Azimuth=" + azimuth);
  double fullDistance = calculator.getOrthodromicDistance();
  System.out.println("distance " + fullDistance);
  System.out.println("% error " + (Math.abs(fullDistance - distance) / fullDistance) * 100);
  calculator = new GeodeticCalculator(crs);
  calculator.setStartingGeographicPoint(startLongitude, startLatitude);
  calculator.setDirection(azimuth, fullDistance);
  Point2D computedEndPoint = calculator.getDestinationGeographicPoint();
  System.out.println("computedEndPoint=" + computedEndPoint);

  calculator = new GeodeticCalculator(crs);
  calculator.setStartingGeographicPoint(endLongitude, endLatitude);
  calculator.setDestinationGeographicPoint(computedEndPoint);
  distance = calculator.getOrthodromicDistance();
  System.out.println("Distance=" + distance);
  System.out.println("% error " + ((distance / fullDistance) * 100));

    gives me:

    Azimuth=44.971973670068415
distance 130893.86215735915
% error 0.04098575880994952
computedEndPoint=Point2D.Double[-4.000599999999989, 41.49682]
Distance=9.64120596409175E-10
% error 7.365666964965247E-13

    As you can see all the error comes in the 1st calculation, the forward/reverse trip gives an identical point and 9e-10m error in distance. That should be fine for any domain.