search

radio-astro-tools / spectral-cube

Library for reading and analyzing astrophysical spectral data cubes
http://spectral-cube.rtfd.org
BSD 3-Clause "New" or "Revised" License
98 stars 65 forks source link

BUNIT value from FITS file header not passed to SpectralCube header #902

Open astro-pablo opened 9 months ago

astro-pablo commented 9 months ago

When reading a FITS file with the following commands:

# read the FITS file and produce Specral Cube
hdul         = fits.open(fits_dir+fits_file+'.fits') 
hdr          = hdul[0].header
hdr['CD3_3'] = hdr['CDELT3'] # to repair value misread from the FITS header. It seems it was looking for 'CD3_3', as it was already using 'CD1_1', and 'CD2_2' instead of CDELT1, CDELT2 which contain wrong values.
w            = WCS(hdr)
data       = hdul[0].data
cube     = SpectralCube(data=hdul[0].data,wcs=w)

the hdr[‘BUNIT’] = 'ergs/s/cm^2/A' string is not copied over to the Specral Cube header (cube.unit is empty). Here is the header of the original FITS file:

NAXIS3          =                12401                                                  
EXTEND        =                    T                                                  
XTENSION    = 'IMAGE   '           / Image extension                                
PCOUNT       =                    0 / number of parameters                           
GCOUNT      =                    1 / number of groups                               
BUNIT           = 'ergs/s/cm^2/A'                                                       
CRPIX1         =                 57.5                                                  
CRPIX2         =                 58.0                                                  
CDELT1         =                5E-11                                                  
CDELT2        =                  1.0                                                  
CUNIT1        = 'deg     '                                                            
CTYPE1        = 'RA---TAN'                                                            
CTYPE2        = 'DEC--TAN'                                                            
CRVAL1        =           134.861708                                                  
CRVAL2       =           -43.737239                                                  
LATPOLE     =                 90.0                                                  
WCSNAME  = 'DEFAULTS'                                                            
MJDREF      =                  0.0                                                  
EXTNAME   = 'FLUX    '           / extension name                                 
CD1_1          = -0.00513888888888888                                                  
CD1_2         =                 -0.0                                                  
CD2_1          =                 -0.0                                                  
CD2_2         = 0.005138888888888889                                                  
CUNIT2       = 'deg     '                                                            
CDELT3       =                  0.5                                                  
CRPIX3        =                  1.0                                                  
CRVAL3       =               3600.0                                                  
CUNIT3       = 'Angstrom'           / Units of coordinate increment and value        
CTYPE3      = 'WAVE    '           / Air wavelength (linear)                        
RADESYS   = 'FK5     '                                                            
OBJSYS     = 'ICRS    '                                                            
EQUINOX   =               2000.0                                                  
IFUCON     = '16209   '           / NFibers

Here is the output of the cube.wcs command, lacking the "unit" variable.

WCS Keywords

Number of WCS axes: 3
CTYPE : 'RA---TAN' 'DEC--TAN' 'WAVE' 
CRVAL : 134.861708 -43.737239 3.6e-07 
CRPIX : 57.5 58.0 1.0 
PC1_1 PC1_2 PC1_3  : 1.0 0.0 0.0 
PC2_1 PC2_2 PC2_3  : 0.0 1.0 0.0 
PC3_1 PC3_2 PC3_3  : 0.0 0.0 1.0 
CDELT : -0.00513888888888888 0.005138888888888889 5e-11 
NAXIS : 115  115  12401

On a side note, the CDELT1 and CDELT2 values in the original FITS header contains incorrect values, so the CD1_1 and CD2_2 values are used for the reconstruction of the spatial axis instead. In that sense, I had to create hdr['CD3_3'] = hdr['CDELT3'], as it seems Spectral Cube assumed that the variable CD3_3 would exist for the spectral reconstruction (?), while it was not in the original header.

keflavich commented 9 months ago

The problem here is that you're stripping the units:

hdul = fits.open(fits_dir+fits_file+'.fits')
hdr = hdul[0].header
hdr['CD3_3'] = hdr['CDELT3'] # to repair value misread from the FITS header. It seems it was looking for 'CD3_3', as it was already using 'CD1_1', and 'CD2_2' instead of CDELT1, CDELT2 which contain wrong values.
w = WCS(hdr)
data = hdul[0].data
cube = SpectralCube(data=hdul[0].data,wcs=w)

When you do w = WCS(hdr), that removes all non-WCS information from the header. What should work instead is [EDITED to fix]:

data = hdul[0].data * u.erg/u.s/u.cm**2/u.AA
cube = SpectralCube(data=data, header=hdr, wcs=w)

or even simpler:

cube = SpectralCube.read(hdul[0])

Note for @e-koch and myself: We don't demo this case on the reading docs; we should.

https://spectral-cube.readthedocs.io/en/latest/creating_reading.html

keflavich commented 9 months ago

Also, on the last part: cube.wcs will never include units, for the same reason: units are not part of the coordinate system.

And, yes, I think you did the right thing with CD3_3; mixing CDELT and the CD matrix is not allowed in the FITS standard

astro-pablo commented 9 months ago

Hi, thanks for the quick answer. I can not remove "wcs=w" from the SpectralCube(data=hdul[0].data, wcs=w) instruction, as I get the following error:

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Cell In[34], line 64
     62 w            = WCS(hdr)
     63 data         = hdul[0].data
---> 64 cube         = SpectralCube(data=hdul[0].data, header=hdr) # Initiate a SpectralCube. The data in the cube lives as an ndarray with shape (n_s, n_y, n_x)
     67 print()
     68 #print('Original Cube:')
     69 #print(repr(hdr))

TypeError: __init__() missing 1 required positional argument: 'wcs'

I also tried:

cube   = SpectralCube(data=hdul[0].data, header=hdr, wcs=w)

but it did not help. Units are still absent.

print(cube)
print('unit:'+str(cube.unit))
print(cube.wcs)

SpectralCube with shape=(12401, 115, 115):
 n_x:    115  type_x: RA---TAN  unit_x: deg    range:   134.450737 deg:  135.265531 deg
 n_y:    115  type_y: DEC--TAN  unit_y: deg    range:   -44.029442 deg:  -43.443603 deg
 n_s:  12401  type_s: WAVE      unit_s: Angstrom  range:     3600.000 Angstrom:    9800.000 Angstrom
unit:
WCS Keywords

Number of WCS axes: 3
CTYPE : 'RA---TAN' 'DEC--TAN' 'WAVE' 
CRVAL : 134.861708 -43.737239 3.6e-07 
CRPIX : 57.5 58.0 1.0 
PC1_1 PC1_2 PC1_3  : 1.0 0.0 0.0 
PC2_1 PC2_2 PC2_3  : 0.0 1.0 0.0 
PC3_1 PC3_2 PC3_3  : 0.0 0.0 1.0 
CDELT : -0.00513888888888888 0.005138888888888889 5e-11 
NAXIS : 115  115  12401
keflavich commented 9 months ago

You're right, you can't use SpectralCube(...) for this unless you initialize the data with a unit yourself.

Use SpectralCube.read instead. That's designed to handle FITS headers.

astro-pablo commented 9 months ago

Great!. Now I got a related warning:

WARNING: Could not parse unit ergs[/s/cm](http://localhost:8888/s/cm)^2[/A.](http://localhost:8888/A.) If you know the correct unit, try u.add_enabled_units(u.def_unit(['ergs[/s/cm](http://localhost:8888/s/cm)^2[/A](http://localhost:8888/A)'], represents=u.<correct_unit>)) [spectral_cube.cube_utils]

so this solved the problem:

u.add_enabled_units(u.def_unit(['ergs/s/cm^2/A'], represents=u.erg / u.s /u.cm**2 / u.Angstrom))

BUT, since I can not create hdr['CD3_3'] = hdr['CDELT3'] in this way, the spectral axis is messed up (it should go up to 9800 A and not 16000 A, this is because in the original FITS file: CDELT3 = 0.9999999999999999 instead of CDELT3 = 0.49999999999999994). How can I either include CD3_3, or modify CDELT3 in the SpectralCube header now?. I tried cube.header.set('CDELT3', 0.49999999999999994) but it did not modify the read-in CDELT3 value.

Here is an example of correct units, but wrong spectral axis:

cube         = SpectralCube.read(fits_dir+fits_file+'.fits')
cube.header.set('CDELT3', 0.49999999999999994)
print(cube)
print(repr(cube.header))

**SpectralCube with shape=(12401, 115, 115) and unit=ergs/s/cm^2/A:**
 n_x:    115  type_x: RA---TAN  unit_x: deg    range:   134.450737 deg:  135.265531 deg
 n_y:    115  type_y: DEC--TAN  unit_y: deg    range:   -44.029442 deg:  -43.443603 deg
 n_s:  12401  type_s: WAVE      unit_s: Angstrom  range:     3600.000 Angstrom:   16000.000 Angstrom

SIMPLE  =                    T / conforms to FITS standard                      
BITPIX  =                  -64 / array data type                                
NAXIS   =                    3                                                  
NAXIS1  =                  115                                                  
NAXIS2  =                  115                                                  
NAXIS3  =                12401                                                  
EXTEND  =                    T                                                  
XTENSION= 'IMAGE   '           / Image extension                                
PCOUNT  =                    0 / number of parameters                           
GCOUNT  =                    1 / number of groups                               
BUNIT   = 'erg Angstrom-1 s-1 cm-2'                                             
EXTNAME = 'FLUX    '           / extension name                                 
OBJSYS  = 'ICRS    '                                                            
IFUCON  = '16209   '           / NFibers                                        
WCSAXES =                    3 / Number of coordinate axes                      
CRPIX1  =                 57.5 / Pixel coordinate of reference point            
CRPIX2  =                 58.0 / Pixel coordinate of reference point            
CRPIX3  =                  1.0 / Pixel coordinate of reference point            
CDELT1  =  -0.0051388888888889 / [deg] Coordinate increment at reference point  
CDELT2  =   0.0051388888888889 / [deg] Coordinate increment at reference point  
**CDELT3  =   0.9999999999999999 / [m] Coordinate increment at reference point**    
CUNIT1  = 'deg'                / Units of coordinate increment and value        
CUNIT2  = 'deg'                / Units of coordinate increment and value        
CUNIT3  = 'Angstrom'           / Units of coordinate increment and value        
CTYPE1  = 'RA---TAN'           / Right ascension, gnomonic projection           
CTYPE2  = 'DEC--TAN'           / Declination, gnomonic projection               
CTYPE3  = 'WAVE'               / Vacuum wavelength (linear)                     
CRVAL1  =           134.861708 / [deg] Coordinate value at reference point      
CRVAL2  =           -43.737239 / [deg] Coordinate value at reference point      
CRVAL3  =    3599.999999999999 / [m] Coordinate value at reference point        
LONPOLE =                180.0 / [deg] Native longitude of celestial pole       
LATPOLE =           -43.737239 / [deg] Native latitude of celestial pole        
WCSNAME = 'DEFAULTS'           / Coordinate system title                        
MJDREF  =                  0.0 / [d] MJD of fiducial time                       
RADESYS = 'FK5'                / Equatorial coordinate system                   
EQUINOX =               2000.0 / [yr] Equinox of equatorial coordinates

Here is an example of incorrect units, but correct spectral axis:

# read the FITS file and produce Specral Cube
hdul         = fits.open(fits_dir+fits_file+'.fits') # for the SpectralCube function
hdr          = hdul[0].header
hdr['CD3_3'] = hdr['CDELT3'] # to repair value misread from the FITS header. It seems it was looking for 'CD3_3', as it was already using 'CD1_1', and 'CD2_2'. instead of CDELT1, CDELT2
w            = WCS(hdr)
data         = hdul[0].data
cube         = SpectralCube(data=hdul[0].data, header=hdr, wcs=w) 
print(cube)
print(repr(cube.header))

**SpectralCube with shape=(12401, 115, 115):**
 n_x:    115  type_x: RA---TAN  unit_x: deg    range:   134.450737 deg:  135.265531 deg
 n_y:    115  type_y: DEC--TAN  unit_y: deg    range:   -44.029442 deg:  -43.443603 deg
 n_s:  12401  type_s: WAVE      unit_s: Angstrom  range:     3600.000 Angstrom:    9800.000 Angstrom

SIMPLE  =                    T / conforms to FITS standard                      
BITPIX  =                  -64 / array data type                                
NAXIS   =                    3                                                  
NAXIS1  =                  115                                                  
NAXIS2  =                  115                                                  
NAXIS3  =                12401                                                  
EXTEND  =                    T                                                  
XTENSION= 'IMAGE   '           / Image extension                                
PCOUNT  =                    0 / number of parameters                           
GCOUNT  =                    1 / number of groups                               
BUNIT   = ''                                                                    
EXTNAME = 'FLUX    '           / extension name                                 
OBJSYS  = 'ICRS    '                                                            
IFUCON  = '16209   '           / NFibers                                        
WCSAXES =                    3 / Number of coordinate axes                      
CRPIX1  =                 57.5 / Pixel coordinate of reference point            
CRPIX2  =                 58.0 / Pixel coordinate of reference point            
CRPIX3  =                  1.0 / Pixel coordinate of reference point            
CDELT1  =  -0.0051388888888889 / [deg] Coordinate increment at reference point  
CDELT2  =   0.0051388888888889 / [deg] Coordinate increment at reference point  
**CDELT3  =  0.49999999999999994 / [m] Coordinate increment at reference point**    
CUNIT1  = 'deg'                / Units of coordinate increment and value        
CUNIT2  = 'deg'                / Units of coordinate increment and value        
CUNIT3  = 'Angstrom'           / Units of coordinate increment and value        
CTYPE1  = 'RA---TAN'           / Right ascension, gnomonic projection           
CTYPE2  = 'DEC--TAN'           / Declination, gnomonic projection               
CTYPE3  = 'WAVE'               / Vacuum wavelength (linear)                     
CRVAL1  =           134.861708 / [deg] Coordinate value at reference point      
CRVAL2  =           -43.737239 / [deg] Coordinate value at reference point      
CRVAL3  =    3599.999999999999 / [m] Coordinate value at reference point        
LONPOLE =                180.0 / [deg] Native longitude of celestial pole       
LATPOLE =           -43.737239 / [deg] Native latitude of celestial pole        
WCSNAME = 'DEFAULTS'           / Coordinate system title                        
MJDREF  =                  0.0 / [d] MJD of fiducial time                       
RADESYS = 'FK5'                / Equatorial coordinate system                   
EQUINOX =               2000.0 / [yr] Equinox of equatorial coordinates
e-koch commented 9 months ago

I would alter the header in the HDU before passing it into spectral-cube. The SpectralCube object isn't expecting to have state changes made after instantiating, so while you could hack in changes to the header, it may cause other unintended breaks.

e.g.,

hdul[0].header['CD3_3'] = hdr['CDELT3']
## Any other changes to the hdu, units, header here
cube = SpectralCube.read(hdul)
keflavich commented 9 months ago

You need to have only CD or only CDELT, so you probably want:

hdul         = fits.open(fits_dir+fits_file+'.fits') # for the SpectralCube function
hdr          = hdul[0].header
hdr['CD3_3'] = hdr['CDELT3'] 
del hdr['CDELT3']
cube         = SpectralCube.read(hdul)

I'll echo what @e-koch said: don't try to modify the cube after you've loaded it, it is designed not to be changed.

astro-pablo commented 9 months ago

You need to have only CD or only CDELT, so you probably want:

hdul         = fits.open(fits_dir+fits_file+'.fits') # for the SpectralCube function
hdr          = hdul[0].header
hdr['CD3_3'] = hdr['CDELT3'] 
del hdr['CDELT3']
cube         = SpectralCube.read(hdul)

I'll echo what @e-koch said: don't try to modify the cube after you've loaded it, it is designed not to be changed.

This worked!. Thanks


SpectralCube with shape=(12401, 115, 115) and unit=ergs/s/cm^2/A:
 n_x:    115  type_x: RA---TAN  unit_x: deg    range:   134.450737 deg:  135.265531 deg
 n_y:    115  type_y: DEC--TAN  unit_y: deg    range:   -44.029442 deg:  -43.443603 deg
 n_s:  12401  type_s: WAVE      unit_s: Angstrom  range:     3600.000 Angstrom:    9800.000 Angstrom
SIMPLE  =                    T / conforms to FITS standard                      
BITPIX  =                  -64 / array data type                                
NAXIS   =                    3                                                  
NAXIS1  =                  115                                                  
NAXIS2  =                  115                                                  
NAXIS3  =                12401                                                  
EXTEND  =                    T                                                  
XTENSION= 'IMAGE   '           / Image extension                                
PCOUNT  =                    0 / number of parameters                           
GCOUNT  =                    1 / number of groups                               
BUNIT   = 'erg Angstrom-1 s-1 cm-2'                                             
EXTNAME = 'FLUX    '           / extension name                                 
OBJSYS  = 'ICRS    '                                                            
IFUCON  = '16209   '           / NFibers                                        
WCSAXES =                    3 / Number of coordinate axes                      
CRPIX1  =                 57.5 / Pixel coordinate of reference point            
CRPIX2  =                 58.0 / Pixel coordinate of reference point            
CRPIX3  =                  1.0 / Pixel coordinate of reference point            
CDELT1  =  -0.0051388888888889 / [deg] Coordinate increment at reference point  
CDELT2  =   0.0051388888888889 / [deg] Coordinate increment at reference point  
CDELT3  =  0.49999999999999994 / [m] Coordinate increment at reference point    
CUNIT1  = 'deg'                / Units of coordinate increment and value        
CUNIT2  = 'deg'                / Units of coordinate increment and value        
CUNIT3  = 'Angstrom'           / Units of coordinate increment and value        
CTYPE1  = 'RA---TAN'           / Right ascension, gnomonic projection           
CTYPE2  = 'DEC--TAN'           / Declination, gnomonic projection               
CTYPE3  = 'WAVE'               / Vacuum wavelength (linear)                     
CRVAL1  =           134.861708 / [deg] Coordinate value at reference point      
CRVAL2  =           -43.737239 / [deg] Coordinate value at reference point      
CRVAL3  =    3599.999999999999 / [m] Coordinate value at reference point        
LONPOLE =                180.0 / [deg] Native longitude of celestial pole       
LATPOLE =           -43.737239 / [deg] Native latitude of celestial pole        
WCSNAME = 'DEFAULTS'           / Coordinate system title                        
MJDREF  =                  0.0 / [d] MJD of fiducial time                       
RADESYS = 'FK5'                / Equatorial coordinate system                   
EQUINOX =               2000.0 / [yr] Equinox of equatorial coordinates   
``