Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fix SwathDefinition get_bbox_lonlats returning counter-clockwise coordinates #389

Merged
merged 5 commits into from
Nov 12, 2021

Conversation

djhoese
Copy link
Member

@djhoese djhoese commented Nov 5, 2021

This is a continuation of #385. That PR originally started as a "my data isn't clockwise so things are incorrect". I then found out that the fake data I was using was actually incorrect and didn't match the real world data I was trying to work with. So switching to 64-bit floats fixed that issue.

Now I'm working with other real world data (MiRS algorithm output) where the lon/lat arrays have their first element (row 0, column 0) in the south-west corner which results in the get_bbox_lonlats producing coordinates in a counter-clockwise path. This PR fixes this issue in a very similar way as my original fix #385 except for a few important points:

  1. I now only use the first two pixels to determine if the first and second side are in increasing in value. This should get around some of the concerns @mraspaud had from Fix SphPolygon producing unexpected results for 32-bit float coordinates #385.
  2. I discovered that the way I was using the utility create_test_longitude was actually not allowed. That is, requesting longitudes in decreasing order resulted in the function thinking we were crossing the anti-meridian. This is now fixed in this PR.

I also discovered that get_edge_lonlats was depending on the non-clockwise path of the data or at least a test it had was checking for that. So I added the force_clockwise keyword argument and set it to False for get_edge_lonlats as this makes sure the lonlats follow the data rather than provide a clockwise bbox.

  • Tests added
  • Tests passed
  • Passes git diff origin/main **/*py | flake8 --diff
  • Fully documented

@djhoese djhoese added the bug label Nov 5, 2021
@djhoese djhoese requested a review from mraspaud November 5, 2021 16:39
@djhoese djhoese self-assigned this Nov 5, 2021
@codecov
Copy link

codecov bot commented Nov 5, 2021

Codecov Report

Merging #389 (4a22787) into main (4389e19) will increase coverage by 0.00%.
The diff coverage is 100.00%.

Impacted file tree graph

@@           Coverage Diff           @@
##             main     #389   +/-   ##
=======================================
  Coverage   93.80%   93.80%           
=======================================
  Files          65       65           
  Lines       11038    11061   +23     
=======================================
+ Hits        10354    10376   +22     
- Misses        684      685    +1     
Flag Coverage Δ
unittests 93.80% <100.00%> (+<0.01%) ⬆️

Flags with carried forward coverage won't be shown. Click here to find out more.

Impacted Files Coverage Δ
pyresample/spherical.py 97.81% <ø> (ø)
pyresample/geometry.py 86.61% <100.00%> (+0.21%) ⬆️
pyresample/test/test_geometry.py 99.30% <100.00%> (+<0.01%) ⬆️
pyresample/test/utils.py 75.70% <100.00%> (-0.94%) ⬇️

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 4389e19...4a22787. Read the comment docs.

@coveralls
Copy link

coveralls commented Nov 5, 2021

Coverage Status

Coverage increased (+0.004%) to 93.609% when pulling 4a22787 on djhoese:bugfix-ccw-swath-bbox into 4389e19 on pytroll:main.

@djhoese
Copy link
Member Author

djhoese commented Nov 6, 2021

I've also added a small patch to the DynamicAreaDefinition. I discovered that dask nanmin/nanmax fails if one of the chunks provided to it is empty. This can happen with the way we were using boolean arrays: np.nanmin(xarr[xarr >= 0])

Copy link
Member

@mraspaud mraspaud left a comment

Choose a reason for hiding this comment

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

Thanks for looking into this! I believe this could be made a bit more generic.

pyresample/geometry.py Outdated Show resolved Hide resolved
Comment on lines 311 to 315
# compare the first two pixels in the right column
lat_is_increasing = lats[1][0] < lats[1][1]
# compare the first two pixels in the "top" column
lon_is_increasing = lons[0][0] < lons[0][1]
is_ccw = (lon_is_increasing and lat_is_increasing) or (not lon_is_increasing and not lat_is_increasing)
Copy link
Member

Choose a reason for hiding this comment

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

I can see two issues with this approach:

  1. This assumes two things: a. the top and right edges of the swath are perpendicular, and b. the swath is never aligned with the equator. For b., it would indeed be a very strange orbit, so I think we can live with this assumption (although it should be mentioned at least in the code), while a. could be a bit problematic with sensors that are not our usual imagers. An alternative test could be to check that the oriented angle around the top-right corner of the swath is not reflex.
  2. pole or antimeridian situations. A solution could be to convert the three points to cartesian and use the angle check described above.

Copy link
Member Author

Choose a reason for hiding this comment

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

Oh good point on the poles. So I could do the math on the first corner to get the angle between the two vectors in 3D space, but won't that always give the smallest angle between them? It's been a while since I've done this type of math (and I have a cold) but looking at https://www.omnicalculator.com/math/angle-between-two-vectors, how do I know if it is clockwise or counter-clockwise in 3D space?

Copy link
Member

Choose a reason for hiding this comment

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

because you are on the surface of a spheroid, you can give it an orientation...

Copy link
Member

Choose a reason for hiding this comment

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

i.e. the normal of the angle has to point away from the centre of the earth

Copy link
Member

Choose a reason for hiding this comment

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

This is important for compatibility with
:class:`pyresample.spherical.AreaBoundary`. This is required in
cases where data is not oriented in the traditional way where
the first element of data (row 0, col 0) is the north-west
Copy link
Member

Choose a reason for hiding this comment

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

An ascending pass and descending pass of a well behaved sensor (clockwise edges) have opposite first element orientations: north-west and south-east (not necessarily respectively, depending on how you look at it).
So I don't think north-west should be mentioned here, since the flip of lons and lats doesn't actually change the position of the first corner of the swath, right?

Copy link
Member Author

Choose a reason for hiding this comment

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

What if I mentioned north-west and south-east?

Copy link
Member

Choose a reason for hiding this comment

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

I think north-east and south-west cases can also happen (I'm thinking with Sentinel1 where I've seen this happen)

Copy link
Member Author

Choose a reason for hiding this comment

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

So you're saying that Sentinel1 data has either north-east or south-west as its first corner and the path followed by the existing logic still produces a clockwise series of coordinates? I suppose anything's possible. I'll see if I can remove the mention of north-west and be clearer about what I'm saying.

Copy link
Member

Choose a reason for hiding this comment

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

yes exactly. The sentinel 1 data we use comes flipped horizontally or vertically depending whether the pass is ascending or descending...

Copy link
Member

@mraspaud mraspaud left a comment

Choose a reason for hiding this comment

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

Good work, just some minor notes

pyresample/geometry.py Outdated Show resolved Hide resolved
pyresample/geometry.py Outdated Show resolved Hide resolved
pyresample/geometry.py Outdated Show resolved Hide resolved
@@ -1050,8 +1078,10 @@ def _compute_bound_centers_dask(self, proj_dict, lons, lats):
y_is_pole = (ymax >= 90 - epsilon) or (ymin <= -90 + epsilon)
if crs.is_geographic and x_passes_antimeridian and not y_is_pole:
Copy link
Member

Choose a reason for hiding this comment

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

If the crs is geographic, I would rather go to cartesian coordinates to do the computations here... but maybe this should be another PR?

Copy link
Member Author

Choose a reason for hiding this comment

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

Yes probably another PR. However, I'm not sure why cartesian coordinates would work better for this math. It is meant to be a simple workaround so my hope was that it wouldn't get too complicated. But if you see something then of course we could do it.

Copy link
Member

Choose a reason for hiding this comment

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

it's just that I don't like the thresholding needed around poles and antimeridian...

pyresample/test/utils.py Outdated Show resolved Hide resolved
pyresample/geometry.py Outdated Show resolved Hide resolved
Copy link
Member

@mraspaud mraspaud left a comment

Choose a reason for hiding this comment

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

LGTM

@mraspaud mraspaud merged commit 568e5d4 into pytroll:main Nov 12, 2021
@djhoese djhoese deleted the bugfix-ccw-swath-bbox branch November 12, 2021 13:04
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
Projects
None yet
Development

Successfully merging this pull request may close these issues.

3 participants