locationtech / spatial4j

LocationTech Spatial4j: A Geospatial Library for Java
https://projects.eclipse.org/projects/locationtech.spatial4j
935 stars 168 forks source link

A question about distance calculations #221

Open StFS opened 2 years ago

StFS commented 2 years ago

I'm quite confused about distance calculations.

I'm using this website as a reference: https://www.omnicalculator.com/other/latitude-longitude-distance

Then I got an implementation of a distance calculation (haversine) from here: https://tutorialspoint.dev/algorithm/geometric-algorithms/program-distance-two-points-earth

And I want to compare this with various different implementations from Spatial4J.

I'm using these coordinates to test with: Point 1: 40.688939, -74.04455 Point 2: 40.746853, -73.985633

And I'm getting pretty big differences between different ways of calculating the distance between these two points.

First, what mainly concerns me is that the website (omnicalculator) and the implementation from tutorialspoint agree completely on the distance: 8.132 km

But none of my Spatial4J calculations agree with that number. The one that comes closest to it is the CartesianDistCalc implementation at 8.262 km. The tutorialspoint demo code claims to be using haversine but the output of Spatial4J haversine DistCalc implementation is quite a ways off at 7.313 km.

But can somebody explain to me where these differences are coming from and what the "correct" one is?

Below is my experimental code:

import org.junit.jupiter.api.Test;
import org.locationtech.spatial4j.context.SpatialContext;
import org.locationtech.spatial4j.distance.CartesianDistCalc;
import org.locationtech.spatial4j.distance.GeodesicSphereDistCalc;

class GeoDesicCalculationTest {
    @Test
    void testGeoDesicCalculations(){
        SpatialContext ctx = SpatialContext.GEO;

        var startPoint = ctx.getShapeFactory().pointLatLon(40.688939, -74.04455);
        var endPoint = ctx.getShapeFactory().pointLatLon(40.746853, -73.985633);

        System.out.println("GEO spatial context:         " + ctx.calcDistance(startPoint, endPoint) * 100);
        System.out.println("Haversine:                   " + new GeodesicSphereDistCalc.Haversine().distance(startPoint, endPoint) * 100);
        System.out.println("Law of cosine:               " + new GeodesicSphereDistCalc.LawOfCosines().distance(startPoint, endPoint) * 100);
        System.out.println("Vincenty:                    " + new GeodesicSphereDistCalc.Vincenty().distance(startPoint, endPoint) * 100);
        System.out.println("Cartesian:                   " + new CartesianDistCalc().distance(startPoint, endPoint) * 100);
        System.out.println("Tutorials Point (haversine): " + distance(startPoint.getLat(), endPoint.getLat(), startPoint.getLon(), endPoint.getLon()));
    }

    public static double distance(double lat1, double lat2, double lon1, double lon2) {
        // The math module contains a function
        // named toRadians which converts from
        // degrees to radians.
        lon1 = Math.toRadians(lon1);
        lon2 = Math.toRadians(lon2);
        lat1 = Math.toRadians(lat1);
        lat2 = Math.toRadians(lat2);

        // Haversine formula
        double dlon = lon2 - lon1;
        double dlat = lat2 - lat1;
        double a = Math.pow(Math.sin(dlat / 2), 2)
                + Math.cos(lat1) * Math.cos(lat2)
                * Math.pow(Math.sin(dlon / 2),2);

        double c = 2 * Math.asin(Math.sqrt(a));

        // Radius of earth in kilometers. Use 3956
        // for miles
        double r = 6371;

        // calculate the result
        return(c * r);
    }
}

And the output of running it:

GEO spatial context:         7.31307025220976
Haversine:                   7.31307025220976
Law of cosine:               7.313070251733588
Vincenty:                    7.3130702522095286
Cartesian:                   8.261503667613857
Tutorials Point (haversine): 8.131763102409689

I'm multiplying the Spatial4J calculations by 100 which is also confusing to me... it doesn't really make sense that Spatial4J is giving me answers as 1/100th of a kilometer???

I realise I must be doing something wrong or completely misunderstanding some premises here. I would really appreciate some help in understanding what I'm doing wrong.