Finding out if two areas on the globe intersect

• Introduction
• Examples of the problem
• The algorithm
• The proof of the algorithm: Forward
• The proof of the algorithm: the other direction
• A sample application: GIS database queries

Introduction

Testing two ``rectangular'' areas on the globe for intersection is markedly more complicated than it may appear. The challenge comes from longitude's wrapping around, that is, abruptly jumping from 180 degrees (i.e., 180E) to -180 degrees when we cross the International Date Line. Latitudes do not wrap around. Therefore, in the body of this article, we concentrate on the one-dimensional problem of checking two segments of a circle for intersection. The last section generalizes to georectangles and geographical information systems.

To identify points on a circle, we pick one point as the origin and refer to the others by their angular distance: ``longitudes''. The longitudes are measured in degrees, counter-clockwise from the point of origin. We will keep the longitudes within the range [-180,180), by adding or subtracting 360 if necessary.

A segment X of a circle is specified by its left end, Lx, and its right end, Rx. In other words, the segment [Lx, Rx] starts at the point Lx and stretches counter-clockwise to the point Rx. The coordinate -- longitude -- of the point Lx may be less, equal or greater than that of the point Rx. We assume that all segments are not degenerate (i.e., the endpoints are distinct) and shorter than the full circle. We call two segments intersecting if one segment contains an internal point of the other.

Examples of the problem

The following examples illustrate how non-obvious the intersection testing is.

    Segment A      Segment B     Intersect?
     [1, 10]        [15, 20]        NO
     [10,1]         [15, 20]        YES
     [5, 6]         [4, 7]          YES
     [6, 5]         [4, 7]        still YES
     [178, -178]    [179, -179]     YES
     [-178, 178]    [179, -179]     NO
     [178, -178]    [-170, -100]    NO
     [-178, 178]    [-170, -100]    YES

It is quite difficult to devise the single inequality involving La, Lb, Ra, and Rb that holds if and only if the segments intersect.

The algorithm

Suppose we are given two segments: segment A [La, Ra] and segment B [Lb, Rb]. We assume that the longitudes, La, Ra, Lb, and Rb, are the numbers within [-180,180).

Then

    diff(Ra-La) + diff(Lb-Ra) + diff(Rb-Lb) + diff(La-Rb) = 360

if and only if the segments A and B do not intersect.

Here, diff x is x mod 360. It normalizes x to be within [0,360). That is,

    diff x  = x       if x >= 0
              360+x   otherwise

The algorithm is that simple: a mere computation of four sets of differences, adding 360 to make all of them positive. If these four numbers add to 360, the segments in question do not intersect.

The proof of the algorithm: Forward

We will now prove that if the segments A and B do not intersect, the sum of the four diff terms above is equal to 360.

Suppose [La, Ra] and [Lb, Rb] do not intersect. If we start at the point La and walk counter-clockwise, we shall encounter the following points: La, Ra, Lb, Rb, and back to La. That means the segments [La, Ra], [Ra, Lb], [Lb, Rb] and [Rb, La] are all disjoint and together make up the full circle. The function diff x is designed so that diff (Rx-Lx) is the length of a segment [Lx, Rx]. The (angular) length is a number within [0,360). The sum of diff(Ra-La), diff(Lb-Ra), diff(Rb-Lb) and diff(La-Rb) is then the length of the full circle.

The proof of the algorithm: the other direction

In the other direction, if the sum of the four diff terms is equal to 360, then the segments do not intersect. We will prove it by contradiction.

Let us assume that the sum is equal to 360 and the segments intersect. There are three major cases to consider: (i) segment B is entirely within segment A (or the other way around); (ii) Segment A contains Lb but not Rb; (iii) Segment A contains both endpoints Lb and Rb, but not the whole [Lb, Rb].

Case 1: Segment B is entirely within Segment A. If we walk from La counter-clockwise we will pass the points Lb, Rb, Ra, and come back to La. Noting that diff(-x) + diff x = 360 for x > 0, we write:

    diff (Lb-Ra) 
       = 360 - diff(Ra-Lb) 
       = 360 - diff(Rb-Lb) - diff(Ra-Rb)
    diff (La-Rb)
       = 360 - diff(Rb-La)
       = 360 - diff(Rb-Lb) - diff(Lb-La)
    diff (Ra-La) 
       = diff(Lb-La) + diff(Rb-Lb) + diff(Ra-Rb)

Thus we have

    diff(Ra-La) + diff(Lb-Ra) + diff(Rb-Lb) + diff(La-Rb) = 720

which is contrary to our assumption that this sum is 360.

