astropy: Incorrect WCS conversion with SIP distortions (error ~1 arcmin!)

I am trying to convert world to pixel coordinates using the following header:

https://gist.githubusercontent.com/astrofrog/57bcb965eab79f99d953d2c1b162af02/raw/29ecd814239f11237ebe3bb5668582977976e9d1/header.hdr

Here’s an example of converting (19500,6500) in pixel coordinates to (ra,dec):

In [1]: from astropy.io import fits

In [2]: from astropy.wcs import WCS

In [3]: header = fits.Header.fromtextfile('header.hdr')

In [4]: wcs = WCS(header)

In [5]: wcs.all_pix2world(19500, 6500, 1)
Out[5]: [array(215.17565987230577), array(52.917709465941066)]

I’m using origin=1 to be able to compare to ds9 results. If I now open ds9, and convert this same position, I get:

ra=215.0174 dec=52.948131

which is almost an arcminute.

If I then try and convert the world coordinates back to pixel coordinates, I get an error, even though the world coordinates are inside the image:

In [8]: wcs.all_world2pix(215.17565987230577, 52.917709465941066, 1)
---------------------------------------------------------------------------
NoConvergence                             Traceback (most recent call last)
<ipython-input-8-e392834de60b> in <module>()
----> 1 wcs.all_world2pix(215.17565987230577, 52.917709465941066, 1)

/Users/tom/miniconda3/envs/dev35/lib/python3.5/site-packages/astropy/wcs/wcs.py in all_world2pix(self, *args, **kwargs)
   1806                 adaptive=adaptive, detect_divergence=detect_div,
   1807                 quiet=quiet),
-> 1808             'input', *args, **kwargs
   1809         )
   1810 

/Users/tom/miniconda3/envs/dev35/lib/python3.5/site-packages/astropy/wcs/wcs.py in _array_converter(self, func, sky, *args, **kwargs)
   1242                     "a 1-D array for each axis, followed by an origin.")
   1243 
-> 1244             return _return_list_of_arrays(axes, origin)
   1245 
   1246         raise TypeError(

/Users/tom/miniconda3/envs/dev35/lib/python3.5/site-packages/astropy/wcs/wcs.py in _return_list_of_arrays(axes, origin)
   1198             if ra_dec_order and sky == 'input':
   1199                 xy = self._denormalize_sky(xy)
-> 1200             output = func(xy, origin)
   1201             if ra_dec_order and sky == 'output':
   1202                 output = self._normalize_sky(output)

/Users/tom/miniconda3/envs/dev35/lib/python3.5/site-packages/astropy/wcs/wcs.py in <lambda>(*args, **kwargs)
   1805                 *args, tolerance=tolerance, maxiter=maxiter,
   1806                 adaptive=adaptive, detect_divergence=detect_div,
-> 1807                 quiet=quiet),
   1808             'input', *args, **kwargs
   1809         )

/Users/tom/miniconda3/envs/dev35/lib/python3.5/site-packages/astropy/wcs/wcs.py in _all_world2pix(self, world, origin, tolerance, maxiter, adaptive, detect_divergence, quiet)
   1786                     .format(k), best_solution=pix,
   1787                     accuracy=np.abs(dpix), niter=k,
-> 1788                     slow_conv=ind, divergent=inddiv)
   1789 
   1790         return pix

NoConvergence: 'WCS.all_world2pix' failed to converge to the requested accuracy.
After 2 iterations, the solution is diverging at least for one input point.

@mcara @nden - do you have any idea what could be going on here?