13

I have a line with two points in latitude and longitude
A: 3.222895, 101.719751
B: 3.227511, 101.724318

and 1 point
C: 3.224972, 101.722932

How can I calculate minimum distance between point C and a line consists of point A and B? It will be convenient if you can provide the calculation and objective-c code too. The distance is around 89 meters (using ruler in Google Earth).

Azam
  • 797
  • 2
  • 8
  • 23
  • possible duplicate of [Distance from point to line on Earth](http://stackoverflow.com/questions/7803004/distance-from-point-to-line-on-earth) – Lior Kogan Nov 29 '13 at 16:28
  • 1
    that question doesn't have any points to calculate. I need calculation example. – Azam Dec 01 '13 at 02:42

5 Answers5

16

Thanks to mimi and this great article http://www.movable-type.co.uk/scripts/latlong.html but they don't give the whole picture. Here is a detail one. All this points are collected using Google Earth using Placemark to mark the locations. Make sure lat/long are set to decimal degrees in Preferences.

lat A = 3.222895  
lon A = 101.719751  
lat B = 3.222895  
lon B = 101.719751  
lat C = 3.224972  
lon C = 101.722932  
Earth radius, R = 6371

1. First you have to find the bearing from A to C and A to B.
Bearing formula 

bearingAC = atan2( sin(Δλ)*cos(φ₂), cos(φ₁)*sin(φ₂) − sin(φ₁)*cos(φ₂)*cos(Δλ) )  
bearingAB = atan2( sin(Δλ)*cos(φ₂), cos(φ₁)*sin(φ₂) − sin(φ₁)*cos(φ₂)*cos(Δλ) )

φ is latitude, λ is longitude, R is earth radius

2. Find A to C distance using spherical law of cosines 

distanceAC = acos( sin(φ₁)*sin(φ₂) + cos(φ₁)*cos(φ₂)*cos(Δλ) )*R

3. Find cross-track distance 

distance = asin(sin(distanceAC/ R) * sin(bearingAC − bearing AB)) * R

Objective-C code 

double lat1 = 3.227511;
double lon1 = 101.724318;
double lat2 = 3.222895;
double lon2 = 101.719751;
double lat3 = 3.224972;
double lon3 = 101.722932;

double y = sin(lon3 - lon1) * cos(lat3);
double x = cos(lat1) * sin(lat3) - sin(lat1) * cos(lat3) * cos(lat3 - lat1);
double bearing1 = radiansToDegrees(atan2(y, x));
bearing1 = 360 - ((bearing1 + 360) % 360);

double y2 = sin(lon2 - lon1) * cos(lat2);
double x2 = cos(lat1) * sin(lat2) - sin(lat1) * cos(lat2) * cos(lat2 - lat1);
double bearing2 = radiansToDegrees(atan2(y2, x2));
bearing2 = 360 - ((bearing2 + 360) % 360);

double lat1Rads = degreesToRadians(lat1);
double lat3Rads = degreesToRadians(lat3);
double dLon = degreesToRadians(lon3 - lon1);

double distanceAC = acos(sin(lat1Rads) * sin(lat3Rads)+cos(lat1Rads)*cos(lat3Rads)*cos(dLon)) * 6371;  
double min_distance = fabs(asin(sin(distanceAC/6371)*sin(degreesToRadians(bearing1)-degreesToRadians(bearing2))) * 6371);

NSLog(@"bearing 1: %g", bearing1);  
NSLog(@"bearing 2: %g", bearing2);  
NSLog(@"distance AC: %g", distanceAC);  
NSLog(@"min distance: %g", min_distance);

Actually there's a library for this. You can find it here https://github.com/100grams/CoreLocationUtils

stuyam
  • 9,561
  • 2
  • 21
  • 40
Azam
  • 797
  • 2
  • 8
  • 23
  • 2
    Thank you! I believe you have minor errors in the code (computing bearing) -- the lines should say `360 - ((bearing + 360) % 360)`, because module has higher priority than addition. – astrowalker Dec 11 '16 at 17:51
  • Maybe I have made some stupid typo, but... For a line/arc the readout should always be the same, no matter if you call `func(point,start,end)` or `func(point,end,start)`. And for `start = { Latitude = 36.496902, Longitude = -93.223531 }`, `end = { Latitude = 36.496902, Longitude = -90.151775 }` and `point = { Latitude = 36.501227, Longitude = -93.368405 }` I have different (vastly) readouts. I found also http://stackoverflow.com/questions/32771458/distance-from-lat-lng-point-to-minor-arc-segment which gives distance to arc **segment** and after modifications it runs fine. – greenoldman Mar 04 '17 at 16:04
  • Once you've calculated the minimum distance, is there a way to get the bearing to travel so that you only need to travel minimum distance to the curve? You have bearing to A and bearing to C, but could you provide bearing to the intersection point of the original curve and the minimum distance curve. – Josh Gafni Mar 25 '17 at 19:56
  • 1
    There seems to be mistakes in the code. Sometimes you calculate the sine of an angle in degrees (for example when you calculate bearings), and sometimes angles in radians (on the last two lines). Both can't be right! Still thanks for your answer, it's nicely detailed! – Dici Nov 13 '19 at 23:47
3

Calculate bearing for each: C to A , and C to B:

var y = Math.sin(dLon) * Math.cos(lat2);
var x = Math.cos(lat1)*Math.sin(lat2) -
        Math.sin(lat1)*Math.cos(lat2)*Math.cos(dLon);
var brng = Math.atan2(y, x).toDeg();

dLon= lon2-lon1;

Calculate cross-track distance:

var dXt = Math.asin(Math.sin(distance_CB/R)*Math.sin(bearing_CA-bearing_CB)) * R;

R is the radius of earth, dXt is the minimum distance you wanted to calculate.

1

Code to carry out this calculation is posted at here. This implements an accurate solution in terms of ellipsoidal geodesics. For the basic geodesic calculations, you can use GeographicLib or the port of these algorithms to C which are included in version 4.9.0 of PROJ.4. This C interface is documented here.

Here's the result of compiling and running intercept.cpp:

$ echo 3.222895 101.719751 3.227511 101.724318 3.224972 101.722932 | ./intercept 
Initial guess 3.225203 101.7220345
Increment 0.0003349040566247297 0.0003313413822354505
Increment -4.440892098500626e-16 0
Increment 0 0
...
Final result 3.225537904056624 101.7223658413822
Azimuth to A1 -135.1593040635131
Azimuth to A2 44.84069593652217
Azimuth to B1 134.8406959363608

Distance to line is 88.743m:

$ echo 3.224972 101.722932 3.225537904056624 101.7223658413822 | GeodSolve -i
-45.15927221 -45.15930407 88.743
cffk
  • 1,839
  • 1
  • 12
  • 17
1

See post here: https://stackoverflow.com/a/33343505/4083623

For distance up to a few thousands meters I would simplify the issue from sphere to plane. Then, the issue is pretty simply as a easy triangle calculation can be used:

We have points A and B and look for a distance X to line AB. Then:

Location a;
Location b;
Location x;

double ax = a.distanceTo(x);
double alfa = (Math.abs(a.bearingTo(b) - a.bearingTo(x))) / 180
            * Math.PI;
double distance = Math.sin(alfa) * ax;
Community
  • 1
  • 1
Martin Koubek
  • 433
  • 4
  • 8
0

If you know how to calculate the distance of two points, get the distances between each two points, you get AB, AC, and BC. You want to know the closest distance between point C and line AB.

First get the value of P

P=（AB+BC+AC)/2

Using P, you need to get S

S=SQRT(（P(P-AC)（P-AB)(P-AC))

SQRT means square root. Then you get what you want by

2*S/AB
David Lin
  • 13,168
  • 5
  • 46
  • 46
  • 1
    it won't be precise because latitude and longitude is not on a flat plane. it's on a big sphere (the earth). – Azam Nov 27 '13 at 01:42
  • You are right, I should point it out. Though I believe the discrepancy caused by the earth's round surface is small if the distances are calculated with this considered. – David Lin Nov 27 '13 at 02:31