Fix off-by-one Error and related WCS changes

pyproject.toml

@@ -90,4 +90,9 @@ filterwarnings = [
     "ignore:Cannot merge meta key.*:astropy.utils.metadata.MergeConflictWarning",
     # Raised when saving fits files, not so important to fix:
     "ignore:.*a HIERARCH card will be created.*:astropy.io.fits.verify.VerifyWarning",
+    # Web-related issues, fix at some point
+    "ignore:'urllib3.contrib.pyopenssl'*:DeprecationWarning",
+    "ignore:Blowfish*:UserWarning",
+    # Not sure what that is but it's everywhere...
+    "ignore:'cgi'*:DeprecationWarning",
 ]

scopesim/effects/detector_list.py

@@ -138,7 +138,8 @@ def apply_to(self, obj, **kwargs):
         if isinstance(obj, FOVSetupBase):
 
             hdr = self.image_plane_header
-            x_mm, y_mm = calc_footprint(hdr, "D")
+            xy = calc_footprint(hdr, "D")
+            x_mm, y_mm = xy[:, 0], xy[:, 1]
             pixel_size = hdr["CDELT1D"]              # mm
             pixel_scale = kwargs.get("pixel_scale", self.meta["pixel_scale"])   # ["]
             pixel_scale = utils.from_currsys(pixel_scale)
@@ -163,7 +164,8 @@ def fov_grid(self, which="edges", **kwargs):
             self.meta = utils.from_currsys(self.meta)
 
             hdr = self.image_plane_header
-            x_mm, y_mm = calc_footprint(hdr, "D")
+            xy = calc_footprint(hdr, "D")
+            x_mm, y_mm = xy[:, 0], xy[:, 1]
             pixel_size = hdr["CDELT1D"]              # mm
             pixel_scale = self.meta["pixel_scale"]   # ["]
             x_sky = x_mm * pixel_scale / pixel_size  # x["] = x[mm] * ["] / [mm]
@@ -265,7 +267,8 @@ def plot(self, axes=None):
             _, axes = figure_factory()
 
         for hdr in self.detector_headers():
-            x_mm, y_mm = calc_footprint(hdr, "D")
+            xy = calc_footprint(hdr, "D")
+            x_mm, y_mm = xy[:, 0], xy[:, 1]
             axes.plot(list(close_loop(x_mm)), list(close_loop(y_mm)))
             axes.text(*np.mean((x_mm, y_mm), axis=1), hdr["ID"],
                       ha="center", va="center")

scopesim/effects/psf_utils.py

@@ -145,9 +145,8 @@ def get_strehl_cutout(fov_header, strehl_imagehdu):
 
     image = np.zeros((fov_header["NAXIS2"], fov_header["NAXIS1"]))
     canvas_hdu = fits.ImageHDU(header=fov_header, data=image)
-    canvas_hdu = imp_utils.add_imagehdu_to_imagehdu(strehl_imagehdu,
-                                                    canvas_hdu, spline_order=0,
-                                                    conserve_flux=False)
+    canvas_hdu = imp_utils.add_imagehdu_to_imagehdu(
+        strehl_imagehdu, canvas_hdu, spline_order=0, conserve_flux=False)
     canvas_hdu.data = canvas_hdu.data.astype(int)
 
     return canvas_hdu

scopesim/optics/fov.py

@@ -71,6 +71,7 @@ def __init__(self, header, waverange, detector_header=None, **kwargs):
         self.header["NAXIS1"] = header["NAXIS1"]
         self.header["NAXIS2"] = header["NAXIS2"]
         self.header.update(header)
+        self._ensure_deg_header()
         self.detector_header = detector_header
         self.hdu = None
 
@@ -599,13 +600,15 @@ def make_cube_hdu(self):
         return canvas_cube_hdu      # [ph s-1 AA-1 (arcsec-2)]
 
     def volume(self, wcs_prefix=""):
