CARTAvis / carta

To CARTA users, this repo holds the CARTA release packages. Please use this repo to log bugs and feature requests. These will be triaged by the development team and prioritised as necessary in the development cycles.
19 stars 0 forks source link

Can't open a particular (otherwise valid) fits file #154

Open kmhess opened 2 years ago

kmhess commented 2 years ago

Describe the bug In going through a selection of cubes I have, I found one that CARTA can't read, and I don't know why. The tip off is that upon trying to load the data, CARTA says "No image hdus found. Select another file from the folder" and then the Load button is grayed out. On the other hand, I can read the file into astropy just fine and plot it with the correct WCS, etc:

In [1]: from astropy.io import fits
In [2]: a=fits.open('hcg44_3axes.fits')
In [3]: a.info()
Filename: hcg44_3axes.fits
No.    Name      Ver    Type      Cards   Dimensions   Format
  0  PRIMARY       1 PrimaryHDU     330   (700, 700, 158)   float32   

The file also opens fine in ds9. I'm happy to upload the file somewhere for someone to look at. It's 295 MB.

To Reproduce Attempt to open the file.

Expected behavior Display header and load file.

Screenshots or videos If applicable, add screenshots/videos to help explain your problem.

Platform info (please complete the following information):

Additional context I found this while seeing if the following issue was still a problem in v3-beta: https://github.com/CARTAvis/carta-frontend/issues/1877

kswang1029 commented 2 years ago

Could you attach the header if you load with ds9? Does it have multiple hdus?

kmhess commented 2 years ago

It has a single hdu. Here is the beginning of the header (minus all the rest of the the HISTORY):

