Geodetic (lon/lat) R-tree with great-circle queries#233
Conversation
urschrei
force-pushed
the
shugel/geodetic_3d
branch
2 times, most recently
from
430e405 to
1acb3a0
Compare
|
My first reaction is not enthusiastic. This is too much code to review at once and the likelihood of it eventually landing would be higher if it were split into more easily digestible parts. Also considering
there is the question of whether this just should be a separate crate altogether. Personally, I am certainly not keen on maintaining 7.500 lines of LLM-generated code which also pulls in a new test framework. (The sheer verbosity of it all is the most impressive property of it IMHO.) I also think the rationale in the crate-level docs is too strongly worded IMHO. While the errors are untenable if one wants to index geometry over the whole globe, the just-treat-geodetic-like-cartesian-coordinates can work and is much simpler/more efficient for reasonably bounded regions like geometry within a single European state, i.e. just having geodetic coordinates at all is not automatically a reason to use an S^2-based index. |
|
Hmm I hadn't considered that it could "just" be a new crate. I'm not overly concerned about the maintenance burden (I wrote the module, so it's ultimately my responsibility). The lack of stacked diffs makes the review process difficult, and I tried to mitigate that as much as possible with logical commits, almost all of which are independently reviewable. That having been said, maybe it makes more sense as a standalone crate. |
I don't think it works like that. This is an ecosystem crate under joint maintenance by multiple people, so at least my expectation is that most of the maintainers understand most of the code. I also think that this is the signal this crate sends via its metadata and documentation. If the goal really is an additive bus-factor-one module, then it should be a separate crate. Otherwise we need to work through the process were code merged into this repository is reviewed and understood by multiple maintainers which given the limited bandwidth most likely requires smaller easier-to-understand increments. It is not just about what the code does or even how efficiently it does it, the main issue is that coding is theory building and FOSS is building shared theory building, so it is primarily about using code to teach people something. |
That's fair enough. In terms of the "core" code that needs to be thoroughly understood (though it will all have to be reviewed in some form), the footprint is a lot smaller (I estimate around 2500 LoC). The rest is almost all tests. I'm open to any potential review strategy and/or re-arrangement of the commit series within reason. But I'm also open to a separate crate – it doesn't have the potential advantage you note above, of course. |
…istance) Add the `geodetic` feature and its embedding foundation: GeodeticCoord (degrees, longitude first), the UnitVec n-vector embedding of (lon, lat) onto the unit sphere, and the great-circle distance helpers (haversine_distance, squared_chord, and the squared-chord/metre conversions), with inline unit and Hegel property tests. The leaf types and the tree build on these; there is no public tree yet.
Add GeodeticPoint (a lon/lat leaf) and GeodeticRTree, a great-circle spatial index built on a stock RTree of unit vectors. Supports nearest-neighbour, nearest-with-distance (metres), radius (locate_within_distance), exact-location, and longitude/latitude window queries, with the antimeridian and poles handled as ordinary interior points. The tree is generic over its leaf type via the open GeodeticObject marker trait, defaulting to GeodeticPoint. Covered by integration, property (Hegel), and embedding tests with a shared oracle module.
Add the great-circle arc primitives that extent leaves and custom leaf types build on: arc_distance_2 (squared-chord distance from a point to a great-circle arc), nearest_point_on_arc, arc_bounding_box (a unit-sphere box covering an arc's bulge), and arc_contains_point. Validated by a property test whose oracle and ARC_GOLDEN rows are anchored to an external 60-digit reference (tests/reference/arc_distance_reference.py, cross-checked against s2sphere). The GeodeticObject documentation now demonstrates a custom arc-based leaf.
Add two extent leaf types for the geodetic R-tree. GeodeticLineString is a great-circle polyline with arc-aware bounding boxes (each box inflated to cover its edges' bulge) and nearest-point distance. GeodeticPolygon (with GeodeticRing) is a filled spherical polygon: membership is a great-circle ray cast (zero distance inside, exact distance outside via robust orient3d predicates), so any simple polygon is supported, including ones larger than a hemisphere. Both index through the existing GeodeticObject leaf trait. Adds the `robust` dependency for the polygon predicates, behind the `geodetic` feature, and a property test suite anchored to the arc oracle.
Add the optional `geodetic-wgs84` feature: an ellipsoidal geodesic refine for the point GeodeticRTree, layered on the spherical index as a filter/refine. The reference ellipsoid is an Ellipsoid value (WGS84, GRS80, or a custom Ellipsoid::new / from_inverse_flattening). nearest_neighbor_on_ellipsoid, nearest_neighbor_with_distance_on_ellipsoid, and locate_within_distance_on_ellipsoid re-rank or filter the spherical candidates by the exact geodesic distance (Karney, via geographiclib-rs), with _wgs84-suffixed shorthands; geodesic_distance is the standalone reference distance. The spherical lower bound keeps the branch-and-bound search sound. The feature requires std; the base geodetic feature stays no_std.
Add the geodetic benchmark suite (point, extent, and WGS84 benches) and update rstar-benches to rand 0.10, dropping the geo-types dependency. Wire the geodetic and geodetic-wgs84 features into CI (test, clippy with -D warnings, and the no_std build), document the feature in both READMEs, and record the new public API in the changelog.
urschrei
force-pushed
the
shugel/geodetic_3d
branch
from
1acb3a0 to
43bf657
Compare
It is your work and hence your choice. If you want to get this out there for people to use with the minimum amount of fuzz, a separate crate is a reasonable choice IMHO. If you want it to become part of this project, you'll probably have to break it down a bit to take along the rest of us and teach us how it works. Personally, I am interested in that and would like to invest some effort into it, but I understand that the amount of time you can spend on breaking things down is limited. So I am not asking for a slide deck, just a bare-bones version first without any of the bells and whistles and extra verification layers. |
2 participants
This PR Adds an optional
geodeticmodule: an R-tree over longitude/latitude data thatanswers nearest-neighbour and radius queries, returning great-circle distance, for
points, linestrings, and polygons. The module is feature-gated and purely additive
What it is
GeodeticRTree<G>indexesGeodeticPoint,GeodeticLineString, andGeodeticPolygonleaves (generic over an openGeodeticObjecttrait, so downstreamcrates can supply their own). Queries take degrees in and return metres out; the
antimeridian and poles need no special handling.
How it works
Each
(lon, lat)is mapped to a unit vector on the sphere (an n-vector) and stored ina stock
RTreewithAABB<UnitVec>envelopes. Because the embedding is continuousover the whole sphere, the ±180° seam and the poles are ordinary interior points requiring no
wrapping or splitting. Squared chord length (3-D Euclidean) is strictly monotonic in
great-circle angle, so the R*-tree's branch-and-bound traversal orders results by true
great-circle distance, and the box
MINDISTis a sound lower bound – the unit-sphereindexing method from a 2013 paper by Schubert, Zimek & Kriegel that proves the
lower-bound and monotonicity properties the pruning relies on. The same embedding also
underlies PostGIS
geography, S2, and H3.with nearest-point-on-arc distance at the leaf; polygon membership is a great-circle
ray cast with pole containment, so any simple polygon works, including ones larger
than a hemisphere. Envelopes are precomputed and cached.
The optional
geodetic-wgs84feature layers an exact ellipsoidal refine (Karney'sgeodesic, via
geographiclib-rs) on top of the spherical index as a soundfilter/refine, selectable per
Ellipsoid(WGS84/GRS80/custom).geodeticisno_std;geodetic-wgs84requiresstd.How it's tested
scan; the node envelope is a true
MINDISTlower bound; arc boxes contain their arcs;point-to-geometry distance matches a cross-track/cap oracle.
mpmath+s2spherereference, with a curated golden set (ARC_GOLDEN). Because thatground truth is computed by external tools rather than by rstar's own code, the index
can't pass a test by sharing the same mistake as the oracle it is checked against.
polygons.
Headline benchmarks (M2 Pro)
Geodetic vs a naive planar
RTree<[f64; 2]>over the same lon/lat data (100k-pointtree, criterion):
Reviewing
This is obviously a large feature. Where possible, I've tried to break up the commits logically, so they can (largely) be reviewed bottom-up, and individual module components have been kept relatively compact for ease of review.
LLM Use
I wrote this module in concert with Claude Opus 4.8, particularly the property tests, embedding, and distance calculations – I'm happy to answer any questions about the code.
rstar/CHANGELOG.mdif knowledge of this change could be valuable to users.