-        xs, ys = imp_utils.calc_footprint(self.header, wcs_suffix=wcs_prefix)
+        xy = imp_utils.calc_footprint(self.header, wcs_suffix=wcs_prefix)
+        xs, ys = xy[:, 0], xy[:, 1]
+        unit = self.header[f"CUNIT1{wcs_prefix}"].lower()
         # FIXME: This is unused!!
         # wave_corners = self.waverange
         self._volume = {"xs": [min(xs), max(xs)],
                         "ys": [min(ys), max(ys)],
                         "waves": self.waverange,
-                        "xy_unit": "mm" if wcs_prefix == "D" else "deg",
+                        "xy_unit": unit,
                         "wave_unit": "um"}
         return self._volume
 
@@ -625,6 +628,10 @@ def data(self):
     @property
     def corners(self):
         """Return sky footprint, image plane footprint."""
+        # Couldn't find where this is ever used, put warning here just in case
+        logging.warning("calc_footprint has been updated, any code that "
+                        "relies on this .corners property must likely be "
+                        "adapted as well.")
         sky_corners = imp_utils.calc_footprint(self.header)
         imp_corners = imp_utils.calc_footprint(self.header, "D")
         return sky_corners, imp_corners
@@ -698,6 +705,17 @@ def background_fields(self):
                 if isinstance(field, fits.ImageHDU)
                 and field.header.get("BG_SRC", False)]
 
+    def _ensure_deg_header(self):
+        cunit = u.Unit(self.header["CUNIT1"].lower())
+        convf = cunit.to(u.deg)
+        self.header["CDELT1"] *= convf
+        self.header["CDELT2"] *= convf
+        self.header["CRVAL1"] *= convf
+        self.header["CRVAL2"] *= convf
+        self.header["CUNIT1"] = "deg"
+        self.header["CUNIT2"] = "deg"
+
+
     def __repr__(self):
         waverange = [self.meta["wave_min"].value, self.meta["wave_max"].value]
         msg = (f"{self.__class__.__name__}({self.header!r}, {waverange!r}, "

scopesim/optics/fov_manager.py

@@ -149,7 +149,8 @@ def generate_fovs_list(self):
                                                 [ys_min, ys_max],
                                                 pixel_scale=pixel_scale / 3600.)
 
-            x_sky, y_sky = ipu.calc_footprint(skyhdr)
+            xy = ipu.calc_footprint(skyhdr)
+            x_sky, y_sky = xy[:, 0], xy[:, 1]
             x_det = x_sky / plate_scale_deg
             y_det = y_sky / plate_scale_deg
             dethdr = ipu.header_from_list_of_xy(x_det, y_det, pixel_size, "D")

scopesim/optics/fov_manager_utils.py

@@ -217,7 +217,8 @@ def _get_sky_hdrs(hdrs):
     pixel_size = pixel_scale / plate_scale
     plate_scale_deg = plate_scale / 3600.   # ["/mm] / 3600 = [deg/mm]
     for skyhdr in sky_hdrs:
-        x_sky, y_sky = imp_utils.calc_footprint(skyhdr)
+        xy = imp_utils.calc_footprint(skyhdr)
+        x_sky, y_sky = xy[:, 0], xy[:, 1]
         x_det = x_sky / plate_scale_deg
         y_det = y_sky / plate_scale_deg
 

scopesim/optics/fov_utils.py

@@ -6,8 +6,8 @@
 from astropy.table import Table, Column
 from synphot import SourceSpectrum, Empirical1D
 
-from scopesim import utils
-from scopesim.optics import image_plane_utils as imp_utils
+from .. import utils
+from . import image_plane_utils as imp_utils
 
 
 def is_field_in_fov(fov_header, field, wcs_suffix=""):
@@ -38,16 +38,21 @@ def is_field_in_fov(fov_header, field, wcs_suffix=""):
             y = list(utils.quantity_from_table("y", field,
                                                u.arcsec).to(u.deg).value)
             s = wcs_suffix
-            cdelt = utils.quantify(fov_header["CDELT1" + s], u.deg).value
+            # cdelt = utils.quantify(fov_header["CDELT1" + s], u.deg).value
+            cdelt = fov_header[f"CDELT1{s}"] * u.Unit(fov_header[f"CUNIT1{s}"]).to(u.deg)
             field_header = imp_utils.header_from_list_of_xy(x, y, cdelt, s)
         elif isinstance(field, (fits.ImageHDU, fits.PrimaryHDU)):
             field_header = field.header
         else:
             logging.warning("Input was neither Table nor ImageHDU: %s", field)
             return False
 
