Suppose I have the following:
How can I determine:
I'm looking for a complete formula/algorithm, rather than a lesson in the math, per-se.
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:
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()
:
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.
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)
Assumptions:
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).
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.
If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With