Geo3D doc-values distance sorts with a small radius treat about half of the nearby documents as missing

Geo3DDocValuesField.newDistanceSort(field, lat, lon, maxRadiusMeters, planetModel) promises ascending distance in metres and +Infinity only for documents missing the field; newOutsideDistanceSort promises 0.0 for points inside the circle. Both comparators decode each doc value (21 bits per coordinate, a grid of about 6 m) to a point a few metres off the planet surface (measured up to 4.48 m) and pass it, unprojected, to the circle's membership test, a single plane that for a 1 km circle lies only about 8 cm below the surface at the centre. So whether the read-back point lands above or below that plane, not its distance, decides in or out: points read back more than about 8 cm below the surface get +Infinity (or a non-zero outside distance) even at the centre, and points read back well above it can be accepted kilometres beyond the rim. With WGS84 and a 1 km radius, 276 of 576 documents inside the circle get +Infinity while Geo3DPoint.newDistanceQuery over the same circle matches all 576; the outside sort gives inside documents up to 988 m. The sampled finite distances were within 4.5 m. For documents well inside the circle the error fades with the radius (125 of 576 at 5 km, 1 at 10 km, 0 at 20 km), but a band near the rim keeps flipping (at 20 km, 83 of 360 documents at 0.98 r get +Infinity), and has been in every release since 6.1.0 (LUCENE-7212, 2016). Projecting the decoded point onto the surface with PlanetModel.createSurfacePoint in the two comparators fixes it; the fix was validated with the repro, a wider independent test and the whole spatial3d module.

• Geo3DDocValuesField (lucene/spatial3d): a per-document location for sorting. It stores x, y and z "21 bits each" in one long (Geo3DDocValuesField.java:57-58); one grid step is about 9.5e-7 planet radii, about 6 m (PlanetModel.java:644-653).
• newDistanceSort: "orders documents by ascending distance from the location", the value is "the distance in meters", and "If a document is missing the field, then by default it is treated as having Double#POSITIVE_INFINITY distance" (Geo3DDocValuesField.java:167-171). maxRadiusMeters "is the maximum radius in meters" (:179).
• newOutsideDistanceSort: "Points within the circle have distance 0.0" (:232-233).
• Membership gate: a shape's distance is +Infinity unless isWithin(x, y, z) (GeoBaseDistanceShape.java:47-53); its outside distance is 0.0 only if isWithin (GeoBaseMembershipShape.java:47-53). For a circle, isWithin is one plane test (GeoStandardCircle.java:174-181).

1. Both comparators decode the stored long into x, y, z with planetModel.getDocValueEncoder() and use that point as is: Geo3DPointDistanceComparator.java:114-116 and :127-131 (compareBottom), :172-177 (computeMinimumDistance, used by copy and compareTop); Geo3DPointOutsideDistanceComparator.java:92-101 and :142-147.
2. Rounding each coordinate to its grid moves the point off the ellipsoid: up to 4.48 m radially and 4.94 m along the surface (measured); half a grid step per axis bounds the move at about 5.3 m.
3. The plane of a cap of radius r lies R(1 - cos r) below the surface at the centre: about 7.8 cm for 1 km, 2.0 m for 5 km, 31 m for 20 km. A point read back e below the surface is inside only within about sqrt(r^2 - 2e) of the centre, which for a 1 km circle and an e of metres is nowhere; a point read back e above the surface passes out to about sqrt(r^2 + 2e), about 7.6 km for e = 4.5 m. (Small-angle, locally spherical estimates in planet radii; on WGS84 the circle is a plane-defined ellipse.)
4. So whether each document's read-back point lands above or below the plane, fixed by its stored cell, decides in or out (a point read back only slightly below the surface can stay inside). Of 10,000 points well inside 1 km circles at 200 centres, 4,805 decoded points were rejected, all of them below the ellipsoid, and 0 after projecting them onto the surface.
5. The geom package requires the opposite. The comparators call the (x, y, z) overload of computeDistance, documented as the distance "to the GeoPoint" (GeoDistance.java:41-46); a GeoPoint "represents a point on the surface of a sphere or ellipsoid" (GeoPoint.java:24), and its (x, y, z) tuple "must be on the desired ellipsoid" (:202-203). The user's points are on it; Lucene's own decode moves them off.

