半正矢公式(两个 GPS 点之间的方位和距离)

问题

我想知道如何得到 两个 GPS 点之间的距离和方位。 我研究过半正矢公式。 有人告诉我，我也可以找到方位使用相同的数据。

剪辑

一切正常，但轴承还没有完全正常工作。轴承输出负，但应在0-360度之间。 设定数据应使水平轴承 96.02166666666666 是:

Start point: 53.32055555555556 , -1.7297222222222221
Bearing:  96.02166666666666
Distance: 2 km
Destination point: 53.31861111111111, -1.6997222222222223
Final bearing: 96.04555555555555

这是我的新代码:

from math import *


Aaltitude = 2000
Oppsite  = 20000


lat1 = 53.32055555555556
lat2 = 53.31861111111111
lon1 = -1.7297222222222221
lon2 = -1.6997222222222223


lon1, lat1, lon2, lat2 = map(radians, [lon1, lat1, lon2, lat2])


dlon = lon2 - lon1
dlat = lat2 - lat1
a = sin(dlat/2)**2 + cos(lat1) * cos(lat2) * sin(dlon/2)**2
c = 2 * atan2(sqrt(a), sqrt(1-a))
Base = 6371 * c




Bearing =atan2(cos(lat1)*sin(lat2)-sin(lat1)*cos(lat2)*cos(lon2-lon1), sin(lon2-lon1)*cos(lat2))


Bearing = degrees(Bearing)
print ""
print ""
print "--------------------"
print "Horizontal Distance:"
print Base
print "--------------------"
print "Bearing:"
print Bearing
print "--------------------"




Base2 = Base * 1000
distance = Base * 2 + Oppsite * 2 / 2
Caltitude = Oppsite - Aaltitude


a = Oppsite/Base
b = atan(a)
c = degrees(b)


distance = distance / 1000


print "The degree of vertical angle is:"
print c
print "--------------------"
print "The distance between the Balloon GPS and the Antenna GPS is:"
print distance
print "--------------------"
205668 次浏览
小开
最佳答案

下面是 Python 的一个版本:

from math import radians, cos, sin, asin, sqrt


def haversine(lon1, lat1, lon2, lat2):
"""
Calculate the great circle distance in kilometers between two points
on the earth (specified in decimal degrees)
"""
# convert decimal degrees to radians
lon1, lat1, lon2, lat2 = map(radians, [lon1, lat1, lon2, lat2])


# haversine formula
dlon = lon2 - lon1
dlat = lat2 - lat1
a = sin(dlat/2)**2 + cos(lat1) * cos(lat2) * sin(dlon/2)**2
c = 2 * asin(sqrt(a))
r = 6371 # Radius of earth in kilometers. Use 3956 for miles. Determines return value units.
return c * r
小开

你可以通过增加360 ° 来解决负向轴承的问题。 不幸的是，这可能导致轴承大于360 ° 的积极轴承。 这是模运算符的一个很好的候选者，所以总的来说，您应该添加这一行

Bearing = (Bearing + 360) % 360

在你的方法结束的时候。

小开

方位计算不正确，需要将输入交换到 atan2。

    bearing = atan2(sin(long2-long1)*cos(lat2), cos(lat1)*sin(lat2)-sin(lat1)*cos(lat2)*cos(long2-long1))
bearing = degrees(bearing)
bearing = (bearing + 360) % 360

这会给你正确的方位。

小开

你可以试试 haversine 软件包: Https://pypi.org/project/haversine/

示例代码:

from haversine import haversine
haversine((45.7597, 4.8422),(48.8567, 2.3508), unit='mi')
243.71209416020253
小开

在缺省情况下，atan2中的 Y 是第一个参数。这是 文件。您将需要切换您的输入，以获得正确的轴承角度。

bearing = atan2(sin(lon2-lon1)*cos(lat2), cos(lat1)*sin(lat2)in(lat1)*cos(lat2)*cos(lon2-lon1))
bearing = degrees(bearing)
bearing = (bearing + 360) % 360
小开

请参考此链接: https://gis.stackexchange.com/questions/84885/whats-the-difference-between-vincenty-and-great-circle-distance-calculations

这实际上给出了两种获得距离的方法。他们是哈弗辛和文森特。根据我的研究，我发现文森特的描述相对准确。还可以使用 import 语句实现。

小开

这里有两个计算距离和方位的函数，它们基于以前消息中的代码和 https://gist.github.com/jeromer/2005586(为了清晰起见，为地理点添加了 lat、 lon 格式的元组类型)。我测试了这两种功能，它们似乎工作正常。

#coding:UTF-8
from math import radians, cos, sin, asin, sqrt, atan2, degrees


def haversine(pointA, pointB):


if (type(pointA) != tuple) or (type(pointB) != tuple):
raise TypeError("Only tuples are supported as arguments")


lat1 = pointA[0]
lon1 = pointA[1]


lat2 = pointB[0]
lon2 = pointB[1]


# convert decimal degrees to radians
lat1, lon1, lat2, lon2 = map(radians, [lat1, lon1, lat2, lon2])


# haversine formula
dlon = lon2 - lon1
dlat = lat2 - lat1
a = sin(dlat/2)**2 + cos(lat1) * cos(lat2) * sin(dlon/2)**2
c = 2 * asin(sqrt(a))
r = 6371 # Radius of earth in kilometers. Use 3956 for miles
return c * r




def initial_bearing(pointA, pointB):


if (type(pointA) != tuple) or (type(pointB) != tuple):
raise TypeError("Only tuples are supported as arguments")


lat1 = radians(pointA[0])
lat2 = radians(pointB[0])


diffLong = radians(pointB[1] - pointA[1])


x = sin(diffLong) * cos(lat2)
y = cos(lat1) * sin(lat2) - (sin(lat1)
* cos(lat2) * cos(diffLong))


initial_bearing = atan2(x, y)


# Now we have the initial bearing but math.atan2 return values
# from -180° to + 180° which is not what we want for a compass bearing
# The solution is to normalize the initial bearing as shown below
initial_bearing = degrees(initial_bearing)
compass_bearing = (initial_bearing + 360) % 360


return compass_bearing


pA = (46.2038,6.1530)
pB = (46.449, 30.690)


print haversine(pA, pB)


print initial_bearing(pA, pB)
小开

这些答案大部分都是围绕地球半径的。如果您将这些函数与其他距离计算器(如 geopy)进行比较，则这些函数将关闭。

这种方法很有效:

from math import radians, cos, sin, asin, sqrt


def haversine(lat1, lon1, lat2, lon2):


R = 3959.87433 # this is in miles.  For Earth radius in kilometers use 6372.8 km


dLat = radians(lat2 - lat1)
dLon = radians(lon2 - lon1)
lat1 = radians(lat1)
lat2 = radians(lat2)


a = sin(dLat/2)**2 + cos(lat1)*cos(lat2)*sin(dLon/2)**2
c = 2*asin(sqrt(a))


return R * c


# Usage
lon1 = -103.548851
lat1 = 32.0004311
lon2 = -103.6041946
lat2 = 33.374939


print(haversine(lat1, lon1, lat2, lon2))
小开

下面是@Michael Dunn 给出的一个简单的矢量化半正矢公式实现，相对于大矢量，它提供了10-50倍的改进。

from numpy import radians, cos, sin, arcsin, sqrt


def haversine(lon1, lat1, lon2, lat2):
"""
Calculate the great circle distance between two points
on the earth (specified in decimal degrees)
"""


#Convert decimal degrees to Radians:
lon1 = np.radians(lon1.values)
lat1 = np.radians(lat1.values)
lon2 = np.radians(lon2.values)
lat2 = np.radians(lat2.values)


#Implementing Haversine Formula:
dlon = np.subtract(lon2, lon1)
dlat = np.subtract(lat2, lat1)


a = np.add(np.power(np.sin(np.divide(dlat, 2)), 2),
np.multiply(np.cos(lat1),
np.multiply(np.cos(lat2),
np.power(np.sin(np.divide(dlon, 2)), 2))))
c = np.multiply(2, np.arcsin(np.sqrt(a)))
r = 6371


return c*r
小开

