DM-12924: SpherePoint.offset should work at the poles and for negative offsets #302

r-owen · 2018-01-02T19:52:39Z

No description provided.

kfindeisen

The code changes look fine, but the behavior of bearingTo and offset at the poles have not been documented as promised. I'd also like to see some test cases testing these two methods' new behavior at the poles.

kfindeisen · 2018-01-03T21:05:23Z

include/lsst/afw/geom/SpherePoint.h

@@ -98,6 +98,9 @@ class SpherePoint final {
     * @throws pex::exceptions::InvalidParameterError
     *         Thrown if `vector` is the zero vector.
     *
+     * @note If the SpherePoint is at a pole then longitude is set to 0.
+     * That provides predictable behavior for methods such as bearingTo and offset.
+     *


Are there any guarantees for what northPole.getLongitude() or southPole.getLongitude() returns, or does the pole's longitude depend on the arguments to SpherePoint(lon, lat)?

IIRC bearingTo will auto-generate a link, but offset will not. Have you checked?

I'll add @ref

kfindeisen · 2018-01-03T21:20:59Z

include/lsst/afw/geom/SpherePoint.h

-     *         Thrown if `this.atPole()`.
-     * @throws pex::exceptions::InvalidParameterError
-     *         Thrown if `amount` is negative.
-     *


If we're making all operations well-defined at the poles, do we still need atPole? IIRC it was originally created as a way to guard against these exceptions.

Please update the description to give the behavior at the poles (depending on how the new code relates to bearingTo, this may be just a matter of removing the phrase "not at a coordinate pole").

I believe atPole may still be useful, as some operations are not practical at the pole. One example: a user may query it to find out if longitude was set to 0 by the Point3D constructor; this could help avoid discontinuities in some situations.

kfindeisen · 2018-01-03T21:27:56Z

src/geom/SpherePoint.cc

-        buffer << "Cannot calculate offset from pole " << *this << ".";
-        throw pexExcept::DomainError(buffer.str());
-    }
-


As I understand it, bearingTo can give different results for the same point (e.g., the north pole) depending on how the longitude has been initialized. If this is the desired behavior, consider adding a warning to the documentation to make this explicit.

I think the documentation is sufficiently clear, now that I have reworked it.

kfindeisen · 2018-01-03T21:30:54Z

include/lsst/afw/geom/SpherePoint.h

