astropy / astropy

Astronomy and astrophysics core library
https://www.astropy.org
BSD 3-Clause "New" or "Revised" License
4.3k stars 1.72k forks source link

u.Magnitude() does not produce a FITS conform output #13457

Open olebole opened 1 year ago

olebole commented 1 year ago

Description

When creating a table column with u.mag() as unit, the FITS writer refuses to take it as a FITS mag standard column but serializes it in a astropy specific manner.

Expected behavior

I'd expect that both u.mag and u.mag() result in a FITS unit mag, as their difference is quite subtle.

Steps to Reproduce

>>> from astropy.table import QTable
>>> import astropy.units as u
>>> from astropy.io import fits
>>> QTable({'fmag': [1.0] * u.mag}).write('mag-units.fits')
>>> fits.open('mag-units.fits')[1].header['T*1']
TTYPE1  = 'fmag    '                                                            
TFORM1  = 'D       '                                                            
TUNIT1  = 'mag     '                                                            
>>> QTable({'fmag': [1.0] * u.mag()}).write('mag-magnitude.fits')
WARNING: The unit 'mag(1)' could not be saved in native FITS format and hence will be lost to non-astropy fits readers. Within astropy, the unit can roundtrip using QTable, though one has to enable the unit before reading. [astropy.io.fits.convenience]
>>> fits.open('mag-magnitude.fits')[1].header['T*1']
TTYPE1  = 'fmag    '                                                            
TFORM1  = 'D       '                                                            

System Details

mhvk commented 1 year ago

This is a tricky one: if we special-case the dimensionless magnitude, then the data will no longer roundtrip. That said, perhaps the best of both worlds would be to properly use TUNIT1 with mag and also serialize the astropy unit (which I think will generally have precedence on reading back in astropy).

olebole commented 1 year ago

I am still wondering what the difference between u.mag() and u.mag from the user perspective is. What is the use case for u.mag at all? If we could make u.mag a u.MagUnit(), wouldn't that resolve the problem?

mhvk commented 1 year ago

Part of the difficulty is in composite units like u.mag / u.day, which are only OK if the magnitude does not have a physical dimension. So, that suggests a mag unit that behaves like a regular unit is still useful. That said, in a way it is a hack, and perhaps a better approach would allow also things like STmag/arcsec^2... Overall, those logarithmic units remain super tricky...

For the problem at hand, I guess one question is whether by default a FITS column with mag unit should become a Magnitude.

olebole commented 1 year ago

I still don't see the difference between u.mag and the dimensionless u.mag(). And to me, STmag/arcsec² looks dangerously sloppy; IMO it should be written as mag(ST/arcsec²).

pllim commented 1 year ago

IMO it should be written as mag(ST/arcsec²).

I don't quite like wrapping arcsec2 inside mag(). The arcsec2 isn't in log unit, just the ST part... right?

olebole commented 1 year ago

@plim Not really, but the other way: ST is a flux unit; mag(ST) converts it to a magnitude. If one wants to convert a surface flux, the natural way is mag(ST/arcsec²). STmag/arcsec² is just lazy writing and not to be taken literally.

pllim commented 1 year ago

Ah, okay. Thanks for clarification and sorry for the noise!