Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

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

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.

like image 652
Nick Johnson Avatar asked Dec 26 '08 19:12

Nick Johnson


4 Answers

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.

like image 135
Jason S Avatar answered Oct 01 '22 02:10

Jason S


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)

like image 35
David Lehavi Avatar answered Oct 01 '22 03:10

David Lehavi


  • 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).

like image 28
user103811 Avatar answered Oct 01 '22 03:10

user103811


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.

like image 24
aib Avatar answered Oct 01 '22 01:10

aib