Skip to content

FIX: boundary attachment arc direction#2658

Open
greglucas wants to merge 1 commit intoSciTools:mainfrom
greglucas:fix/boundary-arc-direction
Open

FIX: boundary attachment arc direction#2658
greglucas wants to merge 1 commit intoSciTools:mainfrom
greglucas:fix/boundary-arc-direction

Conversation

@greglucas
Copy link
Copy Markdown
Contributor

Previously _attach_lines_to_boundary always connected projected LineString endpoints by walking forward around the boundary, even when the shorter path went backward. For cut lines that split large polygons near the middle of a boundary edge (e.g. ObliqueMercator over Alaska/Russia), this caused the ring to trace ~98% of the perimeter instead of ~2%, producing a polygon that is the complement of the intended shape — an inside-out land or ocean feature.

Additionally, it had the possibility of adding a midpoint on the far edge boundary while doing that tracing, even though the two endpoints were not actually on the boundary. So the linestring should really be closed with itself and not attached to the boundary. So I've added in an extra check to make sure that we are only adding boundary points/attaching when we are within the threshold of the boundary, not in the middle of the domain.

This is a pretty major overhaul of the _attach_line_to_boundary function to try and make the intent a little more clear and improve the performance some as well. It is largely the same logic, just refactored with one fix to choose the shortest walk direction. This may fix more cases because it is a pretty low-level thing that is done a lot, but this specifically addresses #2650.

Fixes #2650

@rcomer
Copy link
Copy Markdown
Member

rcomer commented Apr 5, 2026

Cycling to pick up #2659

@rcomer rcomer closed this Apr 5, 2026
@rcomer rcomer reopened this Apr 5, 2026
Copy link
Copy Markdown
Member

@rcomer rcomer left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This certainly makes the code easier to read! I also learned about a new module as I wasn't familiar with bisect.

I think I'm confused about the need to take the shortest route around the edge.

For the basic case with one linestring: if you go forward, you always get an exterior polygon, which may fill up most of the map. If you go backwards, you get the inverse polygon, but it's identified as an interior. So it will get inverted later on, and you end up with the same polygon as if you went forward.

So was the problem case one that should have been an interior contained within a polygon formed from a different linestring? If so, shouldn't the "find the next edge thing" logic have found the start of that outer linestring before going all the way round?

Comment thread lib/cartopy/crs.py
"""
boundary = self.ccw_boundary if is_ccw else self.cw_boundary
perimeter = boundary.length
threshold = self.threshold
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Do we have it documented anywhere what self.threshold is? I have not previously known about it.

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

👀 apparently not! I thought we did, but I don't see a docstring on it immediately. It is a useful tool though when trying to make lines look smoother. I did add an example using it in my other adaptive-resampling PR because I think that illustrates the point quite well when it is symmetric and shows you what different values will do.

https://github.com/SciTools/cartopy/pull/2647/changes#diff-44b7ad8f2f0073da2576021a84484f3e4dcf524edf5ddaad9d12d0225451577b

Comment thread lib/cartopy/crs.py Outdated
Comment on lines +1110 to +1122
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))
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Could this be done within arc_corners? Since it already has the check on which route is shortest there wouldn't be too much to add in there. Also I got myself a bit confused trying to map "start" <==> "to", "end" <==> "from" in my head, so having all that similar logic together might make it easier to follow.

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, let me think about this some more and how to consolidate it.

Comment thread lib/cartopy/crs.py Outdated
Comment on lines +1147 to +1152
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))
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Should this be skipped if the start and end are in the same place, similar to the case where there are no remaining lines?

Comment thread lib/cartopy/crs.py
@rcomer
Copy link
Copy Markdown
Member

rcomer commented Apr 5, 2026

This also gives a big improvement on the left panel from #2584 (comment), though it's still not quite there

image

@greglucas
Copy link
Copy Markdown
Contributor Author

Thanks for adding that other example here too. I can try and investigate that specific case some more and see if there is another condition we can add to capture that as well. It again looks like something where we are closing/adding points in wrong locations.

Your comments actually got me to thinking, I wonder if we could pass some of this information into project_linear() so that we could tell project_linear itself that this was a ring coming in, so all outputs need to be closed by either being attached to a boundary or upon themselves. That would likely mean moving (at least some of) the attach_to_boundary logic into trace.pyx, but I'm wondering if it could simplify some of the dealing with linestrings/rings on the Python caller side because then we could say "project all polygon exteriors, then project all holes and do a difference on the results". We wouldn't have to keep track of rings/linestrings together and then determine if a string should have been a ring etc. That would be encoded in the calling logic.

I'll have a think about some of this and experiment around to see if there is a potential for any simplification, or if it is really just pushing the same complicated logic around to a different place.

@rcomer
Copy link
Copy Markdown
Member

rcomer commented Apr 15, 2026

I tried rebasing this against main, and found that it additionally fixes #2137, so I guess that needs a combination of this and #2665

image

@rcomer
Copy link
Copy Markdown
Member

rcomer commented Apr 15, 2026

I also found that it solves the error from the second case at #2176 (comment), but the plot is probably not what the user wants
image

Background
----------
The old _project_linear_ring / _attach_lines_to_boundary pair treated
all projected line fragments identically regardless of whether the source
geometry was a ring. This caused two classes of correctness failures:

1. Boundary arc direction ambiguity: when a ring projected to two or
   more fragments, the boundary arc connecting consecutive fragments
   could be chosen in either direction around the projection boundary.

2. Inverted ring misclassification: a ring whose projection completely
   surrounds the destination domain (so the "interior" of the ring in the
   source CRS maps to the area outside the visible map) was silently
   discarded rather than inverted.

What changed in crs.py
----------------------
* _project_linear_ring rewritten to use the is_ring=True / as_ring_fragments()
  interface added to trace.pyx. Fragments are returned in traversal order,
  so the correct boundary arc direction is determined by walking the boundary
  from fragment[i].end to fragment[i+1].start in the winding direction
  dictated by is_ccw.

* _attach_lines_to_boundary updated to accept the ordered fragment list from
  _project_linear_ring and insert the boundary arc between consecutive
  fragment pairs using the correct arc direction.

* _rings_to_multi_polygon: added inverted-ring detection. When all exterior
  rings (by area) are smaller than the largest interior ring, the source
  polygon wraps around the projection domain. The small exterior rings are
  reclassified as interior rings and inverted against the domain bounding box.

* _project_polygon: now passes is_ccw (computed once from the source polygon
  exterior) down to each ring projection call so that interior rings use
  the correct opposite winding convention.

* Area guard: near-zero-area rings produced when a source vertex lands
  exactly on the projection boundary are discarded before they reach
  _attach_lines_to_boundary, preventing geometry exceptions.
@greglucas greglucas force-pushed the fix/boundary-arc-direction branch from cbe072d to 8e66f93 Compare April 30, 2026 03:11
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Projects

None yet

2 participants