19

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.

1
  • (-24680.1613, 6708860.65389) <- in which coordinate system is that point given? Commented Jun 26, 2014 at 11:54

3 Answers 3

24

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.
5
  • Hi Jan-Philip, The purpose is to figure out the coordinate for each one and use it to draw a rectangle 1km around the center spot. Commented Jun 26, 2014 at 11:55
  • So you want to get four points based on the origin, each one offset by 1 km into one of the four cardinal directions. You should adjust your question accordingly ;) Commented Jun 26, 2014 at 11:58
  • And then there is a nice writeup here, with reference information and implementations: movable-type.co.uk/scripts/latlong.html Commented Jun 26, 2014 at 14:02
  • 1
    now if only geopy had the opposite, the ability to find the vincenty direction between two points...
    – Michael
    Commented Oct 7, 2014 at 18:00
  • @Jan-PhilipGehrcke I have a similar question, except instead of moving 'north', I'm trying to move between two points. Is this possible? Question is here: stackoverflow.com/questions/38619835/…
    – tempomax
    Commented Jul 27, 2016 at 21:03
20

A 2020 update for this question, based on Dr. Jan-Philip Gehrcke's answer.

VincentyDistance is officially deprecated, and was never fully precise and sometimes inaccurate.

This snippet shows how to use with the latest (and future versions of GeoPy - Vincenty will be deprecated in 2.0)

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 is a new Point object, which when printed, returns 48 51m 43.1721s N, 2 20m 56.4s E

Which as you can see is more accurate than Vincenty, and should maintain better accuracy near the poles.

Hope it helps!

1

I had to deal with adding meters to longitude and latitude.

Here's what I did, inspired by this source :

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
1
  • But it seems your initial_location is not used. I suppose you need to replace location by initial_location or the opposite.
    – mike123
    Commented Apr 18, 2020 at 10:40

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.