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,19 +138,21 @@ def apply_to(self, obj, **kwargs):
         if isinstance(obj, FOVSetupBase):
 
             hdr = self.image_plane_header
-            x_mm, y_mm = calc_footprint(hdr, "D")
+            xy_mm = calc_footprint(hdr, "D")
             pixel_size = hdr["CDELT1D"]              # mm
             pixel_scale = kwargs.get("pixel_scale", self.meta["pixel_scale"])   # ["]
             pixel_scale = utils.from_currsys(pixel_scale)
-            x_sky = x_mm * pixel_scale / pixel_size  # x["] = x[mm] * ["] / [mm]
-            y_sky = y_mm * pixel_scale / pixel_size  # y["] = y[mm] * ["] / [mm]
 
-            obj.shrink(axis=["x", "y"], values=([min(x_sky), max(x_sky)],
-                                                [min(y_sky), max(y_sky)]))
-            obj.detector_limits = {"xd_min": min(x_mm),
-                                   "xd_max": max(x_mm),
-                                   "yd_min": min(y_mm),
-                                   "yd_max": max(y_mm)}
+            # x["] = x[mm] * ["] / [mm]
+            xy_sky = xy_mm * pixel_scale / pixel_size
+
+            obj.shrink(axis=["x", "y"],
+                       values=(tuple(zip(xy_sky.min(axis=0),
+                                         xy_sky.max(axis=0)))))
+
+            lims = np.array((xy_mm.min(axis=0), xy_mm.max(axis=0)))
+            keys = ["xd_min", "xd_max", "yd_min", "yd_max"]
+            obj.detector_limits = dict(zip(keys, lims.T.flatten()))
 
         return obj
 
@@ -163,13 +165,15 @@ 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_mm = calc_footprint(hdr, "D")
             pixel_size = hdr["CDELT1D"]              # mm
             pixel_scale = self.meta["pixel_scale"]   # ["]
-            x_sky = x_mm * pixel_scale / pixel_size  # x["] = x[mm] * ["] / [mm]
-            y_sky = y_mm * pixel_scale / pixel_size  # y["] = y[mm] * ["] / [mm]
 
-            aperture_mask = ApertureMask(array_dict={"x": x_sky, "y": y_sky},
+            # x["] = x[mm] * ["] / [mm]
+            xy_sky = xy_mm * pixel_scale / pixel_size
+
+            aperture_mask = ApertureMask(array_dict={"x": xy_sky[:, 0],
+                                                     "y": xy_sky[:, 1]},
                                          pixel_scale=pixel_scale)
 
         return aperture_mask
@@ -265,9 +269,10 @@ def plot(self, axes=None):
             _, axes = figure_factory()
 
         for hdr in self.detector_headers():
-            x_mm, y_mm = calc_footprint(hdr, "D")
-            axes.plot(list(close_loop(x_mm)), list(close_loop(y_mm)))
-            axes.text(*np.mean((x_mm, y_mm), axis=1), hdr["ID"],
+            xy_mm = calc_footprint(hdr, "D")
+            outline = np.array(list(close_loop(xy_mm)))
+            axes.plot(outline[:, 0], outline[:, 1])
+            axes.text(*xy_mm.mean(axis=0), hdr["ID"],
                       ha="center", va="center")
 
         axes.set_aspect("equal")

scopesim/effects/psf_utils.py

@@ -1,4 +1,6 @@
 import numpy as np
+import logging
+
 from scipy import ndimage as spi
 from scipy.interpolate import RectBivariateSpline, griddata
 from scipy.ndimage import zoom
@@ -95,11 +97,11 @@ def make_strehl_map_from_table(tbl, pixel_scale=1*u.arcsec):
                                          np.arange(-25, 26))).T,
                     method="nearest")
 
-    hdr = imp_utils.header_from_list_of_xy(np.array([-25, 25]) / 3600.,
-                                           np.array([-25, 25]) / 3600.,
-                                           pixel_scale=1/3600)
+    new_wcs, _ = imp_utils.create_wcs_from_points(np.array([[-25, -25],
+                                                            [25, 25]]),
+                                                  pixel_scale=1, arcsec=True)
 
-    map_hdu = fits.ImageHDU(header=hdr, data=smap)
+    map_hdu = fits.ImageHDU(header=new_wcs.to_header(), data=smap)
 
     return map_hdu
 
@@ -113,6 +115,7 @@ def rescale_kernel(image, scale_factor, spline_order=None):
 
     # Re-centre kernel
     im_shape = image.shape
+    # TODO: this might be another off-by-something
     dy, dx = np.divmod(np.argmax(image), im_shape[1]) - np.array(im_shape) // 2
     if dy > 0:
         image = image[2*dy:, :]
@@ -129,14 +132,23 @@ def rescale_kernel(image, scale_factor, spline_order=None):
     return image
 
 
