两点之间的正交距离
与计算地球上两点之间的大圆距离相关的计算算法
引言
计算地球表面两点之间的大圆（正交）距离是地理信息系统 (GIS) 的核心问题之一。这项看似简单的任务需要相当复杂的算法解决方案。 确实，如果问题涉及平面几何，那么毕达哥拉斯定理将提供一种“不费吹灰之力”的解决方案。 但实际的 GIS 计算处理的是 3D 模型，即球形地球的表示，这需要更复杂的解决方案。 另一个复杂层级与更精确的椭球地球模型有关，这对于大多数实际应用来说有些“过度”。 球形模型会导致系统误差范围在 0.3% 以内，这对于大多数商业级应用程序来说是可以接受的。 另一种（即地球的椭球模型）理论上可以将误差范围限制在毫米级，但会大大增加计算复杂度。 以下解决方案基于球形地球模型，描述了 3 种用 C# 编写的实用算法，它们在计算性能/精度方面有所不同 [1-3]。
背景
从数学上讲，下面描述的所有三种算法都计算地球上两点之间的大圆（正交）距离，尽管精度和性能不同。 它们都基于地球的球形模型，并提供合理的近似值，误差范围通常不超过纽约市边界内的几米。 更精确的椭球地球模型和相应的 Vincenty 高精度解决方案 [4] 存在，可以将误差范围降低到毫米级，但也会大大增加计算复杂度，超出合理水平。 因此，基于球形地球模型的 3 种算法已被开发并推荐用于通用商业应用程序，具有良好的计算性能和合理的精度。
使用代码
以下是与计算地球表面两点之间的大圆（正交）距离相关的三种算法解决方案
列表 1. 用于计算地球表面两点之间距离的基本 GIS 函数
/*****************************************************************************************
Module : GIS.cs |Class Lib
Description : Calculate distance between two geo-points on surface
*****************************************************************************************
Author : Alexander Bell
Copyright : 2011-2015 Infosoft International Inc
*****************************************************************************************
DISCLAIMER : This Module is provided on AS IS basis without any warranty
*****************************************************************************************
TERMS OF USE : This module is copyrighted. Please keep the Copyright notice intact.
*****************************************************************************************/
using System;
namespace BusNY
{
internal enum UnitSystem { SI = 0, US = 1 }
internal static class GIS
{
#region internal: properties (read-only)
internal static double EarthRadiusKm { get {return _radiusEarthKM;} }
internal static double EarthRadiusMiles { get { return _radiusEarthMiles; } }
internal static double m2km { get { return _m2km; } }
internal static double Deg2rad { get { return _toRad; } }
#endregion
#region private: const
private const double _radiusEarthMiles = 3959;
private const double _radiusEarthKM = 6371;
private const double _m2km = 1.60934;
private const double _toRad = Math.PI / 180;
#endregion
#region Method 1: Haversine algo
/// <summary>
/// Distance between two geographic points on surface, km/miles
/// Haversine formula to calculate
/// great-circle (orthodromic) distance on Earth
/// High Accuracy, Medium speed
/// re: http://en.wikipedia.org/wiki/Haversine_formula
/// </summary>
/// <param name="Lat1">double: 1st point Latitude</param>
/// <param name="Lon1">double: 1st point Longitude</param>
/// <param name="Lat2">double: 2nd point Latitude</param>
/// <param name="Lon2">double: 2nd point Longitude</param>
/// <returns>double: distance, km/miles</returns>
internal static double DistanceHaversine(double Lat1,
double Lon1,
double Lat2,
double Lon2,
UnitSystem UnitSys ){
try {
double _radLat1 = Lat1 * _toRad;
double _radLat2 = Lat2 * _toRad;
double _dLatHalf = (_radLat2 - _radLat1) / 2;
double _dLonHalf = Math.PI * (Lon2 - Lon1) / 360;
// intermediate result
double _a = Math.Sin(_dLatHalf);
_a *= _a;
// intermediate result
double _b = Math.Sin(_dLonHalf);
_b *= _b * Math.Cos(_radLat1) * Math.Cos(_radLat2);
// central angle, aka arc segment angular distance
double _centralAngle = 2 * Math.Atan2(Math.Sqrt(_a + _b), Math.Sqrt(1 - _a - _b));
// great-circle (orthodromic) distance on Earth between 2 points
if (UnitSys == UnitSystem.SI) { return _radiusEarthKM * _centralAngle; }
else { return _radiusEarthMiles * _centralAngle; }
}
catch { throw; }
}
#endregion
#region Method 2: Spherical Law of Cosines
/// <summary>
/// Distance between two geographic points on surface, km/miles
/// Spherical Law of Cosines formula to calculate
/// great-circle (orthodromic) distance on Earth;
/// High Accuracy, Medium speed
/// re: http://en.wikipedia.org/wiki/Spherical_law_of_cosines
/// </summary>
/// <param name="Lat1">double: 1st point Latitude</param>
/// <param name="Lon1">double: 1st point Longitude</param>
/// <param name="Lat2">double: 2nd point Latitude</param>
/// <param name="Lon2">double: 2nd point Longitude</param>
/// <returns>double: distance, km/miles</returns>
internal static double DistanceSLC(double Lat1,
double Lon1,
double Lat2,
double Lon2,
UnitSystem UnitSys ){
try {
double _radLat1 = Lat1 * _toRad;
double _radLat2 = Lat2 * _toRad;
double _radLon1 = Lon1 * _toRad;
double _radLon2 = Lon2 * _toRad;
// central angle, aka arc segment angular distance
double _centralAngle = Math.Acos(Math.Sin(_radLat1) * Math.Sin(_radLat2) +
Math.Cos(_radLat1) * Math.Cos(_radLat2) * Math.Cos(_radLon2 - _radLon1));
// great-circle (orthodromic) distance on Earth between 2 points
if (UnitSys == UnitSystem.SI) { return _radiusEarthKM * _centralAngle; }
else { return _radiusEarthMiles * _centralAngle; }
}
catch { throw; }
}
#endregion
#region Method 3: Spherical Earth projection
/// <summary>
/// Distance between two geographic points on surface, km/miles
/// Spherical Earth projection to a plane formula (using Pythagorean Theorem)
/// to calculate great-circle (orthodromic) distance on Earth.
/// central angle =
/// Sqrt((_radLat2 - _radLat1)^2 + (Cos((_radLat1 + _radLat2)/2) * (Lon2 - Lon1))^2)
/// Medium Accuracy, Fast,
/// relative error less than 0.1% in search area smaller than 250 miles
/// re: http://en.wikipedia.org/wiki/Geographical_distance
/// </summary>
/// <param name="Lat1">double: 1st point Latitude</param>
/// <param name="Lon1">double: 1st point Longitude</param>
/// <param name="Lat2">double: 2nd point Latitude</param>
/// <param name="Lon2">double: 2nd point Longitude</param>
/// <returns>double: distance, km/miles</returns>
public static double DistanceSEP(double Lat1,
double Lon1,
double Lat2,
double Lon2,
UnitSystem UnitSys ){
try
{
double _radLat1 = Lat1 * _toRad;
double _radLat2 = Lat2 * _toRad;
double _dLat = (_radLat2 - _radLat1);
double _dLon = (Lon2 - Lon1) * _toRad;
double _a = (_dLon) * Math.Cos((_radLat1 + _radLat2) / 2);
// central angle, aka arc segment angular distance
double _centralAngle = Math.Sqrt(_a * _a + _dLat * _dLat);
// great-circle (orthodromic) distance on Earth between 2 points
if (UnitSys == UnitSystem.SI) { return _radiusEarthKM * _centralAngle; }
else { return _radiusEarthMiles * _centralAngle; }
}
catch { throw; }
}
#endregion
}
}
关注点
上述算法已在实时纽约市公交跟踪应用程序的上下文中部分描述，目前可在网上获取。
历史
- 2012 年 10 月：该主题在物联网竞赛提交文章中进行了简要讨论。
- 2015 年 6 月：已发布扩展版本