LUCENE-7147: Improve disjoint check for geo distance query traversal

lucene/CHANGES.txt

@@ -37,6 +37,9 @@ Optimizations
 * LUCENE-7115: Speed up FieldCache.CacheEntry toString by setting initial
   StringBuilder capacity (Gregory Chanan)
 
+* LUCENE-7147: Improve disjoint check for geo distance query traversal
+  (Ryan Ernst, Robert Muir, Mike McCandless)
+
 Bug Fixes
 
 * LUCENE-7127: Fix corner case bugs in GeoPointDistanceQuery. (Robert Muir)

lucene/sandbox/src/java/org/apache/lucene/document/LatLonPointDistanceQuery.java

@@ -108,6 +108,8 @@ final class LatLonPointDistanceQuery extends Query {
       maxPartialDistance = Double.POSITIVE_INFINITY;
     }
 
+    final double axisLat = GeoUtils.axisLat(latitude, radiusMeters);
+
     return new ConstantScoreWeight(this) {
 
       @Override
@@ -171,8 +173,9 @@ final class LatLonPointDistanceQuery extends Query {
                            
                            // algorithm: we create a bounding box (two bounding boxes if we cross the dateline).
                            // 1. check our bounding box(es) first. if the subtree is entirely outside of those, bail.
-                           // 2. see if the subtree is fully contained. if the subtree is enormous along the x axis, wrapping half way around the world, etc: then this can't work, just go to step 3.
+                           // 2. check if the subtree is disjoint. it may cross the bounding box but not intersect with circle
-                           // 3. recurse naively.
+                           // 3. see if the subtree is fully contained. if the subtree is enormous along the x axis, wrapping half way around the world, etc: then this can't work, just go to step 3.
+                           // 4. recurse naively (subtrees crossing over circle edge)
                            @Override
                            public Relation compare(byte[] minPackedValue, byte[] maxPackedValue) {
                              if (StringHelper.compare(Integer.BYTES, minPackedValue, 0, maxLat, 0) > 0 ||
@@ -193,6 +196,17 @@ final class LatLonPointDistanceQuery extends Query {
                              double latMax = LatLonPoint.decodeLatitude(maxPackedValue, 0);
                              double lonMax = LatLonPoint.decodeLongitude(maxPackedValue, Integer.BYTES);
 
+                             if ((longitude < lonMin || longitude > lonMax) && (axisLat+GeoUtils.AXISLAT_ERROR < latMin || axisLat-GeoUtils.AXISLAT_ERROR > latMax)) {
+                               // circle not fully inside / crossing axis
+                               if (SloppyMath.haversinMeters(latitude, longitude, latMin, lonMin) > radiusMeters &&
+                                   SloppyMath.haversinMeters(latitude, longitude, latMin, lonMax) > radiusMeters &&
+                                   SloppyMath.haversinMeters(latitude, longitude, latMax, lonMin) > radiusMeters &&
+                                   SloppyMath.haversinMeters(latitude, longitude, latMax, lonMax) > radiusMeters) {
+                                 // no points inside
+                                 return Relation.CELL_OUTSIDE_QUERY;
+                               }
+                             }
+
                              if (lonMax - longitude < 90 && longitude - lonMin < 90 &&
                                  SloppyMath.haversinMeters(latitude, longitude, latMin, lonMin) <= radiusMeters &&
                                  SloppyMath.haversinMeters(latitude, longitude, latMin, lonMax) <= radiusMeters &&

lucene/spatial/src/java/org/apache/lucene/spatial/geopoint/search/GeoPointDistanceQueryImpl.java

@@ -19,6 +19,7 @@ package org.apache.lucene.spatial.geopoint.search;
 import org.apache.lucene.search.MultiTermQuery;
 import org.apache.lucene.spatial.geopoint.document.GeoPointField.TermEncoding;
 import org.apache.lucene.spatial.util.GeoRect;
+import org.apache.lucene.spatial.util.GeoUtils;
 import org.apache.lucene.util.SloppyMath;
 
 /** Package private implementation for the public facing GeoPointDistanceQuery delegate class.
@@ -32,6 +33,9 @@ final class GeoPointDistanceQueryImpl extends GeoPointInBBoxQueryImpl {
   // optimization, maximum partial haversin needed to be a candidate
   private final double maxPartialDistance;
   
+  // optimization, used for detecting axis cross
+  final double axisLat;
+  
   GeoPointDistanceQueryImpl(final String field, final TermEncoding termEncoding, final GeoPointDistanceQuery q,
                             final double centerLonUnwrapped, final GeoRect bbox) {
     super(field, termEncoding, bbox.minLat, bbox.maxLat, bbox.minLon, bbox.maxLon);
@@ -46,6 +50,7 @@ final class GeoPointDistanceQueryImpl extends GeoPointInBBoxQueryImpl {
     } else {
       maxPartialDistance = Double.POSITIVE_INFINITY;
     }
+    axisLat = GeoUtils.axisLat(distanceQuery.centerLat, distanceQuery.radiusMeters);
   }
 
   @Override
@@ -65,15 +70,22 @@ final class GeoPointDistanceQueryImpl extends GeoPointInBBoxQueryImpl {
 
     @Override
     protected boolean cellCrosses(final double minLat, final double maxLat, final double minLon, final double maxLon) {
+      // bounding box check
       if (maxLat < GeoPointDistanceQueryImpl.this.minLat ||
           maxLon < GeoPointDistanceQueryImpl.this.minLon ||
           minLat > GeoPointDistanceQueryImpl.this.maxLat ||
           minLon > GeoPointDistanceQueryImpl.this.maxLon) {
         return false;
-      } else {
+      } else if ((centerLon < minLon || centerLon > maxLon) && (axisLat+GeoUtils.AXISLAT_ERROR < minLat || axisLat-GeoUtils.AXISLAT_ERROR > maxLat)) {
-        return true;
+        if (SloppyMath.haversinMeters(distanceQuery.centerLat, centerLon, minLat, minLon) > distanceQuery.radiusMeters &&
+            SloppyMath.haversinMeters(distanceQuery.centerLat, centerLon, minLat, maxLon) > distanceQuery.radiusMeters &&
+            SloppyMath.haversinMeters(distanceQuery.centerLat, centerLon, maxLat, minLon) > distanceQuery.radiusMeters &&
+            SloppyMath.haversinMeters(distanceQuery.centerLat, centerLon, maxLat, maxLon) > distanceQuery.radiusMeters) {
+          return false;
         }
       }
+      return true;
+    }
 
     @Override
     protected boolean cellWithin(final double minLat, final double maxLat, final double minLon, final double maxLon) {

lucene/spatial/src/java/org/apache/lucene/spatial/util/GeoUtils.java

@@ -109,7 +109,8 @@ public final class GeoUtils {
   public static GeoRect circleToBBox(final double centerLat, final double centerLon, final double radiusMeters) {
     final double radLat = TO_RADIANS * centerLat;
     final double radLon = TO_RADIANS * centerLon;
-    double radDistance = radiusMeters / SEMIMAJOR_AXIS;
+    // LUCENE-7143
+    double radDistance = (radiusMeters + 7E-2) / SEMIMAJOR_AXIS;
     double minLat = radLat - radDistance;
     double maxLat = radLat + radDistance;
     double minLon;
@@ -176,4 +177,51 @@ public final class GeoUtils {
     return cos(a - PIO2);
   }
 
+  /** maximum error from {@link #axisLat(double, double)}. logic must be prepared to handle this */
+  public static final double AXISLAT_ERROR = 0.1D / SEMIMAJOR_AXIS * TO_DEGREES;
+
+  /**
+   * Calculate the latitude of a circle's intersections with its bbox meridians.
+   * <p>
+   * <b>NOTE:</b> the returned value will be +/- {@link #AXISLAT_ERROR} of the actual value.
+   * @param centerLat The latitude of the circle center
+   * @param radiusMeters The radius of the circle in meters
+   * @return A latitude
+   */
+  public static double axisLat(double centerLat, double radiusMeters) {
+    // A spherical triangle with:
+    // r is the radius of the circle in radians
+    // l1 is the latitude of the circle center
+    // l2 is the latitude of the point at which the circle intersect's its bbox longitudes
+    // We know r is tangent to the bbox meridians at l2, therefore it is a right angle.
+    // So from the law of cosines, with the angle of l1 being 90, we have:
+    // cos(l1) = cos(r) * cos(l2) + sin(r) * sin(l2) * cos(90)
+    // The second part cancels out because cos(90) == 0, so we have:
+    // cos(l1) = cos(r) * cos(l2)
+    // Solving for l2, we get:
+    // l2 = acos( cos(l1) / cos(r) )
+    // We ensure r is in the range (0, PI/2) and l1 in the range (0, PI/2]. This means we
+    // cannot divide by 0, and we will always get a positive value in the range [0, 1) as
+    // the argument to arc cosine, resulting in a range (0, PI/2].
+
+    double l1 = TO_RADIANS * centerLat;
+    double r = (radiusMeters + 7E-2) / SEMIMAJOR_AXIS;
+
+    // if we are within radius range of a pole, the lat is the pole itself
+    if (Math.abs(l1) + r >= MAX_LAT_RADIANS) {
+      return centerLat >= 0 ? MAX_LAT_INCL : MIN_LAT_INCL;
+    }
+
+    // adjust l1 as distance from closest pole, to form a right triangle with bbox meridians
+    // and ensure it is in the range (0, PI/2]
+    l1 = centerLat >= 0 ? PIO2 - l1 : l1 + PIO2;
+
+    double l2 = Math.acos(Math.cos(l1) / Math.cos(r));
+    assert !Double.isNaN(l2);
+
+    // now adjust back to range [-pi/2, pi/2], ie latitude in radians
+    l2 = centerLat >= 0 ? PIO2 - l2 : l2 - PIO2;
+
+    return TO_DEGREES * l2;
+  }
 }

lucene/spatial/src/test/org/apache/lucene/spatial/util/GeoTestUtil.java

@@ -17,8 +17,11 @@
 package org.apache.lucene.spatial.util;
 
 import java.util.ArrayList;
+import java.util.List;
 import java.util.Random;
 
+import org.apache.lucene.util.SloppyMath;
+
 import com.carrotsearch.randomizedtesting.RandomizedContext;
 
 /** static methods for testing geo */
@@ -266,4 +269,232 @@ public class GeoTestUtil {
   private static Random random() {
    return RandomizedContext.current().getRandom();
   }
+  
+  // craziness for plotting stuff :)
+  
+  private static double wrapLat(double lat) {
+    //System.out.println("wrapLat " + lat);
+    if (lat > 90) {
+      //System.out.println("  " + (180 - lat));
+      return 180 - lat;
+    } else if (lat < -90) {
+      //System.out.println("  " + (-180 - lat));
+      return -180 - lat;
+    } else {
+      //System.out.println("  " + lat);
+      return lat;
+    }
+  }
+
+  private static double wrapLon(double lon) {
+    //System.out.println("wrapLon " + lon);
+    if (lon > 180) {
+      //System.out.println("  " + (lon - 360));
+      return lon - 360;
+    } else if (lon < -180) {
+      //System.out.println("  " + (lon + 360));
+      return lon + 360;
+    } else {
+      //System.out.println("  " + lon);
+      return lon;
+    }
+  }
+  
+  private static void drawRectApproximatelyOnEarthSurface(String name, String color, double minLat, double maxLat, double minLon, double maxLon) {
+    int steps = 20;
+    System.out.println("        var " + name + " = WE.polygon([");
+    System.out.println("          // min -> max lat, min lon");
+    for(int i=0;i<steps;i++) {
+      System.out.println("          [" + (minLat + (maxLat - minLat) * i / steps) + ", " + minLon + "],");
+    }
+    System.out.println("          // max lat, min -> max lon");
+    for(int i=0;i<steps;i++) {
+      System.out.println("          [" + (maxLat + ", " + (minLon + (maxLon - minLon) * i / steps)) + "],");
+    }
+    System.out.println("          // max -> min lat, max lon");
+    for(int i=0;i<steps;i++) {
+      System.out.println("          [" + (minLat + (maxLat - minLat) * (steps-i) / steps) + ", " + maxLon + "],");
+    }
+    System.out.println("          // min lat, max -> min lon");
+    for(int i=0;i<steps;i++) {
+      System.out.println("          [" + minLat + ", " + (minLon + (maxLon - minLon) * (steps-i) / steps) + "],");
+    }
+    System.out.println("          // min lat, min lon");
+    System.out.println("          [" + minLat + ", " + minLon + "]");
+    System.out.println("        ], {color: \"" + color + "\", fillColor: \"" + color + "\"});");
+    System.out.println("        " + name + ".addTo(earth);");
+  }
+  
+  private static void plotLatApproximatelyOnEarthSurface(String name, String color, double lat, double minLon, double maxLon) {
+    System.out.println("        var " + name + " = WE.polygon([");
+    double lon;
+    for(lon = minLon;lon<=maxLon;lon += (maxLon-minLon)/36) {
+      System.out.println("          [" + lat + ", " + lon + "],");
+    }
+    System.out.println("          [" + lat + ", " + maxLon + "],");
+    lon -= (maxLon-minLon)/36;
+    for(;lon>=minLon;lon -= (maxLon-minLon)/36) {
+      System.out.println("          [" + lat + ", " + lon + "],");
+    }
+    System.out.println("        ], {color: \"" + color + "\", fillColor: \"#ffffff\", opacity: " + (color.equals("#ffffff") ? "0.3" : "1") + ", fillOpacity: 0.0001});");
+    System.out.println("        " + name + ".addTo(earth);");
+  }
+
+  private static void plotLonApproximatelyOnEarthSurface(String name, String color, double lon, double minLat, double maxLat) {
+    System.out.println("        var " + name + " = WE.polygon([");
+    double lat;
+    for(lat = minLat;lat<=maxLat;lat += (maxLat-minLat)/36) {
+      System.out.println("          [" + lat + ", " + lon + "],");
+    }
+    System.out.println("          [" + maxLat + ", " + lon + "],");
+    lat -= (maxLat-minLat)/36;
+    for(;lat>=minLat;lat -= (maxLat-minLat)/36) {
+      System.out.println("          [" + lat + ", " + lon + "],");
+    }
+    System.out.println("        ], {color: \"" + color + "\", fillColor: \"#ffffff\", opacity: " + (color.equals("#ffffff") ? "0.3" : "1") + ", fillOpacity: 0.0001});");
+    System.out.println("        " + name + ".addTo(earth);");
+  }
+
+  // http://www.webglearth.org has API details:
+  public static void polysToWebGLEarth(List<double[][]> polys) {
+    System.out.println("<!DOCTYPE HTML>");
+    System.out.println("<html>");
+    System.out.println("  <head>");
+    System.out.println("    <script src=\"http://www.webglearth.com/v2/api.js\"></script>");
+    System.out.println("    <script>");
+    System.out.println("      function initialize() {");
+    System.out.println("        var earth = new WE.map('earth_div');");
+
+    int count = 0;
+    for (double[][] poly : polys) {
+      System.out.println("        var poly" + count + " = WE.polygon([");
+      for(int i=0;i<poly[0].length;i++) {
+        double lat = poly[0][i];
+        double lon = poly[1][i];
+        System.out.println("          [" + lat + ", " + lon + "],");
+      }
+      System.out.println("        ], {color: '#00ff00'});");    
+      System.out.println("        poly" + count + ".addTo(earth);");
+    }
+
+    System.out.println("        WE.tileLayer('http://{s}.tile.openstreetmap.org/{z}/{x}/{y}.png',{");
+    System.out.println("          attribution: '© OpenStreetMap contributors'");
+    System.out.println("        }).addTo(earth);");
+    System.out.println("      }");
+    System.out.println("    </script>");
+    System.out.println("    <style>");
+    System.out.println("      html, body{padding: 0; margin: 0;}");
+    System.out.println("      #earth_div{top: 0; right: 0; bottom: 0; left: 0; position: absolute !important;}");
+    System.out.println("    </style>");
+    System.out.println("    <title>WebGL Earth API: Hello World</title>");
+    System.out.println("  </head>");
+    System.out.println("  <body onload=\"initialize()\">");
+    System.out.println("    <div id=\"earth_div\"></div>");
+    System.out.println("  </body>");
+    System.out.println("</html>");
+  }
+
+  // http://www.webglearth.org has API details:
+  public static void toWebGLEarth(double rectMinLatitude, double rectMaxLatitude,
+                                   double rectMinLongitude, double rectMaxLongitude,
+                                   double centerLatitude, double centerLongitude,
+                                   double radiusMeters) {
+    GeoRect box = GeoUtils.circleToBBox(centerLatitude, centerLongitude, radiusMeters);
+    System.out.println("<!DOCTYPE HTML>");
+    System.out.println("<html>");
+    System.out.println("  <head>");
+    System.out.println("    <script src=\"http://www.webglearth.com/v2/api.js\"></script>");
+    System.out.println("    <script>");
+    System.out.println("      function initialize() {");
+    System.out.println("        var earth = new WE.map('earth_div', {center: [" + centerLatitude + ", " + centerLongitude + "]});");
+    System.out.println("        var marker = WE.marker([" + centerLatitude + ", " + centerLongitude + "]).addTo(earth);");
+    drawRectApproximatelyOnEarthSurface("cell", "#ff0000", rectMinLatitude, rectMaxLatitude, rectMinLongitude, rectMaxLongitude);
+    System.out.println("        var polygonB = WE.polygon([");
+    StringBuilder b = new StringBuilder();
+    inverseHaversin(b, centerLatitude, centerLongitude, radiusMeters);
+    System.out.println(b);
+    System.out.println("        ], {color: '#00ff00'});");    
+    System.out.println("        polygonB.addTo(earth);");
+    drawRectApproximatelyOnEarthSurface("bbox", "#00ff00", box.minLat, box.maxLat, box.minLon, box.maxLon);
+    System.out.println("        WE.tileLayer('http://{s}.tile.openstreetmap.org/{z}/{x}/{y}.png',{");
+    System.out.println("          attribution: '© OpenStreetMap contributors'");
+    System.out.println("        }).addTo(earth);");
+    plotLatApproximatelyOnEarthSurface("lat0", "#ffffff", 4.68, 0.0, 360.0);
+    plotLatApproximatelyOnEarthSurface("lat1", "#ffffff", 180-93.09, 0.0, 360.0);
+    plotLatApproximatelyOnEarthSurface("axisLat", "#00ff00", GeoUtils.axisLat(centerLatitude, radiusMeters), box.minLon, box.maxLon);
+    plotLonApproximatelyOnEarthSurface("axisLon", "#00ff00", centerLongitude, box.minLat, box.maxLat);
+    System.out.println("      }");
+    System.out.println("    </script>");
+    System.out.println("    <style>");
+    System.out.println("      html, body{padding: 0; margin: 0;}");
+    System.out.println("      #earth_div{top: 0; right: 0; bottom: 0; left: 0; position: absolute !important;}");
+    System.out.println("    </style>");
+    System.out.println("    <title>WebGL Earth API: Hello World</title>");
+    System.out.println("  </head>");
+    System.out.println("  <body onload=\"initialize()\">");
+    System.out.println("    <div id=\"earth_div\"></div>");
+    System.out.println("  </body>");
+    System.out.println("</html>");
+  }
+
+  private static void inverseHaversin(StringBuilder b, double centerLat, double centerLon, double radiusMeters) {
+    double angle = 0;
+    int steps = 100;
+
+    newAngle:
+    while (angle < 360) {
+      double x = Math.cos(Math.toRadians(angle));
+      double y = Math.sin(Math.toRadians(angle));
+      double factor = 2.0;
+      double step = 1.0;
+      int last = 0;
+      double lastDistanceMeters = 0.0;
+      //System.out.println("angle " + angle + " slope=" + slope);
+      while (true) {
+        double lat = wrapLat(centerLat + y * factor);
+        double lon = wrapLon(centerLon + x * factor);
+        double distanceMeters = SloppyMath.haversinMeters(centerLat, centerLon, lat, lon);
+
+        if (last == 1 && distanceMeters < lastDistanceMeters) {
+          // For large enough circles, some angles are not possible:
+          //System.out.println("  done: give up on angle " + angle);
+          angle += 360./steps;
+          continue newAngle;
+        }
+        if (last == -1 && distanceMeters > lastDistanceMeters) {
+          // For large enough circles, some angles are not possible:
+          //System.out.println("  done: give up on angle " + angle);
+          angle += 360./steps;
+          continue newAngle;
+        }
+        lastDistanceMeters = distanceMeters;
+
+        //System.out.println("  iter lat=" + lat + " lon=" + lon + " distance=" + distanceMeters + " vs " + radiusMeters);
+        if (Math.abs(distanceMeters - radiusMeters) < 0.1) {
+          b.append("          [" + lat + ", " + lon + "],\n");
+          break;
+        }
+        if (distanceMeters > radiusMeters) {
+          // too big
+          //System.out.println("    smaller");
+          factor -= step;
+          if (last == 1) {
+            //System.out.println("      half-step");
+            step /= 2.0;
+          }
+          last = -1;
+        } else if (distanceMeters < radiusMeters) {
+          // too small
+          //System.out.println("    bigger");
+          factor += step;
+          if (last == -1) {
+            //System.out.println("      half-step");
+            step /= 2.0;
+          }
+          last = 1;
+        }
+      }
+      angle += 360./steps;
+    }
+  }
 }

lucene/spatial/src/test/org/apache/lucene/spatial/util/TestGeoUtils.java

@@ -16,6 +16,8 @@
  */
 package org.apache.lucene.spatial.util;
 
+import java.util.Locale;
+
 import org.apache.lucene.util.BytesRefBuilder;
 import org.apache.lucene.util.LuceneTestCase;
 import org.apache.lucene.util.SloppyMath;
@@ -320,4 +322,136 @@ public class TestGeoUtils extends LuceneTestCase {
       }
     }
   }
+
+  public void testAxisLat() {
+    double earthCircumference = 2D * Math.PI * GeoUtils.SEMIMAJOR_AXIS;
+    assertEquals(90, GeoUtils.axisLat(0, earthCircumference / 4), 0.0D);
+
+    for (int i = 0; i < 100; ++i) {
+      boolean reallyBig = random().nextInt(10) == 0;
+      final double maxRadius = reallyBig ? 1.1 * earthCircumference : earthCircumference / 8;
+      final double radius = maxRadius * random().nextDouble();
+      double prevAxisLat = GeoUtils.axisLat(0.0D, radius);
+      for (double lat = 0.1D; lat < 90D; lat += 0.1D) {
+        double nextAxisLat = GeoUtils.axisLat(lat, radius);
+        GeoRect bbox = GeoUtils.circleToBBox(lat, 180D, radius);
+        double dist = SloppyMath.haversinMeters(lat, 180D, nextAxisLat, bbox.maxLon);
+        if (nextAxisLat < GeoUtils.MAX_LAT_INCL) {
+          assertEquals("lat = " + lat, dist, radius, 0.1D);
+        }
+        assertTrue("lat = " + lat, prevAxisLat <= nextAxisLat);
+        prevAxisLat = nextAxisLat;
+      }
+
+      prevAxisLat = GeoUtils.axisLat(-0.0D, radius);
+      for (double lat = -0.1D; lat > -90D; lat -= 0.1D) {
+        double nextAxisLat = GeoUtils.axisLat(lat, radius);
+        GeoRect bbox = GeoUtils.circleToBBox(lat, 180D, radius);
+        double dist = SloppyMath.haversinMeters(lat, 180D, nextAxisLat, bbox.maxLon);
+        if (nextAxisLat > GeoUtils.MIN_LAT_INCL) {
+          assertEquals("lat = " + lat, dist, radius, 0.1D);
+        }
+        assertTrue("lat = " + lat, prevAxisLat >= nextAxisLat);
+        prevAxisLat = nextAxisLat;
+      }
+    }
+  }
+  
+  // TODO: does not really belong here, but we test it like this for now
+  // we can make a fake IndexReader to send boxes directly to Point visitors instead?
+  public void testCircleOpto() throws Exception {
+    for (int i = 0; i < 50; i++) {
+      // circle
+      final double centerLat = -90 + 180.0 * random().nextDouble();
+      final double centerLon = -180 + 360.0 * random().nextDouble();
+      final double radius = 50_000_000D * random().nextDouble();
+      final GeoRect box = GeoUtils.circleToBBox(centerLat, centerLon, radius);
+      // TODO: remove this leniency!
+      if (box.crossesDateline()) {
+        --i; // try again...
+        continue;
+      }
+      final double axisLat = GeoUtils.axisLat(centerLat, radius);
+
+      for (int k = 0; k < 1000; ++k) {
+
+        double[] latBounds = {-90, box.minLat, axisLat, box.maxLat, 90};
+        double[] lonBounds = {-180, box.minLon, centerLon, box.maxLon, 180};
+        // first choose an upper left corner
+        int maxLatRow = random().nextInt(4);
+        double latMax = randomInRange(latBounds[maxLatRow], latBounds[maxLatRow + 1]);
+        int minLonCol = random().nextInt(4);
+        double lonMin = randomInRange(lonBounds[minLonCol], lonBounds[minLonCol + 1]);
+        // now choose a lower right corner
+        int minLatMaxRow = maxLatRow == 3 ? 3 : maxLatRow + 1; // make sure it will at least cross into the bbox
+        int minLatRow = random().nextInt(minLatMaxRow);
+        double latMin = randomInRange(latBounds[minLatRow], Math.min(latBounds[minLatRow + 1], latMax));
+        int maxLonMinCol = Math.max(minLonCol, 1); // make sure it will at least cross into the bbox
+        int maxLonCol = maxLonMinCol + random().nextInt(4 - maxLonMinCol);
+        double lonMax = randomInRange(Math.max(lonBounds[maxLonCol], lonMin), lonBounds[maxLonCol + 1]);
+
+        assert latMax >= latMin;
+        assert lonMax >= lonMin;
+
+        if (isDisjoint(centerLat, centerLon, radius, axisLat, latMin, latMax, lonMin, lonMax)) {
+          // intersects says false: test a ton of points
+          for (int j = 0; j < 200; j++) {
+            double lat = latMin + (latMax - latMin) * random().nextDouble();
+            double lon = lonMin + (lonMax - lonMin) * random().nextDouble();
+
+            if (random().nextBoolean()) {
+              // explicitly test an edge
+              int edge = random().nextInt(4);
+              if (edge == 0) {
+                lat = latMin;
+              } else if (edge == 1) {
+                lat = latMax;
+              } else if (edge == 2) {
+                lon = lonMin;
+              } else if (edge == 3) {
+                lon = lonMax;
+              }
+            }
+            double distance = SloppyMath.haversinMeters(centerLat, centerLon, lat, lon);
+            try {
+            assertTrue(String.format(Locale.ROOT, "\nisDisjoint(\n" +
+                    "centerLat=%s\n" +
+                    "centerLon=%s\n" +
+                    "radius=%s\n" +
+                    "latMin=%s\n" +
+                    "latMax=%s\n" +
+                    "lonMin=%s\n" +
+                    "lonMax=%s) == false BUT\n" +
+                    "haversin(%s, %s, %s, %s) = %s\nbbox=%s",
+                centerLat, centerLon, radius, latMin, latMax, lonMin, lonMax,
+                centerLat, centerLon, lat, lon, distance, GeoUtils.circleToBBox(centerLat, centerLon, radius)),
+                distance > radius);
+            } catch (AssertionError e) {
+              GeoTestUtil.toWebGLEarth(latMin, latMax, lonMin, lonMax, centerLat, centerLon, radius);
+              throw e;
+            }
+          }
+        }
+      }
+    }
+  }
+
+  static double randomInRange(double min, double max) {
+    return min + (max - min) * random().nextDouble();
+  }
+  
+  static boolean isDisjoint(double centerLat, double centerLon, double radius, double axisLat, double latMin, double latMax, double lonMin, double lonMax) {
+    if ((centerLon < lonMin || centerLon > lonMax) && (axisLat+GeoUtils.AXISLAT_ERROR < latMin || axisLat-GeoUtils.AXISLAT_ERROR > latMax)) {
+      // circle not fully inside / crossing axis
+      if (SloppyMath.haversinMeters(centerLat, centerLon, latMin, lonMin) > radius &&
+          SloppyMath.haversinMeters(centerLat, centerLon, latMin, lonMax) > radius &&
+          SloppyMath.haversinMeters(centerLat, centerLon, latMax, lonMin) > radius &&
+          SloppyMath.haversinMeters(centerLat, centerLon, latMax, lonMax) > radius) {
+        // no points inside
+        return true;
+      }
+    }
+    
+    return false;
+  }
 }