Bug fix: Long SMAP shape can't shift if it crosses both anti and prime meridian - #22
Bug fix: Long SMAP shape can't shift if it crosses both anti and prime meridian#22Cameronsplaze wants to merge 5 commits into
Conversation
0634a20 to
2c8d4b1
Compare
|
The write-up behind this bug is here: https://github.com/asfadmin/CIRRUS-ASF/pull/702#discussion_r3680297774 I had to use Claude's help initially, then fixed up the comments. I keep thinking there must be a simpler way for this, but the only other example I can find takes a similar path: |
92b3aec to
bd39def
Compare
bd39def to
8a8b62c
Compare
reweeden
left a comment
There was a problem hiding this comment.
Honestly, I think this approach looks pretty straight forward to me! I was expecting much worse from your comment 😆. I like that it only has to touch the _shift_polygon function and that everything else then can remain the same.
There are some smaller refactors you could do, but I think the overall approach is good! It might be good to add some unit tests for _shift_polygon. I'm not sure if any exist yet.
| # but can't cross BOTH sides. | ||
| revolutions = [0] | ||
|
|
||
| prev_lon = coords[0][0] |
There was a problem hiding this comment.
Can we handle the case of empty polygons too? Would be nice to add a test for that to catch that regression
There was a problem hiding this comment.
So I guess the way you handled this is by making it a precondition of this function that polygons will not be empty? That's probably a valid approach but then I would assert that precondition with an assert statement at the top.
What I was actually trying to get at was a test for something like this assert _shift_polygon(Polygon([])) == Polygon([]) as this is how the code would have behaved before.
| coords = list(polygon.exterior.coords) | ||
|
|
||
| # Which "period" a coord is on around the earth. For a given polygon, | ||
| # each element is expected to be either {-1, 0} or {0, 1}, depending |
There was a problem hiding this comment.
I don't quite understand what the } notation means. Is it supposed to be a set notation? But then why wouldn't you just list it all together like {-1, 0, 1}?
There was a problem hiding this comment.
I was trying to say the revolutions = [0] list could either contain the values of -1 or 0, OR the values of 0, -1. And then the last two lines of this method are what shift the {-1, 0} points to be {0, 1}.
But looking at this again, I think that restriction was only there because _shift_polygon_back doesn't have a mod to deal with any number of rotations. I bet you could have it allow {-2, -1, 0, 1} even, and _shift_polygon could shift them to be {0,1,2,3}, which the mod in _shift_polygon_back would unwrap.
| if delta > 180: | ||
| revolutions.append(revolutions[-1] - 1) | ||
| elif delta < -180: | ||
| revolutions.append(revolutions[-1] + 1) | ||
| else: | ||
| revolutions.append(revolutions[-1]) |
There was a problem hiding this comment.
You could potentially simplify this a tiny bit:
| if delta > 180: | |
| revolutions.append(revolutions[-1] - 1) | |
| elif delta < -180: | |
| revolutions.append(revolutions[-1] + 1) | |
| else: | |
| revolutions.append(revolutions[-1]) | |
| prev_revolution = revolutions[-1] | |
| inc = 0 | |
| if delta > 180: | |
| inc = -1 | |
| elif delta < -180: | |
| inc = 1 | |
| revolutions.append(prev_revolution + inc) |
What I've done here:
- Pull out the parts of if blocks that's the same in each outside (
revolutions[-1]appears in each as does the call torevolutions.append) - Drop the else by moving it to a variable initialization (
inc = 0)
There was a problem hiding this comment.
That's fair, I like how readable revolutions.append(prev_revolution + inc) is. Though I might change inc to step?
lol it doesn't seem as bad anymore. I think trying to understand the math on my Friday night, really drained me. I think I might have to modify |
| def test_split_polygon_on_antimeridian_ccw_multiple_wraps(): | ||
| """A thin ribbon spiralling ~1.8 times around the earth. | ||
|
|
||
| It crosses the antimeridian twice in the same direction, so it | ||
| has to be cut on 180 and 540 once shifted, rather than just once. | ||
| """ | ||
| # Spiral east from longitude 30 to 690, climbing in latitude, 10 tall. | ||
| unwrapped = [30 + 60 * i for i in range(12)] | ||
|
|
||
| def lat(lon): | ||
| return -60 + (lon - 30) / 660 * 100 | ||
|
|
||
| def wrap(lon): | ||
| return ((lon + 180) % 360) - 180 | ||
|
|
||
| coords = [(wrap(lon), lat(lon)) for lon in unwrapped] | ||
| coords += [(wrap(lon), lat(lon) + 10) for lon in reversed(unwrapped)] | ||
| polygon = Polygon(coords) | ||
|
|
||
| polygons = list(split_polygon_on_antimeridian_ccw(polygon)) | ||
|
|
||
| for poly in polygons: | ||
| assert poly.exterior.is_ccw | ||
| assert poly.exterior.is_valid | ||
| # Everything has to land back in range, including the middle piece | ||
| # that spans a whole period. | ||
| min_lon, _, max_lon, _ = poly.bounds | ||
| assert min_lon >= -180 | ||
| assert max_lon <= 180 | ||
|
|
||
| # One piece per period the ribbon passes through. | ||
| assert len(polygons) == 3 | ||
| assert [poly.bounds[0] for poly in polygons] == [30.0, -180.0, -180.0] | ||
| assert [poly.bounds[2] for poly in polygons] == [180.0, 180.0, -30.0] | ||
|
|
||
| # No area is dropped or double counted by the extra cuts. The ribbon is | ||
| # 660 degrees of longitude by 10 of latitude. | ||
| assert sum(poly.area for poly in polygons) == pytest.approx(6600.0) |
There was a problem hiding this comment.
Before:
POLYGON ((30 -60, 90 -50.909091, 150 -41.818182, -150 -32.727273, -90 -23.636364, -30 -14.545455, 30 -5.454545, 90 3.636364, 150 12.727273, -150 21.818182, -90 30.909091, -30 40, -30 50, -90 40.909091, -150 31.818182, 150 22.727273, 90 13.636364, 30 4.545455, -30 -4.545455, -90 -13.636364, -150 -22.727273, 150 -31.818182, 90 -40.909091, 30 -50, 30 -60))
After:
GEOMETRYCOLLECTION ( POLYGON ((180 -37.272727, 180 -27.272727, 150 -31.818182, 90 -40.909091, 30 -50, 30 -60, 90 -50.909091, 150 -41.818182, 180 -37.272727)), POLYGON ((-180 -27.272727, -180 -37.272727, -150 -32.727273, -90 -23.636364, -30 -14.545455, 30 -5.454545, 90 3.636364, 150 12.727273, 180 17.272727, 180 27.272727, 150 22.727273, 90 13.636364, 30 4.545455, -30 -4.545455, -90 -13.636364, -150 -22.727273, -180 -27.272727)), POLYGON ((-180 27.272727, -180 17.272727, -150 21.818182, -90 30.909091, -30 40, -30 50, -90 40.909091, -150 31.818182, -180 27.272727)) )
| last = math.ceil((max_lon - 180) / 360) - 1 | ||
|
|
||
| polygons = [polygon] | ||
| for lon in (180 + 360 * n for n in range(first, last + 1)): |
There was a problem hiding this comment.
You could consider pulling out the generator here to a variable for readability maybe.
| for lon in (180 + 360 * n for n in range(first, last + 1)): | |
| antimeridian_lons = (180 + 360 * n for n in range(first, last + 1)) | |
| for lon in antimeridian_lons: |
| # Only antimeridians strictly inside the span can cut anything. A shape | ||
| # wrapping the globe several times crosses several of them. | ||
| first = int((min_lon - 180) // 360 + 1) | ||
| last = math.ceil((max_lon - 180) / 360) - 1 |
There was a problem hiding this comment.
I also noticed that last is defined as last = (...) - 1 but when it's used you're doing last + 1. Maybe you can just define it as 'one past the end' and drop both the -1 from the definition and the +1 from where it's used?
So basically you would set it up so that range(first, last) makes sense
| # but can't cross BOTH sides. | ||
| revolutions = [0] | ||
|
|
||
| prev_lon = coords[0][0] |
There was a problem hiding this comment.
So I guess the way you handled this is by making it a precondition of this function that polygons will not be empty? That's probably a valid approach but then I would assert that precondition with an assert statement at the top.
What I was actually trying to get at was a test for something like this assert _shift_polygon(Polygon([])) == Polygon([]) as this is how the code would have behaved before.
|
|
||
| # No area is dropped or double counted by the extra cuts. The ribbon is | ||
| # 660 degrees of longitude by 10 of latitude. | ||
| assert sum(poly.area for poly in polygons) == pytest.approx(6600.0) |
There was a problem hiding this comment.
Are you avoiding spelling out the whole point list for a particular reason? I think it's nice to have that even if it does mean you have to update the test every once in a while when the output changes a tiny bit because it means you can see how the code changes affect the output in the tests. For instance, I think it's nice to see in this PR that the rounding of certain points seems to have changed slightly.
reweeden
left a comment
There was a problem hiding this comment.
Oh, actually I just realized that the version number has not been bumped in the pyproject.toml file! The latest release is 0.4.0 and the version number is 0.4.0 so you'll need to bump it appropriately. I think probably 0.4.1 would be right since the API of the external functions is not changing.
Pull Request Checklist
I have: