Merge pull request #74 from AstarVienna/hb/fix_source_from_imagehdu

Fix scaling of spectrum

docs/notebooks/starting.ipynb

@@ -60,8 +60,8 @@
     "hdr = fits.Header(dict(NAXIS=2,\n",
     "                       NAXIS1=img.shape[0]+1,\n",
     "                       NAXIS2=img.shape[1]+1,\n",
-    "                       CRPIX1=0,\n",
-    "                       CRPIX2=0,\n",
+    "                       CRPIX1=img.shape[0] / 2,\n",
+    "                       CRPIX2=img.shape[1] / 2,\n",
     "                       CRVAL1=0,\n",
     "                       CRVAL2=0,\n",
     "                       CDELT1=0.2/3600,\n",
@@ -227,7 +227,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.8.12"
+   "version": "3.11.5"
   }
  },
  "nbformat": 4,

docs/starting.rst

@@ -28,8 +28,8 @@ user
     hdr = fits.Header(dict(NAXIS=2,
                            NAXIS1=img.shape[0]+1,
                            NAXIS2=img.shape[1]+1,
-                           CRPIX1=0,
-                           CRPIX2=0,
+                           CRPIX1=(img.shape[0] + 1) / 2,
+                           CRPIX2=(img.shape[1] + 1) / 2,
                            CRVAL1=0,
                            CRVAL2=0,
                            CDELT1=0.2/3600,

scopesim_templates/micado/spectral_calibrations.py

@@ -75,6 +75,7 @@ def line_list(unit_flux=1*PHOTLAM,
         flux = signal.convolve(flux, kernel, mode="same")
 
     # make spectrum
+    assert len(wave) > 3, "At least 4 wavelength points are required."
     spec = SourceSpectrum(Empirical1D, points=wave * u.nm,
                           lookup_table=flux * unit_flux)
 

scopesim_templates/misc/misc.py

@@ -193,10 +193,12 @@ def source_from_imagehdu(image_hdu, filter_name, pixel_unit_amplitude=None,
 
     if filter_name is not None and isinstance(filter_name, str):
         filter_curve = tu.get_filter(filter_name)
-        waves = filter_curve.waverange
+        waverange = filter_curve.waverange
 
-    spec = SourceSpectrum(Empirical1D, points=waves, lookup_table=[1, 1])
-    tu.scale_spectrum(spec, filter_name=filter_name, amplitude=amp_unit)
+    waves = np.linspace(waverange[0], waverange[1], num=10)
+    lookup_table = [1,] * len(waves)
+    spec = SourceSpectrum(Empirical1D, points=waves, lookup_table=lookup_table)
+    spec = tu.scale_spectrum(spec, filter_name=filter_name, amplitude=amp_unit)
 
     if image_hdu.header.get("SPEC_REF") is None:
         image_hdu.header["SPEC_REF"] = 0
@@ -430,6 +432,7 @@ def poorman_cube_source(filename=None, hdu=None, ext=1, pixel_scale=None,
     zpix = np.arange(specwcs.spectral.array_shape[0])
     # keeping wavelengths in native coordinates
     waves = specwcs.pixel_to_world(zpix)
+    assert len(waves) > 3, "At least 4 wavelength points are required."
     zero_spec = np.zeros(shape=zpix.shape) * bunit
 
     hdus = []

scopesim_templates/tests/pyobjects/source_objects.py

@@ -1,3 +1,5 @@
+import pytest
+
 import numpy as np
 
 from astropy.io import fits
@@ -26,6 +28,34 @@ def _basic_imagehdu(n=11):
     return hdu
 
 
+@pytest.fixture(scope="module")
+def starting_star_imagehdu():
+    """Create the same star as in starting.ipynb"""
+    n = 100
+    sigma = 5
+    x, y = np.meshgrid(np.arange(n), np.arange(n))
+    img = np.exp(-1 * (((x - n / 2) / sigma) ** 2 + ((y - n / 2) / sigma) ** 2))
+
+    # Fits headers of the image. Yes it needs a WCS
+    hdr = fits.Header(dict(NAXIS=2,
+                           NAXIS1=n,
+                           NAXIS2=n,
+                           CRPIX1=(n + 1) / 2,
+                           CRPIX2=(n + 1) / 2,
+                           CRVAL1=0,
+                           CRVAL2=0,
+                           CDELT1=0.2 / 3600,
+                           CDELT2=0.2 / 3600,
+                           CUNIT1="DEG",
+                           CUNIT2="DEG",
+                           CTYPE1='RA---TAN',
+                           CTYPE2='DEC--TAN'))
+
+    # Creating an ImageHDU object
+    hdu = fits.ImageHDU(data=img, header=hdr)
+    return hdu
+
+
 def _make_dummy_cube(scale, wave_unit, ref_wave, wave_step, wave_type, bunit):
 
     if isinstance(scale, u.Quantity) is False:

scopesim_templates/tests/test_misc/test_misc.py

@@ -4,13 +4,19 @@
 from pytest import raises
 from astropy import units as u
 from astropy.io import fits
+import numpy as np
+from matplotlib.colors import LogNorm
+from scopesim.utils import figure_factory
 
 from scopesim_templates.misc import misc
 from scopesim_templates.tests.pyobjects import source_objects as so
-from scopesim_templates.rc import Source
+from scopesim_templates.tests.pyobjects.source_objects import starting_star_imagehdu
+from scopesim_templates.rc import Source, load_example_optical_train
 
 METIS_FILTER_PATH = r"F:/Work/irdb/METIS/filters/TC_filter_H2O-ice.dat"
 
+PLOTS = False
+
 
 class TestSourceFromImageHDU:
     def test_initialises_for_basic_imagehdu(self):
@@ -51,6 +57,92 @@ def test_is_happy_with_metis_filter_and_bunit(self):
 
         assert isinstance(src, Source)
 
+    def test_actually_produces_image(self, starting_star_imagehdu):
+        src1 = misc.source_from_imagehdu(
+            starting_star_imagehdu,
+            filter_name="J",
+            pixel_unit_amplitude=10e12 * u.Jy)
+
+        width_height = 4096
+        opt = load_example_optical_train()
+        opt['psf'].include = False
+        opt.cmds["!OBS.psf_fwhm"] = 0.01
+        opt.cmds["!TEL.area"] = 1000 * u.m**2
+        opt.cmds["!INST.pixel_scale"] = 0.004
+        opt.cmds["!INST.plate_scale"] = 0.4
+        opt.cmds["!DET.width"] = width_height
+        opt.cmds["!DET.height"] = width_height
+        opt.cmds["!DET.dit"] = 30
+        opt.cmds["!DET.ndit"] = 120
+
+        opt["source_fits_keywords"].include = False
+
+        opt.observe(src1)
+        hdul = opt.readout()[0]
+
+        data = hdul[1].data
+        # Is the background okay?
+        assert 500 < np.median(data) < 5000
+        # Is there a bright source?
+        assert data.max() > 10000000
+        # Is the bright source approximately in the center.
+        x_cen, y_cen = np.unravel_index(data.argmax(), data.shape)
+        # TODO: Figure out why source is not in the center!
+        fudge = 100
+        assert width_height / 2 - fudge < x_cen < width_height / 2 + fudge
+        assert width_height / 2 - fudge < y_cen < width_height / 2 + fudge
+
+        if PLOTS:
+            fig, ax = figure_factory(3, 1)
+            ax[0].imshow(src1.fields[0].data)
+            ax[1].imshow(opt.image_planes[0].data, norm=LogNorm())
+            ax[2].imshow(hdul[1].data, norm=LogNorm())
+            fig.show()
+
+    def test_scales(self, starting_star_imagehdu):
+        """Test whether source_from_imagehdu scales w/ pixel_unit_amplitude."""
+        flux1 = 10e12 * u.Jy
+        src1 = misc.source_from_imagehdu(
+            starting_star_imagehdu,
+            filter_name="J",
+            pixel_unit_amplitude=flux1)
+
+        src2 = misc.source_from_imagehdu(
+            starting_star_imagehdu,
+            filter_name="J",
+            pixel_unit_amplitude=flux1 * 10)
+
+        width_height = 4096
+        opt = load_example_optical_train()
+        opt['psf'].include = False
+        opt.cmds["!OBS.psf_fwhm"] = 0.01
+        opt.cmds["!TEL.area"] = 1000 * u.m**2
+        opt.cmds["!INST.pixel_scale"] = 0.004
+        opt.cmds["!INST.plate_scale"] = 0.4
+        opt.cmds["!DET.width"] = width_height
+        opt.cmds["!DET.height"] = width_height
+        opt.cmds["!DET.dit"] = 30
+        opt.cmds["!DET.ndit"] = 120
+
+        opt["source_fits_keywords"].include = False
+
+        opt.observe(src1)
+        hdul1 = opt.readout()[0]
+        data1 = hdul1[1].data
+
+        opt.observe(src2, update=True)
+        hdul2 = opt.readout()[0]
+        data2 = hdul2[1].data
+
+        if PLOTS:
+            fig, ax = figure_factory(1, 1)
+            ax.imshow(data2 / data1)
+            fig.show()
+
+        max1 = data1.max()
+        max2 = data2.max()
+        assert 9 < max2 / max1 < 11
+
 
 class TestPointSource:
     def test_initialize_from_spextra(self):
@@ -59,7 +151,7 @@ def test_initialize_from_spextra(self):
         assert isinstance(src, Source)
 
     def test_initialize_from_synphot(self):
-        sp = synphot.SourceSpectrum(synphot.Empirical1D, points=[1000, 10000], lookup_table=[1, 1])
+        sp = synphot.SourceSpectrum(synphot.Empirical1D, points=[1000] * 4, lookup_table=[1] * 4)
         src = misc.point_source(sed=sp)
         assert isinstance(src, Source)
 
@@ -71,7 +163,7 @@ def test_initialize_from_spextra(self):
         assert isinstance(src, Source)
 
     def test_initialize_from_synphot(self):
-        sp = synphot.SourceSpectrum(synphot.Empirical1D, points=[1000, 10000], lookup_table=[1, 1])
+        sp = synphot.SourceSpectrum(synphot.Empirical1D, points=[1000] * 4, lookup_table=[1] * 4)
         src = misc.uniform_source(sed=sp)
         assert isinstance(src, Source)
 

scopesim_templates/utils/general_utils.py

@@ -24,6 +24,7 @@ def hdu_to_synphot(hdu):
     flux = hdu.data["flux"]
     flux_unit = u.Unit(hdu.header["TUNIT2"])
 
+    assert len(wave) > 3, "At least 4 wavelength points are required."
     spec = SourceSpectrum(Empirical1D, points=wave*wave_unit,
                           lookup_table=flux*flux_unit)