Determining if a lat-long rect and a circle on a sphere overlap

Question

Suppose I have the following:

• A region defined by minimum and maximum latitude and longitude (commonly a 'lat-long rect', though it's not actually rectangular except in certain projections).
• A circle, defined by a center lat/long and a radius

How can I determine:

1. Whether the two shapes overlap?
2. Whether the circle is entirely contained within the rect?

I'm looking for a complete formula/algorithm, rather than a lesson in the math, per-se.

Answer 1

EDIT I retract my answer since it was mathematically incorrect. Sorry, I didn't give the problem enough thought.

Answer 2

warning: this can be tricky if the circles / "rectangles" span large portions of the sphere, e.g.:

"rectangle": min long = -90deg, max long = +90deg, min lat = +70deg, max lat = +80deg

circle: center = lat = +85deg, long = +160deg, radius = 20deg (e.g. if point A is on the circle and point C is the circle's center, and point O is the sphere's center, then angle AOC = 40deg).

These intersect but the math is likely to have several cases to check intersection/containment. The following points lie on the circle described above: P1=(+65deg lat,+160deg long), P2=(+75deg lat, -20deg long). P1 is outside the "rectangle" and P2 is inside the "rectangle" so the circle/"rectangle" intersect in at least 2 points.

OK, here's my shot at an outline of the solution:

---

Let C = circle center with radius R (expressed as a spherical angle as above). C has latitude LATC and longitude LONGC. Since the word "rectangle" is kind of misleading here (lines of constant latitude are not segments of great circles), I'll use the term "bounding box".

• function InsideCircle(P) returns +1,0,or -1: +1 if point P is inside the circle, 0 if point P is on the circle, and -1 if point P is outside the circle: calculation of great-circle distance D (expressed as spherical angle) between C and any point P will tell you whether or not P is inside the circle: InsideCircle(P) = sign(R-D) (As user @Die in Sente mentioned, great circle distances have been asked on this forum elsewhere)

• Define PANG(x) = the principal angle of x = MOD(x+180deg, 360deg)-180deg. PANG(x) is always between -180deg and +180deg, inclusive (+180deg should map to -180deg).

• To define the bounding box, you need to know 4 numbers, but there's a slight issue with longitude. LAT1 and LAT2 represent the bounding latitudes (assuming LAT1 < LAT2); there's no ambiguity there. LONG1 and LONG2 represent the bounding longitudes of a longitude interval, but this gets tricky, and it's easier to rewrite this interval as a center and width, with LONGM = the center of that interval and LONGW = width. NOTE that there are always 2 possibilities for longitude intervals. You have to specify which of these cases it is, whether you are including or excluding the 180deg meridian, e.g. the shortest interval from -179deg to +177deg has LONGM = +179deg and LONGW = 4deg, but the other interval from -179deg to +177deg has LONGM = -1deg and LONGW = 356deg. If you naively try to do "regular" comparisons with the interval [-179,177] you will end up using the larger interval and that's probably not what you want. As an aside, point P, with latitude LATP and longitude LONGP, is inside the bounding box if both of the following are true:

  • LAT1 <= LATP and LATP <= LAT2 (that part is obvious)
  • abs(PANG(LONGP-LONGM)) < LONGW/2

The circle intersects the bounding box if ANY of the following points P in PTEST = union(PCORNER,PLAT,PLONG) as described below, do not all return the same result for InsideCircle():

• PCORNER = the bounding box's 4 corners
• the points PLAT on the bounding box's sides (there are either none or 2) which share the same latitude as the circle's center, if LATC is between LAT1 and LAT2, in which case these points have the latitude LATC and longitude LONG1 and LONG2.
• the points PLONG on the bounding box's sides (there are either none or 2 or 4!) which share the same longitude as the circle's center. These points have EITHER longitude = LONGC OR longitude PANG(LONGC-180). If abs(PANG(LONGC-LONGM)) < LONGW/2 then LONGC is a valid longitude. If abs(PANG(LONGC-180-LONGM)) < LONGW/2 then PANG(LONGC-180) is a valid longitude. Either or both or none of these longitudes may be within the longitude interval of the bounding box. Choose points PLONG with valid longitudes, and latitudes LAT1 and LAT2.

These points PLAT and PLONG as listed above are the points on the bounding box that are "closest" to the circle (if the corners are not; I use "closest" in quotes, in the sense of lat/long distance and not great-circle distance), and cover the cases where the circle's center lies on one side of the bounding box's boundary but points on the circle "sneak across" the bounding box boundary.

If all points P in PTEST return InsideCircle(P) == +1 (all inside the circle) then the circle contains the bounding box in its entirety.

If all points P in PTEST return InsideCircle(P) == -1 (all outside the circle) then the circle is contained entirely within the bounding box.

Otherwise there is at least one point of intersection between the circle and the bounding box. Note that this does not calculate where those points are, although if you take any 2 points P1 and P2 in PTEST where InsideCircle(P1) = -InsideCircle(P2), then you could find a point of intersection (inefficiently) by bisection. (If InsideCircle(P) returns 0 then you have a point of intersection, though equality in floating-point math is generally not to be trusted.)

There's probably a more efficient way to do this but the above should work.

Answer 3

How about this?

Find vector v that connects the center of the rectangle, point Cr, to the center of the circle. Find point i where v intersects the rectangle. If ||i-Cr|| + r > ||v|| then they intersect.

In other words, the length of the segment inside the rectangle plus the length of the segment inside the circle should be greater than the total length (of v, the center-connecting line segment).

Finding point i should be the tricky part, especially if it falls on a longitude edge, but you should be able to come up with something faster than I can.

Edit: This method can't tell if the circle is completely within the rectangle. You'd need to find the distance from its center to all four of the rectangle's edges for that.

Edit: The above is incorrect. There are some cases, as Federico Ramponi has suggested, where it does not work even in Euclidean geometry. I'll post another answer. Please unaccept this and feel free to vote down. I'll delete it shortly.

Answer 4

One more try at this...

I think the solution is to test a set of points, just as Jason S has suggested, but I disagree with his selection of points, which I think is mathematically wrong.

You need to find the points on the sides of the lat/long box where the distance to the center of the circle is a local minimum or maximum. Add those points to the set of corners and then the algorithm above should be correct.

I.e, letting longitude be the x dimension and latitude be the y dimension, let each side of the box be a parametric curve P(t) = P0 + t (P1-P0) for o <= t <= 1.0, where P0 and P1 are two adjacent corners.

Let f(P) = f(P.x, P.y) be the distance from the center of the cirle.

Then f (P0 + t (P1-P0)) is a distance function of t: g(t). Find all the points where the derivative of the distance funtion is zero: g'(t) == 0. (Discarding solutions outsize the domain 0 <= t <= 1.0, of course)

Unfortunately, this needs to find the zero of a transcendental expression, so there's no closed form solution. This type of equation can only solved by Newton-Raphson iteration.

OK, I realize that you wanted code, not the math. But the math is all I've got.

Answer 5

For the Euclidean geometry answer, see: Circle-Rectangle collision detection (intersection)

Answer 6

This should work for any points on earth. If you want to change it to a different size sphere just change the kEarchRadiusKms to whatever radius you want for your sphere.

This method is used to calculate the distance between to lat and lon points.

I got this distance formula from here: http://www.codeproject.com/csharp/distancebetweenlocations.asp

public static double Calc(double Lat1, double Long1, double Lat2, double Long2)
{
    double dDistance = Double.MinValue;
    double dLat1InRad = Lat1 * (Math.PI / 180.0);
    double dLong1InRad = Long1 * (Math.PI / 180.0);
    double dLat2InRad = Lat2 * (Math.PI / 180.0);
    double dLong2InRad = Long2 * (Math.PI / 180.0);

    double dLongitude = dLong2InRad - dLong1InRad;
    double dLatitude = dLat2InRad - dLat1InRad;

    // Intermediate result a.
    double a = Math.Pow(Math.Sin(dLatitude / 2.0), 2.0) +
               Math.Cos(dLat1InRad) * Math.Cos(dLat2InRad) *
               Math.Pow(Math.Sin(dLongitude / 2.0), 2.0);

    // Intermediate result c (great circle distance in Radians).
    double c = 2.0 * Math.Atan2(Math.Sqrt(a), Math.Sqrt(1.0 - a));

    // Distance.
    // const Double kEarthRadiusMiles = 3956.0;
    const Double kEarthRadiusKms = 6376.5;
    dDistance = kEarthRadiusKms * c;

    return dDistance;
}

If the distance between any vertex of the rectangle is less than the distance of the radius of the circle then the circle and rectangle overlap. If the distance between the center of the circle and all of the vertices is greater than the radius of the circle and all of those distances are shorter than the width and height of the rectangle then the circle should be inside of the rectangle.

Feel free to correct my code if you can find a problem with it as I'm sure there some condition that I have not thought of.

Also I'm not sure if this works for a rectangle that spans the ends of the hemispheres as the distance equation might break down.

public string Test(double cLat,
    double cLon,
    double cRadius,
    double rlat1,
    double rlon1,
    double rlat2,
    double rlon2,
    double rlat3,
    double rlon3,
    double rlat4,
    double rlon4)
{
    double d1 = Calc(cLat, cLon, rlat1, rlon1);
    double d2 = Calc(cLat, cLon, rlat2, rlon2);
    double d3 = Calc(cLat, cLon, rlat3, rlon3);
    double d4 = Calc(cLat, cLon, rlat4, rlon4);

    if (d1 <= cRadius ||
        d2 <= cRadius ||
        d3 <= cRadius ||
        d4 <= cRadius)
    {

        return "Circle and Rectangle intersect...";
    }

    double width = Calc(rlat1, rlon1, rlat2, rlon2);
    double height = Calc(rlat1, rlon1, rlat4, rlon4);

    if (d1 >= cRadius &&
        d2 >= cRadius &&
        d3 >= cRadius &&
        d4 >= cRadius &&
        width >= d1 &&
        width >= d2 &&
        width >= d3 &&
        width >= d4 &&
        height >= d1 &&
        height >= d2 &&
        height >= d3 &&
        height >= d4)
    {
        return "Circle is Inside of Rectangle!";
    }



    return "NO!";
}

Answer 7

Use the Stereographic projection. All circles (specifically latitudes, longitudes and your circle) map to circles (or lines) in the plane. Now it's just a question about circles and lines in plane geometry (even better, all the longitues are lines through 0, and all the latitudes are circles around 0)

Answer 8

• Yes, if the box corners contain the circle-center.
• Yes, if any of the box corners are within radius of circle-center.
• Yes, if the box contains the longitude of circle-center and the longitude intersection of the box-latitude closest to circle-center-latitude is within radius of circle-center.
• Yes, if the box contains the latitude of circle-center and the point at radius distance from circle-center on shortest-intersection-bearing is "beyond" the closest box-longitude; where shortest-intersection-bearing is determined by finding the initial bearing from circle-center to a point at latitude zero and a longitude that is pi/2 "beyond" the closest box-longitude.
• No, otherwise.

Assumptions:

• You can find the initial-bearing of a minimum course from point A to point B.
• You can find the distance between two points.

The first check is trivial. The second check just requires finding the four distances. The third check just requires finding the distance from circle-center to (closest-box-latitude, circle-center-longitude).

The fourth check requires finding the longitude line of the bounding box that is closest to the circle-center. Then find the center of the great circle on which that longitude line rests that is furthest from circle-center. Find the initial-bearing from circle-center to the great-circle-center. Find the point circle-radius from circle-center on that bearing. If that point is on the other side of the closest-longitude-line from circle-center, then the circle and bounding box intersect on that side.

It seems to me that there should be a flaw in this, but I haven't been able to find it.

The real problem that I can't seem to solve is to find the bounding-box that perfectly contains the circle (for circles that don't contain a pole). The bearing to the latitude min/max appears to be a function of the latitude of circle-center and circle-radius/(sphere circumference/4). Near the equator, it falls to pi/2 (east) or 3*pi/2 (west). As the center approaches the pole and the radius approaches sphere-circumference/4, the bearing approach zero (north) or pi (south).