还有一个 向量化实现，它允许使用4个数字数组代替坐标的标量值:

def distance(s_lat, s_lng, e_lat, e_lng):


# approximate radius of earth in km
R = 6373.0


s_lat = s_lat*np.pi/180.0
s_lng = np.deg2rad(s_lng)
e_lat = np.deg2rad(e_lat)
e_lng = np.deg2rad(e_lng)


d = np.sin((e_lat - s_lat)/2)**2 + np.cos(s_lat)*np.cos(e_lat) * np.sin((e_lng - s_lng)/2)**2


return 2 * R * np.arcsin(np.sqrt(d))
小开

考虑到你的目标是测量两点之间的距离(由地理坐标表示) ，将留下以下三个选项:

  1. 半正矢公式

  2. 使用 GeoPy测地距离

  3. 使用 GeoPy大圆距离

选择一

半正矢公式将会完成这项工作，但是需要注意的是，这样做的结果是将地球近似为一个球体，并且有一个误差(看看这个答案)——因为地球不是一个球体。

为了使用半正矢公式首先需要确定地球的半径。这本身可能会引起一些争议。考虑到以下三个来源

我将使用值 6371公里作为地球半径的参考。

# Radius of the Earth
r = 6371.0

我们将利用 math模块。

在半径之后，一个移动到坐标，一个开始将坐标转换成弧度，以使用 数学是三角函数。对于这种情况，导入 math.radians(x)并按以下方式使用它们

#Import radians from math module
from math import radians


# Latitude and Longitude for the First Point (let's consider 40.000º and 21.000º)
lat1 = radians(40.000)
lon1 = radians(21.000)


# Latitude and Longitude for the Second Point (let's consider 30.000º and 25.000º)
lat2 = radians(30.000)
lon2 = radians(25.000)

现在人们已经准备好应用半正矢公式了。第一个是用点1的经度减去点2的经度

dlon = lon2 - lon1
dlat = lat2 - lat1

然后，在这里有几个三角函数，一个是要使用，更具体地说，math.sin()，math.cos()，和 math.atan2()。我们也将使用 math.sqrt()

# Import sin, cos, atan2, and sqrt from math module
from math import sin, cos, atan2, sqrt


a = sin(dlat / 2)**2 + cos(lat1) * cos(lat2) * sin(dlon / 2)**2
c = 2 * atan2(sqrt(a), sqrt(1 - a))
d = r * c

然后通过打印 d得到距离。

这可能会有所帮助，让我们将所有内容集中在一个函数中(受到 @ Michael Dunn 的回答的启发)

from math import radians, cos, sin, atan2, sqrt


def haversine(lon1, lat1, lon2, lat2):
"""
Calculate the great-circle distance (in km) between two points
using their longitude and latitude (in degrees).
"""
# Radius of the Earth
r = 6371.0


# Convert degrees to radians
# First point
lat1 = radians(lat1)
lon1 = radians(lon1)


# Second Point
lat2 = radians(lat2)
lon2 = radians(lon2)


# Haversine formula
dlon = lon2 - lon1
dlat = lat2 - lat1
a = sin(dlat / 2)**2 + cos(lat1) * cos(lat2) * sin(dlon / 2)**2
c = 2 * atan2(sqrt(a), sqrt(1 - a))
return r * c

选择二

一个是使用 GeoPy 的距离，更确切地说，是 geodesic。

我们可以在公里或英里(来源)上获得结果

# Import Geopy's distance
from geopy import distance


wellington = (-41.32, 174.81)
salamanca = (40.96, -5.50)
print(distance.distance(wellington, salamanca).km) # If one wants in miles, change `km` to `miles`


[Out]: 19959.6792674

选择三

一个是使用 GeoPy 的距离，更确切地说，是 great-circle。

我们可以在公里或英里(来源)上获得结果

# Import Geopy's distance
from geopy import distance


newport_ri = (41.49008, -71.312796)
cleveland_oh = (41.499498, -81.695391)


print(distance.great_circle(newport_ri, cleveland_oh).miles) # If one wants in km, change `miles` to `km`


[Out]: 536.997990696

提交答案