SIMPLE  =                    T / SIMPLE FITS FORMAT                             
BITPIX  =                  -32 / NUMBER OF BITS PER PIXEL                       
NAXIS   =                    3 / NUMBER OF AXES                                 
NAXIS1  =                  700 / LENGTH OF AXIS                                 
NAXIS2  =                  700 / LENGTH OF AXIS                                 
NAXIS3  =                  158 / LENGTH OF AXIS                                 
BLOCKED =                    T / TAPE MAY BE BLOCKED                            
CDELT1  =  -2.222222306910E-03 / PRIMARY PIXEL SEPARATION                       
CRPIX1  =   3.510000000000E+02 / PRIMARY REFERENCE PIXEL                        
CRVAL1  =   1.547258373260E+02 / PRIMARY REFERENCE VALUE                        
CTYPE1  = 'RA---NCP          ' / PRIMARY AXIS NAME                              
CUNIT1  = 'DEGREE            ' / PRIMARY AXIS UNITS                             
CDELT2  =   2.222222306910E-03 / PRIMARY PIXEL SEPARATION                       
CRPIX2  =   3.510000000000E+02 / PRIMARY REFERENCE PIXEL                        
CRVAL2  =   2.199180661480E+01 / PRIMARY REFERENCE VALUE                        
CTYPE2  = 'DEC--NCP          ' / PRIMARY AXIS NAME                              
CUNIT2  = 'DEGREE            ' / PRIMARY AXIS UNITS                             
CDELT3  =   8.243965186790E+03 / PRIMARY PIXEL SEPARATION                       
CRPIX3  =   0.000000000000E+00 / PRIMARY REFERENCE PIXEL                        
CRVAL3  =   8.803330280140E+05 / PRIMARY REFERENCE VALUE                        
CTYPE3  = 'VELO-HEL          ' / PRIMARY AXIS NAME                              
CUNIT3  = 'M/S               ' / PRIMARY AXIS UNITS                             
CDELT4  =   1.000000000000E+00 / PRIMARY PIXEL SEPARATION                       
CRVAL4  =   1.000000000000E+00 / PRIMARY REFERENCE VALUE                        
CTYPE4  = 'STOKES            ' / PRIMARY AXIS NAME                              
CUNIT4  = 'PAR               ' / PRIMARY AXIS UNITS                             
EPOCH   =   2.000000000000E+03 / EPOCH                                          
FREQ0   =   1.420405758580E+09 / REST FREQUENCY                                 
INSTRUME= 'WSRT              ' / INSTRUMENT                                     
DATAMAX =         9.255453E-03 / MAXIMUM DATA VALUE                             
DATAMIN =        -1.244319E-03 / MINIMUM DATA VALUE                             
ORIGIN  = 'WFITS VERSION 1.3' / VERSION OF THE GIPSY PROGRAM                    
OBJECT  = 'hcg44   '  /                                                         
OBSERVER= 'ProcessDZB v11.'  /                                                  
CELLSCAL= '1/F     '  /                                                         
BUNIT   = 'JY/BEAM '  /                                                         
BMAJ    =    1.47306183353E-02  /                                               
BMIN    =    9.07219853252E-03  /                                               
DATE-OBS= '2010-02-07T18:38:43.0'  /                                            
BTYPE   = 'intensity'  /                                                        
SPECSYS = 'BARYCENT'  /                                                         
BPA     =    6.48472976685E+00  /                                               
NITERS  =                19976  /                                               
LWIDTH  =    2.00000000000E+00  /                                               
LSTEP   =    2.00000000000E+00  /                                               
LSTART  =    4.00000000000E+02  /                                               
VOBS    =   -9.42896080017E+00  /                                               
LTYPE   = 'channel '  /                                                         
HISTORY SORTORD = 'TB      '                                                    
HISTORY AIPS WTSCAL = 1.0                                                       
HISTORY FITS: Miriad Fits: version 1.2 07-Jul-00                                
HISTORY FITS: Executed on: 10OCT04:15:22:50.0                                   
HISTORY FITS: Command line inputs follow:                                       
HISTORY FITS:   in=11000196_S0_T0.UVF                                           
HISTORY FITS:   op=uvin                                                         
HISTORY FITS:   out=11000196_S0_T0.uv                                           
HISTORY FITS:   velocity=optbary,1350.00,513                                    
HISTORY PUTHD: Miriad PutHd: Version 1.0 11-Jan-96                              
HISTORY PUTHD: Executed on: 10OCT04:15:23:29.0                                  
HISTORY PUTHD: Command line inputs follow:                                      
HISTORY PUTHD:   in=11000196_S0_T0.uv/restfreq                                  
HISTORY PUTHD:   value=1.420405752                                              
HISTORY PUTHD:   type=double                                                    
HISTORY UVFLAG: Executed on: 10OCT04:15:23:40.0                                 
HISTORY UVFLAG: Command line inputs follow:                                     
HISTORY UVFLAG:   vis=11000196_S0_T0.uv                                         
HISTORY UVFLAG:   select="shadow(27)"                                           
HISTORY UVFLAG:   flagval=flag                                                  
HISTORY UVFLAG: Overview of flagging on visibility file 11000196_S0_T0.uv       
kswang1029 commented 2 years ago

It has a single hdu. Here is the beginning of the header (minus all the rest of the the HISTORY):

