Enable convertSphericalToCartesian to handle latitudes across poles #68

fmahebert · 2023-04-26T05:23:49Z

The atlas library generates points with latitudes outside the usual [-90°, 90°] interval, for example halo points on a regular Gaussian grid. Passing these points through eckit's convertSphericalToCartesian throws an exception because the code expects latitudes within the usual interval.

This PR preserves the current default behavior, but adds the option to renormalise points with latitudes "across the pole".

(I'd be happy to update the PR to remove the bool and always apply the renormalisation, if that is the preferred behavior. This would enable the change to propagate transparently to atlas without adding the bool, but would affect users who relied on the exception for points outside [-90,90]).

FussyDuck · 2023-04-26T05:23:54Z

All committers have signed the CLA.

fmahebert · 2023-04-26T05:26:14Z

@wdeconinck, pinging you as requested.

wdeconinck · 2023-04-26T09:00:35Z

I think this is a great development, thanks @fmahebert .

Indeed for some numerics in Atlas, there can be stencil points that are placed "across the pole".
As for this algorithm, I believe it is supposed to do what it says, accurately and correctly, whichever the provided latitude is.

Regarding the boolean flag, I would remove it, as well as the exception throwing.
We usually treat exceptions truly as exceptions, meaning the program will (hopefully) exit gracefully when this happens.
Therefore, there should not be programs that rely on the exception to be thrown to follow a different code path.
Removing the exception should therefore be safe to do, and applications that need to guarantee this, could always assert on the application side.

codecov-commenter · 2023-04-26T09:11:00Z

Codecov Report

Patch coverage: 47.88% and project coverage change: -0.02 ⚠️

Comparison is base (d97a4b6) 62.53% compared to head (39adf87) 62.51%.

Additional details and impacted files

@@             Coverage Diff             @@
##           develop      #68      +/-   ##
===========================================
- Coverage    62.53%   62.51%   -0.02%     
===========================================
  Files          777      777              
  Lines        44012    44067      +55     
===========================================
+ Hits         27523    27550      +27     
- Misses       16489    16517      +28

Impacted Files	Coverage Δ
src/eckit/exception/Exceptions.h	`47.36% <ø> (ø)`
src/eckit/filesystem/LocalPathName.cc	`66.33% <ø> (ø)`
src/eckit/filesystem/StdDir.cc	`72.22% <ø> (ø)`
src/eckit/io/AIOHandle.cc	`62.09% <ø> (ø)`
src/eckit/io/FDataSync.cc	`44.44% <ø> (ø)`
src/eckit/io/FOpenDataHandle.cc	`0.00% <ø> (ø)`
src/eckit/io/PeekHandle.cc	`0.00% <0.00%> (ø)`
src/eckit/io/SeekableHandle.cc	`0.00% <0.00%> (ø)`
src/eckit/linalg/LinearAlgebraDense.cc	`87.50% <ø> (ø)`
src/eckit/linalg/LinearAlgebraSparse.cc	`87.50% <ø> (ø)`
... and 26 more

... and 2 files with indirect coverage changes

Help us with your feedback. Take ten seconds to tell us how you rate us. Have a feature suggestion? Share it here.

☔ View full report in Codecov by Sentry.
📢 Do you have feedback about the report comment? Let us know in this issue.

fmahebert · 2023-04-26T17:14:20Z

Regarding the boolean flag, I would remove it, as well as the exception throwing.

Done. The second commit halfway undoes the first commit, so this'll be easier to review as a total diff.

fmahebert · 2023-04-26T17:33:37Z

For increased safety, I could handle cases where latitudes are outside [-270,270] ("across two poles") either by throwing an exception or by generalizing the normalization a bit more. Opinions on this?

wdeconinck · 2023-04-27T09:19:05Z

That would be even more robust yes, preferably by generalising normalisation more.

fmahebert · 2023-04-27T21:28:20Z

I've made the normalisation of the latitudes more robust and added a new test case. I think I'm satisfied, let me know what you think.

MarekWlasak · 2023-04-28T09:19:43Z

src/eckit/geometry/Sphere.cc

+    const auto& normalise_latitude = normalise_longitude;
+    const bool lat_across_pole = (normalise_latitude(Alonlat[1], -90.) > 90.);
+
+    const double lambda_deg = normalise_longitude(Alonlat[0] + (lat_across_pole ? 180. : 0.), -180.);


Hi - Maybe I am missing something but I can't quite see how you are dealing with the South Pole.

const bool lat_across_pole = (normalise_latitude(Alonlat[1], -90.) > 90.);

does not apply normalisation to the used latitude. It just checks if values are between -90 and +90.
normalise_latitude(Alonlat[1], -90.) -> always value between -90. and 270.
So I'm wondering, how is this different from

const bool lat_across_pole=(Alonlat[1] > 90. || Alonlat[1] < -90.);

If the input latitude has sufficiently large values that it "wraps around the poles" a few times, then the two are different. This is because a point can have very large latitudes but still be in the same hemisphere. As an extreme example, consider (lon,lat) = (x,370); this point is also (x,10), in the same hemisphere, and the coordinates can be "fixed" without needing to shift the longitude to x+180.

In this example...

(Alonlat[1] > 90. || Alonlat[1] < -90.) would evaluate to true, even though we do not want to shift the longitude value

normalise_latitude(Alonlat[1], -90.) returns 10, which is in [-90,90], consistent with not editing the longitude

Does that help? Or have I missed something?

@MarekWlasak the statement normalise_latitude(Alonlat[1], -90.) means that points "below the south pole" are re-labeled with a latitude "above the north pole"... e.g. lat=-100 --> lat=260. Using the normalise_latitude function avoids writing complicated branching to achieve the desired behavior.

Do I need to write more documentation? :)

