Fix #490: make rescale_imagehdu more robust against dimension mismatches

[oczoske]

A failure in image_plane_utils.rescale_imagehdu was noted when running scopesim with a cube source (e.g. in notebook LSS_AGN-01_preparation.ipynb in irdb/METIS).
The main change here is that the length of the zoom factor is taken from the data dimensions, not from the WCS. The fix has been tested with the three LSS notebooks for METIS and covered by an additional unit test.

Warnings about a mismatch between data and WCS dimensions occur often, which is likely more annoying than harmful; the origin of that lies elsewhere, in the creation of the ImagePlaneHDU.

[codecov]

Codecov Report

All modified and coverable lines are covered by tests ✅

Project coverage is 77.18%. Comparing base (f1a140f) to head (99f5ac5).
Report is 5 commits behind head on main.

Additional details and impacted files

@@            Coverage Diff             @@
##             main     #503      +/-   ##
==========================================
- Coverage   77.19%   77.18%   -0.02%     
==========================================
  Files          66       66              
  Lines        8209     8205       -4     
==========================================
- Hits         6337     6333       -4     
  Misses       1872     1872              

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[hugobuddel]

I had some trouble understanding where this equation came from. Maybe we can add the derivation as a comment to the code?

Suggested change

	ww.wcs.crpix[:2] = (zoom[:2] + 1) / 2 + (ww.wcs.crpix[:2] - 1) * zoom[:2]
	# In the old coordinate system the following holds (assuming linearity):
	# VAL = CRVAL + (PIX - CRPIX) * CDELT
	# Denoting the zoomed coordinate system with primes we get:
	# VAL' = CRVAL' + (PIX' - CRPIX') * CDELT'
	# Where by definition CDELT' = CDELT / ZOOM, and CRVAL' = CRVAL.
	#
	# There is always a fixed point in such a transformation.
	# All values refer to the center of the pixels, so it is the
	# edge of the first pixel whose VAL is conserved.
	# That is, for PIX = 1/2, PIX' = 1/2, and VAL' = VAL.
	#
	# Filling the above two equations with the values for the
	# edge of the first pixel yields:
	# VAL = CRVAL + (1/2 - CRPIX) * CDELT
	# VAL = CRVAL + (1/2 - CRPIX') * CDELT / ZOOM
	# Equating these equations allows us to solve for CRPIX':
	# (1/2 - CRPIX) * CDELT = (1/2 - CRPIX') * CDELT / ZOOM
	# (1/2 - CRPIX) * ZOOM = 1/2 - CRPIX'
	# CRPIX' = 1/2 - (1/2 - CRPIX) * ZOOM
	# CRPIX' = 1/2 + (CRPIX - 1/2) * ZOOM
	# CRPIX' = 1/2 + CRPIX * ZOOM - ZOOM / 2
	# Which is (effectively) the same as the equation below.
	ww.wcs.crpix[:2] = (zoom[:2] + 1) / 2 + (ww.wcs.crpix[:2] - 1) * zoom[:2]

[oczoske]

I've added the derivation in a somewhat shorter form and a more rational form of the equation (your line 554).

[hugobuddel]

Thanks! I found my derivation a bit too verbose too.

I was wondering whether there is a way to just let astropy.wcs do this for us. But that isn't really relevant now we already have this nicely working.

[teutoburg]

To avoid all those [:2], we could use wcs.celestial, but I think that only works if CTYPE is something like RA---TAN/DEC--TAN, so would not work for the mm-based detector WCS. So at this point, I can't come up something better 👍

I'm in favor of including the explanatory comment suggested by @hugobuddel .

[oczoske]

LINEAR is not wcs.celestial, so unfortunately that doesn't work. I'll merge the current version.

[teutoburg]

LINEAR is not wcs.celestial, so unfortunately that doesn't work.

Hmm, yes that's what I thought.