DM-5159: "Please use angle and Coord where possible"

python/lsst/validate/drp/calcsrd/amx.py

@@ -25,7 +25,7 @@
 import astropy.units as u
 
 from lsst.validate.base import MeasurementBase
-from ..util import averageRaDecFromCat
+from ..util import averageRaDecFromCat, sphDist
 
 
 class AMxMeasurement(MeasurementBase):
@@ -255,36 +255,6 @@ def magInRange(cat):
     return rmsDistances
 
 
-def sphDist(ra1, dec1, ra2, dec2):
-    """Calculate distance on the surface of a unit sphere.
-
-    Input and Output are in radians.
-
-    Notes
-    -----
-    Uses the Haversine formula to preserve accuracy at small angles.
-
-    Law of cosines approach doesn't work well for the typically very small
-    differences that we're looking at here.
-    """
-    # Haversine
-    dra = ra1-ra2
-    ddec = dec1-dec2
-    a = np.square(np.sin(ddec/2)) + \
-        np.cos(dec1)*np.cos(dec2)*np.square(np.sin(dra/2))
-    dist = 2 * np.arcsin(np.sqrt(a))
-
-    # This is what the law of cosines would look like
-#    dist = np.arccos(np.sin(dec1)*np.sin(dec2) + np.cos(dec1)*np.cos(dec2)*np.cos(ra1 - ra2))
-
-    # Could use afwCoord.angularSeparation()
-    #  but (a) that hasn't been made accessible through the Python interface
-    #  and (b) I'm not sure that it would be faster than the numpy interface.
-    #    dist = afwCoord.angularSeparation(ra1-ra2, dec1-dec2, np.cos(dec1), np.cos(dec2))
-
-    return dist
-
-
 def matchVisitComputeDistance(visit_obj1, ra_obj1, dec_obj1,
                               visit_obj2, ra_obj2, dec_obj2):
     """Calculate obj1-obj2 distance for each visit in which both objects are seen.

python/lsst/validate/drp/matchreduce.py

@@ -35,10 +35,10 @@
 from lsst.afw.fits.fitsLib import FitsError
 from lsst.validate.base import BlobBase
 
-from .util import (getCcdKeyName, averageRaDecFromCat)
+from .util import (getCcdKeyName, averageRaDecFromCat, sphDist)
 
 
-__all__ = ['MatchedMultiVisitDataset', 'positionRms']
+__all__ = ['MatchedMultiVisitDataset', 'positionRmsFromCat']
 
 
 class MatchedMultiVisitDataset(BlobBase):
@@ -314,7 +314,7 @@ def safeFilter(cat):
         self.safeMatches = safeMatches
 
 
-def positionRms(cat):
+def positionRmsFromCat(cat):
     """Calculate the RMS for RA, Dec for a set of observations an object.
 
     Parameters
@@ -325,16 +325,7 @@ def positionRms(cat):
     Returns
     -------
     pos_rms -- RMS of positions in milliarcsecond.  Float.
-
-    This routine doesn't handle wrap-around
     """
     ra_avg, dec_avg = averageRaDecFromCat(cat)
     ra, dec = cat.get('coord_ra'), cat.get('coord_dec')
-    # Approximating that the cos(dec) term is roughly the same
-    #   for all observations of this object.
-    ra_var = np.var(ra) * np.cos(dec_avg)**2
-    dec_var = np.var(dec)
-    pos_rms = np.sqrt(ra_var + dec_var)  # radians
-    pos_rms = afwGeom.radToMas(pos_rms)  # milliarcsec
-
-    return pos_rms
+    return positionRms(ra_avg, dec_avg, ra, dec)

python/lsst/validate/drp/util.py

@@ -25,6 +25,7 @@
 
 import os
 
+import numpy as np
 import yaml
 
 import lsst.daf.persistence as dafPersist
@@ -63,6 +64,65 @@ def averageRaDecFromCat(cat):
     return averageRaDec(cat.get('coord_ra'), cat.get('coord_dec'))
 
 
+def positionRms(ra_avg, dec_avg, ra, dec):
+    """Calculate the RMS between RA_avg, Dec_avg and RA, Dec
+
+    Parameters
+    ----------
+    ra_avg -- Average RA  [rad]
+    dec_avg -- Average Dec [rad]
+    ra -- np.array of RA  [rad]
+    dec -- np.array of Dec  [rad]
+
+    Returns
+    -------
+    pos_rms -- RMS of positions in milliarcsecond.  Float.
+
+
+    The RMS of a single-element array will be returned as 0.
+    The RMS of an empty array will be returned as NaN.
+    """
+    separations = sphDist(ra_avg, dec_avg, ra, dec)
+    # Note we don't want `np.std` of separations, which would give us the
+    #   std around the average of separations.
+    # We've already taken out the average,
+    #   so we want the sqrt of the mean of the squares.
+    pos_rms_rad = np.sqrt(np.mean(separations**2))  # radians
+    pos_rms_mas = afwGeom.radToMas(pos_rms_rad)  # milliarcsec
+
+    return pos_rms_mas
+
+
+def sphDist(ra1, dec1, ra2, dec2):
+    """Calculate distance on the surface of a unit sphere.
+
+    Input and Output are in radians.
+
+    Notes
+    -----
+    Uses the Haversine formula to preserve accuracy at small angles.
+
+    Law of cosines approach doesn't work well for the typically very small
+    differences that we're looking at here.
+    """
+    # Haversine
+    dra = ra1-ra2
+    ddec = dec1-dec2
+    a = np.square(np.sin(ddec/2)) + \
+        np.cos(dec1)*np.cos(dec2)*np.square(np.sin(dra/2))
+    dist = 2 * np.arcsin(np.sqrt(a))
+
+    # This is what the law of cosines would look like
+#    dist = np.arccos(np.sin(dec1)*np.sin(dec2) + np.cos(dec1)*np.cos(dec2)*np.cos(ra1 - ra2))
+
+    # Could use afwCoord.angularSeparation()
+    #  but (a) that hasn't been made accessible through the Python interface
+    #  and (b) I'm not sure that it would be faster than the numpy interface.
+    #    dist = afwCoord.angularSeparation(ra1-ra2, dec1-dec2, np.cos(dec1), np.cos(dec2))
+
+    return dist
+
+
 def getCcdKeyName(dataid):
     """Return the key in a dataId that's referring to the CCD or moral equivalent.
 