@fmahebert - I see ... but that creates a discontinuity at the pole. I am not too fussed provided it does not break existing atlas functionality (e.g with structuredcolumns).

@MarekWlasak can you clarify what discontinuity you're referring to? This is just choosing a consistent "branch cut" of the periodic coordinates before converting to cartesian coordinates... I don't expect anything will be discontinuous in the output space...

@fmahebert - I was misunderstanding the use case. All is well.

tlmquintino · 2023-04-28T13:53:38Z

I would like @pmaciel to review this as this may have impact to https://github.com/ecmwf/mir

pmaciel · 2023-05-02T10:39:24Z

Hi @MarekWlasak, this is a very welcome change, because it addresses a fundamental issue on the assumptions regarding spherical coordinates: that we normalise the longitude, but not the latitude -- asserting that latitudes are within [-90, 90].

Let me add that this is an improvement, definitely. I have recently written very similar code to handle exactly this case (for a different project, this code returns normalised coordinates):
` lat = normal(lat, -90.);

    if (lat > 90.) {
        lat = 180. - lat;
        lon += 180.;
    }

    return {lat, lat == -90. || lat == 90. ? 0. : normal(lon, lon_minimum)};

`
However, there are two issues with this PR:

Is keeping the assumption [-90 <= lat <= 90 ] (the one we are removing) useful? My opinion (like yours, and @wdeconinck) is to remove this assumption -- but the consequence is that maybe asserting on [-90 <= lat <= 90 ] has to be moved downstream the ECMWF software stack (Atlas and MIR.) This ihas been very useful in the past, for robustness on user requests (in the millions per day) and improving the software stack. For this, we need the input from @tlmquintino and @wdeconinck (before discussion, I'm supporting this change)
That this PR is not complete: This is only affecting eckit::geometry::Sphere::convertSphericalToCartesian, but not the method (in the same class) centralAngle. And, in addition, EllipsoidOfRevolution::convertSphericalToCartesian and polygon::LonLatPolygon::contains (this specific one being optional).

Whichever direction we do, this can be brought in quite quickly, as this code is really well tested. So discussion should take place first.

fmahebert · 2023-05-02T15:09:23Z

@pmaciel Thanks for your feedback. I'll be happy to address (2), if the conclusion after discussing (1) is to proceed.

fmahebert · 2023-05-10T18:37:40Z

@pmaciel @wdeconinck Is there any update on whether item (1) is satisfactory? I don't want to put in the work for (2) if the PR will be rejected, but if doing (2) would help acceptance, I can do it now, just let me know...

pmaciel · 2023-05-11T07:48:18Z

Hi @fmahebert, sorry to not get back sooner. We decided to proceed with (1), ie. remove these assertions and apply similar changes to the remainder of the code. The other classes/method I listed is not comprehensive, but sufficient for this PR, could you please do these too?

I suggest to rename normalise_longitude to normalise_angle (and remove normalise_latitude).

fmahebert · 2023-05-11T23:20:55Z

@pmaciel Thanks! I've added 3 commits that make the same change to methods Sphere::centralAngle, EllipsoidOfRevolution::convertSphericalToCartesian and LonLatPolygon::contains, and also renames to normalise_angle.

I've also added a 4th commit that pulls all this common logic into a helper function — this may be a tad less efficient (a few more computations and stack variables), but should be easier to maintain. Up to you if you want this commit, if not I'll revert it.

pmaciel · 2023-05-15T19:34:12Z

Hi @fmahebert , @ytremolet, sorry for the late reply. Yes, for me this is done -- thank you for the contribution.

I'm currently working on additional functionality that will supercede this (it will include a PointLonLat type, similar to Atlas/MIR) but in the meantime this is a great addition. We will follow this with changes downstream to protect against wrongly calculated latitudes.

pmaciel · 2023-05-15T19:40:31Z

@fmahebert please note, recently a PR was merged into the develop brach moving the code base to C++17 (as a minimum requirement). This is planned for a release soon. Can you please rebase your PR to an updated develop tree? Thanks! (sorry for the extra work!)

…calToCartesian

fmahebert · 2023-05-15T19:50:05Z

@pmaciel I've rebased; please let me know if there's anything else.

Your described future developments on this front sound good too :)

