根据距离和方向计算点

Question

我想用GeoDjango或GeoPy来根据方向和距离计算一个点。

举个例子，如果我有一个点坐标是(-24680.1613, 6708860.65389)，我想知道在这个点的北边1公里、东边1公里、南边1公里和西边1公里的点，想用Vincenty距离公式来计算。

我找到的最接近的东西是distance.py里的一个“destination”函数（https://code.google.com/p/geopy/source/browse/trunk/geopy/distance.py?r=105）。不过我在任何地方都找不到这个函数的文档，而且我还没弄明白怎么使用它。

任何帮助都非常感谢。

Answer 1

我需要处理在经度和纬度上添加米数的问题。

这是我做的事情，灵感来自于这个来源：

import math
from geopy.distance import vincenty

initial_location = '50.966086,5.502027'
lat, lon = (float(i) for i in location.split(','))
r_earth = 6378000
lat_const = 180 / math.pi
lon_const = lat_const / math.cos(lat * math.pi / 180)

# dx = distance in meters on x axes (longitude)
dx = 1000
new_longitude = lon + (dx / r_earth) * lon_const
new_longitude = round(new_longitude, 6)
new_latitude = lat + (dy / r_earth) * lat_const
new_latitude = round(new_latitude, 6)

# dy = distance on y axes (latitude)
new_latitude = lat + (dy / r_earth) * lat_const
new_latitude = round(new_latitude, 6)
new_location = ','.join([str(y_lat), str(x_lon)])

dist_to_location = vincenty(location, new_location).meters

Answer 2

这是关于这个问题的2020年更新，基于Dr. Jan-Philip Gehrcke的回答。

VincentyDistance这个方法现在已经不推荐使用了，而且它从来没有完全准确，有时候还会出错。

下面的代码片段展示了如何在最新版本（以及未来版本）的GeoPy中使用，因为在2.0版本中，Vincenty将会被淘汰。

import geopy
import geopy.distance

# Define starting point.
start = geopy.Point(48.853, 2.349)

# Define a general distance object, initialized with a distance of 1 km.
d = geopy.distance.distance(kilometers=1)

# Use the `destination` method with a bearing of 0 degrees (which is north)
# in order to go from point `start` 1 km to north.
final = d.destination(point=start, bearing=0)

final是一个新的Point对象，当你打印它的时候，会显示48 51m 43.1721s N, 2 20m 56.4s E。

正如你所看到的，这个结果比Vincenty更准确，并且在靠近极地的地方也能保持更好的准确性。

希望这对你有帮助！

Answer 3

编辑 2

好吧，其实用 geopy 有一个现成的解决方案，只是文档不太完善：

import geopy
import geopy.distance

# Define starting point.
start = geopy.Point(48.853, 2.349)

# Define a general distance object, initialized with a distance of 1 km.
d = geopy.distance.VincentyDistance(kilometers = 1)

# Use the `destination` method with a bearing of 0 degrees (which is north)
# in order to go from point `start` 1 km to north.
print d.destination(point=start, bearing=0)

输出结果是 48 52m 0.0s N, 2 21m 0.0s E（或者 Point(48.861992239749355, 2.349, 0.0)）。

90度表示东，180度表示南，依此类推。

旧的回答：

一个简单的解决方案是：

def get_new_point():
    # After going 1 km North, 1 km East, 1 km South and 1 km West
    # we are back where we were before.
    return (-24680.1613, 6708860.65389)

不过，我不确定这是否能完全满足你的需求。

好吧，认真说，你可以开始使用 geopy。首先，你需要在 geopy 能识别的坐标系统中定义你的起点。乍一看，你似乎不能简单地“添加”一个特定的距离到某个方向。我认为原因在于，计算距离是一个没有简单反向解的问题。或者说，我们怎么能反转 measure 函数呢？这个函数定义在 https://code.google.com/p/geopy/source/browse/trunk/geopy/distance.py#217。

因此，你可能需要采用一种迭代的方法。

正如这里所述： https://stackoverflow.com/a/9078861/145400，你可以这样计算两个给定点之间的距离：

pt1 = geopy.Point(48.853, 2.349)
pt2 = geopy.Point(52.516, 13.378)
# distance.distance() is the  VincentyDistance by default.
dist = geopy.distance.distance(pt1, pt2).km

如果你想向北走一公里，你需要逐步将纬度增加，并检查距离。你可以使用 SciPy 中的简单迭代求解器来自动化这个过程：只需通过列表中的优化器找到 geopy.distance.distance().km - 1 的根，具体可以参考 http://docs.scipy.org/doc/scipy/reference/optimize.html#root-finding。

我想很明显，向南时需要将纬度改为负方向，向西和向东则是改变经度。

我对这种地理计算没有经验，这种迭代方法只有在没有简单直接的“向北”某个距离的方式时才有意义。

编辑：我提议的一个示例实现：

import geopy
import geopy.distance
import scipy.optimize


def north(startpoint, distance_km):
    """Return target function whose argument is a positive latitude
    change (in degrees) relative to `startpoint`, and that has a root
    for a latitude offset that corresponds to a point that is 
    `distance_km` kilometers away from the start point.
    """
    def target(latitude_positive_offset):
        return geopy.distance.distance(
            startpoint, geopy.Point(
                latitude=startpoint.latitude + latitude_positive_offset,
                longitude=startpoint.longitude)
            ).km - distance_km
    return target


start = geopy.Point(48.853, 2.349)
print "Start: %s" % start

# Find the root of the target function, vary the positve latitude offset between
# 0 and 2 degrees (which is for sure enough for finding a 1 km distance, but must
# be adjusted for larger distances).
latitude_positive_offset = scipy.optimize.bisect(north(start, 1),  0, 2)


# Build Point object for identified point in space.
end = geopy.Point(
    latitude=start.latitude + latitude_positive_offset,
    longitude=start.longitude
    )

print "1 km north: %s" % end

# Make the control.
print "Control distance between both points: %.4f km." % (
     geopy.distance.distance(start, end).km)

输出：

$ python test.py 
Start: 48 51m 0.0s N, 2 21m 0.0s E
1 km north: 48 52m 0.0s N, 2 21m 0.0s E
Control distance between both points: 1.0000 km.