Bug in coordinate transform with WCSLib 5

[astrofrog]

The change to WCSLib 5 has led to a bug in the following use case:

import numpy as np
from astropy.wcs import WCS
from astropy.utils.data import get_pkg_data_filename

wcs = WCS(get_pkg_data_filename('galactic_center/gc_msx_e.fits'))

world = np.array(
    [[-0.58995335, -0.5       ],
     [ 0.00664326, -0.5       ],
     [-0.58995335, -0.25      ],
     [ 0.00664326, -0.25      ],
     [-0.58995335,  0.        ],
     [ 0.00664326,  0.        ],
     [-0.58995335,  0.25      ],
     [ 0.00664326,  0.25      ],
     [-0.58995335,  0.5       ],
     [ 0.00664326,  0.5       ]]
)

pixel = wcs.wcs_world2pix(world, 1)

print(pixel)

Before 7bac470

[[  1.64400000e+02  -1.51498185e-01]
 [  7.49105110e+01  -1.51498185e-01]
 [  1.64400000e+02   3.73485009e+01]
 [  7.49105110e+01   3.73485009e+01]
 [  1.64400000e+02   7.48485000e+01]
 [  7.49105110e+01   7.48485000e+01]
 [  1.64400000e+02   1.12348499e+02]
 [  7.49105110e+01   1.12348499e+02]
 [  1.64400000e+02   1.49848498e+02]
 [  7.49105110e+01   1.49848498e+02]]

After 7bac470

[[  1.64400000e+02  -1.51498185e-01]
 [  1.50906998e+02  -1.51498185e-01]
 [  7.49105110e+01   3.73485009e+01]
 [  1.50906998e+02   3.73485009e+01]
 [  1.64400000e+02   7.48485000e+01]
 [  1.13406999e+02   7.48485000e+01]
 [  7.49105110e+01   1.12348499e+02]
 [  1.13406999e+02   1.12348499e+02]
 [  1.64400000e+02   1.49848498e+02]
 [  7.59070000e+01   1.49848498e+02]]

This was discovered using the WCSAxes test suite: astrofrog/wcsaxes#168

cc @mdboom

[mdboom]

Wow. That's sure something. Looking into it now.

[cdeil]

The change to WCSLib 5 also broke a test in Gammapy:
https://travis-ci.org/gammapy/gammapy/jobs/68707705#L1115

Here's the print output corresponding to this test case that shows that the issue is in the WCS transformations, it doesn't round-trip:

import numpy as np
from astropy.wcs import WCS
from gammapy.image import make_header

header = make_header(nxpix=2, nypix=1, binsz=10, xref=0, yref=0, proj='CAR')
# GLON pixel edges: (+10, 0, -10)
# GLAT pixel edges: (-5, +5)

EPS = 0.1
data = [
        ( 5,  5,  1),        # in image[0, 0]
        ( 0,  0 + EPS,  2),  # in image[1, 0]
        ( 5, -5 + EPS,  3),  # in image[0, 0]
        ( 5,  5 + EPS, 99),  # outside image
        (10 + EPS, 0, 99),   # outside image
        ]
lon, lat, weights = np.array(data).T

wcs = WCS(header)
origin = 0
x, y = wcs.wcs_world2pix(lon, lat, origin)
lon2, lat2 = wcs.wcs_pix2world(x, y, origin)
print(wcs)
print()
print('weights : ', weights)
print('lon     : ', lon)
print('lon2    : ', lon2)
print('lat     : ', lat)
print('lat2    : ', lat2)
print('x       : ', x)
print('y       : ', y)

Output:

WCS Keywords

Number of WCS axes: 2
CTYPE : 'GLON-CAR'  'GLAT-CAR'  
CRVAL : 0.0  0.0  
CRPIX : 1.5  1.0  
PC1_1 PC1_2  : 1.0  0.0  
PC2_1 PC2_2  : 0.0  1.0  
CDELT : -10.0  10.0  
NAXIS    : 2 1

weights :  [  1.   2.   3.  99.  99.]
lon     :  [  5.    0.    5.    5.   10.1]
lon2    :  [ 5.   5.   0.   0.1  5. ]
lat     :  [ 5.   0.1 -4.9  5.1  0. ]
lat2    :  [ 5.   0.1 -4.9  5.1  0. ]
x       :  [ 0.    0.    0.5   0.49  0.  ]
y       :  [ 0.5   0.01 -0.49  0.51  0.  ]

I didn't check if #3893 fixes this case also.

[astrofrog]

@cdeil - if you don't want to try out #3893, you can also just check whether converting the coordinates one by one works (the bug here is related to iteration over coordinate arrays)

[cdeil]

#3893 did fix the issue in Gammapy also. Thanks, @mdboom !

[mdboom]

Fixed by #3893.