Repository: AstarVienna/ScopeSim

Issue:
Weird behavior for trivial table header
Considering a minimal header from a table with one point source at (0, 0) and arbitrary CDELT:

```python
from scopesim.optics import image_plane_utils as impu
hdr = impu.header_from_list_of_xy([0], [0], 1)
```

then

```python
impu.calc_footprint(hdr)
```

results in

```python
array([[0.5, 0.5],
       [0.5, 0.5],
       [0.5, 0.5],
       [0.5, 0.5]])
```

This actually occurs frequently in the FOV for super simple sources with one point source in the center, but passes silently because the endresult is still the same. I'm wondering why we need a temporary header from the table fields at all in the particular case I'm currently looking at. Anyway, the example shown here should really print

```python
array([[0., 0.],
       [0., 0.],
       [0., 0.],
       [0., 0.]])
```

in my opinion. The bug seems to be in how the header is created for this edgecase.

But it gets worse. For

```python
hdr = impu.header_from_list_of_xy([-1, 1], [-1, 1], 10)
```

the resulting footprint is

```python
array([[5., 5.],
       [5., 5.],
       [5., 5.],
       [5., 5.]])
```

which seems even worse somehow...

To give an example where it actually works:

```python
hdr = impu.header_from_list_of_xy([-1, 1], [-1, 1], 1)
```

produces

```python
array([[-1., -1.],
       [-1.,  1.],
       [ 1.,  1.],
       [ 1., -1.]])
```

which to me seems correct.

Hints from issue comments:
@hugobuddel: I don't even understand the units. The help text says they can be deg or mm.

The `calc_footprint()` seems correctish, because the axis are empty:

```
from scopesim.optics import image_plane_utils as impu
hdr = impu.header_from_list_of_xy([0], [0], 1)
hdr
```
shows
```
NAXIS   =                    2                                                  
NAXIS1  =                    0                                                  
NAXIS2  =                    0                                                  
WCSAXES =                    2 / Number of coordinate axes                      
CRPIX1  =                  0.5 / Pixel coordinate of reference point            
CRPIX2  =                  0.5 / Pixel coordinate of reference point            
CDELT1  =                  1.0 / [deg] Coordinate increment at reference point  
CDELT2  =                  1.0 / [deg] Coordinate increment at reference point  
CUNIT1  = 'deg'                / Units of coordinate increment and value        
CUNIT2  = 'deg'                / Units of coordinate increment and value        
CTYPE1  = 'LINEAR'             / Coordinate type code                           
CTYPE2  = 'LINEAR'             / Coordinate type code                           
CRVAL1  =                  0.0 / [deg] Coordinate value at reference point      
CRVAL2  =                  0.0 / [deg] Coordinate value at reference point      
LATPOLE =                 90.0 / [deg] Native latitude of celestial pole        
MJDREF  =                  0.0 / [d] MJD of fiducial time                       
```
and indeed, without any pixels, the footprint is zero.

I don't know why the CRPIX is at 0.5 though.

Oh that makes sense, because look at your correct example:

```
hdr2 = impu.header_from_list_of_xy([-1, 1], [-1, 1], 1)
hdr2
```
gives
```
NAXIS   =                    2                                                  
NAXIS1  =                    2                                                  
NAXIS2  =                    2                                                  
WCSAXES =                    2 / Number of coordinate axes                      
CRPIX1  =                  1.5 / Pixel coordinate of reference point            
CRPIX2  =                  1.5 / Pixel coordinate of reference point            
CDELT1  =                  1.0 / [deg] Coordinate increment at reference point  
CDELT2  =                  1.0 / [deg] Coordinate increment at reference point  
CUNIT1  = 'deg'                / Units of coordinate increment and value        
CUNIT2  = 'deg'                / Units of coordinate increment and value        
CTYPE1  = 'LINEAR'             / Coordinate type code                           
CTYPE2  = 'LINEAR'             / Coordinate type code                           
CRVAL1  =                  0.0 / [deg] Coordinate value at reference point      
CRVAL2  =                  0.0 / [deg] Coordinate value at reference point      
LATPOLE =                 90.0 / [deg] Native latitude of celestial pole        
MJDREF  =                  0.0 / [d] MJD of fiducial time                       
```

There are 4 pixels, centered around the origin. The FITS counting is one-based, so the center of the image is at 1.5, 1.5.

Maybe we can resolve this by ensuring that all sources must be at least half a pixel from the edge? That is, in your first example, a source at 0,0 should have 1 pixel in it, which should be in the center. Or is that undesirable?

