Google Maps Two Circles Intersection Points

Question

Is there an easy way to get the lat/lng of the intersection points (if available) of two circles in Google Maps API V3? Or should I go with the hard way?

EDIT : In my problem, circles always have the same radius, in case that makes the solution easier.

Answer 1

Yazanpro, sorry for the late response on this.

You may be interested in a concise variant of MBo's approach, which simplifies in two respects :

• firstly by exploiting some of the built in features of the google.maps API to avoid much of the hard math.
• secondly by using a 2D model for the calculation of the included angle, in place of MBo's spherical model. I was initially uncertain about the validity of this simplification but satisfied myself with tests in a fork of MBo's fiddle that the errors are minor at all but the largest of circles with respect to the size of the Earth (eg at low zoom levels).

Here's the function :

function getIntersections(circleA, circleB) {
    /* 
     * Find the points of intersection of two google maps circles or equal radius
     * circleA: a google.maps.Circle object 
     * circleB: a google.maps.Circle object
     * returns: null if 
     *    the two radii are not equal 
     *    the two circles are coincident
     *    the two circles don't intersect
     * otherwise returns: array containing the two points of intersection of circleA and circleB
     */

    var R, centerA, centerB, D, h, h_;

    try {

        R = circleA.getRadius();
        centerA = circleA.getCenter();
        centerB = circleB.getCenter();

        if(R !== circleB.getRadius()) {
            throw( new Error("Radii are not equal.") );
        }
        if(centerA.equals(centerB)) {
            throw( new Error("Circle centres are coincident.") );
        }

        D = google.maps.geometry.spherical.computeDistanceBetween(centerA, centerB); //Distance between the two centres (in meters)

        // Check that the two circles intersect
        if(D > (2 * R)) {
            throw( new Error("Circles do not intersect.") );
        }

        h = google.maps.geometry.spherical.computeHeading(centerA, centerB); //Heading from centre of circle A to centre of circle B. (in degrees)
        h_ = Math.acos(D / 2 / R) * 180 / Math.PI; //Included angle between the intersections (for either of the two circles) (in degrees). This is trivial only because the two radii are equal.

        //Return an array containing the two points of intersection as google.maps.latLng objects
        return [
            google.maps.geometry.spherical.computeOffset(centerA, R, h + h_),
            google.maps.geometry.spherical.computeOffset(centerA, R, h - h_)
        ];
    }
    catch(e) {
        console.error("getIntersections() :: " + e.message);
        return null;
    }
}

No disrespect to MBo by the way - it's an excellent answer.

Answer 2

Yes, for equal circles rather simple solution could be elaborated:
Let's first circle center is A point, second circle center is F, midpoint is C, and intersection points are B,D. ABC is right-angle spherical triangle with right angle C.

We want to find angle A - this is deviation angle from A-F direction. Spherical trigonometry (Napier's rules for right spherical triangles) gives us formula:

cos(A)= tg(AC) * ctg(AB) where one symbol denote spherical angle, double symbols denote great circle arcs' angles (AB, AC). We can see that AB = circle radius (in radians, of course), AC = half-distance between A and F on the great circle arc. To find AC (and other values) - I'll use code from this excellent page

var R = 6371; // km
var dLat = (lat2-lat1).toRad();
var dLon = (lon2-lon1).toRad();
var lat1 = lat1.toRad();
var lat2 = lat2.toRad();

var a = Math.sin(dLat/2) * Math.sin(dLat/2) +
        Math.sin(dLon/2) * Math.sin(dLon/2) * Math.cos(lat1) * Math.cos(lat2); 
var c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1-a)); 

and our

AC = c/2

If circle radius Rd is given is kilometers, then

AB = Rd / R = Rd / 6371

Now we can find angle

A = arccos(tg(AC) * ctg(AB))

Starting bearing (AF direction):

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);

Intersection points' bearings:

B_bearing = brng - A
D_bearing = brng + A

Intersection points' coordinates:

var latB = Math.asin( Math.sin(lat1)*Math.cos(Rd/R) + 
              Math.cos(lat1)*Math.sin(Rd/R)*Math.cos(B_bearing) );
var lonB = lon1.toRad() + Math.atan2(Math.sin(B_bearing)*Math.sin(Rd/R)*Math.cos(lat1), 
                     Math.cos(Rd/R)-Math.sin(lat1)*Math.sin(lat2));

and the same for D_bearing

latB, lonB are in radians

Answer 3

The computation the "hard" way can be simplified for the case r1 = r2 =: r. We still first have to convert the circle centers P1,P2 from (lat,lng) to Cartesian coordinates (x,y,z).

var DEG2RAD = Math.PI/180;
function LatLng2Cartesian(lat_deg,lng_deg)
{
  var lat_rad = lat_deg*DEG2RAD;
  var lng_rad = lng_deg*DEG2RAD;
  var cos_lat = Math.cos(lat_rad);
  return {x: Math.cos(lng_rad)*cos_lat,
          y: Math.sin(lng_rad)*cos_lat,
          z: Math.sin(lat_rad)};
}

var P1 = LatLng2Cartesian(lat1, lng1);
var P2 = LatLng2Cartesian(lat2, lng2);

But the intersection line of the planes holding the circles can be computed more easily. Let d be the distance of the actual circle center (in the plane) to the corresponding point P1 or P2 on the surface. A simple derivation shows (with R the earth's radius):

var R = 6371; // earth radius in km
var r = 100; // the direct distance (in km) of the given points to the intersections points
// if the value rs for the distance along the surface is known, it has to be converted:
// var r = 2*R*Math.sin(rs/(2*R*Math.PI));
var d = r*r/(2*R);

Now let S1 and S2 be the intersections points and S their mid-point. With s = |OS| and t = |SS1| = |SS2| (where O = (0,0,0) is the earth's center) we get from simple derivations:

var a = Math.acos(P1.x*P2.x + P1.y*P2.y + P1.z*P2.z); // the angle P1OP2
var s = (R-d)/Math.cos(a/2);
var t = Math.sqrt(R*R - s*s);

Now since r1 = r2 the points S, S1, S2 are in the mid-plane between P1 and P2. For v_s = OS we get:

function vecLen(v)
{ return Math.sqrt(v.x*v.x + v.y*v.y + v.z*v.z); }
function vecScale(scale,v)
{ return {x: scale*v.x, y: scale*v.y, z: scale*v.z}; }

var v = {x: P1.x+P2.x, y: P1.y+P2.y, z:P1.z+P2.z}; // P1+P2 is in the middle of OP1 and OP2
var S = vecScale(s/vecLen(v), v);

function crossProd(v1,v2)
{
  return {x: v1.y*v2.z - v1.z*v2.y,
          y: v1.z*v2.x - v1.x*v2.z,
          z: v1.x*v2.y - v1.y*v2.x};
}
var n = crossProd(P1,P2); // normal vector to plane OP1P2 = vector along S1S2
var SS1 = vecScale(t/vecLen(n),n);

var S1 = {x: S.x+SS1.x, y: S.y+SS1.y, z: S.z+SS1.z}; // S + SS1
var S2 = {x: S.x-SS1.x, y: S.y-SS2.y, z: S.z-SS1.z}; // S - SS1

Finally we have to convert back to (lat,lng):

function Cartesian2LatLng(P)
{
  var P_xy = {x: P.x, y:P.y, z:0}
  return {lat: Math.atan2(P.y,P.x)/DEG2RAD, lng: Math.atan2(P.z,vecLen(P_xy))/DEG2RAD};
}
var S1_latlng = Cartesian2LatLng(S1);
var S2_latlng = Cartesian2LatLng(S2);