@@ -265,9 +268,6 @@ class SpherePoint final {
     * @param other the point to which to measure the bearing
     * @returns the direction, as defined above, in the interval [0, 2&pi;).
     *
-     * @throws pex::exceptions::DomainError
-     *         Thrown if `this.atPole()`.
-     *


Please add a paragraph documenting the behavior at the poles. The current description implies that the bearing must always be -π/2 at the north pole and π/2 at the south pole, which I suspect is no longer what happens.

OK. I rewrote the description to one whose behavior is clearer at the poles.

kfindeisen · 2018-01-03T21:50:18Z

tests/test_spherePoint.py

-        """Test if 1-argument form of __init__ handles invalid input.
-        """
+    def testPoint3DConstructor(self):
+        # test poles


Did you mean to remove the docstring?

I meant to not include one. I think the name of the method suffices.

kfindeisen · 2018-01-03T21:54:50Z

tests/test_spherePoint.py

-            northPole.bearingTo(safe)
-        with self.assertRaises(pexEx.DomainError):
-            southPole.bearingTo(safe)
-
    def testBearingToValue(self):
        """Test if bearingTo() returns the expected value.
        """


I would like to see some test coverage at the poles. Does bearingTo produce the output it should (especially as a function of the pole's longitude)?

Fair enough. I had not noticed until now that the existing test are all on the equator.

kfindeisen · 2018-01-03T21:58:31Z

tests/test_spherePoint.py

    def testOffsetValue(self):
        """Test if offset() returns the expected value.
        """
        # This should cover arcs over the meridian, across the pole, etc.
        for lon1, lat1 in self._dataset:
            point1 = SpherePoint(lon1, lat1)
-            if point1.atPole():
-                continue


Again, please add some test coverage at the poles (including testing for consistency between offset and bearingTo/separation).

I will add this to testOffsetValues

r-owen · 2018-01-04T00:35:47Z

Good points. Can you please have another look? I improved the documentation (including finding a few more items that needed updating). I also changed SpherePoint.operator== to compare longitude even at poles, because it makes a difference for bearingTo and offset, and updated the documentation and unit test accordingly.

I also documented that rotated and offset return a SpherePoint with longitude=0 if it is at the pole but I am wondering if this is the best we can do. It is easy to implement (it happens automatically right now). But consider: the only time the result is at the pole is if the offset is due east or west (or zero and we started at the pole), and in these cases the obvious longitude is the input longitude. We could easily add a bit of code to make this happen.

However, if we can offset along a slightly different great circle and still end up at the pole due to roundoff error, then that weakens my argument.

kfindeisen

I have a few comments about the documentation, but looks good overall.

kfindeisen · 2018-01-04T19:18:51Z

include/lsst/afw/geom/SpherePoint.h

@@ -63,6 +63,8 @@ class SpherePoint final {
     * Construct a SpherePoint from a longitude and latitude.
     *
     * @param longitude The longitude of the point.
+     *                  It will be wrapped into the range [0, 2 &pi;) radians,
+     *                  and `+/-inf` is wrapped to `nan`


I can see why you added this, but I'd be a bit careful about separation of specification and implementation here. Maybe something like "getLongitude() will return the angle wrapped into..."? Same for the (double, double, AngleUnit) constructor.

Stating that the angle is wrapped here makes it easier to understand why infinity is converted to nan, but I'll try a gentler rewording.

kfindeisen · 2018-01-04T19:19:55Z

include/lsst/afw/geom/SpherePoint.h

@@ -85,6 +89,9 @@ class SpherePoint final {
     *         Thrown if `latitude` isout of range.
     *
     * @exceptsafe Provides strong exception guarantee.
+     *
+     * @warning If longitude is +/- inf then the resulting longitude is `NaN`
+     * since infinity cannot be wrapped into the range


This looks cut off. Did you mean to explicitly restate the range?

I meant to remove the warning entirely, and have now done so. The same information is contained in the description of the longitude argument.

kfindeisen · 2018-01-04T19:24:20Z

include/lsst/afw/geom/SpherePoint.h

-     * allows multiple values of longitude from a pole, they shall all be
-     * treated as valid representations of the same point.
+     * If this point is at a pole, the longitude is the value provided to the constructor,
+     * if using a longitude, latitude constructor, or 0 if using the Point3D constructor.


I would leave this out entirely, or say something more generic like "points at the pole have a well-defined longitude, and are not equivalent for the purposes of bearingTo/offset". Users can't necessarily tell how a particular SpherePoint object was constructed.

kfindeisen · 2018-01-04T19:30:25Z

include/lsst/afw/geom/SpherePoint.h

@@ -254,20 +260,19 @@ class SpherePoint final {
    bool operator!=(SpherePoint const& other) const noexcept;

    /**
-     * Direction from one point to another.
+     * Orientation at this point of the great circle arc from this point to another point.


This is a bit wordy for a brief description; can you make it more concise?

kfindeisen · 2018-01-04T19:32:56Z

include/lsst/afw/geom/SpherePoint.h

+     * the direction of increasing longitude and 90 degrees along the direction of increasing latitude.
+     * Thus for any point except the pole, the arc is due east at 0 degrees and due north at 90 degrees.
+     * To understand the behavior at the poles it may help to imagine the point has been shifted slightly
+     * by decreasing the absolute value of its latitude.


I think this makes it pretty clear how the pole is handled. 👍

kfindeisen · 2018-01-04T19:49:26Z

include/lsst/afw/geom/SpherePoint.h

+     * For any point `A` not at a coordinate pole, any angle `bearing` and any non-negative angle `amount`,
+     * if `B` = A.offset(bearing, amount)` then `A.bearingTo(B) = amount` and `A.separationTo(B) = amount`.
+     * Negative values of `amount` are supported in the obvious manner:
+     * `A.offset(bearing, delta) = A.offset(bearing + 180*degrees, -delta)`


As I understand it, the new unit tests confirm that these assertions still hold true if A is a pole. Did I miss something?

The new tests check all the points in the standard point list, which include points at the pole. Is that what you are asking?

The text says "For any point A not at a coordinate pole" (emphasis added). I don't understand why we still need that disclaimer.

kfindeisen · 2018-01-04T19:50:04Z

include/lsst/afw/geom/SpherePoint.h

-     * @throws pex::exceptions::InvalidParameterError
-     *         Thrown if `amount` is negative.
+     * @returns a new point created by rotating this point. If the point is at the pole
+     * then its longitude will be 0.


Again, it would be clearer if "the point" -> "the new point".

kfindeisen · 2018-01-04T19:57:06Z

tests/test_spherePoint.py

-        self.assertEqual(SpherePoint(360.0, -42.0, degrees),
-                         SpherePoint(0.0, -42.0, degrees))
-        self.assertEqual(SpherePoint(-90.0, -42.0, degrees),
-                         SpherePoint(270.0, -42.0, degrees))


I think the wrapping test is still useful; why remove it?

Good point. Most of the tests in this test method were obsoleted by my other changes but we still want to test wrapping of longitude. I'll restore this under a new test name.

kfindeisen · 2018-01-04T20:04:15Z

tests/test_spherePoint.py

+                        self.assertAnglesAlmostEqual(bearing, 90 * degrees)
+                    else:
+                        self.assertAnglesAlmostEqual(bearing, -90 * degrees)
+


Test bearing from a pole to an arbitrary point not on the same longitude? As I understand it, the answer should be the difference in longitudes minus 90 degrees for the north pole, and the difference plus 90 degrees for the south pole.

kfindeisen · 2018-01-04T20:09:15Z

tests/test_spherePoint.py

                    self.assertAnglesAlmostEqual(
                        point2.getLatitude(), newPoint.getLatitude())
+                    self.assertAnglesAlmostEqual(
+                        point2.getLatitude(), newPoint2.getLatitude())


These assertions are redundant with your use of assertSpherePointsAlmostEqual, above.

I disagree. assertSpherePointsAlmostEqual tests angular separation, and thus cannot tell the difference between two SpherePoints at (or very near) the pole that have different longitudes.

Then can you add a comment to that effect? I'd assumed the difference was mostly one of style.

r-owen · 2018-01-04T23:14:33Z

I'm done. Please have another look at the docs. I updated the lon/lat constructors and operator== (a brief note that wrap is ignored when comparing longitude)

kfindeisen · 2018-01-04T23:33:38Z

include/lsst/afw/geom/SpherePoint.h

@@ -225,12 +223,12 @@ class SpherePoint final {
     */

    /**
-     * `true` if two points represent the same position.
+     * `true` if two points have the same longitude (ignoring wrap) and latitude.


Sorry, but "ignoring wrap" is kind of confusing (does "ignore" mean that 90 degrees and 450 degrees are considered the same, or that they're not?). It also seems redundant, given that getLongitude() is guaranteed to return a wrapped angle or NaN.

I'm feel that I am in a Catch-22. The way I see it is:

If we assume that the constructor wraps longitude on construction (which you complained about when I documented that), then it is unambiguous what equal longitude means in operator==.

If we do not assume that the constructor wraps longitude on construction, then operator== is ambiguous, in that it is not clear if it compares the actual longitude or the wrapped version returned by getLongitude()

I personally prefer to document that longitude is wrapped on construction. That neatly explains why +/-inf is turned into nan and gives an the user a clear model of what is going on which can be extended to methods such as operator==. But I can see why you felt that exposed the implementation too much. If you still feel that way then please suggest clear and unambiguous documentation for operator==

How about removing the parenthetical note from the title, but writing:

@returns true if this point has exactly the same values of getLongitude()andgetLatitude()asother, false otherwise

and possibly removing the mention of getVector() from the warning, for consistency's sake.

(Also, I just noticed while looking at this that @note is in front of the params/returns sections; it should go next to the warning).

OK. I moved one other note as well

kfindeisen · 2018-01-05T00:15:56Z

include/lsst/afw/geom/SpherePoint.h

     *
     * @param other the point to test for equality
-     * @returns true if this point matches `other` exactly, false otherwise
+     * @returns true if this point has exactly the same values ofgetLongitude() and getLatitude()
+     * asother, false otherwise


Couple of missing spaces (last comment, I promise!)

For a SpherePoint constructed from a 3-vector, set longitude to 0 if at a pole and document the same. This allows bearingTo and offset to be well defined for all cases. Also change operator== to compare longitude even at the poles, since different longitudes result in different behavior for bearingTo and offset even at the poles.

r-owen force-pushed the tickets/DM-12924 branch from 60bb814 to 69470ac Compare January 2, 2018 19:56

kfindeisen requested changes Jan 3, 2018

View reviewed changes

r-owen force-pushed the tickets/DM-12924 branch from 140c187 to 4ca458e Compare January 4, 2018 00:33

r-owen force-pushed the tickets/DM-12924 branch 2 times, most recently from 0404461 to 3b35e19 Compare January 4, 2018 01:06

kfindeisen approved these changes Jan 4, 2018

View reviewed changes

r-owen force-pushed the tickets/DM-12924 branch 2 times, most recently from 6093ca5 to 67d7413 Compare January 4, 2018 23:13

kfindeisen reviewed Jan 4, 2018

View reviewed changes

r-owen force-pushed the tickets/DM-12924 branch 2 times, most recently from 5dd2578 to 420ab38 Compare January 5, 2018 00:14

kfindeisen reviewed Jan 5, 2018

View reviewed changes

r-owen force-pushed the tickets/DM-12924 branch from 420ab38 to dda929d Compare January 5, 2018 00:17

r-owen force-pushed the tickets/DM-12924 branch from dda929d to 1e32488 Compare January 5, 2018 22:08

r-owen merged commit 1e32488 into master Jan 8, 2018

ktlim deleted the tickets/DM-12924 branch August 25, 2018 06:44

DM-12924: SpherePoint.offset should work at the poles and for negative offsets #302

DM-12924: SpherePoint.offset should work at the poles and for negative offsets #302

Conversation

r-owen commented Jan 2, 2018

kfindeisen left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

r-owen Jan 3, 2018 • edited

Choose a reason for hiding this comment

r-owen commented Jan 4, 2018 • edited

kfindeisen left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

kfindeisen Jan 4, 2018 • edited

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

kfindeisen Jan 4, 2018 • edited

Choose a reason for hiding this comment

r-owen commented Jan 4, 2018

kfindeisen Jan 4, 2018 • edited

Choose a reason for hiding this comment

Choose a reason for hiding this comment

kfindeisen Jan 5, 2018 • edited

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

r-owen Jan 3, 2018 •

edited

r-owen commented Jan 4, 2018 •

edited

kfindeisen Jan 4, 2018 •

edited

kfindeisen Jan 4, 2018 •

edited

kfindeisen Jan 4, 2018 •

edited

kfindeisen Jan 5, 2018 •

edited