Updates to use visitInfo Tickets/dm 6625

src/CcdImage.cc

@@ -158,6 +158,7 @@ CcdImage::CcdImage(lsst::afw::table::SortedCatalogT<lsst::afw::table::SourceReco
     double lst_obs;
     double ra;
     double dec;
+    double hourAngle = std::numeric_limits<double>::signaling_NaN();
     if (camera == "megacam" || camera == "CFHT MegaCam") {
         airMass = meta->get<double>("AIRMASS");
         jd = meta->get<double>("MJD-OBS");  // Julian date
@@ -176,6 +177,15 @@ CcdImage::CcdImage(lsst::afw::table::SortedCatalogT<lsst::afw::table::SourceReco
         ra = lsst::afw::geom::degToRad(RaStringToDeg(meta->get<std::string>("RA2000")));
         dec = lsst::afw::geom::degToRad(DecStringToDeg(meta->get<std::string>("DEC2000")));
     }
+    else if (camera == "DECAM" || camera == "decam") {
+        airMass = meta->get<double>("AIRMASS");
+        jd = meta->get<double>("MJD-OBS");  // Julian date
+        expTime = meta->get<double>("EXPTIME");
+        latitude = lsst::afw::geom::degToRad(-30.16606);
+        hourAngle = lsst::afw::geom::degToRad(RaStringToDeg(meta->get<std::string>("HA")));
+        ra = lsst::afw::geom::degToRad(RaStringToDeg(meta->get<std::string>("RA")));
+        dec = lsst::afw::geom::degToRad(DecStringToDeg(meta->get<std::string>("DEC")));
+     }
     // TODO: Massive hack to get my test data to run.
     // TODO: this all needs to go away once DM-5501 is dealt with.
     else if (camera == "monkeySim") {
@@ -189,10 +199,11 @@ CcdImage::CcdImage(lsst::afw::table::SortedCatalogT<lsst::afw::table::SourceReco
     else {
         throw LSST_EXCEPT(pex::exceptions::InvalidParameterError,"jointcal::CcdImage does not understand your camera "+camera);
     }
-    hourAngle = (lst_obs-ra);
-    if  (hourAngle>M_PI) hourAngle -= 2*M_PI;
-    if  (hourAngle<-M_PI) hourAngle += 2*M_PI;
-
+    if (std::isnan(hourAngle) == true) {
+        hourAngle = (lst_obs-ra);
+        if  (hourAngle>M_PI) hourAngle -= 2*M_PI;
+        if  (hourAngle<-M_PI) hourAngle += 2*M_PI;
+    }
     if (airMass==1)
        sineta = coseta = tgz = 0;
     else

tests/jointcalTestBase.py

@@ -171,22 +171,23 @@ def _do_plots(self, dataRefs, visitCatalogs,
         astropy.visualization.quantity_support()
         plt.ion()
 
-        plot_all_wcs_deltas(plt, dataRefs, visitCatalogs, self.oldWcsList)
-
+        name = self.id().strip('__main__.')
         old_rms_relative = rms_per_star(old_relative)
         old_rms_absolute = rms_per_star(old_absolute)
         new_rms_relative = rms_per_star(new_relative)
         new_rms_absolute = rms_per_star(new_absolute)
-        print(len(dataRefs))
+        print("N dataRefs:", len(dataRefs))
         print("relative RMS (old, new):", old_rel_total, new_rel_total)
         print("absolute RMS (old, new):", old_abs_total, new_abs_total)
         plot_rms_histogram(plt, old_rms_relative, old_rms_absolute, new_rms_relative, new_rms_absolute,
-                           old_rel_total, old_abs_total, new_rel_total, new_abs_total)
+                           old_rel_total, old_abs_total, new_rel_total, new_abs_total, name)
+
+        plot_all_wcs_deltas(plt, dataRefs, visitCatalogs, self.oldWcsList, name)
 
         # So one can muck-about with things after plotting...
-        plt.show()
-        import pdb
-        pdb.set_trace()
+        # plt.show()
+        # import pdb
+        # pdb.set_trace()
 
 
 def rms_per_star(data):
@@ -213,11 +214,14 @@ def rms_total(data):
     return (np.sqrt(total/n) * u.radian).to(u.arcsecond)
 
 
-def plot_all_wcs_deltas(plt, dataRefs, visitCatalogs, oldWcsList, perCcdPlots=False):
+def plot_all_wcs_deltas(plt, dataRefs, visitCatalogs, oldWcsList, name,
+                        perCcdPlots=False):
     """Various plots of the difference between old and new Wcs."""
 
-    plot_all_wcs_quivers(plt, dataRefs, visitCatalogs, oldWcsList)
-    plot_wcs_magnitude(plt, dataRefs, visitCatalogs, oldWcsList)
+    print('plotting heat maps')
+    plot_wcs_magnitude(plt, dataRefs, visitCatalogs, oldWcsList, name)
+    print('plotting quivers')
+    plot_all_wcs_quivers(plt, dataRefs, visitCatalogs, oldWcsList, name)
 
     if perCcdPlots:
         for i, ref in enumerate(dataRefs):
@@ -248,12 +252,12 @@ def wcs_convert(xv, yv, wcs):
     return xout, yout
 
 
-def plot_all_wcs_quivers(plt, dataRefs, visitCatalogs, oldWcsList):
+def plot_all_wcs_quivers(plt, dataRefs, visitCatalogs, oldWcsList, name):
     """Make quiver plots of the WCS deltas for each CCD in each visit."""
 
     for cat in visitCatalogs:
         fig = plt.figure()
-        fig.set_tight_layout(True)
+        # fig.set_tight_layout(True)
         ax = fig.add_subplot(111)
         for old_wcs, ref in zip(oldWcsList, dataRefs):
             if ref.dataId['visit'] != cat:
@@ -268,6 +272,7 @@ def plot_all_wcs_quivers(plt, dataRefs, visitCatalogs, oldWcsList):
         plt.xlabel('RA')
         plt.ylabel('Dec')
         plt.title('visit: {}'.format(cat))
+        plt.savefig('.plots/{}-{}-quivers.pdf'.format(name, cat))
 
 
 def plot_wcs_quivers(ax, wcs1, wcs2, dim):
@@ -279,7 +284,7 @@ def plot_wcs_quivers(ax, wcs1, wcs2, dim):
     return ax.quiver(x1, y1, uu, vv, units='x', pivot='tail', scale=1e-3, width=1e-5)
 
 
-def plot_wcs_magnitude(plt, dataRefs, visitCatalogs, oldWcsList):
+def plot_wcs_magnitude(plt, dataRefs, visitCatalogs, oldWcsList, name):
     for cat in visitCatalogs:
         fig = plt.figure()
         fig.set_tight_layout(True)
@@ -302,8 +307,8 @@ def plot_wcs_magnitude(plt, dataRefs, visitCatalogs, oldWcsList):
             xmax = x1.max() if x1.max() > xmax else xmax
             ymax = y1.max() if y1.max() > ymax else ymax
             magnitude = (np.linalg.norm((uu, vv), axis=0)*u.radian).to(u.arcsecond).value
-            img = ax.imshow(magnitude, vmin=0, vmax=0.2,
-                            aspect='auto', extent=extent, cmap=plt.get_cmap('viridis'))
+            img = ax.imshow(magnitude, vmin=0, vmax=0.3,
+                            aspect='auto', extent=extent, cmap=plt.get_cmap('magma'))
             # TODO: add CCD bounding boxes to the plot once DM-5503 is finished.
             # TODO: add a circle for the full focal plane.
 
@@ -316,6 +321,7 @@ def plot_wcs_magnitude(plt, dataRefs, visitCatalogs, oldWcsList):
         plt.xlabel('RA')
         plt.ylabel('Dec')
         plt.title('visit: {}'.format(cat))
+        plt.savefig('.plots/{}-{}-heatmap.pdf'.format(name, cat))
 
 
 def plot_wcs(plt, wcs1, wcs2, dim, center=(0, 0), name=""):
@@ -332,13 +338,14 @@ def plot_wcs(plt, wcs1, wcs2, dim, center=(0, 0), name=""):
 
 def plot_rms_histogram(plt, old_rms_relative, old_rms_absolute,
                        new_rms_relative, new_rms_absolute,
-                       old_rel_total, old_abs_total, new_rel_total, new_abs_total):
+                       old_rel_total, old_abs_total, new_rel_total, new_abs_total,
+                       name):
     """Plot histograms of the star separations and their RMS values."""
     plt.figure()
 
     color_rel = 'black'
     ls_old = 'dotted'
-    color_abs = 'blue'
+    color_abs = 'green'
     ls_new = 'dashed'
     plotOptions = {'lw': 2, 'range': (0, 0.1)*u.arcsecond, 'normed': True,
                    'bins': 30, 'histtype': 'step'}
@@ -359,3 +366,4 @@ def plot_rms_histogram(plt, old_rms_relative, old_rms_absolute,
     plt.xlim(plotOptions['range'])
     plt.xlabel('arcseconds')
     plt.legend(loc='best')
+    plt.savefig('.plots/%s-histogram.pdf'%name)

tests/testJointcal_decam.py

@@ -0,0 +1,66 @@
+# See COPYRIGHT file at the top of the source tree.
+
+from __future__ import division, absolute_import, print_function
+
+import unittest
+import os
+
+from astropy import units as u
+
+import lsst.afw.geom
+import lsst.afw.coord as afwCoord
+import lsst.utils
+import lsst.pex.exceptions
+
+import jointcalTestBase
+
+try:
+    data_dir = lsst.utils.getPackageDir('validation_data_decam')
+    os.environ['ASTROMETRY_NET_DATA_DIR'] = os.path.join(data_dir, 'astrometry_net_data')
+except lsst.pex.exceptions.NotFoundError:
+    data_dir = None
+
+# We don't want the absolute astrometry to become significantly worse
+# than the single-epoch astrometry (about 0.040").
+# This value was empirically determined from the first run of jointcal on
+# this data, and will likely vary from survey to survey.
+absolute_error = 48e-3*u.arcsecond
+# Set to True for a comparison plot and some diagnostic numbers.
+do_plot = False
+
+
+# for MemoryTestCase
+def setup_module(module):
+    lsst.utils.tests.init()
+
+
+class JointcalTestDECAM(jointcalTestBase.JointcalTestBase, lsst.utils.tests.TestCase):
+    def setUp(self):
+        jointcalTestBase.JointcalTestBase.setUp(self)
+        self.do_plot = do_plot
+        self.match_radius = 0.1*lsst.afw.geom.arcseconds
+
+        # position of the validation_data_decam catalog
+        center = afwCoord.IcrsCoord(150.1191666*lsst.afw.geom.degrees, 2.20583333*lsst.afw.geom.degrees)
+        radius = 3*lsst.afw.geom.degrees
+        self._prep_reference_loader(center, radius)
+
+        self.input_dir = os.path.join(data_dir, 'data')
+        self.all_visits = [176837, 176846]
+
+    @unittest.skipIf(data_dir is None, "validation_data_decam not setup")
+    def test_jointcalTask_2_visits(self):
+        # NOTE: The relative RMS limit was empirically determined from the
+        # first run of jointcal on this data. We should always do better than
+        # this in the future!
+        relative_error = 25e-3*u.arcsecond
+        self._testJointCalTask(2, relative_error, absolute_error)
+
+
+# TODO: the memory test cases currently fail in jointcal. Filed as DM-6626.
+# class MyMemoryTestCase(lsst.utils.tests.MemoryTestCase):
+#     pass
+
+if __name__ == "__main__":
+    lsst.utils.tests.init()
+    unittest.main()

tests/testJointcal_hsc.py

@@ -2,6 +2,9 @@
 
 from __future__ import division, absolute_import, print_function
 
+import matplotlib
+matplotlib.use('Agg')
+
 import unittest
 import os
 
@@ -26,7 +29,7 @@
 # this data, and will likely vary from survey to survey.
 absolute_error = 52e-3*u.arcsecond
 # Set to True for a comparison plot and some diagnostic numbers.
-do_plot = False
+do_plot = True
 
 
 # for MemoryTestCase
@@ -45,8 +48,9 @@ def setUp(self):
         radius = 5*lsst.afw.geom.degrees
         self._prep_reference_loader(center, radius)
 
-        self.input_dir = os.path.join(data_dir, 'DATA')
-        self.all_visits = [903334, 903336, 903338, 903342, 903344, 903346]
+        self.input_dir = os.path.join(data_dir, 'DATA', 'rerun', 'validate_drp')
+        # self.all_visits = [903334, 903336, 903338, 903342, 903344, 903346]
+        self.all_visits = [903982, 904006, 904828, 904846]
 
     @unittest.skipIf(data_dir is None, "validation_data_hsc not setup")
     def test_jointcalTask_2_visits(self):
@@ -57,12 +61,20 @@ def test_jointcalTask_2_visits(self):
         self._testJointCalTask(2, relative_error, absolute_error)
 
     @unittest.skipIf(data_dir is None, "validation_data_hsc not setup")
-    def test_jointcalTask_6_visits(self):
+    def test_jointcalTask_4_visits(self):
         # NOTE: The relative RMS limit was empirically determined from the
         # first run of jointcal on this data. We should always do better than
         # this in the future!
-        relative_error = 10e-3*u.arcsecond
-        self._testJointCalTask(6, relative_error, absolute_error)
+        relative_error = 17e-3*u.arcsecond
+        self._testJointCalTask(4, relative_error, absolute_error)
+
+    # @unittest.skipIf(data_dir is None, "validation_data_hsc not setup")
+    # def test_jointcalTask_6_visits(self):
+    #     # NOTE: The relative RMS limit was empirically determined from the
+    #     # first run of jointcal on this data. We should always do better than
+    #     # this in the future!
+    #     relative_error = 10e-3*u.arcsecond
+    #     self._testJointCalTask(6, relative_error, absolute_error)
 
 
 # TODO: the memory test cases currently fail in jointcal. Filed as DM-6626.

tests/testJointcal_lsstSim.py

@@ -2,6 +2,9 @@
 
 from __future__ import division, absolute_import, print_function
 
+import matplotlib
+matplotlib.use('Agg')
+
 import unittest
 import os
 
@@ -27,7 +30,7 @@
 # this data, and will likely vary from survey to survey.
 absolute_error = 42e-3*u.arcsecond
 # Set to True for a comparison plot and some diagnostic numbers.
-do_plot = False
+do_plot = True
 
 
 # for MemoryTestCase