WKB/WKT geometry I/O: ingest + dissolved emit (issue #71)#89
Conversation
Codecov Report❌ Patch coverage is
Additional details and impacted files@@ Coverage Diff @@
## main #89 +/- ##
==========================================
+ Coverage 93.61% 96.72% +3.11%
==========================================
Files 25 27 +2
Lines 3728 4311 +583
==========================================
+ Hits 3490 4170 +680
+ Misses 238 141 -97
Flags with carried forward coverage won't be shown. Click here to find out more.
... and 2 files with indirect coverage changes Continue to review full report in Codecov by Harness.
🚀 New features to boost your workflow:
|
espg
left a comment
espg
left a comment
There was a problem hiding this comment.
🤖 from Claude (review)
Adversarial fresh-context review of Phase 1 only (commit 33c4fef — mortie/geometry.py, its tests, and shapely>=2 added to the test extra). Checked against the locked design in issue #71 and CLAUDE.md §4/§7.
What's solid: runtime stays numpy-only (only import numpy as np at module top; backends imported lazily inside _require_backend, mirroring arrow.py:31), shapely>=2 correctly lands in the test extra only and [project].dependencies is untouched (§7 satisfied), the ImportError message names both shapely and spherely, the (x,y)=(lon,lat) axis flip (line 123) and the EWKT SRID= strip-on-ingest (lines 68-79) are correct, plain-vs-EWKB/EWKT default is right per the sign-off, and .copy() on the sliced coordinate columns avoids handing out non-contiguous views. Module is 169 lines — well under the ~1000 limit. The codec round-trip and SRID-opt-in tests are good.
Findings (4, inline):
- (most important) The spherely path is accepted at the gate but cannot actually drive
decompose/SRID. The code calls the shapely functional API verbatim (get_coordinates,get_exterior_ring,get_interior_ring,get_num_interior_rings,get_parts,set_srid,to_wkb(include_srid=...)), none of which exist in spherely's public API. The plan called for a 2-API shim; Phase 1 has none. Either build the shim or stop advertising spherely. The spherely branch is also 100% untested. - Empty geometries are not rejected despite the docstring saying they are —
POLYGON EMPTYetc. share type ids with non-empty forms and slip through to return a zero-length ring, which only fails (misleadingly) downstream in_prep_rings. - Z coordinates silently dropped — correct behavior for 2D mortie, but undocumented and untested; a real EWKB ingest case.
- Dead code —
_POLYGONAL/_LINEARand three_TYPE_*constants are unused (§4).
None of these are blockers for a draft phase commit, but (1) and (2) are correctness/contract issues I'd fold into the next commit before Phase 2 builds ingest on top of decompose. Review only — not approving or requesting changes.
Generated by Claude Code
| "mortie's WKB/WKT geometry I/O requires a geometry backend; install " | ||
| "`shapely>=2` (preferred) or `spherely` (e.g. `pip install shapely`). " | ||
| "mortie's runtime is numpy-only, so the backend is an optional extra." | ||
| ) |
There was a problem hiding this comment.
🤖 from Claude (review)
The spherely path is accepted at the gate but is not actually a working backend — this is the main Phase-1 finding. The module docstring and this gate advertise "shapely>=2 preferred, spherely accepted," and the locked design (issue #71, 00:34 / 10:10 sign-offs) called for a "small 2-API compat shim" because "their coordinate accessors differ." But the rest of the module calls the shapely functional API verbatim (mod.get_coordinates, mod.get_exterior_ring, mod.get_interior_ring, mod.get_num_interior_rings, mod.get_parts, mod.set_srid, mod.to_wkb(include_srid=...)).
spherely's public API (https://github.com/benbovy/spherely docs/api.rst) exposes only from_wkb/from_wkt/to_wkb/to_wkt/get_type_id/get_x/get_y/get_dimension — it has no get_coordinates, get_exterior_ring, get_interior_ring, get_num_interior_rings, get_parts, or set_srid. So if shapely is absent and spherely is present, _require_backend() succeeds, geometry_from_wkt works, but decompose() raises AttributeError: module 'spherely' has no attribute 'get_coordinates' and geometry_to_wkb(geom, srid=4326) raises on set_srid.
Two honest ways forward: (a) build the shim the plan actually called for, or (b) if shapely is the only backend you intend to support in this PR, drop spherely from the gate and the docstring so the code doesn't promise a backend it can't drive. Right now it's a latent bug that no test catches (the test module hard-skips on shapely and calls shapely.* directly — the spherely branch is 100% untested).
Generated by Claude Code
| raise ValueError( | ||
| f"unsupported geometry type for coverage (type id {type_id}); " | ||
| "expected Polygon, MultiPolygon, LineString, or MultiLineString" | ||
| ) |
There was a problem hiding this comment.
🤖 from Claude (review)
Empty geometries are NOT rejected, contradicting the docstring (lines 148-149: "Points, geometry collections, and empty geometries are rejected"). POLYGON EMPTY / MULTIPOLYGON EMPTY / LINESTRING EMPTY carry the same type ids as their non-empty forms (3/6/1), so they hit the accepting branches above, not this raise. Verified against shapely 2.1.2:
>>> shapely.get_type_id(shapely.from_wkt('POLYGON EMPTY'))
3
>>> np.asarray(shapely.get_coordinates(shapely.get_exterior_ring(shapely.from_wkt('POLYGON EMPTY')))).shape
(0, 2)
So decompose(POLYGON EMPTY) returns ("polygonal", [(array([]), array([]))]) — a zero-length ring — rather than raising. That degenerate ring then flows into phase 2 and only blows up later in _prep_rings (coverage.py:63) with a misleading "each ring needs at least 3 vertices". Either reject empties explicitly here (e.g. mod.is_empty(geom) for shapely, but note that's another method spherely lacks — see the backend comment) or drop the "empty geometries are rejected" sentence. Whichever you choose, add a test pinning the behavior; there's currently none for the empty case.
Generated by Claude Code
| """Extract a ring's vertices as ``(lat, lon)`` float64 arrays (degrees).""" | ||
| coords = np.asarray(mod.get_coordinates(ring_geom), dtype=np.float64) | ||
| # WKB/WKT store (x, y) = (lon, lat). | ||
| return coords[:, 1].copy(), coords[:, 0].copy() |
There was a problem hiding this comment.
🤖 from Claude (review)
3D (Z) coordinates are silently dropped, undocumented and untested. mod.get_coordinates(ring_geom) defaults to include_z=False, so for a POLYGON Z (...) / EWKB with Z this returns an (N, 2) array and the Z band is discarded:
>>> np.asarray(shapely.get_coordinates(shapely.from_wkt('POLYGON Z ((0 0 5,1 0 5,1 1 5,0 1 5,0 0 5))'))).shape
(5, 2)
Dropping Z is the right behavior for mortie (2D lon/lat), so no code change is strictly required — but it's a real ingest case (PostGIS EWKB commonly carries Z) that's neither documented nor pinned by a test. Worth a one-line note in the module/decompose docstring ("Z/M dimensions are dropped; mortie works in 2D lon/lat") plus a test feeding a ... Z ... geometry and asserting the rings come back 2D. The _TYPE_* constants at lines 24-31 imply a deliberate type taxonomy; the dimension axis deserves the same explicitness.
Generated by Claude Code
| _TYPE_GEOMETRYCOLLECTION = 7 | ||
|
|
||
| _POLYGONAL = (_TYPE_POLYGON, _TYPE_MULTIPOLYGON) | ||
| _LINEAR = (_TYPE_LINESTRING, _TYPE_MULTILINESTRING, _TYPE_LINEARRING) |
There was a problem hiding this comment.
🤖 from Claude (review)
Minor (dead code, §4 "no dead code"): _POLYGONAL and _LINEAR are defined here but never referenced — decompose branches on the individual _TYPE_* constants instead, and _TYPE_POINT, _TYPE_MULTIPOINT, _TYPE_GEOMETRYCOLLECTION are likewise unused. flake8's --select=E9,F63,F7,F82 won't flag unused module globals, but the non-blocking style pass / ruff (F841 is locals only, so these slip through too) — regardless, they're speculative surface. Either wire decompose to use if type_id in _POLYGONAL: etc. (cleaner) or drop the unused tuples and the three unused single constants.
Generated by Claude Code
Merging this PR will not alter performance
Comparing 
Footnotes
|
espg
left a comment
espg
left a comment
There was a problem hiding this comment.
🤖 from Claude (review)
Phase 2 review (commit 483a9df, "phase 2 of issue #71") — ingest layer + folded Phase-1 self-review fixes. Checked against issue #71's locked design (shapely≥2-primary / spherely-accepted lazy codec, numpy-only runtime, plain-WKB-default) and CLAUDE.md §4/§7.
Overall this is clean and faithful to the agreed design: routing matches the plan (polygonal → morton_coverage / morton_coverage_moc, linear → linestring_coverage), the single-vs-multi linestring return shape mirrors linestring_coverage's contract, EWKT SRID-strip on ingest is handled, exports/__all__ are consistent, and the _require_shapely NotImplementedError-for-spherely choice is honestly flagged in the docstring + PR rather than silently shimmed — a reasonable in-scope call.
Verified empirically against the built extension that the realistic winding worry is a non-issue: for sub-hemisphere polygons (single, multipart, and donut-with-hole, including a reversed-hole spelling) the even-odd descent is orientation-insensitive, so flat(normalize=True), flat(normalize=False), and the moc path all agree cell-for-cell. So WKT polygons coming out of shapely/GEOS with CW exteriors ingest correctly today. Good.
A few findings below — all low severity; none block the phase.
Findings:
- The moc-path winding gap is real but only at hemisphere-plus scale, and it has no test — see inline on
from_geometry. normalizeis silently ignored on the linear path whilemoc/tolerance/max_cellsare loudly rejected — inconsistent guard; inline.- The moc ingest test is tautological w.r.t. winding (both sides use the same un-normalized array) — inline. Worth one CW-spelled WKT assertion to actually exercise the WKB/WKT→cover winding path.
Nothing here is a correctness bug for the common case; (1) and (3) are about test coverage honesty for the deferred hemisphere-plus case, (2) is a small API-consistency nit.
Generated by Claude Code
espg
left a comment
espg
left a comment
There was a problem hiding this comment.
🤖 from Claude (review)
Inline findings for Phase 2 — three low-severity items, none blocking. See per-line comments. Summary of the verification behind them: I checked routing/return-shape against the issue #71 locked design and ran the built extension to confirm sub-hemisphere polygons (single / multipart / donut) are orientation-insensitive across flat(normalize=True/False) and moc, so the common WKB/WKT ingest case is correct today. The notes are about (a) documenting the moc-path winding limitation at hemisphere-plus scale, (b) a normalize guard inconsistency on the linear path, and (c) making the moc/winding tests non-tautological + adding a WKB moc assertion.
Generated by Claude Code
| if moc: | ||
| return morton_coverage_moc( | ||
| lats, lons, order=order, tolerance=tolerance, max_cells=max_cells | ||
| ) | ||
| return morton_coverage(lats, lons, order=order, normalize=normalize) |
There was a problem hiding this comment.
🤖 from Claude (review)
Winding gap on the moc=True branch — real but narrow, and currently untested. morton_coverage_moc has no normalize parameter, so when moc=True the user-supplied normalize is dropped and the moc descent runs with whatever orientation handling the Rust path bakes in. I verified against the built extension that for sub-hemisphere rings this is harmless: single, multipart, and donut(+reversed-hole) covers all agree cell-for-cell across flat(normalize=True), flat(normalize=False), and moc. So the realistic "shapely emits CW exteriors" case ingests correctly.
The gap is only for hemisphere-plus polygons, where orientation is the only thing that disambiguates interior-vs-complement (per morton_coverage's Ring-winding note, coverage.py:154-162). On the flat path from_geometry lets the caller pass normalize=False and wind the ring themselves; on the moc path there is no such lever, so a WKB/WKT hemisphere-plus polygon whose ring is wound "wrong" will select the complement with no way to correct it through ingest.
This matches the design's deferral of hemisphere-plus winding, so I'm not asking you to build a lever here — but the docstring (line 220-221) says only "Ignored when moc=True," which undersells it. Suggest one extra sentence noting that for moc the ring orientation is taken as-authored (hemisphere-plus rings must already follow the CCW-exterior / CW-hole contract), so the limitation is documented at the API surface rather than discoverable only by surprise.
Generated by Claude Code
| if moc or tolerance is not None or max_cells is not None: | ||
| raise ValueError( | ||
| "moc / tolerance / max_cells apply only to polygonal geometry" | ||
| ) |
There was a problem hiding this comment.
🤖 from Claude (review)
Inconsistent guard on the linear path: moc / tolerance / max_cells are loudly rejected here as polygon-only, but normalize — also polygon-only (linestring_coverage has no such parameter) — is silently ignored. I confirmed from_geometry(linestring, normalize=False) returns a cover with no error, so a caller who passes normalize=False expecting it to mean something on a linestring gets a silent no-op.
Two consistent options: either fold normalize into the rejection check on this branch (treat any non-default normalize as a polygon-only arg, same as the other three), or leave it permissive but say so in the docstring. Folding it in matches the "reject polygon-only args on linear input" intent that this very block establishes; right now normalize is the one polygon-only knob that escapes it.
Generated by Claude Code
| def test_ingest_polygon_moc_matches_array_path(): | ||
| want = mortie.morton_coverage_moc(_LATS, _LONS, order=8) | ||
| wkt = _poly_wkt(_LATS, _LONS) | ||
| assert np.array_equal(mortie.from_wkt(wkt, order=8, moc=True), want) |
There was a problem hiding this comment.
🤖 from Claude (review)
This moc ingest test is tautological with respect to winding (and partly with respect to ingest). want comes from morton_coverage_moc(_LATS, _LONS, ...) and the asserted value comes from from_wkt(_poly_wkt(_LATS, _LONS), ..., moc=True), which decomposes back to the same _LATS/_LONS arrays and calls the same morton_coverage_moc. So both sides traverse the identical un-normalized array path — the test proves the moc routing wires up, but it can never detect a winding/normalize mismatch between the geometry path and the array path, because there is no independent reference.
Two gaps worth one assertion each:
- WKB is not exercised for moc.
test_ingest_polygon_matches_array_pathdoes coverfrom_wkbfor the flat case, but the moc test is WKT-only. Afrom_wkb(geometry_to_wkb(...), moc=True)assertion would close that. - Winding is never exercised through ingest. Note
_LATS/_LONSis actually wound clockwise in (lon, lat) (signed area −200), and it happens to round-trip only because the descent is orientation-insensitive at sub-hemisphere scale. A test that ingests a deliberately CW-spelled WKT polygon and asserts it equals the CCW-spelled cover (flat path, bothnormalizevalues) would actually pin the "WKT polygons from shapely may not follow mortie's CCW contract" behavior the issue thread called out — right now nothing does.
Generated by Claude Code
espg
left a comment
espg
left a comment
There was a problem hiding this comment.
🤖 from Claude (review)
Phase 3 fresh-context adversarial review (commit phase 3 of issue #71, the per-cell emit layer). Scope: _per_cell_polygons, to_geometry/to_wkb/to_wkt, exports, and the new emit tests.
Overall this is in good shape and matches the locked design: dissolve=True default raising NotImplementedError pointing at phase 4 is the right interim for a WIP draft (it does not leave a silently-wrong default — it fails loudly with a clear message), per-cell MultiPolygon as the dissolve=False opt-in is the agreed shape, plain WKB/WKT default with SRID opt-in is plumbed correctly through the phase-1 codec, and step=1 (straight chords) is exposed per the plan. The scalar-vs-array return of mort2polygon is handled correctly (scalar branch wraps in a list; array branch is already one-polygon-per-cell), the group-by-order loop preserves all cells, and the (lat,lon)→(lon,lat) flip is right. Reuse of mort2polygon (incl. its _normalize_antimeridian_polygon) for the boundary is sound. Module is 367 lines, well under the ~1000 limit.
Findings are inline: one real test-strength issue and two missing edge-case suggestions. Nothing blocking phase 4. Left standing for @espg to triage.
Generated by Claude Code
espg
left a comment
espg
left a comment
There was a problem hiding this comment.
🤖 from Claude (review)
Phase 3 inline findings (per-cell emit). See the three per-line comments: one real test-strength issue (the sort-then-compare corner check is too weak to catch a lon/lat flip), missing edge cases (no antimeridian/pole/step>1/single-cell emit tests despite the plan's "antimeridian/pole cells don't self-cross" criterion), and one non-blocking line-length style nit. None block phase 4; left standing for @espg to triage.
Generated by Claude Code
| assert geometry.to_wkt(cov, dissolve=False, srid=4326).startswith("SRID=4326;") | ||
| assert int(shapely.get_srid(shapely.from_wkb( | ||
| geometry.to_wkb(cov, dissolve=False, srid=4326)))) == 4326 |
There was a problem hiding this comment.
🤖 from Claude (review)
(most important) This corner-match assertion is weaker than it reads. np.sort(..., axis=0) sorts the lon column and the lat column independently, so it discards both the vertex order and the (lon, lat) pairing. The check passes as long as the multiset of all lons and the multiset of all lats each agree — it would still pass if lon/lat were swapped per-vertex, or if the polygon were scrambled into a self-crossing ring, as long as the column multisets line up. Since the whole point of this test is to prove the (lat,lon)→(lon,lat) flip in _per_cell_polygons is correct, the assertion can't actually catch a flip bug.
Suggest comparing the coordinate pairs as a set instead, e.g. sort rows lexicographically (ring[np.lexsort(ring.T)] on both sides) or compare set(map(tuple, ...)), so a per-vertex lon/lat swap is caught.
Minor, same test: shapely.get_geometry(back, 0) is matched against cov[0] — that ordering only holds because a flat single-order cover has one depth group, so _per_cell_polygons preserves input order. For a mixed-order MOC the emit order is grouped by np.unique(depths), so geometry index 0 would not be moc[0]. Fine for this test, but worth a one-line comment so a future reader doesn't assume index-stable emit ordering in general.
Generated by Claude Code
| cov = mortie.morton_coverage(_LATS, _LONS, order=6) | ||
| with pytest.raises(NotImplementedError, match="phase 4"): | ||
| geometry.to_wkb(cov) | ||
|
|
There was a problem hiding this comment.
🤖 from Claude (review)
Missing edge cases for per-cell emit — the issue #71 plan's Phase-3 acceptance criterion was explicitly "antimeridian/pole cells don't self-cross," but no emit test exercises a dateline- or pole-touching cell. All four fixtures here use the _LATS=[40,50] / _LONS=[-120,-110] mid-latitude box, which is the easy case. _per_cell_polygons leans on mort2polygon's _normalize_antimeridian_polygon to keep antimeridian-touching cells from self-crossing, and that path is currently unverified through the emit layer.
Suggest adding (each cheap):
- a cover near ±180° lon and a polar-cap cell, asserting
shapely.is_validon every emitted part (catches self-intersection); - a
step>1emit asserting the ring has4*stepvertices (the densification knob is plumbed but untested); - a single-cell cover (
np.array([one_word])) to pin thegrp.size == 1scalar branch of_per_cell_polygons— that branch is distinct code (mort2polygon(int(...))scalar return) and currently only hit incidentally, if at all.
Generated by Claude Code
| See :func:`to_geometry` for ``dissolve`` / ``step``. With ``srid`` set | ||
| (e.g. ``4326``), emit EWKB carrying that SRID; otherwise plain WKB. | ||
| """ | ||
| return geometry_to_wkb(to_geometry(morton, dissolve=dissolve, step=step), srid=srid) |
There was a problem hiding this comment.
🤖 from Claude (review)
Style-only (non-blocking): these two geometry_to_wkb(...)/geometry_to_wkt(...) one-liners are 92 and 88 chars, over the non-blocking --max-line-length=88 style pass CI runs (§7). Line 357 in particular exceeds it. The hard flake8 --select=E9,F63,F7,F82 gate won't flag these, but the style pass will. A wrap (or pulling to_geometry(...) into a local first) keeps the file clean.
Generated by Claude Code
|
🤖 from Claude Status — phases 1–3 landed green; phase 4 (the default dissolved emit) continues next run. This run landed the codec adapter (phase 1), ingest (phase 2), and per-cell emit (phase 3), each as its own commit with a fresh-context self-review whose diff-scoped findings were folded into the following commit:
Local gate green: Why I'm checkpointing before phase 4 rather than rushing it: phase 4 is the dissolved-outline emit via native edge-cancellation (approach (i)) — exact cross-order corner identity across the 12 base faces with seam canonicalization, edge subdivision/cancellation, ring chaining, hole/component classification, then antimeridian/pole splitting, plus the Rust kernel you OK'd. That's the single hardest, most correctness-critical piece, and shipping it unverified would be worse than landing it carefully next run. The phases 1–3 infrastructure it builds on is done and tested, so phase 4 starts from a solid base. PR stays Generated by Claude Code |
espg
left a comment
espg
left a comment
There was a problem hiding this comment.
🤖 from Claude (review)
Fresh-context adversarial review of phase 4 (commit 9809406, dissolved-outline emit). This is a clean, well-commented increment and the test set is a genuine improvement (union oracle, hole, disjoint, mixed-order MOC, area-conservation, antimeridian, pole-raise). Findings below — inline comments carry the specifics.
Blocking correctness (please address or justify before clearing draft):
- Ring chaining at pinch / non-manifold vertices (
_boundary_rings_xyz, theout[cur].pop()chase): arbitrary edge choice at a degree>1 vertex can mis-chain diagonally-touching / bowtie covers into a self-touching or wrongly-merged ring. Needs angular-order edge selection and a pinch-cover test; none of the current tests exercise a single-corner touch. - Antimeridian split for >2 crossings (
_split_at_antimeridian,segments[0] = cur + segments[0]): the wrap-merge only holds for exactly two crossings; ≥4-crossing rings (legal, non-pole) get hemispheres mixed. Only the 2-crossing box is tested.
Correctness / robustness (worth fixing now):
3. Hole-nesting fallback (best if best is not None else 0): a hole whose vertex-mean centroid lands in no exterior is silently force-attached to piece 0 → possibly invalid/wrong geometry. Use a guaranteed-interior point and surface (don't swallow) a genuinely-uncontained hole.
4. Global winding sum(areas) < 0: spherical signed area is mod 4π, so the exterior/hole sign flips for covers near a hemisphere. Document + guard the precondition, or derive winding robustly.
5. _DISSOLVE_SNAP = 1e10 under step>1: exact for step=1, but a real ~5e-11 tolerance for cross-cell sub-edge points; add a seam-cancellation test (area test is insensitive to slivers).
6. Pole detection via odd-crossing count is a proxy, not a direct test — a pole-enclosing ring with even crossings would be mis-split silently rather than raised. Note/strengthen.
Design:
7. Pure-Python vs the locked Rust-helper plan — the whole dissolve hot path is Python; CLAUDE.md §7 puts perf-critical core in Rust. Reasonable as a correctness-first draft, but it's a plan deviation needing an @espg call (land-now-port-later vs Rust-before-merge).
Non-blocking notes: module is at 607 lines (well under the ~1000 soft limit — fine). No new dependencies. Style/terseness read clean and match the surrounding file; nothing flagged for the E9/F63/F7/F82 gate. No dead code spotted.
Items 1 and 2 are the ones I'd gate the “ready for review” transition on.
Generated by Claude Code
| for _ in range(net): | ||
| out[a].append(b) | ||
|
|
||
| # Chain the surviving directed edges into closed rings. |
There was a problem hiding this comment.
🤖 from Claude (review)
Blocking (correctness): ring chaining at pinch / non-manifold vertices is order-dependent and can mis-chain.
out[cur].pop() picks an arbitrary outgoing edge whenever a vertex has out-degree > 1. That happens whenever two parts of the outline meet at a single shared HEALPix corner — e.g. two cell blocks that touch only diagonally, or a checkerboard/“bowtie” cover. At such a pinch vertex there are ≥2 valid continuations and .pop() (LIFO over insertion order) may stitch the wrong pair, producing a self-touching ring or two rings merged into one. Edge-cancellation alone does not guarantee a vertex has degree 1.
The robust fix for a planar straight-line graph is to choose the next edge by angular order at the vertex (turn maximally clockwise/counter-clockwise relative to the incoming edge) rather than arbitrary pop; that resolves pinch points deterministically into the correct face rings.
There is no test exercising a diagonal-touch / pinch cover (the current tests are all simply-connected boxes, one hole, and two well-separated components). Please add one — a cover where two blocks share exactly one corner — and confirm the dissolved outline is two rings, not a bowtie.
Generated by Claude Code
| segments.append(cur) | ||
| cur = [(-boundary, la_x)] | ||
| if segments: | ||
| segments[0] = cur + segments[0] # the wrap-around segment closes the first |
There was a problem hiding this comment.
🤖 from Claude (review)
Blocking (correctness): the wrap-around merge segments[0] = cur + segments[0] only holds for exactly two crossings.
With more than two antimeridian crossings (even count — a perfectly legal non-pole case, e.g. a cover that touches ±180 in two separate latitude bands, or an L-shaped cover straddling the seam), the loop appends one segment per crossing and leaves a trailing cur for the wrap. This line then unconditionally glues the wrap piece onto segments[0] — the first appended segment — but the wrap piece and segments[0] are not generally on the same ±180 side, so the result mixes hemispheres and lons.max()-lons.min() can exceed 180°. The two halves that actually belong together are the wrap cur and the last-by-position open segment on the same side, not index 0.
The current test (test_dissolve_antimeridian_split) only covers the single 2-crossing box, so this path is unverified. Please add a ≥4-crossing case and assert each piece still spans ≤ a hemisphere. (Relatedly, the joints at cur + segments[0] should be deduped if the two endpoints coincide, to avoid a zero-length edge.)
Generated by Claude Code
| best = idx | ||
| hole_groups[best if best is not None else 0].append(hole) | ||
|
|
||
| return [ |
There was a problem hiding this comment.
🤖 from Claude (review)
Correctness: hole-nesting fallback best if best is not None else 0 silently mis-attaches an uncontained hole.
If a hole's centroid lands in no exterior piece (it can: the centroid of a non-convex/antimeridian-affected ring need not lie inside its own region, and the ray-cast is run against the split exterior pieces, not the original ring), the hole is force-attached to ext_pieces[0]. That can attach a hole to a polygon that does not contain it → an invalid geometry, or attach it to the wrong component. Two issues compound it:
- Centroid test is unsafe.
np.meanof the ring vertices is not guaranteed inside a concave ring. Use a guaranteed-interior point (e.g. a point just inside an edge, or test the ring against exteriors by an actual containment-of-ring test) instead of the vertex-mean. - Silent fallback. If genuinely no exterior contains the hole, that's a bug in the dissolve, not something to paper over by dumping it on piece 0 — prefer asserting/raising so it surfaces in tests rather than producing a quietly-wrong polygon.
Also note holes are nested by planar containment after the antimeridian split, but holes themselves are only checked for crossings, not split — see the separate note on _split_at_antimeridian interaction.
Generated by Claude Code
| ext_pieces = [] | ||
| holes = [] | ||
| for ring, area in zip(rings_xyz, areas): | ||
| lat, lon = _xyz_to_latlon(ring) |
There was a problem hiding this comment.
🤖 from Claude (review)
Correctness/robustness: the sum(areas) < 0.0 global winding-normalization is ambiguous for large covers.
Spherical signed area is only defined mod 4π, so sum(areas) distinguishes “covered” from “complement” only while the cover is comfortably less than a hemisphere. For a cover near half the sphere (sum near ±2π) the sign is numerically and topologically ambiguous, and the whole exterior/hole classification flips. The docstring says the net is “always positive,” but that's the area — the signed sum it's tested against isn't, once you're near 2π.
This is fine for the box-sized tests here, but it's an unstated correctness boundary. Please either (a) document the “cover must be well under a hemisphere” precondition and guard it (e.g. raise if abs(sum(areas)) is within a margin of 2π), or (b) derive winding from a robust per-cell orientation invariant rather than the global signed-area sign. Tie-in: this is also why the per-ring area >= 0.0 exterior/hole split can flip under the same condition.
Generated by Claude Code
|
|
||
| # Integer-snap every boundary point to a vertex id; a shared corner/sub-edge | ||
| # point collapses to one id, so adjacent cells reference the same vertex. | ||
| snapped = np.round(flat * _DISSOLVE_SNAP).astype(np.int64) |
There was a problem hiding this comment.
🤖 from Claude (review)
Robustness question: the _DISSOLVE_SNAP = 1e10 integer snap and its interaction with step > 1.
For step == 1 the cancellation is exact-by-construction: both adjacent cells get the shared corner from the same vertices() computation, so the float bits are identical and the snap is a no-op for correctness. Good.
For step > 1, the intermediate sub-edge points are produced by path_along_cell_edge per cell, and the two neighbours traverse the shared edge in opposite directions. Those interpolated points are only approximately equal across the two cells, so cancellation now depends on both rounding to the same round(x*1e10) bucket. That's a real floating tolerance, just spelled as a fixed grid: ~5e-11 in each unit-vector component (~0.3 mm on the sphere). I don't have a concrete failure, but two risks deserve a comment or a test: (1) at high orders the cells are small but the sub-edge spacing is still ≫ the snap, so collisions are unlikely — but please confirm the cross-cell sub-edge points actually agree to <5e-11, otherwise interior sub-edges fail to cancel and you get slivers; (2) np.round to int64 of x*1e10 is fine for unit vectors but assumes inputs are exactly in [-1,1] — they are here, so OK. A step=4 adjacency test that asserts the interior seam fully cancels (no interior vertices in the output ring) would lock this down; the current step test only checks area, which is insensitive to small non-cancelled slivers.
Generated by Claude Code
| """ | ||
| rings_xyz = _boundary_rings_xyz(morton, step) | ||
| if not rings_xyz: | ||
| return [] |
There was a problem hiding this comment.
🤖 from Claude (review)
Pole/odd-crossing detection completeness. Pole-enclosing is inferred from an odd antimeridian-crossing count per ring. This catches a ring that winds around a pole crossing ±180 once, but it is a proxy, not a direct test: a pole-enclosing outline that happens to cross the seam an even number of times (e.g. a polar cap whose boundary dips back and forth across ±180), or one that does not reach the seam at all because the cap is offset, would slip through and then be mis-split as if planar. Given this is an explicit “not yet supported, raise” path, a more direct detector (does the ring enclose ±90° lat — e.g. sum of signed longitude deltas ≈ ±360°, or include-the-pole test) would fail safe instead of silently emitting garbage. At minimum, please note this limitation of the heuristic in the docstring/PR thread, since the consequence of a miss is a silently-wrong polygon rather than a raise.
Generated by Claude Code
| return polys | ||
|
|
||
|
|
||
| # ── emit: dissolved-boundary outline (phase 4) ───────────────────────────── |
There was a problem hiding this comment.
🤖 from Claude (review)
Design deviation (non-blocking for this draft increment, but flag for @espg): the locked issue #71 design called for the dissolve engine to be a Rust helper; this phase implements the entire edge-cancellation + ring-chaining + antimeridian split + hole-nesting in pure Python. Per CLAUDE.md §7, perf-critical core belongs in src_rust/ (mortie._rustie), and the dissolve is O(cells·step) with Python-level loops over every edge and a per-vertex Counter/defaultdict chase — exactly the hot path that motivated the Rust core. This is reasonable as a correctness-first draft, but it's a scope/architecture deviation from the agreed plan, so it needs an explicit @espg decision: land Python now and port to Rust as a follow-up phase, or implement the helper in Rust before merge. Raising rather than deciding (per §4/§6).
Generated by Claude Code
|
🤖 from Claude Phase 4 (dissolved-outline emit) — core landed green; two decisions needed before the pole/antimeridian-split sub-piece. This run implemented the native edge-cancellation dissolve (
Why this is
Until then Generated by Claude Code |
|
rust now is fine (we want this to be fast); you can use python for testing/validation if thats helpful ahead of implementing in rust |
this interacts with above... if we can solve the pole issue with spherely, my lean is to invest in a full spherely ring-walk shim... but again, depends on above
top level name space is fine |
2 participants
Closes #71
Adds WKB/WKT geometry ingest and emit to mortie, per the design locked in the issue thread (see the 05:23 sign-off and the edge-handling sign-off).
Design (settled on the issue thread)
mortie/geometry.py), gated likemortie/arrow.pygates pyarrow —shapely>=2primary,spherelyaccepted as a raw codec.shapely>=2is added to thetestextra only, never[project].dependencies.MultiPolygonis the opt-indissolve=Falsegranularity. Dissolve uses native edge-cancellation / boundary-trace (approach (i)) — no backend spatial predicate.SRID=4326opt-in. EPSG:4326 lon/lat degrees is the documented contract. No PROJ registration.Phases
mortie/geometry.pycodec adapter. Lazyshapely/spherelybackend gate; WKB/WKT decode/encode with EWKB/EWKT(SRID) opt-in;decompose()of Polygon / MultiPolygon / holes / LineString / MultiLineString into(lat, lon)ring arrays (empty geometries rejected, Z dropped). Backend-gated tests.from_wkb/from_wkt/from_geometry→decompose→ existing entry points (morton_coverage,morton_coverage_moc,linestring_coverage). Round-trip tests: WKB and WKT ingest reproduce the array-path cover.dissolve=False).to_wkb/to_wkt/to_geometry: morton word → cell corners (reusingmort2polygon, grouped by order so mixed-order MOC works) → per-cellMultiPolygon; EWKB/EWKT SRID opt-in; empty cover → emptyMultiPolygon.shapely.unary_union(machine-precision symmetric difference over 250 randomized covers) and a backend-independent spherical-area invariant.How it was tested
Local, isolated venv,
maturin develop --release:pytest mortie/tests/test_geometry.py→ 36 passed (codec round-trips, ingest parity, per-cell emit, and the new dissolve suite: union-oracle match, holes, disjoint components, mixed-order MOC, step-densified seam cancellation, antimeridian split with area conservation, corner-touch simple-ring, pole-enclosing raises, WKB/WKT/SRID round-trip, empty cover).pytest -q→ 435 passed, 11 skipped;flake8 mortie --select=E9,F63,F7,F82clean;--max-line-length=88clean on the touched files.geometry.pyis 677 lines (under the ~1000 soft limit). No Rust changed this phase; no new runtime dependency.fold phase 4 self-review on #71(angular ring-chaining for non-manifold corner-touch vertices + a pinch test; rigorous net-longitude-winding pole detection; guard + clear raise for exteriors crossing the antimeridian more than twice; hole nesting now tests a guaranteed-interior hole vertex and raises rather than silently mis-attaching; documented the mod-4π winding assumption; seam-cancellation test forstep>1).Current dissolve scope (well-defined
NotImplementedErroroutside it)dissolve=Truehandles contiguous/disjoint covers, holes, mixed-order MOC,step>1, and antimeridian-crossing covers (split into valid pieces). It raises a clearNotImplementedError(pointing todissolve=False) for the three cases whose spherical→planar split is the remaining sub-piece: a cover whose outline encloses a pole (e.g. a polar cap), an exterior crossing the antimeridian more than twice, and an antimeridian-crossing hole.Questions for review
src_rust/helper the locked plan named. I did this deliberately to de-risk the spherical topology before committing it to Rust, and it's an emit/export path rather than a hot inner loop — but it is a deviation from the plan and from §7's "perf-critical code is Rust", so I'm flagging it for your call: (a) land the Python reference now and port the kernel to Rust in a follow-up, or (b) you'd prefer the Rust port before this merges. I'd like your decision before sinking the pole/antimeridian-split work into either language.shapely.unary_union) breaks exactly there, so it can only be self-verified (spherical-area conservation +is_valid+ ≤180° span). How would you like a polar-cap dissolve represented — explicit ±90° pole vertices in a single split MultiPolygon, or isdissolve=Falsean acceptable answer for polar caps in 1.x? This drives (1).NotImplementedErroron the introspection/emit ops it can't drive, or invest in a full spherely ring-walk shim?from_wkb/from_wkt/to_wkb/to_wktat the top-levelmortie.namespace, or acoverage_-prefixed form?Marked
waitingon questions (1)/(2) — the pole-split sub-piece depends on the Python-vs-Rust call. Phases 1–4 core are green and self-reviewed.