diff --git a/src/BLite.Core/Indexing/RTreeIndex.cs b/src/BLite.Core/Indexing/RTreeIndex.cs index 30a67a4..4940664 100644 --- a/src/BLite.Core/Indexing/RTreeIndex.cs +++ b/src/BLite.Core/Indexing/RTreeIndex.cs @@ -57,16 +57,22 @@ public IEnumerable Search(GeoBox area, ITransaction? transacti /// /// Locations of the indexed points within 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. /// public IEnumerable 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; + } } } diff --git a/src/BLite.Core/Indexing/SpatialMath.cs b/src/BLite.Core/Indexing/SpatialMath.cs index 6ce7c1a..d3dc0d8 100644 --- a/src/BLite.Core/Indexing/SpatialMath.cs +++ b/src/BLite.Core/Indexing/SpatialMath.cs @@ -26,22 +26,45 @@ public static double DistanceKm(double lat1, double lon1, double lat2, double lo } /// - /// Creates a Bounding Box (MBR) centered at a point with a given radius in km. + /// Boxes that together cover every point within of the center, for + /// an R-Tree search whose hits are then checked against the exact distance. + /// + /// 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. + /// 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. + /// /// - 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; diff --git a/tests/BLite.Tests/NearAntimeridianPoleTests.cs b/tests/BLite.Tests/NearAntimeridianPoleTests.cs new file mode 100644 index 0000000..1eac867 --- /dev/null +++ b/tests/BLite.Tests/NearAntimeridianPoleTests.cs @@ -0,0 +1,152 @@ +using BLite.Core; +using BLite.Core.Indexing; +using BLite.Shared; + +namespace BLite.Tests; + +/// +/// Indexed Near 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). +/// +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> 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}"); + } + } + } + } + + /// + /// 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. + /// + [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 */ } + } + } +}