can not parse the fits header with more than 2 axis

[henrysting]

It can not work with the fits file which has more than 2 axis, like in radio data.
For example, a header as below :

SIMPLE  =                    T /Standard FITS                                   
BITPIX  =                  -32 /Floating point (32 bit)                         
NAXIS   =                    4                                                  
NAXIS1  =                 1600                                                  
NAXIS2  =                 1600                                                  
NAXIS3  =                    1                                                  
NAXIS4  =                    1                                                  
EXTEND  =                    T                          
EQUINOX =   2.000000000000E+03                                                  
RADESYS = 'FK5     '                                                            
LONPOLE =   1.800000000000E+02                                                  
LATPOLE =  -3.529166667528E+00                                               
CTYPE1  = 'RA---SIN'                                                            
CRVAL1  =   4.201666666401E+01                                                  
CDELT1  =  -8.333333333333E-05                                                  
CRPIX1  =   8.010000000000E+02                                                  
CUNIT1  = 'deg     '                                                            
CTYPE2  = 'DEC--SIN'                                                            
CRVAL2  =  -3.529166667528E+00                                                  
CDELT2  =   8.333333333333E-05                                                  
CRPIX2  =   8.010000000000E+02                                                  
CUNIT2  = 'deg     '                                                            
CTYPE3  = 'STOKES  '                                                            
CRVAL3  =   1.000000000000E+00                                                  
CDELT3  =   1.000000000000E+00                                                  
CRPIX3  =   1.000000000000E+00                                                  
CUNIT3  = '        '                                                            
CTYPE4  = 'FREQ    '                                                            
CRVAL4  =   1.519499768816E+09                                                  
CDELT4  =   1.024001037953E+09                                                  
CRPIX4  =   1.000000000000E+00                                                  
CUNIT4  = 'Hz      '                                                            
PV2_1   =   0.000000000000E+00                                                  
PV2_2   =   0.000000000000E+00                                                  

It will cause the error:

  File "build/bdist.linux-x86_64/egg/pyregion/__init__.py", line 57, in as_imagecoord
  File "build/bdist.linux-x86_64/egg/pyregion/ds9_region_parser.py", line 184, in sky_to_image
  File "build/bdist.linux-x86_64/egg/pyregion/wcs_helper.py", line 231, in _get_radesys
ValueError: too many values to unpack

[mhardcastle]

This is true, but in a sense ds9's fault: as ds9 does not put Stokes or frequency values into its regions even when applied to a cube, pyregion is necessarily going to have trouble applying them to images.

In my ds9 radio flux measurement plugin http://github.com/mhardcastle/radioflux I go to a lot of trouble to 'flatten' radio maps down to the 2D form expected to allow pyregion to apply regions to them. I'm not sure this is logic that wants to be in pyregion itself -- e.g. what should you do if you actually have a cube, as opposed to the (1,1,y,x) shape of many radio maps?

[bsipocz]

@henrysting @mhardcastle - There are plans to deprecate pyregions when the replacement package, regions, is mature enough. You may want to take a look at it, and raise feature requests there rather than here.

Regarding the original issue, maybe the spectral-cube package is more relevant place for it.

[leejjoon]

I have been knowing this issue, but was not quite sure what the right way to fix is.
As a workaournd, I recommend you to read the header as a wcs object and its 'sub' method to retrieve space part of the axes. Maybe we can safely do this within the pyregion,

wcs = pywcs.WCS(f[0].header)
wcs_im = wcs.sub([1,2]) # index starts at 1
# wcs_im = wcs.sub([pywcs.WCSSUB_LONGITUDE, pywcs.WCSSUB_LATITUDE])
r2 = r.as_imagecoord(wcs_im)

[keflavich]

wcs.sub([wcs.WCSSUB_CELESTIAL]) is probably the 'best' way to do this

EDIT: the shortcut for that is just wcs.celestial

[cdeil]

Does anyone here have time to make a pull request adding this example to the docs?

[mhardcastle]

I will try to test this and document it. But it might be even better if this were the default behaviour in pyregion, since ds9 regions are always spatial, and so always need to be applied to the celestial subset of wcs?

[mhardcastle]

Documented, pull request raised.

[cdeil]

@mhardcastle documented this in #81 .

But it might be even better if this were the default behaviour in pyregion, since ds9 regions are always spatial, and so always need to be applied to the celestial subset of wcs?

I don't know. To me both options seem reasonable:

1. letting the user call wcs.sub([wcs.WCSSUB_CELESTIAL]) before passing it to pyregion
2. calling wcs.sub([wcs.WCSSUB_CELESTIAL]) on input in pyregion methods to automatically do this for users as a convenience

Another option would be to check the wcs on input and give a better errror message.

@leejjoon @keflavich @astrofrog - Thoughts?

[leejjoon]

My inclination is option 2, but only when fits header is given as an argument. I would suggest that no automatic call is made If wcs object is given. Just in case of any odd cases.