ceb8
commented
3 years ago This bug has been temporarily patched by simply not using the affected astropy function for the versions of astropy where it causes problems. This is not a long term solution.
The bug is a result of the changes to the astropy fit_wcs_from_points function PR https://github.com/astropy/astropy/pull/10865. It is not entirely clear how much of that PR needs to be walked back or otherwise updated, versus how much our code needs to be changes.
For sure the pixel shape calculation here https://github.com/astropy/astropy/pull/10865/files#diff-e6ff41e178a73eba0e29f8d018f2bc610617f18d9e653b75e2250017fcab34fcR1066 is incorrect. We cannot assume the minimum pixel value is zero, and it could be less than zero.
The trickier part is the +1 correction for FITS versus numpy arrays. This will require more exploration to determine if it is being treated correctly. It's possible it is being corrected at the wrong juncture, or double corrected, or our code could be wrong.
An example that causes the failure is as follows (I'm sure a more minimal example can be created, sorry about that):
from os import path
from astropy.coordinates import SkyCoord
from astropy.io import fits
import numpy as np
from itertools import product
from astrocut import CubeFactory, CutoutFactory
from astrocut.tests.utils_for_test import create_test_ffis
from astropy.wcs.utils import fit_wcs_from_points
cube_maker = CubeFactory()
img_sz = 10
num_im = 100
ffi_files = create_test_ffis(img_sz, num_im)
cube_file = cube_maker.make_cube(ffi_files, path.join("tmpdir", "test_cube.fits"), verbose=False)
# Setting up
myfactory = CutoutFactory()
hdu = fits.open(cube_file)
cube_table = hdu[2].data
myfactory._parse_table_info(cube_table)
cutout_wcs = myfactory.cube_wcs
cutout_shape = (1, 100)
y, x = cutout_shape
i = 1
while (x/i)*(y/i) > 100:
i += 1
xvals = list(reversed(range(x//2, -1, -i)))[:-1] + list(range(x//2, x, i))
if xvals[-1] != x-1:
xvals += [x-1]
if xvals[0] != 0:
xvals = [0] + xvals
yvals = list(reversed(range(y//2, -1, -i)))[:-1] + list(range(y//2, y, i))
if yvals[-1] != y-1:
yvals += [y-1]
if yvals[0] != 0:
yvals = [0] + yvals
pix_inds = np.array(list(product(xvals, yvals)))
world_pix = SkyCoord(cutout_wcs.all_pix2world(pix_inds, 0), unit='deg')
linear_wcs = fit_wcs_from_points([pix_inds[:, 0], pix_inds[:, 1]], world_pix, proj_point='center')
The Python 3.7 with dev astropy tests are failing, this needs to be fixed before the next astropy release.