WCS in Scopesim, off-by-one error, specifically `image_plane_utils.header_from_list_of_xy`

[oczoske]

The WCS treatment in header_from_list_of_xy highlights the misapplication of FITS conventions in Scopesim. This is just part of a larger problem that presumably leads to the long-standing "off-by-one" issue.
The FITS conventions counts pixels starting from 1, where (1, 1) is the centre of the bottom-left pixel in a 2D image. This convention is violated in lines 195/196. Since function pix2val also uses zero offset, this is not immediatly apparent, but it will become a problem when a header is interpreted with the astropy.wcs module (as any WCS of type other than LINEAR should be).
The centre of an image with (NAXIS1, NAXIS2) pixels is therefore at ((NAXIS1 + 1) / 2, (NAXIS2 + 1) / 2), this is violated in line 198.
Unfortunately, this cannot be fixed in isolation but should be fixed throughout Scopesim and be properly regression tested.

[hugobuddel]

Thanks @oczoske , you know this issue the best.

The lines in question:

ScopeSim/scopesim/optics/image_plane_utils.py

Lines 193 to 203 in add2094

	hdr["CRVAL1"+s] = crval1
	hdr["CRVAL2"+s] = crval2
	hdr["CRPIX1"+s] = 0.
	hdr["CRPIX2"+s] = 0.
	xpcen, ypcen = naxis1 // 2, naxis2 // 2
	xscen, yscen = pix2val(hdr, xpcen, ypcen, s)
	hdr["CRVAL1"+s] = float(xscen)
	hdr["CRVAL2"+s] = float(yscen)
	hdr["CRPIX1"+s] = xpcen
	hdr["CRPIX2"+s] = ypcen

I'm not familiar with this code, but it is weird: the lines 195/196 you mention seem to be pointless because a new value is set in 202/203.

Just for clarity, I'm sure you already know this: the CRPIX values can be anything; the CRPIX does not even have to be within the image (which might be a sensible choice for an instrument with multiple detectors). So the problem must be in combination with other parts of ScopeSim, e.g. with pix2val you mention.

How should we proceed with this?

Would it be a good approach to just always use astropy.wcs? Even for the linear transformations? Because we might have to drop the linearity assumption at some point anyway, because I think the distortions at the edges of the detector will be large enough to be noticeable. So just dropping all the homebrew implementations seems like a good start.

[oczoske]

You're right that the actual problem is in pix2val and val2pix, whereas in header_from_list_of_xy it is largely cosmetic. Since the intention here is clearly to refer to the centre of the image, that should be done correctly as well.
I think we should indeed use astropy.wcs throughout, even though it carries some overhead. Note that while LINEAR and WAVE are linear WCSs, RA---TAN/DEC--TAN is not, and that should be noticeable at the edges of the 12000x12000 field.

[hugobuddel]

Two other reasons for switching to WCS everywhere:

• Currently, ScopeSim only supports sources where PC1_1 is 1.0, not -1.0.
• Currently, ScopeSim does not support diagonals in the PC matrix.

[hugobuddel]

Summary of discussion with @teutoburg :

The cause for the off-by-one error in MICADO seems to be the FieldConstantPSF of the relay optics. The cause might be that the PSF in PSF_SCAO_ConstPSF_0_5off.fits has a CRPIX of (128, 128), but center of the PSF (/the peak off the flux values) is at (129, 129), in 1-based FITS coordinates.

We had some thoughts about tests; e.g. these three simulations should yield identical results:

1. A single star at the center of (e.g.) the test instrument basic_instrument
   • a single detector with an small odd number of detectors, say 11
   • with no other effects except and SeeingPSF
   • should get a single star in the center of the image
   • save the image
   • This ensures that we can properly simulate a source in the center w/o any FITS handling
2. The same as 1. above but
   • with the saved image as source
   • no psf effect
   • should give the same result as simulation 1
   • this ensures our FITS writing and reading is at least consistent
3. The same as 1 above but
   • with instead a FieldConstantPSF using the image from 1 as PSF
   • should also give identical results

[hugobuddel]

I created a notebook that simulates a single star in different ways that highlights the problem, see #273 .

Two pdfs of the notebook with against dev_master and against fh/off-by-one (471ad10):

• simplestars.dev_master.pdf
• simplestars.off-by-one.pdf

All figures should look more or less similar, like they do in the second attachment. So we seem to be close to fixing the problem.