python - Calculate point based on distance and direction

Question

I would like to calculate a point based on direction and distance using GeoDjango or GeoPy.

For example, If I have a point that is (-24680.1613, 6708860.65389) I would like to find out a point 1KM North, 1KM East, 1KM Sourh and 1KM west using Vincenty distance formula.

I closest thing I can find is a "destination" function in distance.py (https://code.google.com/p/geopy/source/browse/trunk/geopy/distance.py?r=105). Although I cannot find this documented anywhere and I'm yet to figure out how to use it.

Any help is much appreciated.

Answer 1

Edit 2

Okay, there is an out-of-the-box solution with geopy, it is just not well-documented:

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)

The output is 48 52m 0.0s N, 2 21m 0.0s E (or Point(48.861992239749355, 2.349, 0.0)).

A bearing of 90 degrees corresponds to East, 180 degrees is South, and so on.

Older answers:

A simple solution would be:

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)

However, I am not sure if that serves your purposes in all generality.

Okay, seriously, you can get started using geopy. First of all, you need to define your starting point in a coordinate system known to geopy. At the first glance, it seems that you cannot just "add" a certain distance into a certain direction. The reason, I think, is that calculation of the distance is a problem without simple inverse solution. Or how would we invert the measure function defined in https://code.google.com/p/geopy/source/browse/trunk/geopy/distance.py#217?

Hence, you might want to take an iterative approach.

As stated here: https://stackoverflow.com/a/9078861/145400 you can calculate the distance between two given points like that:

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

For going north by one kilometer you would iteratively change the latitude into a positive direction and check against the distance. You can automate this approach using a simple iterative solver from e.g. SciPy: just find the root of geopy.distance.distance().km - 1 via one of the optimizers listed in http://docs.scipy.org/doc/scipy/reference/optimize.html#root-finding.

I think it is clear that you go south by changing the latitude into a negative direction, and west and east by changing the longitude.

I have no experience with such geo calculations, this iterative approach only makes sense if there is no simple direct way for "going north" by a certain distance.

Edit: an example implementation of my proposal:

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)

Output:

$ 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.