使用 NetTopologySuite 计算点之间的地理距离
Calculate geographic distance between points using NetTopologySuite
我正在尝试使用 NetTopologySuite 计算两点之间的距离。由于我参考了 Microsoft 文档，因此我想出了以下 GeometryExtension 和 GeometryHelper 类:
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);
}
}
}
为了测试我的距离计算是否正确，我打开了 Google 地图并选择了计算 3 对不同点之间的距离。然后我复制粘贴测试中的值以查看它们是否匹配。我从我的代码中获得的值比我从 Google 地图中获得的值大 2 倍。我究竟做错了什么？这是我的测试：
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);
}
}
您需要计算great circle distance。
NetTopologySuite Point.Distance 方法 returns cartesian distance.
尝试以下操作：
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;
}
您还可以使用内置的 GeoCoordinate.GetDistanceTo() which implements the Haversine 公式，该公式对于小距离更准确。
尽管 saeedkazemi 的答案是正确的并且给出的结果更接近真实值，但我的代码给出的结果也非常接近，但我发现我必须翻转 Google 地图给出的坐标。所以如果我得到 55.6765, 12.5672 我需要给公式 12.5672,55.6765。话虽如此，saeedkazemi 的回答提供的结果更接近真实值。
我正在尝试使用 NetTopologySuite 计算两点之间的距离。由于我参考了 Microsoft 文档，因此我想出了以下 GeometryExtension 和 GeometryHelper 类:
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);
}
}
}
为了测试我的距离计算是否正确，我打开了 Google 地图并选择了计算 3 对不同点之间的距离。然后我复制粘贴测试中的值以查看它们是否匹配。我从我的代码中获得的值比我从 Google 地图中获得的值大 2 倍。我究竟做错了什么？这是我的测试：
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);
}
}
您需要计算great circle distance。
NetTopologySuite Point.Distance 方法 returns cartesian distance.
尝试以下操作：
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;
}
您还可以使用内置的 GeoCoordinate.GetDistanceTo() which implements the Haversine 公式，该公式对于小距离更准确。
尽管 saeedkazemi 的答案是正确的并且给出的结果更接近真实值，但我的代码给出的结果也非常接近，但我发现我必须翻转 Google 地图给出的坐标。所以如果我得到 55.6765, 12.5672 我需要给公式 12.5672,55.6765。话虽如此，saeedkazemi 的回答提供的结果更接近真实值。