DM-14809: RingsSkyMap incorrect at south pole and RA=360 #23

PaulPrice · 2018-06-20T23:30:32Z

No description provided.

r-owen

Basically looks good. A few requested changes

r-owen · 2018-06-21T20:17:12Z

examples/plotSkyMap.py

@@ -165,17 +165,15 @@ def plotSkyMap2d(skyMap):
        plt.clf()
        axes = plt.axes()

-        colorCycle = matplotlib.rcParams['axes.color_cycle']
-        for i, tract in enumerate(skyMap):
+        for tract, color in zip(skyMap, cycle(list("rgbcmyk"))):


Note that list is not needed here since strings are iterables. Unless you feel it makes the intent clearer I suggest removing it.

r-owen · 2018-06-21T20:23:06Z

python/lsst/skymap/ringsSkyMap.py

+        verisons. ``version=0`` covers the period from first implementation
+        until DM-14809, at which point bugs were identified in the numbering
+        of tracts (affecting only tracts at RA=0). ``version=1`` uses the
+        post-DM-14809 tract numbering.


Please explain the numbering for both version=0 and version=1 (perhaps in a Notes section)

Thanks. That's very helpful.

r-owen · 2018-06-21T20:25:54Z

python/lsst/skymap/ringsSkyMap.py

-        @param[in] version: software version of this class, to retain compatibility with old instances
-        """
+    def __init__(self, config, version=1):
+        """Constructor"""
        # We count rings from south to north
        # Note: pole caps together count for one additional ring


This comment (pole caps together...) suggests that the south cap and north cap are part of the same tract, which is surely not true. Please clarify this comment, or remove it if it is not useful.

r-owen · 2018-06-21T20:32:40Z

python/lsst/skymap/ringsSkyMap.py

@@ -84,9 +92,17 @@ def getRingIndices(self, index):
            return self.config.numRings, 0


Pre-existing, but... I strongly suggest raising an exception if the input index is out of range. Otherwise it's hard to see what the behavior would be, but I doubt it would be good.

I realize the valid range may be slightly different for version 0 and version 1.

r-owen · 2018-06-21T20:32:49Z

python/lsst/skymap/ringsSkyMap.py

+        else:
+            while ring < self.config.numRings and index >= self._ringNums[ring]:
+                index -= self._ringNums[ring]
+                ring += 1
        return ring, index


I realize this is a pre-existing issue, but this code is very confusing because it uses the variable index to refer to two entirely different things. Please fix that by using two distinct names.

r-owen · 2018-06-21T20:40:30Z

python/lsst/skymap/ringsSkyMap.py

@@ -129,16 +146,18 @@ def findTract(self, coord):
            return self[0]
        elif dec > firstRingStart*-1:
            # Northern cap
-            return self[-1]
+            return self[self._numTracts - 1]


Ouch. Is it then true that relative indexing does not work and so -1 is not the last tract?

r-owen · 2018-06-21T20:44:21Z

python/lsst/skymap/ringsSkyMap.py

@@ -98,7 +114,8 @@ def generateTract(self, index):
            ra, dec = 0, 0.5*math.pi
        else:
            dec = self._ringSize*(ringNum + 1) - 0.5*math.pi
-            ra = math.fmod(self.config.raStart + 2*math.pi*tractNum/self._ringNums[ringNum], 2*math.pi)
+            ra = (self.config.raStart*afwGeom.degrees +


I see several places where you have to to get the starting RA as an Angle; please consider adding an attribute self._raStart that has this as an Angle (either cache it in __init__ or make it a property). I think it would make the resulting code clearer (especially when it mixes degrees and radians).

Also, pre-existing, but please consider having ra and dec be Angle instead of float in radians. The code would be easier to understand (no need to wonder what the units are). Or if you are keen to keep them as floats, please consider renaming to raRad, decRad to clarify the units.

r-owen · 2018-06-21T20:49:31Z

python/lsst/skymap/ringsSkyMap.py

@@ -163,7 +182,7 @@ def findAllTracts(self, coord):
        for r in [ringNum - 1, ringNum, ringNum + 1]:
            if r < 0 or r > self.config.numRings - 1:
                continue
-            tractNum = int(math.fmod(ra - self.config.raStart, 2*math.pi) /
+            tractNum = int((ra - self.config.raStart*afwGeom.degrees).wrap().asRadians() /
                           (2*math.pi/self._ringNums[r]) + 0.5)


Please perform this computation in a method (perhaps a private method, though I'd be tempted to unit test it unless existing tests suffice). I see it at least twice in this code and it is rather complicated and required updating so it seems a very good candidate.

r-owen · 2018-06-21T21:00:00Z

tests/helper/skyMapTestCase.py

+        skymap = self.getSkyMap()
+        centers = set([(coord.getRa().asRadians(), coord.getDec().asRadians()) for
+                       coord in (tract.getCtrCoord() for tract in skymap)])
+        self.assertEqual(len(centers), len(skymap))


Thank you very much for the new tests. However I worry that this test is not reliable in that the tiniest roundoff error in computing tract centers would allow the original problem to be undetected.

One possible solution is to round the recorded RA and Dec to some reasonable amount. A value based on neighborAngularSeparation would be ideal, but a fixed value such as 1e-4 radians would probably suffice.

This isn't testing the centers, so I don't understand your concern. So long as the centers are contained exclusively within the tract (which I think is a strongly desirable property of all skymaps), this test works for what it's intended to prevent.

r-owen · 2018-06-21T21:12:04Z

tests/test_ringsSkyMap.py

+        # Check that the first tract in the last ring exists
+        coord = self.getFirstTractLastRingCoord()
+        tract = self.skymap.findTract(coord)
+        self.assertTrue(tract.contains(coord))


Please eliminate the code duplication with the new methods in SkyMapTestCase.

Preferably please have this test case and the next one inherit from SkyMapTestCase (but the next one, for version 0 maps, would need to override one or both of the base class test methods).

This test case is for testing the specific instance of the skymap used by HSC. It doesn't need the more general tests contained in the SkyMapTestCase. In fact, since it is has many more tracts than the ones usually used with SkyMapTestCase, inheriting would increase the test runtime substantially without providing much (if any) additional safety.

PaulPrice · 2018-06-23T20:44:42Z

Your comments were very helpful in improving this submission and the code, thanks!

r-owen

Thanks for the improvements. I have a few more really minor suggestions (and returned to the suggestion to round RA and Dec in that one test, with more explanation)

r-owen · 2018-06-25T15:13:16Z

python/lsst/skymap/ringsSkyMap.py

+        verisons. ``version=0`` covers the period from first implementation
+        until DM-14809, at which point bugs were identified in the numbering
+        of tracts (affecting only tracts at RA=0). ``version=1`` uses the
+        post-DM-14809 tract numbering.


Thanks. That's very helpful.

r-owen · 2018-06-25T15:17:57Z

python/lsst/skymap/ringsSkyMap.py

+        -------
+        tractNum : `int`
+            Tract number within the ring (starts at 0 for the tract at raStart).
+        """