…riate (API change, to preserve operational behaviour), ./format-sources.sh, cleanup

pmaciel · 2023-05-16T09:19:03Z

From internal discussions, so that we can smoothly transition the behaviour, I've added an extra argument (bool normalise_angle = false) where appropriate. This is an API change, to make optional the assertion -90<=lat<=90 (which is still the default ebhaviour), so you might need minor adaptations.

I've made the necesary changes, including updating the new unit tests. @fmahebert the changes we're proposing are pending your approval.

Feature/normalise lats across poles

fmahebert · 2023-05-17T02:17:30Z

@pmaciel Sorry for the slower response this time; I've merged your commit and added a tiny test fix as well.

tlmquintino · 2023-05-17T05:34:00Z

For what I can see there is a missing test for the assert / throw, when normalise_angle = false, ie. EXPECT_THROWS(...)

fmahebert · 2023-05-17T06:05:09Z

Added some new test statements for the throw

DJDavies2 · 2023-05-17T12:19:23Z

I am getting a compile failure with this:

/home/d03/frwd/cylc-run/eckit-68/share/mo-bundle/eckit/tests/geometry/test_polygon.cc:388:0: error: unterminated argument list invoking macro "EXPECT_THROWS_AS"
}

/home/d03/frwd/cylc-run/eckit-68/share/mo-bundle/eckit/tests/geometry/test_polygon.cc: In function 'void eckit::test::test_88(std::__cxx11::string&, int&, int)':
/home/d03/frwd/cylc-run/eckit-68/share/mo-bundle/eckit/tests/geometry/test_polygon.cc:270:9: error: 'EXPECT_THROWS_AS' was not declared in this scope
EXPECT_THROWS_AS(poly.contains({lonmid, 180. - latmid}, eckit::BadValue);
^~~~~~~~~~~~~~~~
/home/d03/frwd/cylc-run/eckit-68/share/mo-bundle/eckit/tests/geometry/test_polygon.cc:270:9: error: expected '}' at end of input
/home/d03/frwd/cylc-run/eckit-68/share/mo-bundle/eckit/tests/geometry/test_polygon.cc:270:9: error: expected '}' at end of input
/home/d03/frwd/cylc-run/eckit-68/share/mo-bundle/eckit/tests/geometry/test_polygon.cc: At global scope:
/home/d03/frwd/cylc-run/eckit-68/share/mo-bundle/eckit/tests/geometry/test_polygon.cc:270:9: error: expected '}' at end of input
gmake[2]: *** [eckit/tests/geometry/CMakeFiles/eckit_test_geometry_polygon.dir/test_polygon.cc.o] Error 1
gmake[1]: *** [eckit/tests/geometry/CMakeFiles/eckit_test_geometry_polygon.dir/all] Error 2

This is with gnu 7.3.0.

tlmquintino · 2023-05-17T13:01:00Z

I am getting a compile failure with this:

/home/d03/frwd/cylc-run/eckit-68/share/mo-bundle/eckit/tests/geometry/test_polygon.cc:388:0: error: unterminated argument list invoking macro "EXPECT_THROWS_AS" }

