[AstroPy] Running wcs_pix2world for all image pixels
David Berry
d.berry at jach.hawaii.edu
Fri Oct 18 08:36:22 EDT 2013
On 18 October 2013 12:49, Michael Droettboom <mdroe at stsci.edu> wrote:
> On 10/18/2013 04:54 AM, David Berry wrote:
>
>
>
>
> On 17 October 2013 16:28, Maik Riechert <maik.riechert at arcor.de> wrote:
>
> I didn't use tranGrid here, just tran! So no tolerance involved.
>>
>
> Doh! Sorry - didn't spot that.
>
>
>>
>> NAXIS1 = 1468 / image
>> width
>> NAXIS2 = 916 / image
>> height
>> CTYPE1 = 'RA---TAN-SIP' / TAN (gnomic) projection + SIP distortions
>> CTYPE2 = 'DEC--TAN-SIP' / TAN (gnomic) projection + SIP distortions
>> EQUINOX = 2000.0 / Equatorial coordinates definition (yr)
>> LONPOLE = 180.0 / no
>> comment
>> LATPOLE = 0.0 / no
>> comment
>> CRVAL1 = 302.88771538 / RA of reference
>> point
>> CRVAL2 = 41.8941029213 / DEC of reference
>> point
>> CRPIX1 = 735 / X reference
>> pixel
>> CRPIX2 = 459 / Y reference
>> pixel
>> CUNIT1 = 'deg ' / X pixel scale
>> units
>> CUNIT2 = 'deg ' / Y pixel scale
>> units
>> CD1_1 = -0.0157184611387 / Transformation
>> matrix
>> CD1_2 = 0.00022032143241 / no
>> comment
>> CD2_1 = -0.000284856284919 / no
>> comment
>> CD2_2 = -0.0159943489114 / no
>> comment
>> IMAGEW = 1468 / Image width, in
>> pixels.
>> IMAGEH = 916 / Image height, in
>> pixels.
>> A_ORDER = 2 / Polynomial order, axis
>> 1
>> A_0_2 = -5.10167355829E-05 / no
>> comment
>> A_1_1 = -5.54267287337E-05 / no
>> comment
>> A_2_0 = 3.68095680847E-05 / no
>> comment
>> B_ORDER = 2 / Polynomial order, axis
>> 2
>> B_0_2 = 2.49853669312E-05 / no
>> comment
>> B_1_1 = 6.07756179822E-05 / no
>> comment
>> B_2_0 = -1.70350503136E-05 / no
>> comment
>> AP_ORDER= 2 / Inv polynomial order, axis
>> 1
>> AP_0_1 = -0.0014471978634 / no
>> comment
>> AP_0_2 = 5.0935282476E-05 / no
>> comment
>> AP_1_0 = 0.000641040653934 / no
>> comment
>> AP_1_1 = 5.54267224252E-05 / no
>> comment
>> AP_2_0 = -3.674588665E-05 / no
>> comment
>> BP_ORDER= 2 / Inv polynomial order, axis
>> 2
>> BP_0_1 = 0.0010860421788 / no
>> comment
>> BP_0_2 = -2.49999638422E-05 / no
>> comment
>> BP_1_0 = -0.000569392740582 / no
>> comment
>> BP_1_1 = -6.07473585161E-05 / no comment
>> BP_2_0 = 1.7000827442E-05 / no comment
>>
>
> OK - so I suspect the problems may actually be in astropy rather than
> AST. Maybe the maintainers of astropy.wcs could comment.
>
> The above header includes a polynomial distortion correction using the
> Spitzer SIP conventions. The "A_" and "B_" headers give the coefficients of
> the forward ( (x,y) to (ra,dec) ) polynomial, and the "AP_" and "BP_"
> headers give the coefficients of the inverse ( (ra,dec) to (x,y) )
> polynomial.
>
> Using the header as supplied to transform the following pixel positions:
>
> x = -196, 735, 4059
> y = -504, 459, 2328
>
> I get, using AST (converted to degrees):
>
> ra = 328.090902, 302.887715, 265.839543
> dec = 53.570723, 41.894103, 8.257963
>
> whereas astropy gives:
>
> ra = 326.992376, 302.866900, 263.545072
> dec = 54.749009, 41.877821, 10.589047
>
> For an independent check I also used SAO WCSTools which gave:
>
> ra = 328.09090, 302.88772, 265.83954
> dec = 53.57072, 41.89410, 8.25796
>
> So AST and WCSTools agree, but astropy is different. This is not
> conclusive in itself, but I then modified the header by removing all the
> SIP keywords and removing "-SIP" from the CTYPE1/2 headers. With these
> changes the header represents an undistorted TAN projection. I then
> transformed the same pixel positions again using the changed header:
>
> AST gave:
>
> ra = 327.022569, 302.887715, 263.548641
> dec = 54.758313, 41.894103, 10.599558
>
> astropy gave:
>
> ra = 326.992376, 302.866900, 263.545072
> dec = 54.749009, 41.877822, 10.589047
>
> and WCSTools gave:
>
> ra = 327.02257, 302.88772, 263.54864
> dec = 54.75831, 41.89410, 10.59956
>
> So again AST and WCSTools agree, with astropy as the odd one out. But
> more importantly, the astropy results are exactly the same as when using
> the header that includes the SIP distortion. This seems to indicate that in
> fact astropy is ignoring the SIP distortion. You can see that the astropy
> values are much closer to the AST/WCSTools values when no distortion is
> included.
>
>
> You need to use wcs.all_pix2world to include the SIP distortion. SIP is
> not a part of the WCS standard, so we only perform it when explicitly asked
> for. See:
>
>
> https://astropy.readthedocs.org/en/v0.2.4/_generated/astropy.wcs.wcs.WCS.html#astropy.wcs.wcs.WCS.wcs_pix2world
>
>
Sorry - I didn't read the docs carefully enough. So does application code
need to decide for itself whether to use "wcs_*" or "all_*" functions,
based on its own interpretation of the header? Or is it safe just to use
"all_*" all the time, for both standard and non-standard headers? If that's
the case, then some confusion may be avoided by removing "wcs_*"
altogether.
Using all_pix2world, I now get:
ra = 328.059328, 302.866902, 265.838882
dec = 53.564017, 41.877821, 8.245622
which are much closer to the AST/WCSTools values:
ra = 328.090902, 302.887715, 265.839543
dec = 53.570723, 41.894103, 8.257963
but there is still a discrepancy of about a pixel. Also, it is still the
case that the (ra,dec) values given by astropy at the reference pixel
(735,459) are not equal to (crval1,crval2), which they should be since the
forward polynomial in the header contains no offset terms (A_0_0 or B_0_0).
Answering my own question from the start of this email, I tried
using all_world2pix on a standard header formed by stripping out all the
SIP stuff, and I got the appended error.
David
/usr/lib64/python2.7/site-packages/scipy/optimize/nonlin.py:942:
RuntimeWarning: invalid value encountered in divide
d = v / vdot(df, v)
Traceback (most recent call last):
File "./aa.py", line 80, in <module>
x, y = mywcs3.all_world2pix(ra, dec, 0)
File
"/home/dsb/.local/lib/python2.7/site-packages/astropy-0.3.dev5914-py2.7-linux-x86_64.egg/astropy/wcs/wcs.py",
line 1177, in all_world2pix
**kwargs)
File
"/home/dsb/.local/lib/python2.7/site-packages/astropy-0.3.dev5914-py2.7-linux-x86_64.egg/astropy/wcs/wcs.py",
line 990, in _array_converter
return _return_list_of_arrays(axes, origin)
File
"/home/dsb/.local/lib/python2.7/site-packages/astropy-0.3.dev5914-py2.7-linux-x86_64.egg/astropy/wcs/wcs.py",
line 946, in _return_list_of_arrays
output = func(xy, origin)
File
"/home/dsb/.local/lib/python2.7/site-packages/astropy-0.3.dev5914-py2.7-linux-x86_64.egg/astropy/wcs/wcs.py",
line 1175, in <lambda>
self._all_world2pix(*args, tolerance=tolerance, **kwargs),
File
"/home/dsb/.local/lib/python2.7/site-packages/astropy-0.3.dev5914-py2.7-linux-x86_64.egg/astropy/wcs/wcs.py",
line 1166, in _all_world2pix
soln = scipy.optimize.broyden1(func, x0, x_tol=tolerance)
File "<string>", line 8, in broyden1
File "/usr/lib64/python2.7/site-packages/scipy/optimize/nonlin.py", line
330, in nonlin_solve
raise NoConvergence(_array_like(x, x0))
scipy.optimize.nonlin.NoConvergence: [-196. -504.]
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/astropy/attachments/20131018/99ab238b/attachment.html>
More information about the AstroPy
mailing list