Weird behavior for trivial table header

[teutoburg]

Considering a minimal header from a table with one point source at (0, 0) and arbitrary CDELT:

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

then

impu.calc_footprint(hdr)

results in

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

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

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

the resulting footprint is

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

which seems even worse somehow...

To give an example where it actually works:

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

produces

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

which to me seems correct.

[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.