Merge pull request #284 from AstarVienna/dev_master

Release 0.7.0

.github/workflows/tests.yml

@@ -11,6 +11,9 @@ on:
       - master
       - dev_master
       - dev_spectroscopy
+  schedule:
+    - # Run every day at 5:00 UTC
+    - cron: "0 5 * * *"
 
   # Allows you to run this workflow manually from the Actions tab
   workflow_dispatch:

pyproject.toml

@@ -1,6 +1,6 @@
 [project]
 name = "ScopeSim"
-version = "0.6.2"
+version = "0.7.0"
 description = "Generalised telescope observation simulator"
 readme = "README.md"
 requires-python = ">=3.8"
@@ -20,10 +20,11 @@ classifiers=[
     "Topic :: Scientific/Engineering :: Astronomy",
 ]
 dependencies = [
-    "numpy>=1.19",
-    "scipy>=1.4.0",
+    "numpy>=1.19.5",
+    "scipy>=1.10.0",
     "astropy>=5.0",
     "matplotlib>=3.2.0",
+    "pooch>=1.7.0",  # for scipy.datasets
 
     "docutils>=0.15",
     "requests>=2.28.2",
@@ -35,7 +36,7 @@ dependencies = [
     "requests-cache>1.0",
 
     "synphot>=1.1.0",
-    "skycalc_ipy>=0.1.3",
+    "skycalc_ipy>=0.1.5",
     "anisocado>=0.3.0",
 ]
 
@@ -77,3 +78,21 @@ addopts = "--strict-markers"
 markers = [
     "webtest: marks tests as requiring network (deselect with '-m \"not webtest\"')",
 ]
+filterwarnings = [
+    "error",
+    # Should probably be fixed:
+    "ignore::ResourceWarning",
+    "ignore:The fit may be poorly conditioned:astropy.utils.exceptions.AstropyUserWarning",
+    # Perhaps fix?
+    "ignore:divide by zero encountered in scalar divide:RuntimeWarning",
+    "ignore:divide by zero encountered in double_scalars:RuntimeWarning",
+    "ignore:invalid value encountered in multiply:RuntimeWarning",
+    "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
@@ -283,7 +294,7 @@ def get_bkg_level(obj, bg_w):
         if bg_w == 0:
             bkg_level = 0
         else:
-            mask = np.zeros_like(obj, dtype=np.bool8)
+            mask = np.zeros_like(obj, dtype=bool)
             if bg_w > 0:
                 mask[bg_w:-bg_w, bg_w:-bg_w] = True
             bkg_level = np.ma.median(np.ma.masked_array(obj, mask=mask))
@@ -292,7 +303,7 @@ def get_bkg_level(obj, bg_w):
         if bg_w == 0:
             bkg_level = np.array([0] * obj.shape[0])
         else:
-            mask = np.zeros_like(obj, dtype=np.bool8)
+            mask = np.zeros_like(obj, dtype=bool)
             if bg_w > 0:
                 mask[:, bg_w:-bg_w, bg_w:-bg_w] = True
             bkg_level = np.ma.median(np.ma.masked_array(obj, mask=mask),

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
 
@@ -513,7 +514,7 @@ def make_cube_hdu(self):
         3. Find Image fields (see ``FieldOfView._make_cube_imagefields()``).
         4. Find Table fields (see ``FieldOfView._make_cube_tablefields()``).
 
-        ``PHOTLAM = ph/s/m2/um``.
+        ``PHOTLAM = ph / (s * m2 * um)``.
         Original source fields are in units of:
 
         - tables: (PHOTLAM in spectrum)
@@ -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