From cbe072d15c513ba37e47c6ad775f7349e20cb026 Mon Sep 17 00:00:00 2001 From: Greg Lucas Date: Thu, 2 Apr 2026 21:20:39 -0600 Subject: [PATCH] FIX: boundary attachment arc direction Make the boundary attachment method more robust by iterating through and checking the forward/backward directions for the shortest path. Making sure to check closing upon the ring itself and iterate through other line attachment options. --- lib/cartopy/crs.py | 395 ++++++++++++++---------------- lib/cartopy/tests/test_polygon.py | 51 ++++ 2 files changed, 229 insertions(+), 217 deletions(-) diff --git a/lib/cartopy/crs.py b/lib/cartopy/crs.py index 500057ade..f7231e82a 100644 --- a/lib/cartopy/crs.py +++ b/lib/cartopy/crs.py @@ -10,6 +10,7 @@ """ from abc import ABCMeta +import bisect from collections import OrderedDict from functools import lru_cache import io @@ -988,198 +989,186 @@ def _project_polygon(self, polygon, src_crs): def _attach_lines_to_boundary(self, multi_line_strings, is_ccw): """ - Return a list of LinearRings by attaching the ends of the given lines - to the boundary, paying attention to the traversal directions of the - lines and boundary. + Return a list of LinearRings by closing projected line segments with + boundary arcs. - """ - debug = False - debug_plot_edges = False + Each projected segment has both its start and end point lying on the + projection boundary. The algorithm connects them into closed rings by + inserting boundary arcs between consecutive endpoints. - # Accumulate all the boundary and segment end points, along with - # their distance along the boundary. - edge_things = [] + How the forward walk and shorter arc work together: + Segments are sorted by their start-point distance along the boundary. + To close a ring, we walk *forward* (in the boundary's parameterisation + direction) from each segment's end to find the next segment whose start + arrives first, or the ring's own start if it comes around sooner. + This forward walk determines which segment is next and whether to + close. - # Get the boundary as a LineString of the correct orientation - # so we can compute distances along it. - if is_ccw: - boundary = self.ccw_boundary - else: - boundary = self.cw_boundary + However, the boundary-corner vertices inserted between two consecutive + endpoints use the *shorter* of the two possible arcs (forward or + backward). These decisions are independent: the forward walk gives + correct ring topology, while the shorter-arc rule correctly handles the + edge case where two crossings straddle a boundary corner on opposite + sides of the perimeter seam. - def boundary_distance(xy): - return boundary.project(sgeom.Point(*xy)) + Note: the shorter-arc rule assumes a convex projection boundary. - # Squash all the LineStrings into a single list. - line_strings = [] - for multi_line_string in multi_line_strings: - line_strings.extend(multi_line_string.geoms) - - # Record the positions of all the segment ends - for i, line_string in enumerate(line_strings): - first_dist = boundary_distance(line_string.coords[0]) - thing = _BoundaryPoint(first_dist, False, - (i, 'first', line_string.coords[0])) - edge_things.append(thing) - last_dist = boundary_distance(line_string.coords[-1]) - thing = _BoundaryPoint(last_dist, False, - (i, 'last', line_string.coords[-1])) - edge_things.append(thing) - - # Record the positions of all the boundary vertices - for xy in boundary.coords[:-1]: - point = sgeom.Point(*xy) - dist = boundary.project(point) - thing = _BoundaryPoint(dist, True, point) - edge_things.append(thing) - - if debug_plot_edges: - import matplotlib.pyplot as plt - current_fig = plt.gcf() - fig = plt.figure() - # Reset the current figure so we don't upset anything. - plt.figure(current_fig.number) - ax = fig.add_subplot(1, 1, 1) - - # Order everything as if walking around the boundary. - # NB. We make line end-points take precedence over boundary points - # to ensure that end-points are still found and followed when they - # coincide. - edge_things.sort(key=lambda thing: (thing.distance, thing.kind)) - remaining_ls = dict(enumerate(line_strings)) - - prev_thing = None - for edge_thing in edge_things[:]: - if (prev_thing is not None and - not edge_thing.kind and - not prev_thing.kind and - edge_thing.data[0] == prev_thing.data[0]): - j = edge_thing.data[0] - # Insert a edge boundary point in between this geometry. - mid_dist = (edge_thing.distance + prev_thing.distance) * 0.5 - mid_point = boundary.interpolate(mid_dist) - new_thing = _BoundaryPoint(mid_dist, True, mid_point) - if debug: - print(f'Artificially insert boundary: {new_thing}') - ind = edge_things.index(edge_thing) - edge_things.insert(ind, new_thing) - prev_thing = None + Parameters + ---------- + multi_line_strings : list of MultiLineString + One MultiLineString per source ring (exterior or interior). + is_ccw : bool + Expected winding of the source polygon exterior ring. + """ + boundary = self.ccw_boundary if is_ccw else self.cw_boundary + perimeter = boundary.length + threshold = self.threshold + + # Pre-compute boundary corner (distance, coord) pairs once, sorted + # for O(log N) bisect queries inside fwd_arc_corners(). + # Explicitly extract (x, y) to strip any z value that Shapely may + # include in coords tuples, keeping corner_coords homogeneous 2-tuples. + corner_pairs = sorted( + (boundary.project(sgeom.Point(c)), (c[0], c[1])) + for c in boundary.coords[:-1] + ) + corner_dists = [d for d, _ in corner_pairs] + corner_coords = [c for _, c in corner_pairs] + + def fwd_arc_corners(d_from, arc_len): + """Boundary corners on the forward arc of *arc_len* from *d_from*.""" + if arc_len == 0: + return [] + if d_from + arc_len <= perimeter: + lo = bisect.bisect_right(corner_dists, d_from) + hi = bisect.bisect_left(corner_dists, d_from + arc_len) + return corner_coords[lo:hi] else: - prev_thing = edge_thing - - if debug: - print() - print('Edge things') - for thing in edge_things: - print(' ', thing) - if debug_plot_edges: - for thing in edge_things: - if isinstance(thing.data, sgeom.Point): - ax.plot(*thing.data.xy, marker='o') - else: - ax.plot(*thing.data[2], marker='o') - ls = line_strings[thing.data[0]] - coords = np.array(ls.coords) - ax.plot(coords[:, 0], coords[:, 1]) - ax.text(coords[0, 0], coords[0, 1], thing.data[0]) - ax.text(coords[-1, 0], coords[-1, 1], - f'{thing.data[0]}.') - - def filter_last(t): - return t.kind or t.data[1] == 'first' - - edge_things = list(filter(filter_last, edge_things)) - - processed_ls = [] - while remaining_ls: - # Rename line_string to current_ls - i, current_ls = remaining_ls.popitem() - - if debug: - import sys - sys.stdout.write('+') - sys.stdout.flush() - print() - print(f'Processing: {i}, {current_ls}') + wrap_end = d_from + arc_len - perimeter + lo = bisect.bisect_right(corner_dists, d_from) + hi = bisect.bisect_left(corner_dists, wrap_end) + return corner_coords[lo:] + corner_coords[:hi] + + def arc_corners(d_from, d_to): + """Boundary corners on the shorter arc from *d_from* to *d_to*.""" + d_fwd = (d_to - d_from) % perimeter + if d_fwd <= perimeter / 2: + return fwd_arc_corners(d_from, d_fwd) + else: + return list(reversed(fwd_arc_corners(d_to, perimeter - d_fwd))) + + def endpoint_info(xy): + """Return (boundary_distance, is_on_boundary) for a segment endpoint.""" + pt = sgeom.Point(xy) + d = boundary.project(pt) + bd_pt = boundary.interpolate(d) + return d, math.hypot(xy[0] - bd_pt.x, xy[1] - bd_pt.y) < threshold + + # Flatten all source ring segments into one pool. + line_strings = [ls for mls in multi_line_strings for ls in mls.geoms] + if not line_strings: + return [] + + completed = [] + remaining_segs = [] # (d_start, d_end, seg_idx) for boundary segments + + # Partition segments by whether both endpoints are on the boundary. + # Interior segments (e.g. a coastline running entirely inside the + # domain that trace.pyx emits because a nearby cut-line was crossed) + # have endpoints far from the boundary. Close them directly as rings + # rather than running them through the topology loop, which could + # incorrectly add a midpoint along the boundary and change the shape. + for k, ls in enumerate(line_strings): + d_start, start_on = endpoint_info(ls.coords[0]) + d_end, end_on = endpoint_info(ls.coords[-1]) + if start_on and end_on: + remaining_segs.append((d_start, d_end, k)) + else: + completed.append(list(ls.coords)) + + # Sort by d_start so the forward-walk can always find the next segment + # with a single bisect query. + remaining_segs.sort(key=lambda t: t[0]) + # Parallel list of just the d_start values for bisect lookups. + remaining_segs_d = [t[0] for t in remaining_segs] + + while remaining_segs: + # Start a new ring from the segment with the smallest d_start. + # This choice ensures each ring consumes exactly the segments + # on its forward arc before reaching any other ring's start. + d_ring_start, d_end, i = remaining_segs.pop(0) + remaining_segs_d.pop(0) + ring_coords = list(line_strings[i].coords) + # d_end tracks the boundary distance of ring_coords[-1] without + # requiring an extra boundary.project() call each inner iteration. - added_linestring = set() while True: - # Find out how far around this linestring's last - # point is on the boundary. We will use this to find - # the next point on the boundary. - d_last = boundary_distance(current_ls.coords[-1]) - if debug: - print(f' d_last: {d_last!r}') - next_thing = _find_first_ge(edge_things, d_last) - # Remove this boundary point from the edge. - edge_things.remove(next_thing) - if debug: - print(' next_thing:', next_thing) - if next_thing.kind: - # We've just got a boundary point, add it, and keep going. - if debug: - print(' adding boundary point') - boundary_point = next_thing.data - combined_coords = (list(current_ls.coords) + - [(boundary_point.x, boundary_point.y)]) - current_ls = sgeom.LineString(combined_coords) - - elif next_thing.data[0] == i: - # We've gone all the way around and are now back at the - # first boundary thing. - if debug: - print(' close loop') - processed_ls.append(current_ls) - if debug_plot_edges: - coords = np.array(current_ls.coords) - ax.plot(coords[:, 0], coords[:, 1], color='black', - linestyle='--') + if not remaining_segs: + # No more segments; close this ring. + corners = arc_corners(d_end, d_ring_start) + ring_coords += corners + if not corners: + # No boundary corners on the closing arc. Insert a + # midpoint on the shorter arc so the ring has a point + # on the boundary (needed for correct winding). Skip + # when the arc has zero length (start == end). + d_fwd = (d_ring_start - d_end) % perimeter + if d_fwd > 0: + if d_fwd <= perimeter / 2: + d_mid = (d_end + d_fwd / 2) % perimeter + else: + d_mid = (d_end - (perimeter - d_fwd) / 2) % perimeter + mid_pt = boundary.interpolate(d_mid) + ring_coords.append((mid_pt.x, mid_pt.y)) + completed.append(ring_coords) break - else: - if debug: - print(' adding line') - j = next_thing.data[0] - line_to_append = line_strings[j] - if j in remaining_ls: - remaining_ls.pop(j) - coords_to_append = list(line_to_append.coords) - - # Build up the linestring. - current_ls = sgeom.LineString(list(current_ls.coords) + - coords_to_append) - - # Catch getting stuck in an infinite loop by checking that - # linestring only added once. - if j not in added_linestring: - added_linestring.add(j) - else: - if debug_plot_edges: - plt.show() - raise RuntimeError('Unidentified problem with ' - 'geometry, linestring being ' - 're-added. Please raise an issue.') - - # filter out any non-valid linear rings - def makes_valid_ring(line_string): - if len(line_string.coords) == 3: - # When sgeom.LinearRing is passed a LineString of length 3, - # if the first and last coordinate are equal, a LinearRing - # with 3 coordinates will be created. This object will cause - # a segfault when evaluated. - coords = list(line_string.coords) - return coords[0] != coords[-1] and line_string.is_valid - else: - return len(line_string.coords) > 3 and line_string.is_valid - - linear_rings = [ - sgeom.LinearRing(line_string) - for line_string in processed_ls - if makes_valid_ring(line_string)] - - if debug: - print(' DONE') - return linear_rings + # Find the next segment start going forward from d_end. + # bisect_right gives the index of the first remaining_segs_d + # entry > d_end. The modulo handles wrap-around (when all + # remaining starts are <= d_end we cycle back to the smallest). + pos = bisect.bisect_right(remaining_segs_d, d_end) % len(remaining_segs) + d_next, d_next_end, j = remaining_segs[pos] + + # Compare forward distance to the ring's own start vs. d_next. + # On a tie the ring closes without consuming the next segment; + # that segment will become the start of the following ring. + d_fwd_close = (d_ring_start - d_end) % perimeter + d_fwd_next = (d_next - d_end) % perimeter + + if d_fwd_close <= d_fwd_next: + # Ring start arrives first: close the ring. + corners = arc_corners(d_end, d_ring_start) + ring_coords += corners + if not corners and d_fwd_close > 0: + # No boundary corners on the closing arc. Insert a + # midpoint on the shorter arc so the ring has a point + # on the boundary (needed for correct winding). + if d_fwd_close <= perimeter / 2: + mid_d = (d_end + d_fwd_close / 2) % perimeter + else: + mid_d = (d_end - (perimeter - d_fwd_close) / 2) % perimeter + mid_pt = boundary.interpolate(mid_d) + ring_coords.append((mid_pt.x, mid_pt.y)) + completed.append(ring_coords) + break + else: + # Connect to segment j, then continue building the ring. + ring_coords += arc_corners(d_end, d_next) + ring_coords += list(line_strings[j].coords) + remaining_segs.pop(pos) + remaining_segs_d.pop(pos) + d_end = d_next_end + + # Build LinearRings, skipping degenerate coordinate lists. + # The minimum useful input is 3 distinct points. A 3-coord list + # where coords[0] == coords[-1] represents only 2 unique points and + # would produce a degenerate ring that isn't valid. + return [ + sgeom.LinearRing(coords) + for coords in completed + if len(coords) > 3 or (len(coords) == 3 and coords[0] != coords[-1]) + ] def _rings_to_multi_polygon(self, rings, is_ccw): exterior_rings = [] @@ -1267,8 +1256,14 @@ def _rings_to_multi_polygon(self, rings, is_ccw): if not polygon.is_empty: if isinstance(polygon, sgeom.MultiPolygon): polygon_bits.extend(polygon.geoms) - else: + elif isinstance(polygon, sgeom.Polygon): polygon_bits.append(polygon) + elif isinstance(polygon, sgeom.GeometryCollection): + for geom in polygon.geoms: + if isinstance(geom, sgeom.Polygon): + polygon_bits.append(geom) + elif isinstance(geom, sgeom.MultiPolygon): + polygon_bits.extend(geom.geoms) return sgeom.MultiPolygon(polygon_bits) @@ -3270,40 +3265,6 @@ def __init__(self, rotation=45, false_easting=0.0, false_northing=0.0, globe=Non ] -class _BoundaryPoint: - def __init__(self, distance, kind, data): - """ - A representation for a geometric object which is - connected to the boundary. - - Parameters - ---------- - distance: float - The distance along the boundary that this object - can be found. - kind: bool - Whether this object represents a point from the - pre-computed boundary. - data: point or namedtuple - The actual data that this boundary object represents. - - """ - self.distance = distance - self.kind = kind - self.data = data - - def __repr__(self): - return f'_BoundaryPoint({self.distance!r}, {self.kind!r}, {self.data})' - - -def _find_first_ge(a, x): - for v in a: - if v.distance >= x: - return v - # We've gone all the way around, so pick the first point again. - return a[0] - - def epsg(code): """ Return the projection which corresponds to the given EPSG code. diff --git a/lib/cartopy/tests/test_polygon.py b/lib/cartopy/tests/test_polygon.py index 0b5ecb25b..a3170f958 100644 --- a/lib/cartopy/tests/test_polygon.py +++ b/lib/cartopy/tests/test_polygon.py @@ -281,6 +281,57 @@ def test_project_degenerate_poly(self): polygons = target.project_geometry(polygon, source) assert isinstance(polygons, sgeom.MultiPolygon) + def test_shorter_boundary_arc(self): + # Verify that arc_corners() picks the *shorter* boundary arc and not + # always the forward (CCW) one. + # + # PlateCarree ccw_boundary parameterisation with distances (d) + # BL (-180,-90) d=0 -> BR (180,-90) d=360 -> + # TR (180, 90) d=540 -> TL (-180, 90) d=900 -> BL d=1080 + # + # Seg A: (180, 0) [d=450, right edge] -> (-180, 80) [d=910, left edge] + # Seg B: (-170, 90) [d=890, top edge] -> (180, 0) [d=450, right edge] + proj = ccrs.PlateCarree() + mls = sgeom.MultiLineString([ + [(180, 0), (-180, 80)], # Seg A: right edge d=450, left edge d=910 + [(-170, 90), (180, 0)], # Seg B: top edge d=890, right edge d=450 + ]) + rings = proj._attach_lines_to_boundary([mls], is_ccw=True) + assert rings, "Expected at least one ring to be produced" + domain_area = proj.domain.area + for ring in rings: + ring_area = sgeom.Polygon(ring).area + assert ring_area < 0.5 * domain_area, ( + f"Ring covers {ring_area / domain_area:.0%} of the domain; " + "expected the shorter boundary arc to be chosen." + ) + + def test_interior_segment_no_spurious_midpoint(self): + # A line-string whose endpoints lie entirely within the domain + # (not touching the boundary) shouldn't add a midpoint on the boundary + # + # ObliqueMercator along the Seattle–Tokyo route: right boundary at + # x ≈ 20 038 km. We feed a 4-point open ring near x ≈ 10 000 km + # (firmly interior). Both endpoints happen to project to the right + # boundary edge, but the actual coordinates are ~10 000 km away. + proj = ccrs.ObliqueMercator( + central_longitude=-161.07, + central_latitude=54.55, + azimuth=90.0, + scale_factor=1, + ) + x0 = 9_950_000.0 + x1 = 10_050_000.0 + coords = [ + [(x0, -50_000), (x1, -50_000), (x1, 50_000), (x0, 50_000)], + ] + mls = sgeom.MultiLineString(coords) + rings = proj._attach_lines_to_boundary([mls], is_ccw=True) + + # The ring should just close itself and not attach to the + # right boundary / add any points + assert len(rings) == 1 + assert rings[0] == sgeom.LinearRing(mls.geoms[0].coords) class TestQuality: def setup_class(self):