[AstroPy] Running wcs_pix2world for all image pixels

Fri Oct 18 10:01:52 EDT 2013

Tom, by "if it is set to zero, the pixel coordinates have to be zero 
based" do you include in that the reference pixel coordinate in the FITS 
header? If yes, then there's not really a choice here to say it's 
zero-based, as the FITS standard is 1-based?

Cheers
Maik

On 18/10/13 15:01, Thomas Robitaille wrote:
> In the Astropy routine, you need to take care to set the origin
> correctly (the last argument in the call to all_pix2world) - if it is
> set to zero, the pixel coordinates have to be zero based (the bottom
> left corner of the FITS image is then (0,0)). Tools such as ds9 assume
> the bottom left corner is (1,1). Since you mention a 1 pixel offset, I
> thought it could be related to this. Does AST always assume the pixel
> values are 1-based?
>
> Cheers,
> Tom
>
>
> On 18 October 2013 14:36, David Berry <d.berry at jach.hawaii.edu> wrote:
>>
>>
>> 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.]
>>
>>
>>
>> _______________________________________________
>> AstroPy mailing list
>> AstroPy at scipy.org
>> http://mail.scipy.org/mailman/listinfo/astropy
>>
> _______________________________________________
> AstroPy mailing list
> AstroPy at scipy.org
> http://mail.scipy.org/mailman/listinfo/astropy