-        ext_xsky, ext_ysky = imp_utils.calc_footprint(field_header, wcs_suffix)
-        fov_xsky, fov_ysky = imp_utils.calc_footprint(fov_header, wcs_suffix)
+        xy = imp_utils.calc_footprint(field_header, wcs_suffix)
+        ext_xsky, ext_ysky = xy[:, 0], xy[:, 1]
+        xy = imp_utils.calc_footprint(fov_header, wcs_suffix)
+        fov_xsky, fov_ysky = xy[:, 0], xy[:, 1]
+        fov_xsky *= u.Unit(fov_header["CUNIT1"].lower()).to(u.deg)
+        fov_ysky *= u.Unit(fov_header["CUNIT2"].lower()).to(u.deg)
 
         is_inside_fov = (min(ext_xsky) < max(fov_xsky) and
                          max(ext_xsky) > min(fov_xsky) and
@@ -91,7 +96,8 @@ def combine_table_fields(fov_header, src, field_indexes):
     tbl : Table
 
     """
-    fov_xsky, fov_ysky = imp_utils.calc_footprint(fov_header)
+    xy = imp_utils.calc_footprint(fov_header)
+    fov_xsky, fov_ysky = xy[:, 0], xy[:, 1]
 
     x, y, ref, weight = [], [], [], []
 
@@ -282,8 +288,16 @@ def extract_area_from_imagehdu(imagehdu, fov_volume):
     hdr = imagehdu.header
     new_hdr = {}
     naxis1, naxis2 = hdr["NAXIS1"], hdr["NAXIS2"]
-    x_hdu, y_hdu = imp_utils.calc_footprint(imagehdu)  # field edges in "deg"
-    x_fov, y_fov = fov_volume["xs"], fov_volume["ys"]
+    xy = imp_utils.calc_footprint(imagehdu)  # field edges in "deg"
+    x_hdu, y_hdu = xy[:, 0], xy[:, 1]
+    x_fov, y_fov = np.array(fov_volume["xs"]), np.array(fov_volume["ys"])
+
+    if hdr["CUNIT1"] == "arcsec":
+        x_hdu *= u.arcsec.to(u.deg)
+        y_hdu *= u.arcsec.to(u.deg)
+    if fov_volume["xy_unit"] == "arcsec":
+        x_fov *= u.arcsec.to(u.deg)
+        y_fov *= u.arcsec.to(u.deg)
 
     x0s, x1s = max(min(x_hdu), min(x_fov)), min(max(x_hdu), max(x_fov))
     y0s, y1s = max(min(y_hdu), min(y_fov)), min(max(y_hdu), max(y_fov))

scopesim/optics/image_plane.py

@@ -3,6 +3,7 @@
 
 from astropy.io import fits
 from astropy.table import Table
+from astropy.wcs import WCS
 
 from .image_plane_utils import add_table_to_imagehdu, add_imagehdu_to_imagehdu
 
@@ -55,9 +56,13 @@ def __init__(self, header, **kwargs):
             raise ValueError(f"header must have a valid image-plane WCS: "
                              f"{dict(header)}")
 
-        image = np.zeros((header["NAXIS2"]+1, header["NAXIS1"]+1))
+        # image = np.zeros((header["NAXIS2"]+1, header["NAXIS1"]+1))
+        image = np.zeros((header["NAXIS2"], header["NAXIS1"]))
         self.hdu = fits.ImageHDU(data=image, header=header)
 
+        self._det_wcs = self._get_wcs(header, "D")
+        self._sky_wcs = self._get_wcs(header, " ")
+
     def add(self, hdus_or_tables, sub_pixel=None, spline_order=None,
             wcs_suffix=""):
         """
@@ -137,3 +142,17 @@ def image(self):
     def view(self, sub_pixel):
         # for consistency with FieldOfView
         return self.data
+
+    @staticmethod
+    def _get_wcs(header: fits.Header, key: str) -> WCS:
+        sky_alias = {" ", "S"}
+        try:
+            wcs = WCS(header, key=key)
+        except KeyError:
+            # retry with alias
+            sky_alias.discard(key)
+            try:
+                wcs = WCS(header, key=sky_alias.pop())
+            except KeyError:
+                wcs = None
+        return wcs