fix(spatial): cover the antimeridian and the poles in Near queries (#166)

src/BLite.Core/Indexing/RTreeIndex.cs

@@ -57,16 +57,22 @@ public IEnumerable<DocumentLocation> Search(GeoBox area, ITransaction? transacti
 
     /// <summary>
     /// Locations of the indexed points within <paramref name="radiusKm"/> of the center. The tree
-    /// is searched with the circle's bounding box, then each hit is checked against the
-    /// great-circle distance, so points in the box corners are left out.
+    /// is searched with boxes covering the circle (split at the ±180° meridian, widened to all
+    /// longitudes around a pole), then each hit is checked against the great-circle distance,
+    /// so points in the box corners are left out.
     /// </summary>
     public IEnumerable<DocumentLocation> SearchNear(double latitude, double longitude, double radiusKm, ITransaction? transaction = null)
     {
-        var box = SpatialMath.BoundingBox(latitude, longitude, radiusKm);
-        foreach (var (mbr, location) in SearchEntries(box, transaction))
+        // The boxes never overlap, so no point is returned twice: a split happens only when no
+        // pole is inside the circle, and then the longitude half-width asin(sin r / cos lat) is
+        // below 90°, so the two halves together span less than 360°.
+        foreach (var box in SpatialMath.NearQueryBoxes(latitude, longitude, radiusKm))
         {
-            if (SpatialMath.DistanceKm(mbr.MinLat, mbr.MinLon, latitude, longitude) <= radiusKm)
-                yield return location;
+            foreach (var (mbr, location) in SearchEntries(box, transaction))
+            {
+                if (SpatialMath.DistanceKm(mbr.MinLat, mbr.MinLon, latitude, longitude) <= radiusKm)
+                    yield return location;
+            }
         }
     }
 

src/BLite.Core/Indexing/SpatialMath.cs

@@ -26,22 +26,45 @@ public static double DistanceKm(double lat1, double lon1, double lat2, double lo
     }
 
     /// <summary>
-    /// Creates a Bounding Box (MBR) centered at a point with a given radius in km.
+    /// Boxes that together cover every point within <paramref name="radiusKm"/> of the center, for
+    /// an R-Tree search whose hits are then checked against the exact distance.
+    /// <list type="bullet">
+    ///   <item>A circle that contains a pole reaches every longitude, so its box spans all
+    ///   longitudes from the circle's far edge up to the pole.</item>
+    ///   <item>Otherwise the longitude half-width is the circle's widest point,
+    ///   asin(sin r / cos lat), which lies poleward of the center; a box that crosses the ±180°
+    ///   meridian is split into one box on each side.</item>
+    /// </list>
     /// </summary>
-    internal static GeoBox BoundingBox(GeoPoint center, double radiusKm) => BoundingBox(center.Latitude, center.Longitude, radiusKm);
+    internal static GeoBox[] NearQueryBoxes(double lat, double lon, double radiusKm)
+    {
+        double angularRadius = radiusKm / EarthRadiusKm;
+        if (angularRadius >= Math.PI)
+            return [new GeoBox(-90, -180, 90, 180)]; // the circle covers the whole globe
 
-    internal static GeoBox InternalBoundingBox(double lat, double lon, double radiusKm) => BoundingBox(lat, lon, radiusKm);
+        // Keep the centre in [-180, 180) so the meridian split below yields tight boxes.
+        lon = ((lon + 180) % 360 + 360) % 360 - 180;
 
-    internal static GeoBox BoundingBox(double lat, double lon, double radiusKm)
-    {
-        double dLat = ToDegrees(radiusKm / EarthRadiusKm);
-        double dLon = ToDegrees(radiusKm / (EarthRadiusKm * Math.Cos(ToRadians(lat))));
-
-        return new GeoBox(
-            lat - dLat,
-            lon - dLon,
-            lat + dLat,
-            lon + dLon);
+        double dLat = ToDegrees(angularRadius);
+        double minLat = lat - dLat;
+        double maxLat = lat + dLat;
+
+        // The circle contains a pole exactly when sin(r) >= cos(lat). Testing that ratio, rather
+        // than lat ± dLat in degrees, also keeps the asin argument below 1 when rounding puts
+        // |lat| + r within a few ulp of 90°.
+        double ratio = Math.Sin(angularRadius) / Math.Cos(ToRadians(lat));
+        if (maxLat >= 90 || minLat <= -90 || ratio >= 1)
+            return [new GeoBox(Math.Max(minLat, -90), -180, Math.Min(maxLat, 90), 180)];
+
+        double dLon = ToDegrees(Math.Asin(ratio));
+        double minLon = lon - dLon;
+        double maxLon = lon + dLon;
+
+        if (minLon < -180)
+            return [new GeoBox(minLat, minLon + 360, maxLat, 180), new GeoBox(minLat, -180, maxLat, maxLon)];
+        if (maxLon > 180)
+            return [new GeoBox(minLat, minLon, maxLat, 180), new GeoBox(minLat, -180, maxLat, maxLon - 360)];
+        return [new GeoBox(minLat, minLon, maxLat, maxLon)];
     }
 
     private static double ToRadians(double degrees) => degrees * Math.PI / 180.0;

tests/BLite.Tests/NearAntimeridianPoleTests.cs

@@ -0,0 +1,152 @@
+using BLite.Core;
+using BLite.Core.Indexing;
+using BLite.Shared;
+
+namespace BLite.Tests;
+
+/// <summary>
+/// Indexed <c>Near</c> must find every point within the radius, including points across the
+/// ±180° meridian, on the far side of a pole, and at the circle's widest longitude (#166).
+/// </summary>
+public class NearAntimeridianPoleTests : IDisposable
+{
+    private readonly string _dbPath;
+    private readonly TestDbContext _db;
+
+    public NearAntimeridianPoleTests()
+    {
+        _dbPath = Path.Combine(Path.GetTempPath(), $"near_edges_{Guid.NewGuid():N}.db");
+        _db = new TestDbContext(_dbPath);
+    }
+
+    public void Dispose()
+    {
+        _db.Dispose();
+        foreach (var f in new[] { _dbPath, Path.ChangeExtension(_dbPath, ".wal") })
+            try { if (File.Exists(f)) File.Delete(f); } catch { /* best-effort */ }
+    }
+
+    private async Task<List<string>> NearAsync((double, double) center, double radiusKm)
+    {
+        var viaIndex = (await _db.GeoItems.NearAsync("idx_spatial", center, radiusKm).ToListAsync()).Select(e => e.Name).OrderBy(n => n).ToList();
+        var viaLinq = _db.GeoItems.AsQueryable().Where(p => p.Location.Near(center, radiusKm)).ToList().Select(e => e.Name).OrderBy(n => n).ToList();
+        Assert.Equal(viaIndex, viaLinq);
+        return viaIndex;
+    }
+
+    [Fact]
+    public async Task Near_AcrossTheAntimeridian_FindsPointsOnBothSides()
+    {
+        // 22.2 km apart across ±180°.
+        await _db.GeoItems.InsertAsync(new GeoEntity { Name = "east", Location = (0.0, 179.9) });
+        await _db.GeoItems.InsertAsync(new GeoEntity { Name = "west", Location = (0.0, -179.9) });
+
+        Assert.Equal(new[] { "east", "west" }, await NearAsync((0.0, -179.9), 50.0));
+        Assert.Equal(new[] { "east", "west" }, await NearAsync((0.0, 179.9), 50.0));
+    }
+
+    [Fact]
+    public async Task Near_OverTheNorthPole_FindsPointsOnTheFarSide()
+    {
+        // 161.8 km apart; the shortest path crosses the pole.
+        await _db.GeoItems.InsertAsync(new GeoEntity { Name = "over the pole", Location = (89.5, 150.0) });
+        await _db.GeoItems.InsertAsync(new GeoEntity { Name = "too far", Location = (85.0, 180.0) });
+
+        Assert.Equal(new[] { "over the pole" }, await NearAsync((89.0, 0.0), 200.0));
+    }
+
+    [Fact]
+    public async Task Near_AroundTheSouthPole_FindsPointsAtEveryLongitude()
+    {
+        foreach (var lon in new[] { -180.0, -90.0, 0.0, 90.0, 179.0 })
+            await _db.GeoItems.InsertAsync(new GeoEntity { Name = $"lon {lon}", Location = (-89.9, lon) });
+
+        Assert.Equal(5, (await NearAsync((-89.9, 45.0), 30.0)).Count);
+    }
+
+    [Fact]
+    public async Task Near_AtHighLatitude_FindsPointsAtTheWidestLongitude()
+    {
+        // The circle around (80, 0) with a 556 km radius is widest near 81.3°N, where it reaches
+        // ±30.1° of longitude, more than radius / cos(80°) = ±28.8° at the center latitude.
+        (double, double) center = (80.0, 0.0), widest = (81.3, 29.5);
+        Assert.True(SpatialMath.DistanceKm(80.0, 0.0, widest.Item1, widest.Item2) < 556.0);
+        await _db.GeoItems.InsertAsync(new GeoEntity { Name = "widest", Location = widest });
+
+        Assert.Equal(new[] { "widest" }, await NearAsync(center, 556.0));
+    }
+
+    [Fact]
+    public void NearQueryBoxes_AtThePoleBoundary_NeverContainNaN()
+    {
+        // Radii at and just under the exact distance to the pole, where sin(r)/cos(lat) can round to 1 + ε.
+        for (double lat = 60; lat < 90; lat += 0.37)
+        {
+            double toPoleKm = (90 - lat) * Math.PI / 180 * 6371.0;
+            foreach (var radius in new[] { toPoleKm, toPoleKm * (1 - 1e-15), toPoleKm * (1 - 1e-12) })
+            {
+                foreach (var box in SpatialMath.NearQueryBoxes(lat, 10.0, radius))
+                {
+                    Assert.False(double.IsNaN(box.MinLat) || double.IsNaN(box.MaxLat) || double.IsNaN(box.MinLon) || double.IsNaN(box.MaxLon),
+                        $"NaN box for lat {lat}, radius {radius}");
+                }
+            }
+        }
+    }
+
+    /// <summary>
+    /// Radii equal to the distance to the pole where lat ± r in degrees rounds just short of the
+    /// pole while sin(r)/cos(lat) rounds to 1 + ε, so asin returned NaN before the ratio check.
+    /// </summary>
+    [Theory]
+    [InlineData(-2.1702, 9766.228168206266)]
+    [InlineData(2.1702, 9766.228168206266)]
+    public void NearQueryBoxes_AtThePoleBoundary_KnownRoundingCase_HasNoNaN(double lat, double radiusKm)
+    {
+        foreach (var box in SpatialMath.NearQueryBoxes(lat, 10.0, radiusKm))
+            Assert.False(double.IsNaN(box.MinLat) || double.IsNaN(box.MaxLat) || double.IsNaN(box.MinLon) || double.IsNaN(box.MaxLon),
+                $"NaN box {box}");
+    }
+
+    [Fact]
+    public async Task Near_WithAnUnnormalisedCentreLongitude_MatchesTheNormalisedOne()
+    {
+        Assert.Equal(SpatialMath.NearQueryBoxes(0.0, -170.0, 50.0), SpatialMath.NearQueryBoxes(0.0, 190.0, 50.0));
+
+        await _db.GeoItems.InsertAsync(new GeoEntity { Name = "near", Location = (0.0, -170.1) });
+        await _db.GeoItems.InsertAsync(new GeoEntity { Name = "far", Location = (0.0, 10.0) });
+
+        Assert.Equal(new[] { "near" }, await NearAsync((0.0, 190.0), 50.0));
+    }
+
+    [Fact]
+    public void NearQueryBoxes_LargerThanHalfTheEarth_CoverTheWholeGlobe()
+    {
+        var box = Assert.Single(SpatialMath.NearQueryBoxes(10.0, 20.0, 25_000));
+        Assert.Equal((-90.0, -180.0, 90.0, 180.0), (box.MinLat, box.MinLon, box.MaxLat, box.MaxLon));
+    }
+
+    [Fact]
+    public async Task Dynamic_Near_AcrossTheAntimeridian_FindsPointsOnBothSides()
+    {
+        var path = Path.ChangeExtension(_dbPath, ".dyn.db");
+        try
+        {
+            using var engine = new BLiteEngine(path);
+            var col = engine.GetOrCreateCollection("places");
+            await col.CreateSpatialIndexAsync("pos");
+            await col.InsertAsync(col.CreateDocument(["name", "pos"], b => b
+                .AddString("name", "east").AddCoordinates("pos", (0.0, 179.9))));
+
+            var results = await col.NearAsync("idx_spatial_pos", (0.0, -179.9), 50.0).ToListAsync();
+
+            Assert.True(Assert.Single(results).TryGetString("name", out var name));
+            Assert.Equal("east", name);
+        }
+        finally
+        {
+            foreach (var f in new[] { path, Path.ChangeExtension(path, ".wal") })
+                try { if (File.Exists(f)) File.Delete(f); } catch { /* best-effort */ }
+        }
+    }
+}