Merge pull request #242 from mwcraig/temp-fix-notebook-2

Fix bug in catalog search functionality and make comparison notebook compatible with `CatalogData`

stellarphot/core.py

@@ -1089,7 +1089,7 @@ def _tidy_vizier_catalog(data, mag_column_regex, color_column_regex):
     @classmethod
     def from_vizier(
         cls,
-        header_or_center,
+        field_center,
         desired_catalog,
         radius=0.5 * u.degree,
         clip_by_frame=False,
@@ -1106,10 +1106,11 @@ def from_vizier(
         Parameters
         ----------
 
-        header_or_center : FITS header or `astropy.coordinates.SkyCoord`
-            Either a FITS header with WCS information or a `SkyCoord` object.
-            The center of the frame or the input coordinate is the center
-            of the cone search.
+        field_center : `astropy.coordinates.SkyCoord`, `astropy.wcs.WCS`, or FITS header
+            Either a `~astropy.coordinates.SkyCoord` object, a `~astropy.wcs.WCS` object
+            or a FITS header with WCS information. The input coordinate should be the
+            center of the frame; if a header or WCS is the input then the center of the
+            frame will be determined from the WCS.
 
         desired_catalog : str
             Vizier name of the catalog to be searched.
@@ -1164,18 +1165,28 @@ def from_vizier(
         followed by a letter or letters.
         """
 
-        if isinstance(header_or_center, SkyCoord):
+        if isinstance(field_center, SkyCoord):
             # Center was passed in, just use it.
-            center = header_or_center
+            center = field_center
             if clip_by_frame:
                 raise ValueError(
                     "To clip entries by frame you must use "
                     "a WCS as the first argument."
                 )
+        elif isinstance(field_center, WCS):
+            center = SkyCoord(*field_center.wcs.crval, unit="deg")
         else:
-            # Find the center of the frame
-            shape = (header_or_center["NAXIS2"], header_or_center["NAXIS1"])
-            center = WCS(header_or_center).pixel_to_world(shape[1] / 2, shape[0] / 2)
+            wcs = WCS(field_center)
+            # The header may not have contained WCS information. In that case
+            # the WCS CTYPE will be empty strings and we need to raise an
+            # error.
+            if wcs.wcs.ctype[0] == "" and wcs.wcs.ctype[1] == "":
+                raise ValueError(
+                    f"Invalid coordinates in input {field_center}. Make sure the "
+                    "header contains valid WCS information or pass in a WCS or "
+                    "coordinate."
+                )
+            center = SkyCoord(*wcs.wcs.crval, unit="deg")
 
         # Get catalog via cone search
         Vizier.ROW_LIMIT = -1  # Set row_limit to have no limit
@@ -1213,7 +1224,7 @@ def from_vizier(
         # desired.
         if clip_by_frame:
             cat_coords = SkyCoord(ra=cat["ra"], dec=cat["dec"])
-            wcs = WCS(header_or_center)
+            wcs = WCS(field_center)
             x, y = wcs.all_world2pix(cat_coords.ra, cat_coords.dec, 0)
             in_x = (x >= padding) & (x <= wcs.pixel_shape[0] - padding)
             in_y = (y >= padding) & (y <= wcs.pixel_shape[1] - padding)
@@ -1223,17 +1234,18 @@ def from_vizier(
         return cat
 
 
-def apass_dr9(header_or_center, radius=1 * u.degree, clip_by_frame=False, padding=100):
+def apass_dr9(field_center, radius=1 * u.degree, clip_by_frame=False, padding=100):
     """
     Return the items from APASS DR9 that are within the search radius and
     (optionally) within the field of view of a frame.
 
     Parameters
     ----------
-    header_or_center : FITS header or `astropy.coordinates.SkyCoord`
-        Either a FITS header with WCS information or a `SkyCoord` object.
-        The center of the frame or the input coordinate is the center
-        of the cone search.
+    field_center : `astropy.coordinates.SkyCoord`, `astropy.wcs.WCS`, or FITS header
+        Either a `~astropy.coordinates.SkyCoord` object, a `~astropy.wcs.WCS` object
+        or a FITS header with WCS information. The input coordinate should be the
+        center of the frame; if a header or WCS is the input then the center of the
+        frame will be determined from the WCS.
 
     radius : `astropy.units.Quantity`, optional
         Radius around which to search.
@@ -1247,14 +1259,19 @@ def apass_dr9(header_or_center, radius=1 * u.degree, clip_by_frame=False, paddin
         of the frame to be considered in the field of view. Default value
         is 100.
 
+    Returns
+    -------
+
+    `stellarphot.CatalogData`
+        Table of catalog information.
     """
     apass_colnames = {
         "recno": "id",  # There is no APASS ID, this is the one generated by Vizier
         "RAJ2000": "ra",
         "DEJ2000": "dec",
     }
     return CatalogData.from_vizier(
-        header_or_center,
+        field_center,
         "II/336/apass9",
         radius=radius,
         clip_by_frame=clip_by_frame,
@@ -1263,17 +1280,18 @@ def apass_dr9(header_or_center, radius=1 * u.degree, clip_by_frame=False, paddin
     )
 
 
-def vsx_vizier(header_or_center, radius=1 * u.degree, clip_by_frame=False, padding=100):
+def vsx_vizier(field_center, radius=1 * u.degree, clip_by_frame=False, padding=100):
     """
     Return the items from the copy of VSX on Vizier that are within the search
     radius and (optionally) within the field of view of a frame.
 
     Parameters
     ----------
-    header_or_center : FITS header or `astropy.coordinates.SkyCoord`
-        Either a FITS header with WCS information or a `SkyCoord` object.
-        The center of the frame or the input coordinate is the center
-        of the cone search.
+    field_center : `astropy.coordinates.SkyCoord`, `astropy.wcs.WCS`, or FITS header
+        Either a `~astropy.coordinates.SkyCoord` object, a `~astropy.wcs.WCS` object
+        or a FITS header with WCS information. The input coordinate should be the
+        center of the frame; if a header or WCS is the input then the center of the
+        frame will be determined from the WCS.
 
     radius : `astropy.units.Quantity`, optional
         Radius around which to search.
@@ -1287,6 +1305,10 @@ def vsx_vizier(header_or_center, radius=1 * u.degree, clip_by_frame=False, paddi
         of the frame to be considered in the field of view. Default value
         is 100.
 
+    Returns
+    -------
+    `stellarphot.CatalogData`
+        Table of catalog information.
     """
     vsx_map = dict(
         Name="id",
@@ -1302,7 +1324,7 @@ def prepare_cat(cat):
         return cat
 
     return CatalogData.from_vizier(
-        header_or_center,
+        field_center,
         "B/vsx/vsx",
         radius=radius,
         clip_by_frame=clip_by_frame,

stellarphot/tests/test_core.py

@@ -957,10 +957,36 @@ def test_catalog_from_vizier_search_apass():
 
 
 @pytest.mark.remote_data
-def test_catalog_from_vizier_search_vsx():
+@pytest.mark.parametrize("location_method", ["coord", "header", "wcs"])
+def test_catalog_from_vizier_search_vsx(location_method):
     # Do a cone search with a small enough radius to return exaclty one star,
     # DQ Psc, which happens to already be in the test data.
     coordinate = SkyCoord(ra=359.94371 * u.deg, dec=-0.2801 * u.deg)
+
+    if location_method == "coord":
+        center = coordinate
+    else:
+        # We need a WCS for each of these methods, so make that first.
+        # The projection doesn't matter because we only need the center
+        # coordinate. The code below is more or less copied from
+        # https://docs.astropy.org/en/stable/wcs/example_create_imaging.html
+        wcs = WCS(naxis=2)
+        # Put the center in the right place
+        wcs.wcs.crpix = [100, 100]
+        wcs.wcs.crval = [coordinate.ra.degree, coordinate.dec.degree]
+        wcs.array_shape = [200, 200]
+
+        # The rest of these values shouldn't matter for this test
+        wcs.wcs.ctype = ["RA---AIR", "DEC--AIR"]
+        wcs.wcs.set_pv([(2, 1, 45.0)])
+        wcs.wcs.cdelt = np.array([-0.066667, 0.066667])
+        if location_method == "header":
+            center = wcs.to_header()
+            center["NAXIS1"] = wcs.array_shape[0]
+            center["NAXIS2"] = wcs.array_shape[1]
+        else:
+            center = wcs
+
     vsx_map = dict(
         Name="id",
         RAJ2000="ra",
@@ -975,7 +1001,7 @@ def prepare_cat(cat):
         return cat
 
     my_cat = CatalogData.from_vizier(
-        coordinate,
+        center,
         "B/vsx/vsx",
         radius=0.1 * u.arcmin,
         clip_by_frame=False,
@@ -1009,6 +1035,16 @@ def test_from_vizier_with_coord_and_frame_clip_fails():
         _ = CatalogData.from_vizier(cen_coord, "B/vsx/vsx", clip_by_frame=True)
 
 
+def test_from_vizier_with_header_no_wcs_raise_error():
+    # Check that calling from_vizier with a header that has no WCS
+    # information generates the appropriate error.
+    ccd = CCDData(data=np.ones((10, 10)), unit="adu")
+    header = ccd.to_hdu()[0].header
+
+    with pytest.raises(ValueError, match="Invalid coordinates in input"):
+        CatalogData.from_vizier(header, "B/vsx/vsx", clip_by_frame=True)
+
+
 @pytest.mark.remote_data
 @pytest.mark.parametrize(
     "clip, data_file",

stellarphot/utils/comparison_utils.py

@@ -72,13 +72,13 @@ def set_up(sample_image_for_finding_stars, directory_with_images="."):
 
     ccd = CCDData.read(path)
     try:
-        vsx = vsx_vizier(ccd.header, radius=0.5 * u.degree)
+        vsx = vsx_vizier(ccd.wcs, radius=0.5 * u.degree)
     except RuntimeError:
         vsx = []
     else:
-        ra = vsx["RAJ2000"]
-        dec = vsx["DEJ2000"]
-        vsx["coords"] = SkyCoord(ra=ra, dec=dec, unit=(u.hour, u.degree))
+        ra = vsx["ra"]
+        dec = vsx["dec"]
+        vsx["coords"] = SkyCoord(ra=ra, dec=dec, unit=u.degree)
 
     return ccd, vsx
 
@@ -103,7 +103,7 @@ def crossmatch_APASS2VSX(CCD, RD, vsx):
     Returns
     -------
 
-    apass : `astropy.table.Table`
+    apass : `stellarphot.CatalogData`
         Table with APASS stars in the field of view.
 
     v_angle : `astropy.units.Quantity`
@@ -112,10 +112,12 @@ def crossmatch_APASS2VSX(CCD, RD, vsx):
     RD_angle : `astropy.units.Quantity`
         Angular separation between APASS stars and input targets.
     """
-    apass = apass_dr9(CCD)
+    apass = apass_dr9(CCD.wcs)
+    # Use the standard names we have introduced for ra and dec
     ra = apass["ra"]
     dec = apass["dec"]
-    apass["coords"] = SkyCoord(ra=ra, dec=dec, unit=(u.hour, u.degree))
+    # Coordinate units are in degrees now
+    apass["coords"] = SkyCoord(ra=ra, dec=dec, unit=u.degree)
     apass_coord = apass["coords"]
 
     if vsx:
@@ -131,7 +133,9 @@ def crossmatch_APASS2VSX(CCD, RD, vsx):
     return apass, v_angle, RD_angle
 
 
-def mag_scale(cmag, apass, v_angle, RD_angle, brighter_dmag=0.44, dimmer_dmag=0.75):
+def mag_scale(
+    cmag, apass, v_angle, RD_angle, brighter_dmag=0.44, dimmer_dmag=0.75, passband="r"
+):
     """
     Select comparison stars that are 1) not close the VSX stars or to other
     target stars and 2) fall within a particular magnitude range.
@@ -157,6 +161,9 @@ def mag_scale(cmag, apass, v_angle, RD_angle, brighter_dmag=0.44, dimmer_dmag=0.
     dimmer_dmag : float, optional
         Minimum difference in magnitude between the target and comparison stars.
 
+    passband : str, optional
+        Passband to use for selecting the comparison stars.
+
     Returns
     -------
 
@@ -166,8 +173,9 @@ def mag_scale(cmag, apass, v_angle, RD_angle, brighter_dmag=0.44, dimmer_dmag=0.
     good_stars : `astropy.table.Table`
         Table with the comparison stars.
     """
-    high_mag = apass["r_mag"] < cmag + dimmer_dmag
-    low_mag = apass["r_mag"] > cmag - brighter_dmag
+    good_filter = apass["passband"] == passband
+    high_mag = apass["mag"] < cmag + dimmer_dmag
+    low_mag = apass["mag"] > cmag - brighter_dmag
     if len(v_angle) > 0:
         good_v_angle = v_angle > 1.0 * u.arcsec
     else:
@@ -178,7 +186,7 @@ def mag_scale(cmag, apass, v_angle, RD_angle, brighter_dmag=0.44, dimmer_dmag=0.
     else:
         good_RD_angle = True
 
-    good_stars = high_mag & low_mag & good_RD_angle & good_v_angle
+    good_stars = good_filter & high_mag & low_mag & good_RD_angle & good_v_angle
     good_apass = apass[good_stars]
     apass_good_coord = good_apass["coords"]
     return apass_good_coord, good_stars