Case 2: Segment A contains Lb but not Rb. This implies that Rb must be distinct from La and Lb must be distinct from Ra. The same walk as in Case 1 gives us the following sequence of points: La, Lb, Ra, Rb, La. Therefore

    diff (Ra-La)
       = diff(Lb-La) + diff(Ra-Lb)
    diff (Lb-Ra) 
       = 360 - diff(Ra-Lb)
    diff (Rb-Lb)
       = diff(Ra-Lb) + diff(Rb-Ra)
    diff (La-Rb) 
       = 360 - diff(Rb-La)
       = 360 - diff(Lb-La) - diff(Ra-Lb) - diff(Rb-Ra)

which means that the sum

    diff(Ra-La) + diff(Lb-Ra) + diff(Rb-Lb) + diff(La-Rb) = 720

and not 360 as we assumed initially.

Case 3: Segment A contains both points Lb and Rb, but not all the points in between; for example: [6, 5] and [4, 7]. The walk counter-clockwise starting from La goes through the points La, Rb, Lb, Ra, and back to La. Consider the equality

    diff(Rb-La) + diff(Lb-Rb) + diff(Ra-Lb) + diff(La-Ra) = 360

whose right-hand side looks almost the same as the sum we want to evaluate:

    diff(Ra-La) + diff(Lb-Ra) + diff(Rb-Lb) + diff(La-Rb)

However, all the differences are in ``reverse'': diff (Ra-La) rather than diff (La-Ra), etc. Keeping in mind that diff x = 360 - diff(-x) for x > 0, we obtain

    diff(Ra-La) + diff(Lb-Ra) + diff(Rb-Lb) + diff(La-Rb) = 1080

if all four end points are distinct. If La and Rb are the same, the sum will be equal to 720. Ditto if Lb and Ra are the same. In any event, the result contradicts our initial assumption. Since we assumed that segments A and B intersect, one segment must contain an internal point of the other. This fact implies that La and Rb cannot be the same when Lb and Ra are the same.

QED.

A sample application: GIS database queries

The intersection algorithm has been used in a production geographical information system. One of the databases stores satellite imagery. A satellite image is characterized, among other attributes, by its georectangular bounding box. The user can select an area of interest (again, as a georectangle) and find out relevant satellite images, i.e., the images whose bounding box intersects the area of interest. If the area of interest is a circle or a complex polygon, we can inscribe it into a georectangle and use the algorithm above. The test may in this case yield a false positive, but it will never return a false negative. Therefore, the georectangle intersection approach can be used as the first, winnowing stage in checking for intersection of complex areas on the globe.

The following is a query template excerpted from the production code. The query returns the bounding boxes and the descriptions of all satellite images that intersect the user-defined georectangle. The latter is specified by its two corner-points: (NLAT, WLON) and (SLAT,ELON).

    SELECT lat_n, lat_s, lon_w, lon_e, descr
    FROM ImageDescr
    WHERE
      (SELECT image_descr FROM Timestamps) >= '2002-06-07 16:47'
      AND lat_s <= 'NLAT' AND lat_n >= 'SLAT'
      AND q_lon_intersect(lon_w,lon_e,WLON,ELON) <> 360
    ORDER BY lat_n, lat_s, lon_w, lon_e
    FOR READ ONLY;

The query relies on the following stored procedure q_lon_intersect(), which implements the segment intersection algorithm verbatim:

    -- This q_lon_intersect(lon_w,lon_e,Query_WLON,Query_ELON)
    -- procedure determines if two circle segments [lon_w, lon_e] and
    -- [Query_WLON, Query_ELON]  intersect. The procedure works correctly
    -- for lon_w < lon_e and lon_w >= lon_e (and similarly for Query_WLON
    -- and Query_ELON).
    -- The two circle segments intersect when the result of this procedure
    -- is NOT equal to 360.
    
    CREATE PROCEDURE q_lon_intersect
            (La SMALLFLOAT, Ra SMALLFLOAT, Lb SMALLFLOAT, Rb SMALLFLOAT)
            returning SMALLFLOAT;
    DEFINE d1, d2, d3, d4 SMALLFLOAT;
    LET d1, d2, d3, d4 = Ra-La, Lb-Ra, Rb-Lb, La-Rb;
    RETURN  CASE WHEN d1 < 0 THEN 360 + d1 ELSE d1 END +
            CASE WHEN d2 < 0 THEN 360 + d2 ELSE d2 END +
            CASE WHEN d3 < 0 THEN 360 + d3 ELSE d3 END +
            CASE WHEN d4 < 0 THEN 360 + d4 ELSE d4 END;
    END PROCEDURE
    DOCUMENT 'Compute a SMALLFLOAT that can tell if two circle segments',
             'intersect. If they do, the result of this procedure is',
             'NOT equal to 360.'
    WITH LISTING IN '/tmp/sql.complile.log';