Please consider testing for "out of range". I realize it's a private method, but still...it's best to catch errors early.

r-owen · 2018-06-25T15:25:23Z

tests/helper/skyMapTestCase.py

+    got = [skymap.findTract(tract.getCtrCoord()).getId() for tract in skymap]
+    testcase.assertListEqual(got, expect)
+
+    # Check that the tract found for central coordinate of a tract is that tract


I think this comment applies to the test above, not the one below, which seems to assert that no two tracts have the same (numerically identical) centers, with no "find" involved.

I'm not sure this 2nd test is actually required, but I think it would be a lot more useful if it tested that no two tracts had nearly identical centers. One efficient way to do that is to round the RA and Dec values to some reasonable small value before putting the tuple in the set.

r-owen · 2018-06-25T15:26:38Z

tests/helper/skyMapTestCase.py

+    # Check that the tract ID is that used to index the skymap
+    expect = [tract.getId() for tract in skymap]
+    got = [skymap.findTract(tract.getCtrCoord()).getId() for tract in skymap]
+    testcase.assertListEqual(got, expect)


I think this test does what the comment below says. I am not sure it does what its own comment says, I would have expected a test for the comment above to look like this:

expect = range(len(skmap)) got = [tract.getId() for tract in skymap] testcase.assertEqual(expect, got)