-def cutout_kernel(image, fov_header):
+def cutout_kernel(image, fov_header, kernel_header=None):
+    from astropy.wcs import WCS
+
+    wk = WCS(kernel_header)
     h, w = image.shape
     xcen, ycen = 0.5 * w, 0.5 * h
+    xcen_w, ycen_w = wk.wcs_world2pix(np.array([[0., 0.]]), 0).squeeze().round(7)
+    if xcen != xcen_w or ycen != ycen_w:
+        logging.warning("PSF center off")
+
     dx = 0.5 * fov_header["NAXIS1"]
     dy = 0.5 * fov_header["NAXIS2"]
-    x0, x1 = max(0, int(xcen-dx)), min(w, int(xcen+dx))
-    y0, y1 = max(0, int(ycen-dy)), min(w, int(ycen+dy))
-    image_cutout = image[y0:y1, x0:x1]
+
+    # TODO: this is WET with imp_utils, somehow, I think
+    x0, x1 = max(0, np.floor(xcen-dx).astype(int)), min(w, np.ceil(xcen+dx).astype(int))
+    y0, y1 = max(0, np.floor(ycen-dy).astype(int)), min(w, np.ceil(ycen+dy).astype(int))
+    image_cutout = image[y0:y1+1, x0:x1+1]
 
     return image_cutout
 
@@ -145,9 +157,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/effects/psfs.py

@@ -23,6 +23,7 @@ def __init__(self, recursion_call=False):
         pixel_scale = utils.from_currsys("!INST.pixel_scale")
         self.header = {"CDELT1": pixel_scale / 3600.,
                        "CDELT2": pixel_scale / 3600.,
+                       "NAXIS": 2,
                        "NAXIS1": 128,
                        "NAXIS2": 128,
                        }
@@ -634,7 +635,8 @@ def get_kernel(self, fov):
 
                 if ((fov.header["NAXIS1"] < hdr["NAXIS1"]) or
                     (fov.header["NAXIS2"] < hdr["NAXIS2"])):
-                    self.kernel = pu.cutout_kernel(self.kernel, fov.header)
+                    self.kernel = pu.cutout_kernel(self.kernel, fov.header,
+                                                   kernel_header=hdr)
 
         return self.kernel
 
@@ -672,14 +674,16 @@ def make_psf_cube(self, fov):
 
         # We need linear interpolation to preserve positivity. Might think of
         # more elaborate positivity-preserving schemes.
+        # Note: According to some basic profiling, this line is one of the
+        #       single largest hits on performance.
         ipsf = RectBivariateSpline(np.arange(nyin), np.arange(nxin), psf,
                                    kx=1, ky=1)
 
         xcube, ycube = np.meshgrid(np.arange(nxpsf), np.arange(nypsf))
         cubewcs = WCS(naxis=2)
         cubewcs.wcs.ctype = ["LINEAR", "LINEAR"]
         cubewcs.wcs.crval = [0., 0.]
-        cubewcs.wcs.crpix = [nxpsf // 2, nypsf // 2]
+        cubewcs.wcs.crpix = [(nxpsf + 1) / 2, (nypsf + 1) / 2]
         cubewcs.wcs.cdelt = [fov_pixel_scale, fov_pixel_scale]
         cubewcs.wcs.cunit = [fov_pixel_unit, fov_pixel_unit]
 

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)
+        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)],
+        minmax = np.array((xy.min(axis=0), xy.max(axis=0)))
+        self._volume = {"xs": minmax[:, 0],
+                        "ys": minmax[:, 1],
                         "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,10 +149,10 @@ def generate_fovs_list(self):
                                                 [ys_min, ys_max],
                                                 pixel_scale=pixel_scale / 3600.)
 
-            x_sky, y_sky = ipu.calc_footprint(skyhdr)
-            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")
+            xy_sky = ipu.calc_footprint(skyhdr)
+            xy_det = xy_sky / plate_scale_deg
+            dethdr = ipu.header_from_list_of_xy(xy_det[:, 0], xy_det[:, 1],
+                                                pixel_size, "D")
             skyhdr.update(dethdr)
 
             # useful for spectroscopy mode where slit dimensions is not the same

scopesim/optics/fov_manager_utils.py

@@ -217,11 +217,11 @@ 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)
-        x_det = x_sky / plate_scale_deg
-        y_det = y_sky / plate_scale_deg
+        xy_sky = imp_utils.calc_footprint(skyhdr)
+        xy_det = xy_sky / plate_scale_deg
 
-        dethdr = imp_utils.header_from_list_of_xy(x_det, y_det, pixel_size, "D")
+        dethdr = imp_utils.header_from_list_of_xy(xy_det[:, 0], xy_det[:, 1],
+                                                  pixel_size, "D")
         skyhdr.update(dethdr)
         yield skyhdr