/home/d03/frwd/cylc-run/eckit-68/share/mo-bundle/eckit/tests/geometry/test_polygon.cc: In function 'void eckit::test::test_88(std::__cxx11::string&, int&, int)': /home/d03/frwd/cylc-run/eckit-68/share/mo-bundle/eckit/tests/geometry/test_polygon.cc:270:9: error: 'EXPECT_THROWS_AS' was not declared in this scope EXPECT_THROWS_AS(poly.contains({lonmid, 180. - latmid}, eckit::BadValue); ^~~~~~~~~~~~~~~~ /home/d03/frwd/cylc-run/eckit-68/share/mo-bundle/eckit/tests/geometry/test_polygon.cc:270:9: error: expected '}' at end of input /home/d03/frwd/cylc-run/eckit-68/share/mo-bundle/eckit/tests/geometry/test_polygon.cc:270:9: error: expected '}' at end of input /home/d03/frwd/cylc-run/eckit-68/share/mo-bundle/eckit/tests/geometry/test_polygon.cc: At global scope: /home/d03/frwd/cylc-run/eckit-68/share/mo-bundle/eckit/tests/geometry/test_polygon.cc:270:9: error: expected '}' at end of input gmake[2]: *** [eckit/tests/geometry/CMakeFiles/eckit_test_geometry_polygon.dir/test_polygon.cc.o] Error 1 gmake[1]: *** [eckit/tests/geometry/CMakeFiles/eckit_test_geometry_polygon.dir/all] Error 2

This is with gnu 7.3.0.

As per error message "unterminated argument list", there is a parenthesis missing to be closed.

fmahebert · 2023-05-17T14:04:08Z

there is a parenthesis missing to be closed

... serves me right for pushing an "easy" fix without compiling! 🙃

Fixed.

tlmquintino

Looks ok now

fmahebert · 2023-05-19T15:50:31Z

Thanks @pmaciel @tlmquintino and @wdeconinck for iterating on this!

wdeconinck · 2023-06-05T07:53:30Z

@tlmquintino, @simondsmart, please could we create a tagged release with this API change, so I can ensure this API exists in Atlas with a version check?

…oles Enable convertSphericalToCartesian to handle latitudes across poles

danovaro · 2023-06-05T14:28:33Z

@wdeconinck just released with tag 1.24.0

github-actions bot added the contributor label Apr 26, 2023

wdeconinck self-requested a review April 26, 2023 08:48

wdeconinck added the approved-for-ci Approved for CI run label Apr 26, 2023

github-actions bot removed the approved-for-ci Approved for CI run label Apr 26, 2023

wdeconinck added the approved-for-ci Approved for CI run label Apr 27, 2023

github-actions bot removed the approved-for-ci Approved for CI run label Apr 27, 2023

MarekWlasak reviewed Apr 28, 2023

View reviewed changes

tlmquintino requested a review from pmaciel April 28, 2023 13:52

fmahebert force-pushed the feature/normalise_lats_across_poles branch from 367690a to 9ec49f1 Compare May 11, 2023 23:17

fmahebert added 5 commits May 15, 2023 13:46

Enable SphericalToCartesian to handle latitudes across poles

fefda3c

Remove option; always renormalise lats across pole

d3fa921

Improve normalisation robustness for large angles

55e3674

Handle latitudes across poles in Sphere::centralAngle

999ea39

Handle latitudes across poles in EllipsoidOfRevolution::convertSpheri…

a88188e

…calToCartesian

fmahebert added 2 commits May 15, 2023 13:46

Handle latitudes across poles in LonLatPolygon::contains

82c2f45

Factor out coordinate operations into helper

03347f5

fmahebert force-pushed the feature/normalise_lats_across_poles branch from 9ec49f1 to 03347f5 Compare May 15, 2023 19:48

C++17

d21ac4a

pmaciel approved these changes May 15, 2023

View reviewed changes

tlmquintino self-requested a review May 15, 2023 20:25

Added additional argument "bool normalise_angle = false" where approp…

3df502a

…riate (API change, to preserve operational behaviour), ./format-sources.sh, cleanup

fmahebert added 2 commits May 16, 2023 20:13

Merge pull request #1 from ecmwf/feature/normalise_lats_across_poles

f06075e

Feature/normalise lats across poles

Fix Polygon test for new bool argument

1deb275

fmahebert mentioned this pull request May 17, 2023

Enable latitude normalisation in KDTree coordinate transform ecmwf/atlas#140

Merged

Code review: test throwing of exceptions

581440a

Fix paren

e995d0f

tlmquintino approved these changes May 19, 2023

View reviewed changes

tlmquintino merged commit 046ffb1 into ecmwf:develop May 19, 2023

pmaciel pushed a commit that referenced this pull request Jun 5, 2023

Merge pull request #68 from fmahebert/feature/normalise_lats_across_p…

2160a04

…oles Enable convertSphericalToCartesian to handle latitudes across poles

Enable convertSphericalToCartesian to handle latitudes across poles #68

Enable convertSphericalToCartesian to handle latitudes across poles #68

Conversation

fmahebert commented Apr 26, 2023 • edited Loading

FussyDuck commented Apr 26, 2023 • edited Loading

fmahebert commented Apr 26, 2023

wdeconinck commented Apr 26, 2023

codecov-commenter commented Apr 26, 2023

Codecov Report

fmahebert commented Apr 26, 2023

fmahebert commented Apr 26, 2023

wdeconinck commented Apr 27, 2023 • edited Loading

fmahebert commented Apr 27, 2023 • edited Loading

MarekWlasak Apr 28, 2023

Choose a reason for hiding this comment

wdeconinck Apr 28, 2023

Choose a reason for hiding this comment

fmahebert Apr 28, 2023 • edited Loading

Choose a reason for hiding this comment

fmahebert Apr 28, 2023

Choose a reason for hiding this comment

fmahebert Apr 28, 2023

Choose a reason for hiding this comment

MarekWlasak May 2, 2023

Choose a reason for hiding this comment

fmahebert May 2, 2023

Choose a reason for hiding this comment

MarekWlasak May 2, 2023

Choose a reason for hiding this comment

tlmquintino commented Apr 28, 2023

pmaciel commented May 2, 2023

fmahebert commented May 2, 2023

fmahebert commented May 10, 2023

pmaciel commented May 11, 2023

fmahebert commented May 11, 2023

pmaciel commented May 15, 2023

pmaciel commented May 15, 2023

fmahebert commented May 15, 2023

pmaciel commented May 16, 2023 • edited Loading

fmahebert commented May 17, 2023

tlmquintino commented May 17, 2023

fmahebert commented May 17, 2023

DJDavies2 commented May 17, 2023

tlmquintino commented May 17, 2023 • edited Loading

fmahebert commented May 17, 2023

tlmquintino left a comment

Choose a reason for hiding this comment

fmahebert commented May 19, 2023

wdeconinck commented Jun 5, 2023

danovaro commented Jun 5, 2023

fmahebert commented Apr 26, 2023 •

edited

Loading

FussyDuck commented Apr 26, 2023 •

edited

Loading

wdeconinck commented Apr 27, 2023 •

edited

Loading

fmahebert commented Apr 27, 2023 •

edited

Loading

fmahebert Apr 28, 2023 •

edited

Loading

pmaciel commented May 16, 2023 •

edited

Loading

tlmquintino commented May 17, 2023 •

edited

Loading