Calculate geographic distance between points using NetTopologySuite

Viewed 1137

I am trying to calculate distance between two points using NetTopologySuite. Since I am referencing Microsoft documentation, I Came up with following GeometryExtension and GeometryHelper classes:

        public static class GeometryExtensions
    {
        private static readonly CoordinateSystemServices _coordinateSystemServices = new CoordinateSystemServices(new Dictionary<int, string>
            {
                // Coordinate systems:

                [4326] = GeographicCoordinateSystem.WGS84.WKT,

                // CRS for Denmark ESPG 25832. Source: https://epsg.io/25832 and https://sdfe.dk/
                [25832] = @"PROJCS[""ETRS89 / UTM zone 32N"",
                                GEOGCS[""ETRS89"",
                                    DATUM[""European_Terrestrial_Reference_System_1989"",
                                        SPHEROID[""GRS 1980"",6378137,298.257222101,
                                            AUTHORITY[""EPSG"",""7019""]],
                                        TOWGS84[0,0,0,0,0,0,0],
                                        AUTHORITY[""EPSG"",""6258""]],
                                    PRIMEM[""Greenwich"",0,
                                        AUTHORITY[""EPSG"",""8901""]],
                                    UNIT[""degree"",0.0174532925199433,
                                        AUTHORITY[""EPSG"",""9122""]],
                                    AUTHORITY[""EPSG"",""4258""]],
                                PROJECTION[""Transverse_Mercator""],
                                PARAMETER[""latitude_of_origin"",0],
                                PARAMETER[""central_meridian"",9],
                                PARAMETER[""scale_factor"",0.9996],
                                PARAMETER[""false_easting"",500000],
                                PARAMETER[""false_northing"",0],
                                UNIT[""metre"",1,
                                    AUTHORITY[""EPSG"",""9001""]],
                                AXIS[""Easting"",EAST],
                                AXIS[""Northing"",NORTH],
                                AUTHORITY[""EPSG"",""25832""]]"
            }
        );
        


        /// <summary>
        /// Projects a geometry to another SRID
        /// </summary>
        /// <param name="geometry"></param>
        /// <param name="targetSrid"></param>
        /// <param name="defaultSourceSrid">If the geometry SRID has not been specified (i.e. equals 0) defaultSourceSrid is used</param>
        /// <returns></returns>
        public static Geometry ProjectTo(this Geometry geometry, int targetSrid, int? defaultSourceSrid = null)
        {
            if (geometry == null)
                throw new Exception("Geometry is null, cannot project");

            var sourceSrid = geometry.SRID == 0 && defaultSourceSrid.HasValue ? defaultSourceSrid.Value : geometry.SRID;
            var transformation = _coordinateSystemServices.CreateTransformation(sourceSrid, targetSrid);

            var result = geometry.Copy();
            result.Apply(new MathTransformFilter(transformation.MathTransform));

            return result;
        }
}




     public  class GeometryHelper
    {
        private static readonly int _longLatSRID = 4326;
        private static readonly int _targetSRID = 25832;

        /// <summary>
        /// In order to get the distance in meters, we need to project to an appropriate
        /// coordinate system. If no SRID is provided 25832, which covers Denmark is used.
        /// If the provided Points have no SRID, 4326 (longitude/latitude) is assumed.
        /// </summary>
        /// <param name="a"></param>
        /// <param name="b"></param>
        /// <param name="targetSrid"></param>
        /// <returns></returns>
        public static double DistanceInMeters(Point a, Point b, int? targetSrid = null)
        {
            targetSrid ??= _targetSRID;

            try
            {
                //If SRID is not set, assume long/lat, ie. 4326
                return a.ProjectTo(targetSrid.Value, _longLatSRID).Distance(b.ProjectTo(targetSrid.Value, _longLatSRID));
            }
            catch (Exception e)
            {
                throw new Exception("Failed to calculate distance", e);
            }
        }
}

To test if my distance is being calculated properly, I opened Google Maps and selected the calculate distance between 3 different pairs of points. Then I copy pasted the values in my tests to see if they match. The values I am getting from my code is 2 times larger than what I am getting from Google Maps. What am I doing wrong? Here are my tests:

public class GeometryHelperTests
{
    [Theory]
    [InlineData(55.676518126293466, 12.567203858066554, 55.67645068023813, 12.56698376863065, 15.68)]//I am getting 36.1
    [InlineData(55.659368700924475, 12.546625254609248, 55.65940085355421, 12.546601114728679, 3.77)]//I am getting 6.2
    [InlineData(55.65896978705746, 12.546674114514795, 55.6596855501795, 12.547258269821455, 87.09)]//I am getting 173.7
    public void DistanceInMeters(double x1, double y1, double x2, double y2, double expected)
    {
        // setup
        Point point1 = new Point(x1, y1) { SRID = 4326 };
        Point point2 = new Point(x2, y2) { SRID = 4326 };

        //CoordinateSystemFactory csFact = new CoordinateSystemFactory();
        

        // act
        var result = GeometryHelper.DistanceInMeters(point1, point2);

        

        // assert
        result.Should().Be(expected);
    }

}
2 Answers

You need to calculate great circle distance. NetTopologySuite Point.Distance method returns the cartesian distance.

Try the following:

public static double Radians(double x)
{
    return x * Math.PI / 180;
}

public static double GreatCircleDistance(double lon1, double lat1, double lon2, double lat2)
{
    double R = 6371e3; // m

    double sLat1 = Math.Sin(Radians(lat1));
    double sLat2 = Math.Sin(Radians(lat2));
    double cLat1 = Math.Cos(Radians(lat1));
    double cLat2 = Math.Cos(Radians(lat2));
    double cLon = Math.Cos(Radians(lon1) - Radians(lon2));

    double cosD = sLat1*sLat2 + cLat1*cLat2*cLon;

    double d = Math.Acos(cosD);

    double dist = R * d;

    return dist;
}

You can also use the built-in GeoCoordinate.GetDistanceTo() which implements the Haversine formula that is more accurate for small distances.

Even though saeedkazemi's answer is correct and gives results closer to the real value, results given by my code were also quite close, but I figured out that I had to flip the coordinates given by Google maps. So if I was given 55.6765, 12.5672 I need to give 12.5672,55.6765 to the formula. That being said, saeedkazemi's answer provides results closer to the real value.

Related