From a510ec22749f43abe60ec8634208551f4c6d55e6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ricardo=20Gro=C3=9F?= Date: Mon, 12 Oct 2026 14:14:57 +0200 Subject: [PATCH 1/3] fix(spatial): cover the antimeridian and the poles in Near queries (#166) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The R-Tree query box for Near was plain lat ± dLat, lon ± dLon. A box around longitude -179.9 never reached +179.9, and a circle that contains a pole reaches every longitude while its box only spanned lon ± dLon. The longitude half-width radius / cos(lat) was also too narrow at high latitudes: the circle is widest poleward of its center. SpatialMath.NearQueryBoxes replaces BoundingBox: - a circle that contains a pole gets one box over all longitudes, from its far edge up to the pole; - otherwise the half-width is asin(sin r / cos lat), and a box that crosses ±180° is split into one box on each side. RTreeIndex.SearchNear searches each box and keeps the exact great-circle distance check, so the extra area is trimmed. Co-Authored-By: Claude Opus 5.5 --- src/BLite.Core/Indexing/RTreeIndex.cs | 16 +-- src/BLite.Core/Indexing/SpatialMath.cs | 41 ++++--- .../BLite.Tests/NearAntimeridianPoleTests.cs | 102 ++++++++++++++++++ 3 files changed, 139 insertions(+), 20 deletions(-) create mode 100644 tests/BLite.Tests/NearAntimeridianPoleTests.cs diff --git a/src/BLite.Core/Indexing/RTreeIndex.cs b/src/BLite.Core/Indexing/RTreeIndex.cs index 30a67a47..0152771b 100644 --- a/src/BLite.Core/Indexing/RTreeIndex.cs +++ b/src/BLite.Core/Indexing/RTreeIndex.cs @@ -57,16 +57,20 @@ 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. + 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 6ce7c1ae..d0c791b7 100644 --- a/src/BLite.Core/Indexing/SpatialMath.cs +++ b/src/BLite.Core/Indexing/SpatialMath.cs @@ -26,22 +26,35 @@ 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 InternalBoundingBox(double lat, double lon, double radiusKm) => BoundingBox(lat, lon, radiusKm); - - internal static GeoBox BoundingBox(double lat, double lon, double radiusKm) + internal static GeoBox[] NearQueryBoxes(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 angularRadius = radiusKm / EarthRadiusKm; + double dLat = ToDegrees(angularRadius); + double minLat = lat - dLat; + double maxLat = lat + dLat; + + if (maxLat >= 90 || minLat <= -90) + return [new GeoBox(Math.Max(minLat, -90), -180, Math.Min(maxLat, 90), 180)]; + + double dLon = ToDegrees(Math.Asin(Math.Sin(angularRadius) / Math.Cos(ToRadians(lat)))); + 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 00000000..59979c2c --- /dev/null +++ b/tests/BLite.Tests/NearAntimeridianPoleTests.cs @@ -0,0 +1,102 @@ +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 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 */ } + } + } +} From b64355add0436657144237ab18c29c7ec6ef0e8e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ricardo=20Gro=C3=9F?= Date: Mon, 12 Oct 2026 15:46:49 +0200 Subject: [PATCH 2/3] fix(spatial): address review of antimeridian and pole handling (#166) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - The pole case is decided on the same ratio that feeds asin (sin r >= cos lat), so rounding can no longer produce an asin argument above 1 and NaN bounds that match every index entry. - The centre longitude is normalised into [-180, 180), so an unnormalised centre (e.g. 190) yields the same tight boxes as its normalised form. - A radius of at least half the Earth's circumference short-circuits to a whole-globe box. - The non-overlap comment in SearchNear states why (dLon < 90°). Co-Authored-By: Claude Opus 5.5 --- src/BLite.Core/Indexing/RTreeIndex.cs | 4 ++- src/BLite.Core/Indexing/SpatialMath.cs | 14 ++++++-- .../BLite.Tests/NearAntimeridianPoleTests.cs | 36 +++++++++++++++++++ 3 files changed, 51 insertions(+), 3 deletions(-) diff --git a/src/BLite.Core/Indexing/RTreeIndex.cs b/src/BLite.Core/Indexing/RTreeIndex.cs index 0152771b..4940664b 100644 --- a/src/BLite.Core/Indexing/RTreeIndex.cs +++ b/src/BLite.Core/Indexing/RTreeIndex.cs @@ -63,7 +63,9 @@ public IEnumerable Search(GeoBox area, ITransaction? transacti /// public IEnumerable SearchNear(double latitude, double longitude, double radiusKm, ITransaction? transaction = null) { - // The boxes never overlap, so no point is returned twice. + // 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)) { foreach (var (mbr, location) in SearchEntries(box, transaction)) diff --git a/src/BLite.Core/Indexing/SpatialMath.cs b/src/BLite.Core/Indexing/SpatialMath.cs index d0c791b7..d3dc0d87 100644 --- a/src/BLite.Core/Indexing/SpatialMath.cs +++ b/src/BLite.Core/Indexing/SpatialMath.cs @@ -39,14 +39,24 @@ public static double DistanceKm(double lat1, double lon1, double lat2, double lo 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 + + // Keep the centre in [-180, 180) so the meridian split below yields tight boxes. + lon = ((lon + 180) % 360 + 360) % 360 - 180; + double dLat = ToDegrees(angularRadius); double minLat = lat - dLat; double maxLat = lat + dLat; - if (maxLat >= 90 || minLat <= -90) + // 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(Math.Sin(angularRadius) / Math.Cos(ToRadians(lat)))); + double dLon = ToDegrees(Math.Asin(ratio)); double minLon = lon - dLon; double maxLon = lon + dLon; diff --git a/tests/BLite.Tests/NearAntimeridianPoleTests.cs b/tests/BLite.Tests/NearAntimeridianPoleTests.cs index 59979c2c..854a837c 100644 --- a/tests/BLite.Tests/NearAntimeridianPoleTests.cs +++ b/tests/BLite.Tests/NearAntimeridianPoleTests.cs @@ -76,6 +76,42 @@ public async Task Near_AtHighLatitude_FindsPointsAtTheWidestLongitude() 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}"); + } + } + } + } + + [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() { From c82db8283423e2c0538ff2a11bfa985f323cfb2a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ricardo=20Gro=C3=9F?= Date: Mon, 12 Oct 2026 18:57:44 +0200 Subject: [PATCH 3/3] test(spatial): pin pole-boundary inputs that produced NaN boxes (#166) The 0.37 degree sweep in NearQueryBoxes_AtThePoleBoundary_NeverContainNaN never hits the rounding case, so it passed with the fix reverted. The new cases (lat +/-2.1702, radius 9766.228168206266 km, the distance to the pole) give NaN boxes without the sin(r) >= cos(lat) check. Co-Authored-By: Claude Opus 5.5 --- tests/BLite.Tests/NearAntimeridianPoleTests.cs | 14 ++++++++++++++ 1 file changed, 14 insertions(+) diff --git a/tests/BLite.Tests/NearAntimeridianPoleTests.cs b/tests/BLite.Tests/NearAntimeridianPoleTests.cs index 854a837c..1eac867f 100644 --- a/tests/BLite.Tests/NearAntimeridianPoleTests.cs +++ b/tests/BLite.Tests/NearAntimeridianPoleTests.cs @@ -94,6 +94,20 @@ public void NearQueryBoxes_AtThePoleBoundary_NeverContainNaN() } } + /// + /// 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() {