tests/test_positionRms.py

@@ -0,0 +1,110 @@
+#!/usr/bin/env python
+
+#
+# LSST Data Management System
+# Copyright 2012-2017 LSST Corporation.
+#
+# This product includes software developed by the
+# LSST Project (http://www.lsst.org/).
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the LSST License Statement and
+# the GNU General Public License along with this program.  If not,
+# see <http://www.lsstcorp.org/LegalNotices/>.
+#
+
+from __future__ import print_function
+
+import unittest
+
+import numpy as np
+
+from numpy.testing import assert_almost_equal, assert_array_max_ulp
+
+import lsst.utils
+from lsst.validate.drp import util
+
+
+def test_empty_positionRms():
+    ra = np.deg2rad(np.array([]))
+    dec = np.deg2rad(np.array([]))
+    ra_avg = np.mean(ra)
+    dec_avg = np.mean(dec)
+
+    obs = util.positionRms(ra_avg, dec_avg, ra, dec)
+
+    assert np.isnan(obs)
+
+
+def test_single_positionRms():
+    ra = np.deg2rad(np.array([10.0010]))
+    dec = np.deg2rad(np.array([20.001]))
+    ra_avg = np.mean(ra)
+    dec_avg = np.mean(dec)
+
+    exp = 0
+    obs = util.positionRms(ra_avg, dec_avg, ra, dec)
+
+    assert_array_max_ulp(obs, exp)
+
+
+def test_basic_positionRms():
+    ra = np.deg2rad(np.array([10.0010, 10.0005, 10.0000, 10.0005]))
+    dec = np.deg2rad(np.array([20.001, 20.006, 20.002, 20.004]))
+    ra_avg = np.mean(ra)
+    dec_avg = np.mean(dec)
+
+    exp = 0.0019488135998463505 * 3600 * 1000  # degrees * arcsec/degree * milliarcsec/arcsec
+    obs = util.positionRms(ra_avg, dec_avg, ra, dec)
+
+    assert_almost_equal(obs, exp, decimal=7)  # 1e-7 rad == 5.73e-6 deg == 36 milliarcsec
+
+
+def test_north_pole_positionRms():
+    ra = np.deg2rad(np.array([10.0010, 10.0005, 190.0000, 190.0005]))
+    dec = np.deg2rad(np.array([89.999, 89.998, 89.999, 89.998]))
+    ra_avg = np.mean(ra)
+    dec_avg = np.mean(dec)
+
+    exp = 0.0021794464685958273 * 3600 * 1000  # degrees * arcsec/degree * milliarcsec/arcsec
+    obs = util.positionRms(ra_avg, dec_avg, ra, dec)
+
+    assert_almost_equal(obs, exp, decimal=7)  # 1e-7 rad == 5.73e-6 deg == 36 milliarcsec
+
+
+def test_south_pole_positionRms():
+    ra = np.deg2rad(np.array([10.0010, 10.0005, 190.0000, 190.0005]))
+    dec = -np.deg2rad(np.array([89.999, 89.998, 89.999, 89.998]))
+    ra_avg = np.mean(ra)
+    dec_avg = np.mean(dec)
+
+    exp = 0.0021794464685958273 * 3600 * 1000  # degrees * arcsec/degree * milliarcsec/arcsec
+    obs = util.positionRms(ra_avg, dec_avg, ra, dec)
+
+    assert_almost_equal(obs, exp, decimal=7)  # 1e-7 rad == 5.73e-6 deg == 36 milliarcsec
+
+
+def test_wrap_positionRms():
+    ra = np.deg2rad(np.array([10.0010, 10.0005, 370.0000, 370.0005]))
+    dec = np.deg2rad(np.array([20.001, 20.006, 20.002, 20.004]))
+    ra_avg = np.mean(ra % (2*np.pi))
+    dec_avg = np.mean(dec)
+
+    exp = 0.0019488135998463505 * 3600 * 1000  # degrees * arcsec/degree * milliarcsec/arcsec
+    obs = util.positionRms(ra_avg, dec_avg, ra, dec)
+
+    assert_almost_equal(obs, exp, decimal=7)  # 1e-7 rad == 5.73e-6 deg == 36 milliarcsec
+
+
+if __name__ == "__main__":
+    lsst.utils.tests.init()
+    unittest.main()