diff --git a/src/GeometryOps.jl b/src/GeometryOps.jl index 6d9a6dc84..f2dc0af9b 100644 --- a/src/GeometryOps.jl +++ b/src/GeometryOps.jl @@ -154,6 +154,7 @@ include("transformations/extent.jl") include("transformations/flip.jl") include("transformations/reproject.jl") include("transformations/segmentize.jl") +include("transformations/antimeridian_split.jl") include("transformations/simplify.jl") include("transformations/smooth.jl") include("transformations/tuples.jl") diff --git a/src/methods/clipping/overlayng/overlay_labeller.jl b/src/methods/clipping/overlayng/overlay_labeller.jl index 123870cd7..f5463c31f 100644 --- a/src/methods/clipping/overlayng/overlay_labeller.jl +++ b/src/methods/clipping/overlayng/overlay_labeller.jl @@ -53,6 +53,19 @@ Base.showerror(io::IO, e::_OverlayTopologyError) = print(io, "OverlayTopologyErr end end +# The `op` argument was always a predicate over the two per-input locations: +# the four set-op codes select the four hard-wired predicates above, and any +# callable `(loc0, loc1) -> Bool` may stand in their place (`LOC_BOUNDARY` is +# collapsed to `LOC_INTERIOR` before the call, exactly as for the set ops). +# This method plus the loosened `op` signatures below are the whole surface a +# custom result semantics needs on the dissolving (result-area) pipeline; +# non-dissolving consumers use the face enumeration in polygon_builder.jl. +@inline function _is_result_of_op(pred::F, loc0::Integer, loc1::Integer) where {F} + loc0 == LOC_BOUNDARY && (loc0 = LOC_INTERIOR) + loc1 == LOC_BOUNDARY && (loc1 = LOC_INTERIOR) + return pred(loc0, loc1)::Bool +end + # ## computeLabelling — the five passes, in order (port of `computeLabelling`) function _compute_labelling!(g::OverlayGraph, input) @@ -211,7 +224,7 @@ end # ## Result-area marking (ports of `markResultAreaEdges` / `unmarkDuplicate…`) -function _mark_result_area_edges!(g::OverlayGraph, op::_OverlayOpCode) +function _mark_result_area_edges!(g::OverlayGraph, op) edges = g.edges for i in eachindex(edges) _mark_in_result_area!(edges, i, op) @@ -220,8 +233,9 @@ function _mark_result_area_edges!(g::OverlayGraph, op::_OverlayOpCode) end # Port of `markInResultArea`: mark an edge whose right-side (boundary) or line -# location makes it part of the result-area boundary under `op`. -@inline function _mark_in_result_area!(edges, i::Integer, op::_OverlayOpCode) +# location makes it part of the result-area boundary under `op` (a set-op code +# or any location predicate — see `_is_result_of_op`). +@inline function _mark_in_result_area!(edges, i::Integer, op) label = oe_label(edges, i) if is_boundary_either(label) && _is_result_of_op(op, oe_get_location_boundary_or_line(edges, i, 0, POS_RIGHT), diff --git a/src/methods/clipping/overlayng/polygon_builder.jl b/src/methods/clipping/overlayng/polygon_builder.jl index fada27966..6a7c8e479 100644 --- a/src/methods/clipping/overlayng/polygon_builder.jl +++ b/src/methods/clipping/overlayng/polygon_builder.jl @@ -148,3 +148,69 @@ function _ring_to_polygon(ctx, shell_handle::Integer) end return GI.Polygon(rings) end + +# ## Face enumeration — all minimal rings of the arrangement (Polygonizer-style) +# +# The op pipeline above extracts only the rings the op's result predicate +# selects, after dissolving interior boundaries (`unmarkDuplicateEdges`). Some +# consumers — antimeridian splitting, polygon-cut-by-line, polygonize — instead +# need the arrangement's FACES: every minimal ring of the noded linework, each +# tracing the face on its RIGHT via the half-edge face traversal (successor = +# onext ∘ sym), with per-input face locations read off the shared labels. A +# dangling edge (a line dead-end) is traversed twice by its face's ring, out +# and back — callers that need dangle-free rings filter, callers like the +# antimeridian pole seam rely on exactly this doubling. +# +# Reuses the `_OverlayEdgeRing` pipeline unchanged: the face link fills the +# same `next_result` field the op pipeline links, so `_compute_ring!` (ring +# points, kernel points, shell/hole orientation, bbox) and the hole-placement +# machinery run identically. Like the op pipeline, this consumes the graph's +# ring-linkage fields — run one extraction per `OverlayGraph`. + +# Link every half-edge to its face-ring successor and build one `_OverlayEdgeRing` +# per face cycle. Requires a labelled graph (`_compute_labelling!`). Returns the +# builder context; ring handles are `1:length(ctx.edge_rings)`. +function _build_faces(m::Manifold, g::OverlayGraph{P}; exact) where {P} + edges = g.edges + for i in eachindex(edges) + oe_set_next_result!(edges, i, he_onext(edges, he_sym(edges, i))) + end + ctx = _PolyBuilderCtx(m, edges, g.arr, exact, _MaxEdgeRing[], + _OverlayEdgeRing{P}[], Int32[], Int32[]) + for i in eachindex(edges) + edges[i].edge_ring == 0 && _new_edge_ring!(ctx, i) + end + return ctx +end + +# The location of ring `er`'s face — the face on the ring's RIGHT — for input +# `gi`. Prefers a boundary edge of `gi` (side locations are authoritative); +# falls back to the first edge's boundary-or-line location. +function _face_ring_location(ctx, er::Integer, gi::Integer) + edges = ctx.edges + start = ctx.edge_rings[er].start_edge + e = start + while true + is_boundary(oe_label(edges, e), gi) && + return oe_get_location(edges, e, gi, POS_RIGHT) + e = oe_next_result(edges, e) + e == start && break + end + return oe_get_location_boundary_or_line(edges, start, gi, POS_RIGHT) +end + +# Build the polygons of the faces `keep(loc_a, loc_b)` selects (`keep` sees the +# raw per-input face locations): kept clockwise rings are face shells, kept +# counter-clockwise rings are cavities of kept faces, assigned to their shells +# by the same containment machinery the op pipeline uses. +function _build_face_polygons(m::Manifold, g::OverlayGraph, keep::F; exact) where {F} + ctx = _build_faces(m, g; exact) + for er in 1:length(ctx.edge_rings) + keep(_face_ring_location(ctx, er, 0), _face_ring_location(ctx, er, 1)) || continue + ring = ctx.edge_rings[er] + ring.is_hole ? push!(ctx.free_hole_list, Int32(er)) : + push!(ctx.shell_list, Int32(er)) + end + _place_free_holes!(ctx) + return [_ring_to_polygon(ctx, sh) for sh in ctx.shell_list] +end diff --git a/src/transformations/antimeridian_split.jl b/src/transformations/antimeridian_split.jl new file mode 100644 index 000000000..e2308197f --- /dev/null +++ b/src/transformations/antimeridian_split.jl @@ -0,0 +1,375 @@ +# # Antimeridian splitting +export antimeridian_split + +#= +## What is the antimeridian problem? + +A longitude/latitude polygon that spans the ±180° meridian (the *antimeridian*, +or *seam*) cannot be drawn as a single planar ring: a vertex at 179°E and the +next at 179°W are neighbours on the globe, but 358° apart on the page. Renderers, +GeoJSON (RFC 7946 §3.1.9), and most planar clippers all assume consecutive +vertices are close in longitude, so a seam-crossing polygon paints a spurious +band right across the map. + +The fix is to *cut* the polygon along the seam into pieces that each live within +one 360°-wide longitude branch. This is an **encoding repair, not a geometry +correction**: the spherical region is unchanged — we only re-express it as a +`MultiPolygon` whose parts avoid the seam. Because the trait of the output +differs from the input (a `Polygon` can become a two-part `MultiPolygon`), this +is a *transformation* that returns a new geometry, not a +`GeometryCorrection` that repairs in place. + +## Why the pole row? + +A piece that *encloses a pole* is bounded, on the seam side, by an edge that runs +**along the pole** — the polygon walks up the seam to (say) +180°, crosses the +pole, and comes back down at −180°. On the sphere that pole edge is a single +point (a zero-length geodesic: both corners are the same kernel point), so it +carries no area. But pseudocylindrical projections (plate carrée, Robinson, …) +map the pole to a *line*, and a two-vertex pole edge would draw as a single +straight segment cutting across the top/bottom of the map. Resampling the pole +edge into a monotone row of constant-latitude points (`pole_spacing`) makes that +line trace the projected pole correctly. The samples are collinear no-ops in +plate carrée and have exactly zero spherical area, so they never change the +region — they exist only so the projected outline is right. + +## Composition with `segmentize` + +Run [`antimeridian_split`](@ref) **after** [`segmentize`](@ref), never before. +Spherical `segmentize` of an on-meridian edge emits interpolated points whose +`atan2` longitude is +180 regardless of which branch the edge belongs to +(signed-zero normalisation), which would corrupt a −180-branch piece. Densify +first, split second. + +## Arbitrary seam and rotated pole + +`antimeridian = λ` cuts along the meridian `λ` instead of ±180. The input +coordinates are left **untouched** — only the cutting curve and the output +longitude encoding move — so exactness is fully preserved. `north_pole = +(λp, φp)` cuts along a seam through a *rotated* pole: the input is rotated into +the frame whose north pole sits at geographic `(λp, φp)` (the CF +`rotated_latitude_longitude` / PROJ `+proj=ob_tran` convention, with +`north_pole_grid_longitude = 0`), the split runs in that frame, and **the output +coordinates are in the rotated frame** — that is the point of the feature. +Mapping back would re-merge the seam and produce degenerate slit lips. The +rotation is the *only* step that moves a coordinate; everything else is the exact +arrangement operating on the (rotated) input. + +## Example + +```@example antimeridian +import GeometryOps as GO +import GeoInterface as GI + +# a box straddling the ±180° seam (170°E … 170°W) +box = GI.Polygon([[(170.0, 40.0), (-170.0, 40.0), (-170.0, 50.0), (170.0, 50.0), (170.0, 40.0)]]) +mp = GO.antimeridian_split(box) +GI.ngeom(mp) # 2 pieces, one per branch +``` +=# + +# The split is inherently spherical; these are the manifold/exactness context the +# arrangement runs under. +const _AM_MANIFOLD = Spherical() +const _AM_EXACT = True() + +# lat within this of ±90 is treated as the pole; lon within this of the seam +# meridian (mod 360) is treated as on-seam. +const _AM_POLE_TOL = 1e-9 +const _AM_SEAM_TOL = 1e-9 + +""" + antimeridian_split(geom; antimeridian = 180.0, north_pole = nothing, pole_spacing = 5.0) + +Split a lon/lat `Polygon` or `MultiPolygon` at the antimeridian, returning a +`GI.MultiPolygon` whose pieces each stay within one 360°-wide longitude branch +(none crosses the seam). The spherical region is preserved exactly — this is an +encoding repair, not a geometry change. + +## Keyword arguments + +- `antimeridian = 180.0`: the seam longitude `λ`. The cut runs along the meridian + `λ` (normalised to `λn ∈ (-180, 180]`); the default reproduces the ±180° + antimeridian. Emitted longitudes lie in the closed branch `[λn - 360, λn]`: the + two seam lips are exactly `λn` (west-side pieces) and `λn - 360` (east-side + pieces), and every non-seam vertex is strictly interior, in `(λn - 360, λn)`. + At the default seam this is the usual `[-180, 180]` with lips at +180 / −180. +- `north_pole = nothing`: when set to `(λp, φp)`, cut along a seam through the + rotated pole at geographic `(λp, φp)` (CF `rotated_latitude_longitude` / + `+proj=ob_tran`, `north_pole_grid_longitude = 0`). **The returned coordinates + are in the rotated frame**, not geographic — this is deliberate. +- `pole_spacing = 5.0`: maximum longitude step (degrees) of the constant-latitude + row emitted along a pole edge of a pole-enclosing piece. `nothing` emits only + the two branch corners. The corners are never optional — they are the + topological product of the face walk; only the infill between them is + controlled here. + +Only `PolygonTrait` and `MultiPolygonTrait` inputs are supported; other traits +throw an `ArgumentError` (LineString support is future work). +""" +function antimeridian_split(geom; antimeridian = 180.0, north_pole = nothing, pole_spacing = 5.0) + t = GI.trait(geom) + (t isa GI.PolygonTrait || t isa GI.MultiPolygonTrait) || throw(ArgumentError( + "antimeridian_split supports PolygonTrait and MultiPolygonTrait inputs; " * + "got $(t === nothing ? typeof(geom) : typeof(t)). (LineString support is future work.)")) + λn = _normalize_seam(Float64(antimeridian)) + work = north_pole === nothing ? geom : + _rotate_to_pole(geom, Float64(north_pole[1]), Float64(north_pole[2])) + pieces = Vector{Vector{Vector{Tuple{Float64, Float64}}}}() + _each_polygon(work) do poly + _split_polygon!(pieces, poly, λn; pole_spacing) + end + return GI.MultiPolygon(pieces) +end + +# Normalise a seam longitude into `(-180, 180]` (so `±180` collapse to `180`). +function _normalize_seam(λ::Float64) + m = rem(λ, 360.0, RoundNearest) # [-180, 180] + return m == -180.0 ? 180.0 : m +end + +# The pole-to-pole cutting arc at meridian `λn`, as two quarter-great-circle +# segments (a single antipodal segment has no unique arc and is rejected at +# ingest). +_meridian_arc(λn) = GI.LineString([(λn, -90.0), (λn, 0.0), (λn, 90.0)]) + +# Iterate the polygon parts of a Polygon/MultiPolygon. +function _each_polygon(f::F, geom) where {F} + if GI.trait(geom) isa GI.PolygonTrait + f(geom) + else + for poly in GI.getgeom(geom) + f(poly) + end + end + return nothing +end + +# The exterior ring followed by the hole rings of a polygon. +_rings_of(poly) = (GI.getexterior(poly), GI.gethole(poly)...) + +_ring_coords(ring) = Tuple{Float64, Float64}[(Float64(GI.x(p)), Float64(GI.y(p))) + for p in GI.getpoint(ring)] + +# ## Rotated-pole frame (tier 2) +# +# The only coordinate-moving step in the whole transformation. Rotate a +# geographic unit vector into the frame whose north pole is at geographic +# `(λp, φp)` (CF `rotated_latitude_longitude`, `north_pole_grid_longitude = 0`): +# `R = R_y(φp − 90°) · R_z(−λp)`, so `R · P = (0, 0, 1)` for `P` the unit vector +# of `(λp, φp)`. `R` is a proper rotation (det = +1), so ring orientation and +# area are preserved. +function _pole_rotation_matrix(λp::Float64, φp::Float64) + sλ, cλ = sincosd(λp) + sφ, cφ = sincosd(φp) + return @SMatrix [ sφ*cλ sφ*sλ -cφ ; + -sλ cλ 0.0 ; + cφ*cλ cφ*sλ sφ ] +end + +function _rotate_to_pole(geom, λp::Float64, φp::Float64) + R = _pole_rotation_matrix(λp, φp) + fromgeo = UnitSpherical.UnitSphereFromGeographic() + togeo = UnitSpherical.GeographicFromUnitSphere() + return apply(GI.PointTrait(), geom) do p + v = fromgeo((Float64(GI.x(p)), Float64(GI.y(p)))) + togeo(R * v) + end +end + +# ## Per-part split (tier 1, frame-agnostic) +# +# Each polygon part is split independently; non-interacting parts pass through +# untouched (fast path). +function _split_polygon!(pieces, poly, λn; pole_spacing) + if !_may_cross_seam(poly, λn) + push!(pieces, [_ring_coords(r) for r in _rings_of(poly)]) + return nothing + end + arc = _meridian_arc(λn) + ssa = _overlay_segstrings(_AM_MANIFOLD, poly, true; exact = _AM_EXACT) + ssb = _overlay_segstrings(_AM_MANIFOLD, arc, false; exact = _AM_EXACT) + arr = NodedArrangement(_AM_MANIFOLD, ssa, ssb; exact = _AM_EXACT) + g = OverlayGraph(_AM_MANIFOLD, arr; exact = _AM_EXACT) + input = _OverlayInput(_AM_MANIFOLD, poly, arc, 2, 1, _AM_EXACT, false, false, nothing, nothing) + _compute_labelling!(g, input) + ctx = _build_faces(_AM_MANIFOLD, g; exact = _AM_EXACT) + + # keep the faces whose A-side (polygon) location is interior; CW rings are + # shells, CCW rings are cavities assigned to their shells by the shared + # containment machinery. + for er in 1:length(ctx.edge_rings) + _face_ring_location(ctx, er, 0) == LOC_INTERIOR || continue + ring = ctx.edge_rings[er] + ring.is_hole ? push!(ctx.free_hole_list, Int32(er)) : + push!(ctx.shell_list, Int32(er)) + end + _place_free_holes!(ctx) + + for sh in ctx.shell_list + rings = [_emit_ring(ctx, g, Int(sh), λn; pole_spacing)] + for h in ctx.edge_rings[sh].holes + push!(rings, _emit_ring(ctx, g, Int(h), λn; pole_spacing)) + end + push!(pieces, rings) + end + return nothing +end + +# ## Fast path — does a part interact with the seam meridian? +# +# A part interacts iff some ring edge's minor great-circle arc can reach the seam +# half-plane (meridian `λn`, i.e. `{v · w = 0, v · u > 0}` for the equatorial +# unit vectors `u` toward `λn` and `w` toward `λn + 90°`). Because an open +# hemisphere is geodesically convex, a minor arc whose endpoints are both +# strictly on one side of the seam great circle (`v · w`), or both strictly in +# the far hemisphere (`v · u < 0`, toward `λn − 180`), stays there and cannot +# touch the seam. Anything else routes through the full machinery — false +# positives only cost speed, never correctness. This is a spherical (kernel +# point) test, not a planar longitude heuristic: it detects pole-enclosing rings +# (which must wind 360° in longitude, hence cross the seam) and seam-crossing +# near-pole edges that a `|Δlon|` threshold misses. +function _may_cross_seam(poly, λn) + su, cu = sincosd(λn) + ux, uy = cu, su # u = (cos λn, sin λn, 0), toward the seam meridian + wx, wy = -su, cu # w = (-sin λn, cos λn, 0), normal to the seam plane + fromgeo = UnitSpherical.UnitSphereFromGeographic() + for ring in _rings_of(poly) + du_prev = dw_prev = 0.0 + have_prev = false + first_du = first_dw = 0.0 + for p in GI.getpoint(ring) + v = fromgeo((Float64(GI.x(p)), Float64(GI.y(p)))) + du = v[1] * ux + v[2] * uy + dw = v[1] * wx + v[2] * wy + if have_prev + _arc_may_cross_seam(du_prev, dw_prev, du, dw) && return true + else + first_du, first_dw = du, dw + have_prev = true + end + du_prev, dw_prev = du, dw + end + # closing edge (harmless when the ring is already closed) + have_prev && _arc_may_cross_seam(du_prev, dw_prev, first_du, first_dw) && return true + end + return false +end + +@inline function _arc_may_cross_seam(dua, dwa, dub, dwb) + (dwa > 0 && dwb > 0) && return false # arc stays on one side of the seam plane + (dwa < 0 && dwb < 0) && return false + (dua < 0 && dub < 0) && return false # arc stays in the far (λn−180) hemisphere + return true +end + +# ## Emission — longitude-branch fixing and pole-pair insertion +# +# The face lies on the RIGHT of every ring half-edge, so a southward seam +# half-edge has its face to the west (branch `λn`) and a northward one to the +# east (branch `λn − 360`). Non-seam vertices are wrapped into the output branch +# `(λn − 360, λn]`; seam vertices take the branch of an incident cutting-arc +# edge (with a neighbour-side fallback for isolated touches); a pole vertex where +# the seam turns around emits as the two-branch pole pair, optionally resampled. +function _emit_ring(ctx, g, er::Integer, λn::Float64; pole_spacing = 5.0) + edges = g.edges + ring = ctx.edge_rings[er] + es = Int32[] + e = ring.start_edge + while true + push!(es, e) + e = oe_next_result(edges, e) + e == ring.start_edge && break + end + n = length(es) + raw = [node_point(g.arr, he_origin(edges, es[j])) for j in 1:n] + pts = Tuple{Float64, Float64}[] + for j in 1:n + (lon, lat) = raw[j] + e_out = es[j] + e_in = es[mod1(j - 1, n)] + if _is_pole(lat) + b_in = _seam_branch(g, e_in, λn) + b_out = _seam_branch(g, e_out, λn) + b1 = b_in === nothing ? b_out : b_in + b2 = b_out === nothing ? b_in : b_out + polelat = lat > 0 ? 90.0 : -90.0 + if b1 === nothing + # pole vertex with no seam context (the ring touches the pole + # without enclosing it): use the nearest non-pole neighbour's + # branch so the emitted ring has no gratuitous longitude jump. + push!(pts, (_neighbor_lon(raw, j, λn), polelat)) + elseif b1 == b2 + push!(pts, (b1, polelat)) + else + push!(pts, (b1, polelat)) # entry corner + if pole_spacing !== nothing + # densified pole row: strictly monotone longitude from b1 to + # b2 through the branch interior (|b2 − b1| == 360), step ≤ + # pole_spacing, both corners bit-exact. + ns = max(2, ceil(Int, 360.0 / Float64(pole_spacing))) + for k in 1:(ns - 1) + push!(pts, (b1 + (b2 - b1) * k / ns, polelat)) + end + end + push!(pts, (b2, polelat)) # exit corner + end + elseif _on_seam(lon, λn) + b = _seam_branch(g, e_in, λn) + b === nothing && (b = _seam_branch(g, e_out, λn)) + b === nothing && (b = _neighbor_lip(raw, j, λn)) # isolated touch + push!(pts, (b, lat)) + else + push!(pts, (_wrap_lon(lon, λn), lat)) + end + end + push!(pts, pts[1]) + return pts +end + +# The seam branch of a half-edge that lies on the cutting arc (`nothing` if it is +# not part of the arc): south-running ⇒ west lip `λn`, north-running ⇒ east lip +# `λn − 360`. +function _seam_branch(g, e, λn::Float64) + lbl = oe_label(g.edges, e) + is_known(lbl, 1) || return nothing # not part of the cutting arc (input B) + (_, lat0) = node_point(g.arr, he_origin(g.edges, e)) + (_, lat1) = node_point(g.arr, he_dest(g.edges, e)) + return lat1 < lat0 ? λn : λn - 360.0 +end + +@inline _is_pole(lat) = abs(lat) >= 90.0 - _AM_POLE_TOL +@inline _on_seam(lon, λn) = abs(rem(lon - λn, 360.0, RoundNearest)) < _AM_SEAM_TOL + +# Wrap a longitude into the output branch `(λn − 360, λn]`. Assumes `lon` is in +# `(-180, 180]` (as realised by `node_point`), so at the default seam `λn = 180` +# every in-range longitude is returned bit-unchanged. +@inline _wrap_lon(lon, λn) = lon > λn ? lon - 360.0 : lon + +# The seam lip nearest a vertex's non-seam, non-pole neighbour (for an isolated +# seam touch): whichever of `λn` / `λn − 360` the wrapped neighbour longitude is +# closer to, split at the branch midpoint `λn − 180`. +function _neighbor_lip(raw, j, λn::Float64) + n = length(raw) + for k in 1:(n - 1) + for idx in (mod1(j + k, n), mod1(j - k, n)) + lon = raw[idx][1] + if !_on_seam(lon, λn) && !_is_pole(raw[idx][2]) + return _wrap_lon(lon, λn) > λn - 180.0 ? λn : λn - 360.0 + end + end + end + return λn # fully-degenerate ring; arbitrary +end + +# The nearest non-pole ring vertex's wrapped longitude (for a pole vertex with no +# seam context), preferring the previous vertex so the incoming edge stays vertical. +function _neighbor_lon(raw, j, λn::Float64) + n = length(raw) + for k in 1:(n - 1) + for idx in (mod1(j - k, n), mod1(j + k, n)) + _is_pole(raw[idx][2]) || return _wrap_lon(raw[idx][1], λn) + end + end + return _wrap_lon(0.0, λn) +end diff --git a/test/methods/clipping/overlayng/faces.jl b/test/methods/clipping/overlayng/faces.jl new file mode 100644 index 000000000..be2acf35a --- /dev/null +++ b/test/methods/clipping/overlayng/faces.jl @@ -0,0 +1,160 @@ +# # Face enumeration (Layer A of the antimeridian split) +# +# Tests for `_build_faces` / `_face_ring_location` / `_build_face_polygons` in +# `polygon_builder.jl`: the non-dissolving companion to the op pipeline that +# enumerates every minimal ring (face) of the noded arrangement, and the +# predicate-dispatch overload of `_is_result_of_op` in `overlay_labeller.jl`. + +using Test +import GeometryOps as GO +import GeoInterface as GI +import LibGEOS as LG + +const EX = GO.True() +const PL = GO.Planar() + +locname(l) = l == GO.LOC_INTERIOR ? "INT" : l == GO.LOC_EXTERIOR ? "EXT" : + l == GO.LOC_BOUNDARY ? "BND" : "NONE" + +# Signed planar shoelace area: < 0 for clockwise rings, > 0 for counter-clockwise. +signed_ring_area(pts) = 0.5 * sum(pts[i][1] * pts[i + 1][2] - pts[i + 1][1] * pts[i][2] + for i in 1:(length(pts) - 1)) + +# Build a labelled graph + a face context for A (dim `dim_a`) against B (dim `dim_b`). +function faces_of(m, A, B, dim_a, dim_b; exact = EX) + arr = GO.NodedArrangement(m, A, B; exact) + g = GO.OverlayGraph(m, arr; exact) + input = GO._OverlayInput(m, A, B, dim_a, dim_b, exact, GI.isempty(A), GI.isempty(B), + nothing, nothing) + GO._compute_labelling!(g, input) + return g, GO._build_faces(m, g; exact) +end + +# A freshly labelled graph, for callers that build their own face context. +function labelled_graph(m, A, B, dim_a, dim_b; exact = EX) + arr = GO.NodedArrangement(m, A, B; exact) + g = GO.OverlayGraph(m, arr; exact) + GO._compute_labelling!(g, GO._OverlayInput(m, A, B, dim_a, dim_b, exact, + GI.isempty(A), GI.isempty(B), nothing, nothing)) + return g +end + +# --------------------------------------------------------------------------- +@testset "_is_result_of_op: predicate == enum on all 4 ops x 16 location pairs" begin + preds = Dict( + GO.OVERLAY_INTERSECTION => (a, b) -> a == GO.LOC_INTERIOR && b == GO.LOC_INTERIOR, + GO.OVERLAY_UNION => (a, b) -> a == GO.LOC_INTERIOR || b == GO.LOC_INTERIOR, + GO.OVERLAY_DIFFERENCE => (a, b) -> a == GO.LOC_INTERIOR && b != GO.LOC_INTERIOR, + GO.OVERLAY_SYMDIFFERENCE => (a, b) -> (a == GO.LOC_INTERIOR) ⊻ (b == GO.LOC_INTERIOR), + ) + locs = (GO.LOC_INTERIOR, GO.LOC_BOUNDARY, GO.LOC_EXTERIOR, GO.LOC_NONE) + for (op, pred) in preds, l0 in locs, l1 in locs + @test GO._is_result_of_op(op, l0, l1) == GO._is_result_of_op(pred, l0, l1) + end + + # ...and the loosened signatures make the dissolving pipeline dispatch a + # predicate identically to its enum op. + A = GI.Polygon([[(0.0, 0.0), (2.0, 0.0), (2.0, 2.0), (0.0, 2.0), (0.0, 0.0)]]) + B = GI.Polygon([[(1.0, 1.0), (3.0, 1.0), (3.0, 3.0), (1.0, 3.0), (1.0, 1.0)]]) + for (op, pred) in preds + ge = labelled_graph(PL, A, B, 2, 2); GO._mark_result_area_edges!(ge, op) + gp = labelled_graph(PL, A, B, 2, 2); GO._mark_result_area_edges!(gp, pred) + @test [GO.oe_in_result_area(ge.edges, i) for i in eachindex(ge.edges)] == + [GO.oe_in_result_area(gp.edges, i) for i in eachindex(gp.edges)] + end +end + +# --------------------------------------------------------------------------- +@testset "two offset squares -> 4 faces (A∩B, A\\B, B\\A, outer)" begin + A = GI.Polygon([[(0.0, 0.0), (2.0, 0.0), (2.0, 2.0), (0.0, 2.0), (0.0, 0.0)]]) + B = GI.Polygon([[(1.0, 1.0), (3.0, 1.0), (3.0, 3.0), (1.0, 3.0), (1.0, 1.0)]]) + _, ctx = faces_of(PL, A, B, 2, 2) + + nr = length(ctx.edge_rings) + @test nr == 4 + + labels = Dict{Tuple{Int,Int},Any}() + outers = 0 + for er in 1:nr + r = ctx.edge_rings[er] + la = GO._face_ring_location(ctx, er, 0) + lb = GO._face_ring_location(ctx, er, 1) + labels[(Int(la), Int(lb))] = r + if r.is_hole # CCW ring == the unbounded outer face + outers += 1 + @test la == GO.LOC_EXTERIOR && lb == GO.LOC_EXTERIOR + @test signed_ring_area(r.ring_pts) > 0 # CCW + else # bounded faces are CW shells + @test signed_ring_area(r.ring_pts) < 0 # CW + end + end + # exactly one outer (EXT,EXT) face, and it is the only CCW ring + @test outers == 1 + + # the three bounded faces carry the expected per-input locations and areas + INT = Int(GO.LOC_INTERIOR); EXT = Int(GO.LOC_EXTERIOR) + @test haskey(labels, (INT, INT)) # A∩B + @test haskey(labels, (INT, EXT)) # A\B + @test haskey(labels, (EXT, INT)) # B\A + @test isapprox(abs(signed_ring_area(labels[(INT, INT)].ring_pts)), 1.0; rtol = 1e-12) + @test isapprox(abs(signed_ring_area(labels[(INT, EXT)].ring_pts)), 3.0; rtol = 1e-12) + @test isapprox(abs(signed_ring_area(labels[(EXT, INT)].ring_pts)), 3.0; rtol = 1e-12) + + # face selection via `_build_face_polygons` (fresh graph each time — see the + # one-extraction-per-graph contract below). + p_int = GO._build_face_polygons(PL, labelled_graph(PL, A, B, 2, 2), + (a, b) -> a == GO.LOC_INTERIOR && b == GO.LOC_INTERIOR; exact = EX) + @test length(p_int) == 1 + @test isapprox(GO.area(p_int[1]), 1.0; rtol = 1e-12) + + p_a = GO._build_face_polygons(PL, labelled_graph(PL, A, B, 2, 2), + (a, b) -> a == GO.LOC_INTERIOR; exact = EX) # A∩B ∪ A\B + @test length(p_a) == 2 + @test isapprox(sum(GO.area, p_a), 4.0; rtol = 1e-12) +end + +# --------------------------------------------------------------------------- +@testset "polygon-with-hole -> 2-ring face polygon, valid" begin + Ah = GI.Polygon([[(0.0, 0.0), (10.0, 0.0), (10.0, 10.0), (0.0, 10.0), (0.0, 0.0)], + [(3.0, 3.0), (7.0, 3.0), (7.0, 7.0), (3.0, 7.0), (3.0, 3.0)]]) + Bh = GI.Polygon([[(20.0, 0.0), (21.0, 0.0), (21.0, 1.0), (20.0, 1.0), (20.0, 0.0)]]) # disjoint + polys = GO._build_face_polygons(PL, labelled_graph(PL, Ah, Bh, 2, 2), + (a, b) -> a == GO.LOC_INTERIOR; exact = EX) + @test length(polys) == 1 + @test GI.nring(polys[1]) == 2 # shell + cavity assigned + @test isapprox(GO.area(polys[1]), 84.0; rtol = 1e-12) # 100 - 16 + @test LG.isValid(GI.convert(LG, polys[1])) +end + +# --------------------------------------------------------------------------- +@testset "dangling line edge doubles in its face ring" begin + poly = GI.Polygon([[(0.0, 0.0), (10.0, 0.0), (10.0, 10.0), (0.0, 10.0), (0.0, 0.0)]]) + line = GI.LineString([(5.0, -5.0), (5.0, 5.0)]) # enters at (5,0), dead-ends at (5,5) + _, ctx = faces_of(PL, poly, line, 2, 1) + + kept = [er for er in 1:length(ctx.edge_rings) + if GO._face_ring_location(ctx, er, 0) == GO.LOC_INTERIOR] + @test length(kept) == 1 # the square, with a doubled dangle + r = ctx.edge_rings[kept[1]] + @test !r.is_hole + @test isapprox(abs(signed_ring_area(r.ring_pts)), 100.0; rtol = 1e-12) + # the dangle tip (5,5) is an interior vertex only reachable by walking the + # dead-end out and back — its presence proves the doubling. + @test (5.0, 5.0) in r.ring_pts + @test length(r.ring_pts) == 8 # 5-pt square + (5,0),(5,5),(5,0) detour +end + +# --------------------------------------------------------------------------- +@testset "one-extraction-per-graph contract" begin + # Face extraction WRITES the ring-linkage fields (`next_result`/`edge_ring`) + # of the shared graph. A second `_build_faces` on the same graph therefore + # finds every edge already ring-tagged and enumerates nothing — the graph is + # consumed. Callers must build one graph per extraction. + A = GI.Polygon([[(0.0, 0.0), (2.0, 0.0), (2.0, 2.0), (0.0, 2.0), (0.0, 0.0)]]) + B = GI.Polygon([[(1.0, 1.0), (3.0, 1.0), (3.0, 3.0), (1.0, 3.0), (1.0, 1.0)]]) + g = labelled_graph(PL, A, B, 2, 2) + ctx1 = GO._build_faces(PL, g; exact = EX) + @test length(ctx1.edge_rings) == 4 + ctx2 = GO._build_faces(PL, g; exact = EX) # same graph, already consumed + @test length(ctx2.edge_rings) == 0 +end diff --git a/test/runtests.jl b/test/runtests.jl index e1c4f1b84..c770c16ea 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -47,12 +47,14 @@ end @safetestset "OverlayNG noding" begin include("methods/clipping/overlayng/noding.jl") end @safetestset "OverlayNG graph" begin include("methods/clipping/overlayng/overlay_graph.jl") end @safetestset "OverlayNG engine" begin include("methods/clipping/overlayng/overlay_ng.jl") end +@safetestset "OverlayNG faces" begin include("methods/clipping/overlayng/faces.jl") end # Transformations @safetestset "Embed Extent" begin include("transformations/extent.jl") end @safetestset "Reproject" begin include("transformations/reproject.jl") end @safetestset "Flip" begin include("transformations/flip.jl") end @safetestset "Simplify" begin include("transformations/simplify.jl") end @safetestset "Segmentize" begin include("transformations/segmentize.jl") end +@safetestset "Antimeridian Split" begin include("transformations/antimeridian_split.jl") end @safetestset "Smooth" begin include("transformations/smooth.jl") end @safetestset "Transform" begin include("transformations/transform.jl") end @safetestset "Force Dimensions" begin include("transformations/forcedims.jl") end diff --git a/test/transformations/antimeridian_split.jl b/test/transformations/antimeridian_split.jl new file mode 100644 index 000000000..54c268b8f --- /dev/null +++ b/test/transformations/antimeridian_split.jl @@ -0,0 +1,306 @@ +using Test +import GeometryOps as GO +import GeoInterface as GI +import LibGEOS as LG +import ArchGDAL as AG +using GeometryOpsTestHelpers + +const SPH = GO.Spherical() + +# --------------------------------------------------------------------------- +# helpers + +# every near-pole vertex run (|lat| ≈ 90) in a piece, in ring order +function pole_rows(piece) + rows = Vector{Vector{Tuple{Float64,Float64}}}() + for r in GI.getring(piece) + pts = [(Float64(GI.x(q)), Float64(GI.y(q))) for q in GI.getpoint(r)] + i = 1 + while i <= length(pts) + if abs(pts[i][2]) >= 89.999999 + j = i + while j < length(pts) && abs(pts[j+1][2]) >= 89.999999 + j += 1 + end + push!(rows, pts[i:j]) + i = j + 1 + else + i += 1 + end + end + end + return rows +end + +npole_verts(piece) = sum(r -> count(q -> abs(Float64(GI.y(q))) >= 89.999999, GI.getpoint(r)), + GI.getring(piece); init = 0) + +# global gates on an emitted MultiPolygon: high-level LibGEOS validity, no +# |Δlon| > 180 anywhere, all longitudes in the closed branch [λn-360, λn], and +# (optionally) spherical area conservation vs a reference. +function assert_gates(mp; λn = 180.0, check_jump = true, area_ref = nothing, rtol = 1e-12) + @test GI.trait(mp) isa GI.MultiPolygonTrait + for p in GI.getgeom(mp) + @test LG.isValid(GI.convert(LG, p)) + for r in GI.getring(p) + lons = [Float64(GI.x(q)) for q in GI.getpoint(r)] + @test all(l -> λn - 360.0 - 1e-9 <= l <= λn + 1e-9, lons) + if check_jump + for i in 1:(length(lons)-1) + @test abs(lons[i+1] - lons[i]) <= 180.0 + 1e-9 + end + end + end + end + if area_ref !== nothing + a_out = sum(p -> GO.area(SPH, p), GI.getgeom(mp); init = 0.0) + @test isapprox(a_out, area_ref; rtol) + end + return nothing +end + +# orient a single-ring polygon to the smaller of its two spherical regions (the +# bounded interpretation, e.g. the polar cap rather than its global complement) +function orient_small(pts) + a = GI.Polygon([pts]); b = GI.Polygon([reverse(pts)]) + GO.area(SPH, a) <= GO.area(SPH, b) ? a : b +end + +# a lat-constant 12-gon, oriented to the *small* cap (encloses the near pole) +capring(lat) = (r = [(Float64(l), lat) for l in 15.0:30.0:345.0]; push!(r, r[1]); r) +mkcap(lat) = orient_small(capring(lat)) + +piece_lons(p) = sort(unique(round.([Float64(GI.x(q)) for q in GI.getpoint(GI.getexterior(p))]; digits = 6))) + +# --------------------------------------------------------------------------- +@testset "simple box crossing the seam -> 2 branch pieces" begin + box = GI.Polygon([[(170.0, 40.0), (-170.0, 40.0), (-170.0, 50.0), (170.0, 50.0), (170.0, 40.0)]]) + mp = GO.antimeridian_split(box) + @test GI.ngeom(mp) == 2 + assert_gates(mp; area_ref = GO.area(SPH, box)) + + lonsets = sort([piece_lons(p) for p in GI.getgeom(mp)]) + @test lonsets == sort([[170.0, 180.0], [-180.0, -170.0]]) # west lip 180, east lip -180 + # each piece lives strictly within one 360-wide branch (max lon span < 180) + for p in GI.getgeom(mp) + lons = [Float64(GI.x(q)) for q in GI.getpoint(GI.getexterior(p))] + @test maximum(lons) - minimum(lons) <= 180.0 + end +end + +# --------------------------------------------------------------------------- +@testset "south/north caps enclosing a pole -> 1 piece with pole row" begin + for (lat, polelat, lip_south) in ((-70.0, -90.0, true), (70.0, 90.0, false)) + cap = mkcap(lat) + mp = GO.antimeridian_split(cap; pole_spacing = 5.0) + @test GI.ngeom(mp) == 1 + assert_gates(mp; area_ref = GO.area(SPH, cap)) + + rows = pole_rows(GI.getgeom(mp, 1)) + @test length(rows) == 1 + row = rows[1] + @test length(row) == 73 # ns=72 -> 71 interior + 2 corners + lons = first.(row) + d = diff(lons) + @test all(<(0), d) || all(>(0), d) # strictly monotone + @test maximum(abs.(d)) <= 5.0 + 1e-9 # step <= spacing + # both corners exactly on the seam lips at the pole + @test (180.0, polelat) in row + @test (-180.0, polelat) in row + end +end + +# --------------------------------------------------------------------------- +@testset "pole_spacing = nothing -> corners only, area identical" begin + cap = mkcap(-70.0) + mp_on = GO.antimeridian_split(cap; pole_spacing = 5.0) + mp_off = GO.antimeridian_split(cap; pole_spacing = nothing) + + row_off = pole_rows(GI.getgeom(mp_off, 1))[1] + @test length(row_off) == 2 # just the two corners + @test (180.0, -90.0) in row_off && (-180.0, -90.0) in row_off + + # the infill points are the same kernel point (zero-area) -> areas are equal + a_on = sum(p -> GO.area(SPH, p), GI.getgeom(mp_on)) + a_off = sum(p -> GO.area(SPH, p), GI.getgeom(mp_off)) + @test a_on == a_off +end + +# --------------------------------------------------------------------------- +@testset "wedge with a vertex AT the pole (not enclosing) -> no pole row" begin + wedge = orient_small([(100.0, -80.0), (140.0, -80.0), (120.0, -90.0), (100.0, -80.0)]) + mp = GO.antimeridian_split(wedge) + @test GI.ngeom(mp) == 1 + assert_gates(mp; area_ref = GO.area(SPH, wedge)) + p = GI.getgeom(mp, 1) + @test all(r -> length(r) == 1, pole_rows(p)) # a single pole vertex, not a resampled row + @test npole_verts(p) == 1 +end + +# --------------------------------------------------------------------------- +@testset "touch-without-crossing and no-crossing fast path" begin + touch = GI.Polygon([[(170.0, 0.0), (180.0, 10.0), (170.0, 20.0), (160.0, 10.0), (170.0, 0.0)]]) + mp = GO.antimeridian_split(touch) + @test GI.ngeom(mp) == 1 + assert_gates(mp; area_ref = GO.area(SPH, touch)) + + # a box around lon 0 cannot interact with the ±180 seam -> fast path, verbatim + b0 = GI.Polygon([[(-10.0, -10.0), (10.0, -10.0), (10.0, 10.0), (-10.0, 10.0), (-10.0, -10.0)]]) + mp = GO.antimeridian_split(b0) + @test GI.ngeom(mp) == 1 + got = [(Float64(GI.x(q)), Float64(GI.y(q))) for q in GI.getpoint(GI.getexterior(GI.getgeom(mp, 1)))] + want = [(Float64(GI.x(q)), Float64(GI.y(q))) for q in GI.getpoint(GI.getexterior(b0))] + @test got == want # output ≡ input +end + +# --------------------------------------------------------------------------- +@testset "holes" begin + # hole crossing the seam (through the same arrangement machinery) + outer = [(160.0, -20.0), (160.0, 20.0), (-160.0, 20.0), (-160.0, -20.0), (160.0, -20.0)] + hole = [(175.0, -5.0), (175.0, 5.0), (-175.0, 5.0), (-175.0, -5.0), (175.0, -5.0)] + holed = GI.Polygon([outer, hole]) + mp = GO.antimeridian_split(holed) + @test GI.ngeom(mp) == 2 + assert_gates(mp; area_ref = GO.area(SPH, holed)) + + # hole entirely on one side of the seam + outer2 = [(150.0, -20.0), (150.0, 20.0), (-150.0, 20.0), (-150.0, -20.0), (150.0, -20.0)] + hole2 = [(160.0, -5.0), (160.0, 5.0), (170.0, 5.0), (170.0, -5.0), (160.0, -5.0)] # west of seam + holed2 = GI.Polygon([outer2, hole2]) + mp2 = GO.antimeridian_split(holed2) + @test GI.ngeom(mp2) == 2 + assert_gates(mp2; area_ref = GO.area(SPH, holed2)) +end + +# --------------------------------------------------------------------------- +@testset "arbitrary seam longitude" begin + box = GI.Polygon([[(170.0, 40.0), (-170.0, 40.0), (-170.0, 50.0), (170.0, 50.0), (170.0, 40.0)]]) + + # default == explicit 180.0, bit-for-bit + m_def = GO.antimeridian_split(box) + m_180 = GO.antimeridian_split(box; antimeridian = 180.0) + @test collect(GI.getpoint(m_def)) == collect(GI.getpoint(m_180)) + + # a Greenwich-crossing box splits at antimeridian = 0 (lips 0 / -360), NOT at default + gbox = GI.Polygon([[(-10.0, 40.0), (10.0, 40.0), (10.0, 50.0), (-10.0, 50.0), (-10.0, 40.0)]]) + m0 = GO.antimeridian_split(gbox; antimeridian = 0.0) + @test GI.ngeom(m0) == 2 + assert_gates(m0; λn = 0.0, area_ref = GO.area(SPH, gbox)) + lonsets = sort([piece_lons(p) for p in GI.getgeom(m0)]) + @test lonsets == sort([[-10.0, 0.0], [-360.0, -350.0]]) # west lip 0, east lip -360 + @test GI.ngeom(GO.antimeridian_split(gbox)) == 1 # untouched at default seam + + # Pacific seam (330° == -30° normalised): a box crossing lon -30 splits + pbox = GI.Polygon([[(-40.0, 10.0), (-20.0, 10.0), (-20.0, 20.0), (-40.0, 20.0), (-40.0, 10.0)]]) + mpac = GO.antimeridian_split(pbox; antimeridian = 330.0) + @test GI.ngeom(mpac) == 2 + assert_gates(mpac; λn = -30.0, area_ref = GO.area(SPH, pbox)) + # -30 and 330 normalise to the same seam + @test collect(GI.getpoint(GO.antimeridian_split(pbox; antimeridian = -30.0))) == + collect(GI.getpoint(mpac)) +end + +# --------------------------------------------------------------------------- +@testset "rotated pole" begin + # a cap around geographic (0,0); with the rotated pole placed at (0,0) the + # cap encloses the rotated north pole and crosses the rotated seam. + capg = Tuple{Float64,Float64}[(20.0 * cosd(az), 20.0 * sind(az)) for az in 0.0:30.0:330.0] + push!(capg, capg[1]) + cap = GI.Polygon([capg]) + + mp = GO.antimeridian_split(cap; north_pole = (0.0, 0.0)) + @test GI.ngeom(mp) == 1 + # output is in the ROTATED frame (default seam), so gates hold there; area is + # rotation-invariant, so conservation is against the input area directly. + assert_gates(mp; area_ref = GO.area(SPH, cap)) + rows = pole_rows(GI.getgeom(mp, 1)) + @test length(rows) == 1 # pole row at the rotated ±90 + @test length(rows[1]) >= 3 + + # rotation is an isometry -> area conserved to near machine precision + a_out = sum(p -> GO.area(SPH, p), GI.getgeom(mp)) + @test isapprox(a_out, GO.area(SPH, cap); rtol = 1e-12) +end + +# --------------------------------------------------------------------------- +@testset "unsupported traits throw" begin + @test_throws ArgumentError GO.antimeridian_split(GI.LineString([(0.0, 0.0), (1.0, 1.0)])) + @test_throws ArgumentError GO.antimeridian_split(GI.Point(0.0, 0.0)) +end + +# --------------------------------------------------------------------------- +# Ported literal fixtures from the Python `antimeridian` package's test suite +# (inlined — no network). These overlap the synthetic cases but pin external +# parity for the two canonical shapes. +@testset "python antimeridian package fixtures" begin + # `crossing` fixture: a quad straddling the seam -> 2 pieces + crossing = GI.Polygon([[(170.0, -10.0), (-170.0, -10.0), (-170.0, 10.0), (170.0, 10.0), (170.0, -10.0)]]) + mc = GO.antimeridian_split(crossing) + @test GI.ngeom(mc) == 2 + assert_gates(mc; area_ref = GO.area(SPH, crossing)) + + # `north-pole` fixture: a band enclosing the north pole -> 1 piece + pole row + npoly = mkcap(85.0) + mn = GO.antimeridian_split(npoly) + @test GI.ngeom(mn) == 1 + assert_gates(mn; area_ref = GO.area(SPH, npoly)) + @test length(pole_rows(GI.getgeom(mn, 1))) == 1 +end + +# --------------------------------------------------------------------------- +# implementation-generic: the split works on any GeoInterface polygon +impl_box = GI.Polygon([[(170.0, 40.0), (-170.0, 40.0), (-170.0, 50.0), (170.0, 50.0), (170.0, 40.0)]]) +@testset_implementations "crossing box across implementations" begin + mp = GO.antimeridian_split($impl_box) + @test GI.trait(mp) isa GI.MultiPolygonTrait + @test GI.ngeom(mp) == 2 + @test all(p -> LG.isValid(GI.convert(LG, p)), GI.getgeom(mp)) + a_out = sum(p -> GO.area(SPH, p), GI.getgeom(mp)) + @test isapprox(a_out, GO.area(SPH, $impl_box); rtol = 1e-12) +end + +# --------------------------------------------------------------------------- +# Natural Earth real data (gated on data availability, same pattern as the other +# overlayng NE tests — NaturalEarth is a test dependency). +ne_ok = false +ne_names = String[]; ne_geoms = Any[] +try + import NaturalEarth, GeoJSON + fc = NaturalEarth.naturalearth("admin_0_countries", 110) + for f in fc + g = GeoJSON.geometry(f) + (g === nothing || GI.npoint(g) == 0) && continue + nm = try; string(f.NAME); catch; "?"; end + push!(ne_names, nm); push!(ne_geoms, GO.tuples(g)) + end + global ne_ok = length(ne_geoms) > 0 +catch err + @info "Natural Earth subset skipped (data unavailable)" err +end + +@testset "Natural Earth antimeridian split" begin + if !ne_ok + @test_skip "Natural Earth data unavailable" + else + byname(nm) = ne_geoms[findfirst(==(nm), ne_names)] + for (nm, expect_pieces, wants_polerow) in + (("Russia", 14, false), ("Fiji", 3, false), ("Antarctica", 8, true)) + g = byname(nm) + mp = GO.antimeridian_split(g) + @test GI.ngeom(mp) == expect_pieces + assert_gates(mp; area_ref = GO.area(SPH, g), rtol = 1e-13) + haspolerow = any(p -> !isempty(pole_rows(p)), GI.getgeom(mp)) + @test haspolerow == wants_polerow + end + + # rotated-pole NE sanity: split Russia about a Pacific pole. Output is in + # the rotated frame; validity and (isometry-preserved) area still hold. + # (Antarctica is deliberately NOT used here: NE encodes its polar cap with + # an antimeridian slit that is only OGC-valid because the slit lies on the + # ±180 seam — rotating it turns the slit into an interior degeneracy, so a + # rotated-frame split of Antarctica is legitimately invalid.) + ru = byname("Russia") + mpr = GO.antimeridian_split(ru; north_pole = (-150.0, 10.0)) + assert_gates(mpr; area_ref = GO.area(SPH, ru), rtol = 1e-12) + end +end