SIMPLE  =                    T / SIMPLE FITS FORMAT                             
BITPIX  =                  -32 / NUMBER OF BITS PER PIXEL                       
NAXIS   =                    3 / NUMBER OF AXES                                 
NAXIS1  =                  700 / LENGTH OF AXIS                                 
NAXIS2  =                  700 / LENGTH OF AXIS                                 
NAXIS3  =                  158 / LENGTH OF AXIS                                 
BLOCKED =                    T / TAPE MAY BE BLOCKED                            
CDELT1  =  -2.222222306910E-03 / PRIMARY PIXEL SEPARATION                       
CRPIX1  =   3.510000000000E+02 / PRIMARY REFERENCE PIXEL                        
CRVAL1  =   1.547258373260E+02 / PRIMARY REFERENCE VALUE                        
CTYPE1  = 'RA---NCP          ' / PRIMARY AXIS NAME                              
CUNIT1  = 'DEGREE            ' / PRIMARY AXIS UNITS                             
CDELT2  =   2.222222306910E-03 / PRIMARY PIXEL SEPARATION                       
CRPIX2  =   3.510000000000E+02 / PRIMARY REFERENCE PIXEL                        
CRVAL2  =   2.199180661480E+01 / PRIMARY REFERENCE VALUE                        
CTYPE2  = 'DEC--NCP          ' / PRIMARY AXIS NAME                              
CUNIT2  = 'DEGREE            ' / PRIMARY AXIS UNITS                             
CDELT3  =   8.243965186790E+03 / PRIMARY PIXEL SEPARATION                       
CRPIX3  =   0.000000000000E+00 / PRIMARY REFERENCE PIXEL                        
CRVAL3  =   8.803330280140E+05 / PRIMARY REFERENCE VALUE                        
CTYPE3  = 'VELO-HEL          ' / PRIMARY AXIS NAME                              
CUNIT3  = 'M/S               ' / PRIMARY AXIS UNITS                             
CDELT4  =   1.000000000000E+00 / PRIMARY PIXEL SEPARATION                       
CRVAL4  =   1.000000000000E+00 / PRIMARY REFERENCE VALUE                        
CTYPE4  = 'STOKES            ' / PRIMARY AXIS NAME                              
CUNIT4  = 'PAR               ' / PRIMARY AXIS UNITS                             
EPOCH   =   2.000000000000E+03 / EPOCH                                          
FREQ0   =   1.420405758580E+09 / REST FREQUENCY                                 
INSTRUME= 'WSRT              ' / INSTRUMENT                                     
DATAMAX =         9.255453E-03 / MAXIMUM DATA VALUE                             
DATAMIN =        -1.244319E-03 / MINIMUM DATA VALUE                             
ORIGIN  = 'WFITS VERSION 1.3' / VERSION OF THE GIPSY PROGRAM                    
OBJECT  = 'hcg44   '  /                                                         
OBSERVER= 'ProcessDZB v11.'  /                                                  
CELLSCAL= '1/F     '  /                                                         
BUNIT   = 'JY/BEAM '  /                                                         
BMAJ    =    1.47306183353E-02  /                                               
BMIN    =    9.07219853252E-03  /                                               
DATE-OBS= '2010-02-07T18:38:43.0'  /                                            
BTYPE   = 'intensity'  /                                                        
SPECSYS = 'BARYCENT'  /                                                         
BPA     =    6.48472976685E+00  /                                               
NITERS  =                19976  /                                               
LWIDTH  =    2.00000000000E+00  /                                               
LSTEP   =    2.00000000000E+00  /                                               
LSTART  =    4.00000000000E+02  /                                               
VOBS    =   -9.42896080017E+00  /                                               
LTYPE   = 'channel '  /                                                         
HISTORY SORTORD = 'TB      '                                                    
HISTORY AIPS WTSCAL = 1.0                                                       
HISTORY FITS: Miriad Fits: version 1.2 07-Jul-00                                
HISTORY FITS: Executed on: 10OCT04:15:22:50.0                                   
HISTORY FITS: Command line inputs follow:                                       
HISTORY FITS:   in=11000196_S0_T0.UVF                                           
HISTORY FITS:   op=uvin                                                         
HISTORY FITS:   out=11000196_S0_T0.uv                                           
HISTORY FITS:   velocity=optbary,1350.00,513                                    
HISTORY PUTHD: Miriad PutHd: Version 1.0 11-Jan-96                              
HISTORY PUTHD: Executed on: 10OCT04:15:23:29.0                                  
HISTORY PUTHD: Command line inputs follow:                                      
HISTORY PUTHD:   in=11000196_S0_T0.uv/restfreq                                  
HISTORY PUTHD:   value=1.420405752                                              
HISTORY PUTHD:   type=double                                                    
HISTORY UVFLAG: Executed on: 10OCT04:15:23:40.0                                 
HISTORY UVFLAG: Command line inputs follow:                                     
HISTORY UVFLAG:   vis=11000196_S0_T0.uv                                         
HISTORY UVFLAG:   select="shadow(27)"                                           
HISTORY UVFLAG:   flagval=flag                                                  
HISTORY UVFLAG: Overview of flagging on visibility file 11000196_S0_T0.uv       