Repro (WGS84, centre 40N 74W, 1 km): of 36 documents at 500 m, 19 sort before every 3 km document and 17 get +Infinity and sort after the 14 of 36 documents at 3 km that get finite distances (with ties broken by doc ID, 22 of the first 36 hits are 500 m documents); the outside sort gives the same 17 an outside distance above 1 m, the worst 722 m. 4 centres, 576 documents at most 700 m from the centre: 276 get +Infinity, 123 documents 2-3 km away get finite values, the worst outside distance is 988 m. The sampled finite values were within 4.5 m of the true distance (per-radius maxima 3.9-4.5 m). Near the rim the flips persist at every radius tested: at 20 km, 11, 83 and 120 of 360 documents at 0.95, 0.98 and 0.99 r get +Infinity and 88 of 360 at 1.02 r get finite distances; at 50 km, 92 of 360 at 0.998 r; with the fix, every ring tested is right. Documents exactly at the centre get +Infinity at 1 of 6, 1 of 6 and 2 of 6 centres on a true sphere, WGS84 and CLARKE_1866.

• Who is exposed. Users of the Geo3D doc-values sorts with an ordinary "nearest within a few km" radius. Defaults are enough: WGS84, the default codec and IndexWriterConfig, no special data and no index-size threshold (each document's flip depends only on its own stored cell; its rank also depends on doc IDs and the other hits). newPathSort end caps probably fail the same way (not executed).
• What they see. No error. At 1 km, about half of the documents inside the radius sort as if they had no location; documents outside the radius take their places with finite distances; the outside sort reports inside documents as hundreds of metres outside. In one probe a 1 km top-50 page missed 27 of the 50 nearest documents (it held no document 1.5 km or more away: the order inside the radius degrades).
• Scale by radius (576 documents at most 0.7 r from 4 centres, WGS84; documents near the rim are not in this sample):

• Near the rim, at any radius. A rim probe (5 centres x 72 bearings per ring) found +Infinity for 11, 83 and 120 of 360 documents at 0.95, 0.98 and 0.99 r of a 20 km circle and finite distances for 88 of 360 at 1.02 r; at 50 km 92 of 360 at 0.998 r, at 100 km 5 of 360 at 0.998 r. The band is about e R / r wide (e the decode error, R the planet radius): about 1.4 km at 20 km.
• Workaround until fixed. A maximum radius of 20 km or more with the application cutting by distance at no more than 0.9 x that radius (rim probe at 20 km: 0 of 360 wrong at 0.9 r, 11 of 360 at 0.95 r), or LatLonDocValuesField sorts.
• Not affected. The hit set (a sort only orders the hits the query chose), stored data (the fix is query-time only, no reindex), LatLonDocValuesField sorts (all 36 near documents first on the same points).
• How widely used. Rarely: no upstream test sorts with Geo3DDocValuesField, and a read-only code search found no use in Elasticsearch, OpenSearch or Solr. No report was found in the searched sources, so the likelihood is rated plausible.

# Lucene main 645d4fe14a20, assertions on (-ea), seed DEADBEEF: 1 test(s) run, 1 failed
    => java.lang.AssertionError: expected a finite distance near 0 m, but was Infinity
       org.junit.Assert.fail(Assert.java:89)
       org.junit.Assert.assertTrue(Assert.java:42)
       org.apache.lucene.spatial3d.TestGeo3DDistanceSortCentreDocument.testDocumentAtCentreGetsFiniteDistance(TestGeo3DDistanceSortCentreDocument.java:51)

• Contract. Geo3DDocValuesField.java:167-171: "This sort orders documents by ascending distance from the location. The value returned in FieldDoc for the hits contains a Double instance with the distance in meters." and "If a document is missing the field, then by default it is treated as having Double#POSITIVE_INFINITY distance (missing values sort last)." :179: "maxRadiusMeters is the maximum radius in meters." :232-233: "Points within the circle have distance 0.0." GeoDistance.java:28-35: distances are for points "expected to be INSIDE the shape", +Infinity "for points outside of the shape". GeoPoint.java:24 and :202-203: a point "on the surface", whose tuple "must be on the desired ellipsoid".
• The sort disagrees with the matching query. The class javadoc pairs the field with a Geo3DPoint for queries (:43); over the same circle the query matches 576 of 576 inside documents and the sort gives 276 of them +Infinity.
• No test. No upstream test calls newDistanceSort, newOutsideDistanceSort or newPathSort on Geo3DDocValuesField; TestGeo3DDocValues checks only the encoding, with a 0.02 tolerance (TestGeo3DDocValues.java:56).
• History. In #8267 (LUCENE-7212, 2016-04-27) mikemccand asked: "Will the geo3d distance checking APIs be upset about the heavy quantization (21 bits per value, instead of 32 bits per value points is using)? I know we've had issues with "the point is not on the earth's surface" already with 32 bits per x/y/z." The author answered: "the resolution is quite important for determination if in-set/out-of-set for shapes like small circles. However, for distance measurement/ordering, I can't see any need to be that precise, so if the docvalues is only used for distance computation then 64 bits should be OK." The comparators committed three weeks later (2016-05-17/18) gate every distance on that in/out test. The same author proposed projecting quantized points onto the surface for the query tests (#8223, #8280).

The checkers built the strongest "works as designed" case: the class WARNING admits "some loss of precision" (Geo3DDocValuesField.java:45-47); +Infinity is the documented result outside the shape and maxRadiusMeters bounds the work; for Geo3DPoint queries the quantized point is deliberately the truth (LUCENE-7168, LUCENE-7225); and LUCENE-7212 is still open with "tests" listed as still needed.

Why it fails: the WARNING's figures (4.19e-8 and 8.38e-8, millimetres if read as degrees) were copied from the 32-bit lat/lon fields, while the real error is about 5 m, and the failure is not imprecision but an in/out flip for documents hundreds of metres inside the circle, even at its centre, with finite values still right to a few metres. +Infinity is meant for documents beyond the radius, not inside it. The quantized-point rule keeps BKD cell relations consistent for queries; a per-document sort has no such constraint, and it disagrees with that very query. An open feature issue explains why nothing caught it, not a decision to accept wrong results in an API shipped since 6.1.0. Verdict: not intended, a defect.

Option A was built and validated in a scratch worktree of a9779ef3c193; B, C and D come from the brief and the checkers and were not built.

Add a private decodeSurfacePoint(long) to Geo3DPointDistanceComparator and Geo3DPointOutsideDistanceComparator and use it in compareBottom (before the bounds check and the distance call) and in computeMinimumDistance. PlanetModel.createSurfacePoint (PlanetModel.java:417-431) scales the point along its ray onto the ellipsoid, which is what the geom package expects.

// compareBottom, before the bounds check:
final GeoPoint point = decodeSurfacePoint(encoded);
// computeMinimumDistance (likewise computeOutsideDistance):
distanceShape.computeDistance(DistanceStyle.ARC, decodeSurfacePoint(encoded));

private GeoPoint decodeSurfacePoint(long encoded) {
  final PlanetModel.DocValueEncoder encoder = planetModel.getDocValueEncoder();
  return planetModel.createSurfacePoint(
      encoder.decodeXValue(encoded),
      encoder.decodeYValue(encoded),
      encoder.decodeZValue(encoded));
}

Trade-offs

• One square root, a few multiplications and one GeoPoint allocation per doc value read, on the sort's hot path: once the queue is full, compareBottom runs for every collected hit (TopFieldCollector.java:92-93, :203-209) and projects before its bounding-box rejection, and a competitive hit projects again in copy. Not benchmarked; an allocation-free variant scales x, y, z inline by the same factor (not built).
• Finite distances move by at most a few metres; documents beyond the radius now get +Infinity as intended.
• Fixes the doc-values sorts only; the same kind of error on the Geo3DPoint query path (32-bit points, tracked separately, likely behind #13056) needs its own fix, and near the rim the query and the sort can still disagree.
• The two comparators also serve newPathSort and the outside box, polygon, large-polygon and path sorts (Geo3DDocValuesField.java:216-224, :281-373); only the circle sorts were tested directly. It does not fix the separate outside-distance bug for points outside a small circle on WGS84 (LFV-274).

What the checks validated

• Scratch worktree of a9779ef3c193, gradle, JDK 25: stock + the new test, 2 of 3 fail with assertions on and off; + the fix, the new test with TestGeo3DDocValues and TestGeo3DPoint (3 suites, 21 tests, 2 skipped) pass, and with assertions off the new test with TestGeo3DDocValues (2 suites, 4 tests), the whole lucene/spatial3d module (20 suites, 184 tests, 3 skipped) passes with assertions on and off, and :lucene:spatial3d:check -x test (google-java-format, ecjLint, forbiddenApis) passes.
• Standalone runner: the new test 3 of 3 with assertions on and off; the independent 4-centre test 5 of 5, with 0 wrong documents at every radius from 250 m to 50 km; the finder's test passes.
• The finder's equivalent fix passed in every recorded run that applied it: the independent reproducer's run, the adjudication's run, the security probe and 6.1.0. A separate probe showed that projecting the points removes all 17 +Infinity values; the devil's advocate did not run a fix.

Project inside GeoStandardCircle.isWithin (or a shared helper), or give the plane test a tolerance scaled to the point's height above or below the surface.

Trade-offs

Could also fix the query path, but changes query semantics and must stay consistent with the shapes' bounds and getRelationship used by the BKD visitor. A change in GeoStandardCircle.isWithin reaches every circle user, queries included; a shared helper or a general plane tolerance reaches every shape, and a tolerance can admit points beyond the intended rim. Needs the full spatial3d random tests and review. Projecting at index time instead does not work: the per-axis quantization that follows moves the point off the surface again.

Skip the plane test for circles in the comparators.

Trade-offs

The comparators accept any GeoDistanceShape (paths too), so it needs a special case per shape, and on WGS84 a plain arc-distance test would not keep the plane-defined ellipse GeoStandardCircle uses; A is simpler.

Say that a small maxRadiusMeters gives unreliable results on WGS84 and that results near the rim are unreliable at any radius (for example: use at least 20 km and cut at no more than 0.9 x maxRadius).

Trade-offs

Not a fix; the radius is the user's choice and the javadoc gives no hint today.

Method: a hypothesis from a model, then witnessed by a tool (semi-formal).

• In the spatial3d hunt unit, a sibling candidate on the Geo3DPoint query path had an exact-rational solver model of a small cap tested on a quantized point; it bounds the effective rim at about sqrt(r^2 + 2e) for a radial error e.
• The finder applied that bound to the doc-values grid (step about 9.5e-7 planet radii, radial error up to about 8e-7), which predicts a rim shift up to about sqrt(2 x 8e-7) = 1.3e-3 rad, about 8 km. The 21-bit prediction is a closed-form extrapolation of that model, not a separate solver run.
• A radius sweep on the stock main jar (100 m to 20 km) and a public-API JUnit then showed the effect: 57-58 of 108 inside documents with +Infinity at 100 m and 1 km, 23 at 5 km, 0 at 20 km.