7

I have a line segments defined with a start and an end point:

A: 
x1 = 10.7196405787775
y1 = 59.9050401935882

B: 
x2 = 10.7109989561813
y2 = 59.9018650448204

where x defines longitude and y defines latitude.

I also have a point:

P:
x0 = 10.6542116666667
y0 = 59.429105

How do I compute the shortest distance between the line segment and the point? I know how to do this in Cartesian coordinates, but not in long/lat coordinates.

davidism
  • 121,510
  • 29
  • 395
  • 339
bjornasm
  • 2,211
  • 7
  • 37
  • 62
  • 2
    google 'haversine': http://en.wikipedia.org/wiki/Haversine_formula – William Pursell Dec 13 '14 at 17:34
  • 1
    For GIS type questions you might consider posting to [Stack Exchange - Geographic Information Systems](http://gis.stackexchange.com/). Then post actual coding problems here. This may be of interest - [What tools in Python are available for doing great circle distance + line creation?](http://gis.stackexchange.com/q/47/23174) – wwii Dec 13 '14 at 18:24
  • 1
    http://stackoverflow.com/a/865080/948550 – Reut Sharabani Dec 13 '14 at 18:50
  • @ReutSharabani, not only the other post deals with Cartesian coordinates, not spherical, but also the very linked answer answers different question (distance to a line, not to line segment). – greenoldman Jan 13 '19 at 13:23

2 Answers2

0

Here is an implementation of a formula off Wikipedia:

def distance(p0, p1, p2): # p3 is the point
    x0, y0 = p0
    x1, y1 = p1
    x2, y2 = p2
    nom = abs((y2 - y1) * x0 - (x2 - x1) * y0 + x2 * y1 - y2 * x1)
    denom = ((y2 - y1)**2 + (x2 - x1) ** 2) ** 0.5
    result = nom / denom
    return result

print distance((0, 0), (3, 4), (5, 6))

# should probably test less obvious cases
assert 1 == distance((0, 0), (1, 1), (2, 1))
# change 0.001 to whatever accuracy you demand on testing.
# Notice that it's never perfect...
assert 0.5 * (2 ** 0.5) - distance((0, 0), (1, 0), (0, 1)) < 0.001
Reut Sharabani
  • 30,449
  • 6
  • 70
  • 88
-1

Using the helpful Python geocoding library geopy, and the formula for the midpoint of a great circle from Chris Veness's geodesy formulae, we can find the distance between a great circle arc and a given point:

from math import sin, cos, atan2, sqrt, degrees, radians, pi
from geopy.distance import great_circle as distance
from geopy.point import Point


def midpoint(a, b):
    a_lat, a_lon = radians(a.latitude), radians(a.longitude)
    b_lat, b_lon = radians(b.latitude), radians(b.longitude)
    delta_lon = b_lon - a_lon
    B_x = cos(b_lat) * cos(delta_lon)
    B_y = cos(b_lat) * sin(delta_lon)
    mid_lat = atan2(
        sin(a_lat) + sin(b_lat),
        sqrt(((cos(a_lat) + B_x)**2 + B_y**2))
    )
    mid_lon = a_lon + atan2(B_y, cos(a_lat) + B_x)
    # Normalise
    mid_lon = (mid_lon + 3*pi) % (2*pi) - pi
    return Point(latitude=degrees(mid_lat), longitude=degrees(mid_lon))

Which in this example gives:

# Example:
a = Point(latitude=59.9050401935882, longitude=10.7196405787775)
b = Point(latitude=59.9018650448204, longitude=10.7109989561813)
p = Point(latitude=59.429105, longitude=10.6542116666667)

d = distance(midpoint(a, b), p)
print d.km
# 52.8714586903
nelfin
  • 1,019
  • 7
  • 14
  • 14
    This gives the distance to the midpoint. That is not the same as the distance to the line segment - for example, if p is close to a, d should be far less than half the distance from a to b. – Skyler Oct 18 '16 at 03:44