Thanks. I will see if I can reproduce the issue with a fake image with your example header. If not I will need a copy of your test image if possible. Will update you again if I found something.

veggiesaurus commented 2 years ago

BLOCKED = T / TAPE MAY BE BLOCKED

This is a weird header...

kswang1029 commented 2 years ago

if you remove

CDELT4 = 1.000000000000E+00 / PRIMARY PIXEL SEPARATION
CRVAL4 = 1.000000000000E+00 / PRIMARY REFERENCE VALUE
CTYPE4 = 'STOKES ' / PRIMARY AXIS NAME
CUNIT4 = 'PAR ' / PRIMARY AXIS UNITS

the image would load. Note that in this header, only three axes are claimed.

NAXIS = 3 / NUMBER OF AXES
NAXIS1 = 700 / LENGTH OF AXIS
NAXIS2 = 700 / LENGTH OF AXIS
NAXIS3 = 158 / LENGTH OF AXIS

Screen Shot 2022-06-13 at 19 23 59

kmhess commented 2 years ago

Hi @kswang1029,

Thanks for finding that. I don't think that should prevent CARTA from working though. Through conversation with the astropy folks, they suggest that the WCS in the header and the dimensions of the data don't need to match: https://github.com/astropy/astropy/issues/13307#issuecomment-1147500220

Degenerate stokes axes are common, and in fact I have an example of a FITS file that has almost the exact same features which CARTA recognizes just fine. In this case the only difference I see is the CUNIT4 = 'par'. Could this be what is causing CARTA to trip up?

Screen Shot 2022-06-13 at 14 12 12
kswang1029 commented 2 years ago

Hi @kswang1029,

Thanks for finding that. I don't think that should prevent CARTA from working though. Through conversation with the astropy folks, they suggest that the WCS in the header and the dimensions of the data don't need to match: astropy/astropy#13307 (comment)

Degenerate stokes axes are common, and in fact I have an example of a FITS file that has almost the exact same features which CARTA recognizes just fine. In this case the only difference I see is the CUNIT4 = 'par'. Could this be what is causing CARTA to trip up? Screen Shot 2022-06-13 at 14 12 12

The backend log shows the reason as

2022-06-13 12:31:59 SEVERE Cannot set coordinate system: coords.nPixelAxes() == 4, image.ndim() == 3

which I think it is due to mis-matching dimensions between the header and image array. However your new example is a counter-example. Could you also attach the header as seen via ds9 too? Thanks.

veggiesaurus commented 2 years ago

Hi @kswang1029,

Thanks for finding that. I don't think that should prevent CARTA from working though.

I agree, CARTA should be able to deal with this. Worst case, it should fail to parse the WCS, but it should still work, even in pixel coordinates.

kmhess commented 2 years ago

Here you go:

SIMPLE  =                    T / file does conform to FITS standard             
BITPIX  =                  -32 / number of bits per data pixel                  
NAXIS   =                    3 / number of data axes                            
NAXIS1  =                   84 / length of data axis 1                          
NAXIS2  =                   61 / length of data axis 2                          
NAXIS3  =                  100 / length of data axis 3                          
EXTEND  =                    T / FITS dataset may contain extensions            
COMMENT   FITS (Flexible Image Transport System) format is defined in 'Astronomy
COMMENT   and Astrophysics', volume 376, page 359; bibcode: 2001A&A...376..359H 
BLANK   =                   -1  /                                               
BUNIT   = 'JY/BEAM '  /                                                         
DATE-OBS= '2013-12-24T08:38:01.2'  /                                            
OBSRA   =    5.49280416806E+01  /                                               
OBSDEC  =   -3.51602777699E+01  /                                               
TELESCOP= 'ATCA    '  /                                                         
CRPIX1  = -1.330000000000000E+02                                                
CDELT1  =   -5.55555597572E-03  /                                               
CRVAL1  =    5.46208333333E+01  /                                               
CTYPE1  = 'RA---NCP'  /                                                         
CRPIX2  = 1.580000000000000E+02                                                 
CDELT2  =    5.55555597572E-03  /                                               
CRVAL2  =   -3.54502777778E+01  /                                               
CTYPE2  = 'DEC--NCP'  /                                                         
CRPIX3  = -1.750000000000000E+02                                                
CDELT3  =    6.59000000505E+03  /                                               
CRVAL3  =    1.66450003394E+05  /                                               
CTYPE3  = 'VELO-HEL'  /                                                         
CRPIX4  =    1.00000000000E+00  /                                               
CDELT4  =    1.00000000000E+00  /                                               
CRVAL4  =    1.00000000000E+00  /                                               
CTYPE4  = 'STOKES  '  /                                                         
CELLSCAL= '1/F     '  /                                                         
RESTFREQ=    1.42040571183E+09  /                                               
SPECSYS = 'BARYCENT'  /                                                         
VOBS    =    1.45609636307E+01  /                                               
BMAJ    =    2.63608861715E-02  /                                               
BMIN    =    1.84838455170E-02  /                                               
BPA     =    4.37762737274E-01  /                                               
BTYPE   = 'intensity'  /                                                        
EPOCH   =    2.00000000000E+03  /                                               
LTYPE   = 'velocity'  /                                                         
LSTART  =    1.66449996948E+02  /                                               
LSTEP   =    6.59000015259E+00  /                                               
LWIDTH  =    6.59000015259E+00  /                                               
OBJECT  = 'afs_014 '  /                                                         
OBSERVER= 'Serra   '  /                                                         
HISTORY ATLOD: Miriad atlod: Revision 1.53, 2016/08/25 23:02:21 UTC             
HISTORY ATLOD: Executed on: 2019-01-09T23:28:25.0                               
HISTORY ATLOD: Command line inputs follow:                                      
HISTORY ATLOD:   in=../2013-12-11*                                              
HISTORY ATLOD:   out=all.uv                                                     
HISTORY ATLOD:   ifsel=3                                                        
HISTORY ATLOD:   restfreq=1.420405752                                        
kswang1029 commented 2 years ago

Hi @kswang1029, Thanks for finding that. I don't think that should prevent CARTA from working though.

I agree, CARTA should be able to deal with this. Worst case, it should fail to parse the WCS, but it should still work, even in pixel coordinates.

We need to be very careful when making a loosing FITS header validation routine in CARTA. FITS headers can be very "CREATIVE" in general...

kswang1029 commented 2 years ago

Here you go:

SIMPLE  =                    T / file does conform to FITS standard             
BITPIX  =                  -32 / number of bits per data pixel                  
NAXIS   =                    3 / number of data axes                            
NAXIS1  =                   84 / length of data axis 1                          
NAXIS2  =                   61 / length of data axis 2                          
NAXIS3  =                  100 / length of data axis 3                          
EXTEND  =                    T / FITS dataset may contain extensions            
COMMENT   FITS (Flexible Image Transport System) format is defined in 'Astronomy
COMMENT   and Astrophysics', volume 376, page 359; bibcode: 2001A&A...376..359H 
BLANK   =                   -1  /                                               
BUNIT   = 'JY/BEAM '  /                                                         
DATE-OBS= '2013-12-24T08:38:01.2'  /                                            
OBSRA   =    5.49280416806E+01  /                                               
OBSDEC  =   -3.51602777699E+01  /                                               
TELESCOP= 'ATCA    '  /                                                         
CRPIX1  = -1.330000000000000E+02                                                
CDELT1  =   -5.55555597572E-03  /                                               
CRVAL1  =    5.46208333333E+01  /                                               
CTYPE1  = 'RA---NCP'  /                                                         
CRPIX2  = 1.580000000000000E+02                                                 
CDELT2  =    5.55555597572E-03  /                                               
CRVAL2  =   -3.54502777778E+01  /                                               
CTYPE2  = 'DEC--NCP'  /                                                         
CRPIX3  = -1.750000000000000E+02                                                
CDELT3  =    6.59000000505E+03  /                                               
CRVAL3  =    1.66450003394E+05  /                                               
CTYPE3  = 'VELO-HEL'  /                                                         
CRPIX4  =    1.00000000000E+00  /                                               
CDELT4  =    1.00000000000E+00  /                                               
CRVAL4  =    1.00000000000E+00  /                                               
CTYPE4  = 'STOKES  '  /                                                         
CELLSCAL= '1/F     '  /                                                         
RESTFREQ=    1.42040571183E+09  /                                               
SPECSYS = 'BARYCENT'  /                                                         
VOBS    =    1.45609636307E+01  /                                               
BMAJ    =    2.63608861715E-02  /                                               
BMIN    =    1.84838455170E-02  /                                               
BPA     =    4.37762737274E-01  /                                               
BTYPE   = 'intensity'  /                                                        
EPOCH   =    2.00000000000E+03  /                                               
LTYPE   = 'velocity'  /                                                         
LSTART  =    1.66449996948E+02  /                                               
LSTEP   =    6.59000015259E+00  /                                               
LWIDTH  =    6.59000015259E+00  /                                               
OBJECT  = 'afs_014 '  /                                                         
OBSERVER= 'Serra   '  /                                                         
HISTORY ATLOD: Miriad atlod: Revision 1.53, 2016/08/25 23:02:21 UTC             
HISTORY ATLOD: Executed on: 2019-01-09T23:28:25.0                               
HISTORY ATLOD: Command line inputs follow:                                      
HISTORY ATLOD:   in=../2013-12-11*                                              
HISTORY ATLOD:   out=all.uv                                                     
HISTORY ATLOD:   ifsel=3                                                        
HISTORY ATLOD:   restfreq=1.420405752                                        

ok. I see why. In the first header, rest frequency is defined as FREQ0, so CARTA does not "upgrade" the image dimension from 3 to 4. In the second header, rest frequency is defined as RESTFREQ, so CARTA "upgrades" the image dimension from 3 to 4. The backend log says

2022-06-13 12:55:53 INFO FITSCoordinateUtil::fromFITSHeader The WCS for this image contains 1 degenerate axes. 2022-06-13 12:55:53 INFO FITSImage::getImageAttributes Image dimension appears to be less than number of pixel axes in CoordinateSystem 2022-06-13 12:55:53 INFO FITSImage::getImageAttributes + Adding 1 degenerate trailing axes

So if I manually change FREQ0 to RESTFREQ in the first header, the image loads.

And if I manually change RESTFREQ in the 2nd header to FREQ0, the image does not load.

2022-06-13 13:00:35 INFO FITSCoordinateUtil::fromFITSHeader The WCS for this image contains 1 degenerate axes. 2022-06-13 13:00:35 WARN FITSCoordinateUtil::fromFITSHeader Zero or no rest frequency provided for velocity axis. 2022-06-13 13:00:35 SEVERE Cannot set coordinate system: coords.nPixelAxes() == 4, image.ndim() == 3 [2022-06-13 21:00:35.183] [CARTA] [error] No image hdus found.

kmhess commented 2 years ago

Ah, weird. Good catch. I haven't noticed that before...I've messaged the person I got the troublesome cube from. Maybe he can weigh in here.