coordinatescoordinate-systemsef-core-6.0sqlgeographynettopologysuite

How to measure distance between two Points in metres with NetTopologySuite / ProjNet (2022)


This question has been posted multiple times before but the answers no longer appear to match the method signatures available in the latest NuGet packages.

After spending days trying to piece together an answer from other places (1, 2, 3, 4, 5, 6, blogs and even Microsoft's official spiel et al...) this is as far as I've got, but it won't run:

@using NetTopologySuite
@using NetTopologySuite.Geometries
@using ProjNet
@using ProjNet.CoordinateSystems

@using MyApp.Helpers

@{
var windsorCastle = Geography.CreatePoint(51.483884, -0.604455);
var buckinghamPalace = Geography.CreatePoint(51.501576, -0.141208);

var csWgs84 = ProjNet.CoordinateSystems.GeographicCoordinateSystem.WGS84;
const string epsg27700 = @"
                        PROJCS[""OSGB36 / British National Grid"",
                        GEOGCS[""OSGB36"",
                            DATUM[""Ordnance_Survey_of_Great_Britain_1936"",
                                SPHEROID[""Airy 1830"",6377563.396,299.3249646],
                                EXTENSION[""PROJ4_GRIDS"",""OSTN15_NTv2_OSGBtoETRS.gsb""]],
                            PRIMEM[""Greenwich"",0,
                                AUTHORITY[""EPSG"",""8901""]],
                            UNIT[""degree"",0.0174532925199433,
                                AUTHORITY[""EPSG"",""9122""]],
                            AUTHORITY[""EPSG"",""4277""]],
                        PROJECTION[""Transverse_Mercator""],
                        PARAMETER[""latitude_of_origin"",49],
                        PARAMETER[""central_meridian"",-2],
                        PARAMETER[""scale_factor"",0.9996012717],
                        PARAMETER[""false_easting"",400000],
                        PARAMETER[""false_northing"",-100000],
                        UNIT[""metre"",1,
                            AUTHORITY[""EPSG"",""9001""]],
                        AXIS[""Easting"",EAST],
                        AXIS[""Northing"",NORTH],
                        AUTHORITY[""EPSG"",""27700""]]
                        "; // see http://epsg.io/27700

var cs27700 = (CoordinateSystem) ProjNet.IO.CoordinateSystems.CoordinateSystemWktReader.Parse(epsg27700);
var ctFactory = new ProjNet.CoordinateSystems.Transformations.CoordinateTransformationFactory();

var ct = ctFactory.CreateFromCoordinateSystems(csWgs84, cs27700);
var mt = ct.MathTransform;

var gf = new NetTopologySuite.Geometries.GeometryFactory(new PrecisionModel(PrecisionModels.Floating), 27700);

Console.WriteLine(windsorCastle.Distance(buckinghamPalace)); // Need metres, not degrees...
}

It blows up with this:

ArgumentException: Expecting (',') but got a '[' at line 6 column 47.
    ProjNet.IO.CoordinateSystems.WktStreamTokenizer.ReadToken(string expectedToken)
    ProjNet.IO.CoordinateSystems.CoordinateSystemWktReader.ReadGeographicCoordinateSystem(WktStreamTokenizer tokenizer)
    ProjNet.IO.CoordinateSystems.CoordinateSystemWktReader.ReadProjectedCoordinateSystem(WktStreamTokenizer tokenizer)
    ProjNet.IO.CoordinateSystems.CoordinateSystemWktReader.ReadCoordinateSystem(string coordinateSystem, WktStreamTokenizer tokenizer)
    ProjNet.IO.CoordinateSystems.CoordinateSystemWktReader.Parse(string wkt)
    MyApp.Pages.Pages_DistanceTest.ExecuteAsync() in DistanceTest.cshtml

var cs27700 = (CoordinateSystem) ProjNet.IO.CoordinateSystems.CoordinateSystemWktReader.Parse(epsg27700);

Any pointers on a simpler/correct solution would be welcome. I am new to the coordinate system so please be patient. Previously I've done this using TSQL with geography::Point().STDistance but the goal is an EF Core-only end result.


Solution

  • The CoordinateSystem WKT contains an EXTENSION for the GEOGCS|DATUM|SPHEROID definition. You have to remove that, as ProjNET can't handle the EXTENSION tag.

    For the coordinate transformation you need the following instructions:

    var cwc = mt.Transform(windsorCastle.X, windsorCastle.Y);
    var cbp = mt.Transform(buckinghamPalace.X, buckinghamPalace.Y);
    
    Console.WriteLine(new Coordinate(cwc.x, cwc.y).Distance(new Coordinate(cbp.x, cbp.y)))