Also I recommend using assertEqual instead of assertListEqual; it works fine on lists and is the technique recommended by the docs for the unittest module (though I believe it does delegate to assertListEqual for lists).

getRingIndices, findTract and findTractAll had off-by-one bugs. These affect the setting of primary flags (e.g. detect_isPrimary) in pipe_tasks for tracts near the RA=360 line. However, this involves renumbering some of the tracts (specifically, those on the RA=360 line). We use the existing version mechanism to support the old tract numbering. version=0 is valid up until DM-14809, version=1 is valid afterwards. The version is preserved in the skymap pickle, so this Just Works. For version=0, we leave the tract numbering as it is (broken in two places: there is one duplicated tract (the first tract on the first ring) and one missing tract (the first tract on the last ring)), but fix the lookup in the findTract method to work with the broken numbering. Included in these changes is a fix for non-zero raStart: simply using fmod doesn't work to wrap angles because it leaves negative values as negative. Using Angle.wrap works. Added specific tests for RingsSkyMap and generic tests for all skymaps to catch these sorts of bugs. Thanks to Sogo Mineo (NAOJ) for finding and diagnosing these bugs!

matplotlib's axes.color_cycle parameter is deprecated/gone, so use different method of cycling colors. Plot tract number instead of dot.

Requested during review.

r-owen requested changes Jun 21, 2018

View reviewed changes

PaulPrice force-pushed the tickets/DM-14809 branch from b1a6bb0 to 17693bc Compare June 23, 2018 20:41

r-owen approved these changes Jun 25, 2018

View reviewed changes

PaulPrice force-pushed the tickets/DM-14809 branch from 17693bc to ec6d491 Compare June 25, 2018 20:38

PaulPrice added 4 commits June 25, 2018 23:17

git: ignore coverage cache file

f77da49

plotSkyMap: update for modern matplotlib, py3

813c631

matplotlib's axes.color_cycle parameter is deprecated/gone, so use different method of cycling colors. Plot tract number instead of dot.

RingsSkyMap: clarify some docs and comments

9a352e9

Requested during review.

PaulPrice force-pushed the tickets/DM-14809 branch from 08c8846 to 9a352e9 Compare June 26, 2018 03:17

PaulPrice merged commit 9a352e9 into master Jun 27, 2018

ktlim deleted the tickets/DM-14809 branch August 25, 2018 05:56

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

DM-14809: RingsSkyMap incorrect at south pole and RA=360 #23

DM-14809: RingsSkyMap incorrect at south pole and RA=360 #23

PaulPrice commented Jun 20, 2018

r-owen left a comment

r-owen Jun 21, 2018

r-owen Jun 21, 2018

r-owen Jun 25, 2018

r-owen Jun 21, 2018

r-owen Jun 21, 2018

r-owen Jun 21, 2018

r-owen Jun 21, 2018

r-owen Jun 21, 2018

r-owen Jun 21, 2018

r-owen Jun 21, 2018

r-owen Jun 21, 2018

PaulPrice Jun 23, 2018

r-owen Jun 21, 2018

PaulPrice Jun 23, 2018

PaulPrice commented Jun 23, 2018

r-owen left a comment

r-owen Jun 25, 2018

r-owen Jun 25, 2018

r-owen Jun 25, 2018

r-owen Jun 25, 2018

		@@ -84,9 +92,17 @@ def getRingIndices(self, index):
		return self.config.numRings, 0

DM-14809: RingsSkyMap incorrect at south pole and RA=360 #23

DM-14809: RingsSkyMap incorrect at south pole and RA=360 #23

Conversation

PaulPrice commented Jun 20, 2018

r-owen 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

PaulPrice commented Jun 23, 2018

r-owen 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