I don't know@teutoburg: > Maybe we can resolve this by ensuring that all sources must be at least half a pixel from the edge? That is, in your first example, a source at 0,0 should have 1 pixel in it, which should be in the center. Or is that undesirable?

I talked to @oczoske about this shortly after writing the issue. The formula used to calculate the `CRPIXn` is indeed correct, but breaks down for `NAXISn == 0`, so this case needs to be handled appropriately. Either by ensuring at least one full pixel is present (like you say), or by throwing an Exception (because this kind of thing shouldn't really happen) and in the cases where it's expected that this might occur, handle the exception accordingly. I haven't yet decided which option I prefer, or if there are any problems we might create by going with the "at least 1 px" approach.

Solution patch:
diff --git a/.github/dependabot.yml b/.github/dependabot.yml
index a3273c74..f8f74fbf 100644
--- a/.github/dependabot.yml
+++ b/.github/dependabot.yml
@@ -17,16 +17,26 @@ updates:
   - package-ecosystem: "pip"
     directory: "/"
     schedule:
-      interval: "weekly"
+      interval: "monthly"
     labels:
       - "dependencies"
+    ignore:
+      # For all packages, ignore all patch updates, the "version-update" only targets non-security related updates
+      - dependency-name: "*"
+        update-types: ["version-update:semver-patch"]
     groups:
-      minor-versions:
-        # The "version-update" only targets non-security related updates
+      internal-dependencies:
         applies-to: version-updates
         patterns:
-          - "*"
-        # Group minor and patch level update, keep separate PRs for major versions
+          - "astar*"
+          - "scopesim*"
+          - "spextra"
+          - "skycalc-ipy"
+          - "anisocado*"
+      main-dependencies:
+        # Applies to all dependencies defined in the main section of pyproject.toml
+        applies-to: version-updates
+        dependency-type: "production"
         update-types:
+          - "major"
           - "minor"
-          - "patch"
diff --git a/docs/source/5_liners/source_from_images.md b/docs/source/5_liners/source_from_images.md
index 01ab640e..0de22ff1 100644
--- a/docs/source/5_liners/source_from_images.md
+++ b/docs/source/5_liners/source_from_images.md
@@ -16,17 +16,17 @@ kernelspec:
 We can use a FITS image as the Source object for a ScopeSim Simulation
 
 **Warning: The simulation output is only as good as the input**
-   
+
    If the pixel scale of the input (`CDELTn`) is bigger than the pixel scale of the instrument, ScopeSim will simply interpolate the image.
-   
+
    Please don't expect wonders if the input image WCS information is not appropriate for the instrument you are using.
-   
-ScopeSim Source objects can be generated from fits.ImageHDU object in the following ways:  
-  
-1. Just an Image and scaling flux value  
+
+ScopeSim Source objects can be generated from fits.ImageHDU object in the following ways:
+
+1. Just an Image and scaling flux value
 2. An Image and the associated synphot Spectrum
 3. An Image and arrays for wavelength and flux
-  
+
 - ``image_hdu=<fits.ImageHDU>`` + ``flux=<astropy.Quantity>``
 - ``image_hdu=<fits.ImageHDU>`` + ``spectra=<list of synphot.SourceSpectrum>``
 - ``image_hdu=<fits.ImageHDU>`` + ``lam=<array>`` + ``spectra=<list of arrays>``
@@ -42,7 +42,7 @@ import matplotlib.pyplot as plt
 import scopesim
 
 # Make an ImageHDU with some pixel data
-hdu = fits.ImageHDU(data=scipy.misc.face(gray=True))
+hdu = fits.ImageHDU(data=scipy.datasets.face(gray=True))
 
 # Give the header some proper WCS info
 hdu.header.update({"CDELT1": 1, "CUNIT1": "arcsec", "CRPIX1": 0, "CRVAL1": 0,
@@ -61,7 +61,7 @@ def plot(src):
     plt.imshow(src.fields[0].data)
 ```
 
-## 1. Just an Image and scaling flux value 
+## 1. Just an Image and scaling flux value
 
 ``image_hdu=<fits.ImageHDU>`` + ``flux=<astropy.Quantity>``
 
diff --git a/scopesim/optics/image_plane_utils.py b/scopesim/optics/image_plane_utils.py
index 5f6e77ed..83d7356c 100644
--- a/scopesim/optics/image_plane_utils.py
+++ b/scopesim/optics/image_plane_utils.py
@@ -202,24 +202,31 @@ def create_wcs_from_points(points: np.ndarray,
                                         f"pix, got '{naxis.unit}' instead.")
         naxis = naxis.value
 
+    if (naxis == 0).any():
+        # Ensure at least one pixel in each axis.
+        logger.debug("NAXISn == 0, using minimum of 1.")
+        naxis[naxis == 0] = 1
+
     crpix = (naxis + 1) / 2
     crval = (points.min(axis=0) + points.max(axis=0)) / 2
 
-    ctype = "LINEAR" if wcs_suffix in "DX" else "RA---TAN"
+    # Cannot do `in "DX"` here because that would also match the empty string.
+    linsuff = {"D", "X"}
+    lintype = ["LINEAR", "LINEAR"]
+    # skytype = ["RA---TAN", "DEC--TAN"]  # ScopeSim can't handle the truth yet
+    skytype = ["LINEAR", "LINEAR"]
+    ctype = lintype if wcs_suffix in linsuff else skytype
 
     if isinstance(points, u.Quantity):
         cunit = points.unit
     else:
-        if wcs_suffix == "D":
-            cunit = "mm"
-        else:
-            cunit = "deg"
+        cunit = "mm" if wcs_suffix == "D" else "deg"
 
     if isinstance(pixel_scale, u.Quantity):
         pixel_scale = pixel_scale.value
 
     new_wcs = WCS(key=wcs_suffix)
-    new_wcs.wcs.ctype = 2 * [ctype]
+    new_wcs.wcs.ctype = ctype
     new_wcs.wcs.cunit = 2 * [cunit]
     new_wcs.wcs.cdelt = np.array(2 * [pixel_scale])
     new_wcs.wcs.crval = crval
@@ -916,6 +923,8 @@ def calc_footprint(header, wcs_suffix="", new_unit: str = None):
         xy0 = np.array([[0, 0], [0, y_ext], [x_ext, y_ext], [x_ext, 0]])
         xy1 = coords.wcs_pix2world(xy0, 0)
 
+    # FIXME: Catch case of unit mismatch in CUNIT1 and CUNIT2, there should be
+    #        a function for this somewhere but I forgot.
     if (cunit := coords.wcs.cunit[0]) == "deg":
         xy1 = _fix_360(xy1)
 


Test patch:
diff --git a/scopesim/tests/tests_optics/test_ImagePlaneUtils.py b/scopesim/tests/tests_optics/test_ImagePlaneUtils.py
index 0149b084..787442e8 100644
--- a/scopesim/tests/tests_optics/test_ImagePlaneUtils.py
+++ b/scopesim/tests/tests_optics/test_ImagePlaneUtils.py
@@ -347,3 +347,26 @@ class TestSkyDetWCS:
         # Compare header representations because of missing NAXIS info in new
         assert new_wcs.to_header() == sky_wcs.to_header()
         assert all(nax == (hdu.header["NAXIS1"], hdu.header["NAXIS2"]))
+
+
+class TestWCSFromMinimalPoints:
+    def test_all_zero_points(self):
+        pnts = np.array([[0, 0]])
+        wcs, naxis = imp_utils.create_wcs_from_points(pnts, 5)
+        np.testing.assert_array_equal(naxis, [1, 1])
+        np.testing.assert_array_equal(wcs.wcs.crpix, [1, 1])
+        np.testing.assert_array_equal(wcs.wcs.crval, [0, 0])
+
+    def test_all_within_one_pixel(self):
+        pnts = np.array([[-1, -1], [-1, 1], [1, -1], [1, 1]])
+        wcs, naxis = imp_utils.create_wcs_from_points(pnts, 5)
+        np.testing.assert_array_equal(naxis, [1, 1])
+        np.testing.assert_array_equal(wcs.wcs.crpix, [1, 1])
+        np.testing.assert_array_equal(wcs.wcs.crval, [0, 0])
+
+    def test_all_within_one_pixel_along_one_axis(self):
+        pnts = np.array([[-3, 0], [-2, 0], [0, 0], [2, 0]])
+        wcs, naxis = imp_utils.create_wcs_from_points(pnts, 1)
+        np.testing.assert_array_equal(naxis, [5, 1])
+        np.testing.assert_array_equal(wcs.wcs.crpix, [3, 1])
+        np.testing.assert_array_equal(wcs.wcs.crval, [-0.5, 0])