How to compute line segment intersections on a map

Question

I am working with Latitude / Longitude coordinates in a google map.

I have two lines :

• Line A : 48.31508162629726, -2.591741396838972 to 48.40216156645915, -2.2218462112093404
• Line B : 48.383816077371215, -2.274292940053768 to 48.66103546935337, -1.7066197241571377

I then use the following formula to find the point where they cross.

var XAsum = A.LngStart - A.LngEnd;
var XBsum = B.LngStart - B.LngEnd;
var YAsum = A.LatStart - A.LatEnd;
var YBsum = B.LatStart - B.LatEnd;

var LineDenominator = XAsum * YBsum - YAsum * XBsum;
if(LineDenominator == 0.0)
    return false;

var a = A.LngStart * A.LatEnd - A.LatStart * A.LngEnd;
var b = B.LngStart * B.LatEnd - B.LatStart * B.LngEnd;

var x = (a * XBsum - b * XAsum) / LineDenominator;
var y = (a * YBsum - b * YAsum) / LineDenominator;

This tells me that the lines do indeed cross and returns the x and y values.

However, when I plot the returned point, it is offset (not much) from the real intersection.

Is there a better and just as fast algorithm I could use which will return me the correct intersection point ?

It needs to be fast as I am iterating over a large number of lines (~1000).

EDIT : Note this is giving me an error offset of around 7.5 meters

Equal

What language are you programming in? Are there any libraries you can use (or are using) to help you solve GIS/spatial-analysis problems? – Pelvic

Did you solve your problem? I have the same challenge to find the intersection of two lines (arcs) on google maps. The Cartesian intersection produces errors on the map. For each arc I have a pair of latlons and I need to find their intersection on google maps. – Fabri

No, I had to live with the error. Sorry – Equal

Try with negative lat (south part or earth). It don't work at all. I'll have a look at this (I'm in Brasil) – Loam

It seems when LineDenominator = 0, line are stricly paralell. But even if they are not, the question is: do you want to get intersection of lines or of segment? Because here, you seem to het intersection of the line and not of the segment. So in some case you get an intersection coord which is far from the segments. – Loam

Answer 1

I'm assuming the algorithm you're using is the one for finding line intersections on a Cartesian coordinate system (i.e. a plane). Unfortunately, the earth is not a plane (or even a sphere) so using that algorithm will introduce error. Google Maps uses an ellipsoid (specifically WGS84) to approximate the surface of the earth, so you'll need an algorithm for finding the intersections of arcs on an ellipsoid.

This page might contains some helpful information: http://mathhelpforum.com/calculus/90196-point-intersection-two-lines.html

Pelvic

thank you for taking the trouble to answer. You are right in your assumptions. I looked at the page you mention, but I am not sufficiently advanced (!) in maths to be able to translate what is indicated to my use. Would you have any further suggestions ? – Equal

Any particular language you care to use? – Pelvic

Didn't have time to do much research, but this might help? geotools.org Or maybe this? geospatialmethods.org/spheres – Pelvic

Answer 2

After searching and trying to use your code, I discover the main problem came from the fact you're not checking for intersection of segment. In many case the coord of the intersection are outside of the segment, meanig you're checking for intersection of line.

Here is a piece of PHP code. Assuming you have 2 lines in GPS Coord (x1y1-x2y2 for the first one and x3y3-x4y4 for the second one)

$x1 = deg2rad($x1);
$y1 = deg2rad($y1);

$x2 = deg2rad($x2);
$y2 = deg2rad($y2); 

$x3 = deg2rad($x3);
$y3 = deg2rad($y3);

$x4 = deg2rad($x4);
$y4 = deg2rad($y4);     

$denom = (($y4 - $y3) * ($x2 - $x1)) - (($x4 - $x3) * ($y2 - $y1));

// Id denom = 0 -> parallele line
if ($denom != 0)
{
    $ua = (($x4 - $x3) * ($y1 - $y3) - ($y4 - $y3) * ($x1 - $x3))/$denom;
    $ub = (($x2 - $x1) * ($y1 - $y3) - ($y2 - $y1) * ($x1 - $x3))/$denom;

    if ($ua >= 0 && $ua <= 1 && $ub >= 0 && $ub <= 1)
    {
        $intersection_x = ($x1 + $ua*($x2 - $x1));
        $intersection_y = ($y1 + $ua*($y2 - $y1));
        echo "<br><b>Intersection: ".round(rad2deg($intersection_x),6).",".round(rad2deg($intersection_y),6)."</b>";
    }
    else
    {
        echo "<br>No intersection";
    }
} 

Loam