New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Change map FITS unit header key to standard "BUNIT" #1755

Merged
merged 2 commits into from Sep 2, 2018

Conversation

1 participant
@cdeil
Member

cdeil commented Aug 31, 2018

This is a reminder issue for me to look check if we should support BUNIT on Map.read.

the CTA 1DC diffuse model was given like this:

$ ftlist $CTADATA/models/cube_iem.fits.gz K | grep UNIT
CUNIT1  = 'deg'                / Units of coordinate increment and value
CUNIT2  = 'deg'                / Units of coordinate increment and value
BUNIT   = 'photon/cm2/s/MeV/sr' / Photon flux
TUNIT1  = 'MeV     '

and currently Map.read doesn't read BUNIT:

>>> from gammapy.maps import Map
>>> m = Map.read('$CTADATA/models/cube_iem.fits.gz')
>>> m.unit
Unit(dimensionless)

@robertazanin @facero @registerrier - For now, when using the diffuse model, please always load it via SkyDiffuseCube.read - there we set the default unit and the model should come out OK to use with MapEvaluator and MapFit.

m.unit = 'cm-2 s-1 MeV-1 sr-1'

@cdeil cdeil added the question label Aug 29, 2018

@cdeil cdeil added this to the 0.8 milestone Aug 29, 2018

@cdeil cdeil self-assigned this Aug 29, 2018

@cdeil

This comment has been minimized.

Show comment
Hide comment
@cdeil

cdeil Aug 30, 2018

Member

It looks like gammapy.maps uses a FITS header key of UNIT:

>>> from gammapy.maps import Map
>>> Map.create(npix=(10, 5), unit='cm-2 s-1').to_hdulist()[0].header
SIMPLE  =                    T / conforms to FITS standard                      
BITPIX  =                  -32 / array data type                                
NAXIS   =                    2 / number of array dimensions                     
NAXIS1  =                   10                                                  
NAXIS2  =                    5                                                  
WCSAXES =                    2 / Number of coordinate axes                      
CRPIX1  =                  5.5 / Pixel coordinate of reference point            
CRPIX2  =                  3.0 / Pixel coordinate of reference point            
CDELT1  =                 -0.1 / [deg] Coordinate increment at reference point  
CDELT2  =                  0.1 / [deg] Coordinate increment at reference point  
CUNIT1  = 'deg'                / Units of coordinate increment and value        
CUNIT2  = 'deg'                / Units of coordinate increment and value        
CTYPE1  = 'RA---CAR'           / Right ascension, plate caree projection        
CTYPE2  = 'DEC--CAR'           / Declination, plate caree projection            
CRVAL1  =                  0.0 / [deg] Coordinate value at reference point      
CRVAL2  =                  0.0 / [deg] Coordinate value at reference point      
LONPOLE =                  0.0 / [deg] Native longitude of celestial pole       
LATPOLE =                 90.0 / [deg] Native latitude of celestial pole        
RADESYS = 'ICRS'               / Equatorial coordinate system                   
WCSSHAPE= '(10,5)  '                                                            
META    = '{}      '                                                            
UNIT    = 'cm-2 s-1'                                                            

Reading back works of course:

>>> from gammapy.maps import Map
>>> h = Map.create(npix=(10, 5), unit='cm-2 s-1').to_hdulist()
>>> Map.from_hdulist(h).unit
Unit("1 / (cm2 s)")

It looks like we did support BUNIT in SkyImage:
https://github.com/gammapy/gammapy/blob/v0.7/gammapy/image/core.py#L213

And we did support SPECUNIT in SkyCube:
https://github.com/gammapy/gammapy/blob/v0.7/gammapy/cube/core.py#L759

There is no mention of BUNIT or SPECUNIT currently in the Gammapy code base.

Looking at the FITS standard, as far as I can see it only mentions BUNIT and neither UNIT or SPECUNIT as valid keywords:
https://fits.gsfc.nasa.gov/standard40/fits_standard40aa-le.pdf

BUNIT keyword. The value field shall contain a character string describing the physical units in which the quantities in the array ...

So I'm inclined to change gammapy.maps to just use BUNIT instead of UNIT.

Of course we could also support multiple keys, e.g. {BUNIT, UNIT, SPECUNIT} on read, or even have the write configurable, with BUNIT as default. But unless there's a strong reason this is needed, I think it will just create problems and we shouldn't.

@cboisson @woodmd - Please comment!

Member

cdeil commented Aug 30, 2018

It looks like gammapy.maps uses a FITS header key of UNIT:

>>> from gammapy.maps import Map
>>> Map.create(npix=(10, 5), unit='cm-2 s-1').to_hdulist()[0].header
SIMPLE  =                    T / conforms to FITS standard                      
BITPIX  =                  -32 / array data type                                
NAXIS   =                    2 / number of array dimensions                     
NAXIS1  =                   10                                                  
NAXIS2  =                    5                                                  
WCSAXES =                    2 / Number of coordinate axes                      
CRPIX1  =                  5.5 / Pixel coordinate of reference point            
CRPIX2  =                  3.0 / Pixel coordinate of reference point            
CDELT1  =                 -0.1 / [deg] Coordinate increment at reference point  
CDELT2  =                  0.1 / [deg] Coordinate increment at reference point  
CUNIT1  = 'deg'                / Units of coordinate increment and value        
CUNIT2  = 'deg'                / Units of coordinate increment and value        
CTYPE1  = 'RA---CAR'           / Right ascension, plate caree projection        
CTYPE2  = 'DEC--CAR'           / Declination, plate caree projection            
CRVAL1  =                  0.0 / [deg] Coordinate value at reference point      
CRVAL2  =                  0.0 / [deg] Coordinate value at reference point      
LONPOLE =                  0.0 / [deg] Native longitude of celestial pole       
LATPOLE =                 90.0 / [deg] Native latitude of celestial pole        
RADESYS = 'ICRS'               / Equatorial coordinate system                   
WCSSHAPE= '(10,5)  '                                                            
META    = '{}      '                                                            
UNIT    = 'cm-2 s-1'                                                            

Reading back works of course:

>>> from gammapy.maps import Map
>>> h = Map.create(npix=(10, 5), unit='cm-2 s-1').to_hdulist()
>>> Map.from_hdulist(h).unit
Unit("1 / (cm2 s)")

It looks like we did support BUNIT in SkyImage:
https://github.com/gammapy/gammapy/blob/v0.7/gammapy/image/core.py#L213

And we did support SPECUNIT in SkyCube:
https://github.com/gammapy/gammapy/blob/v0.7/gammapy/cube/core.py#L759

There is no mention of BUNIT or SPECUNIT currently in the Gammapy code base.

Looking at the FITS standard, as far as I can see it only mentions BUNIT and neither UNIT or SPECUNIT as valid keywords:
https://fits.gsfc.nasa.gov/standard40/fits_standard40aa-le.pdf

BUNIT keyword. The value field shall contain a character string describing the physical units in which the quantities in the array ...

So I'm inclined to change gammapy.maps to just use BUNIT instead of UNIT.

Of course we could also support multiple keys, e.g. {BUNIT, UNIT, SPECUNIT} on read, or even have the write configurable, with BUNIT as default. But unless there's a strong reason this is needed, I think it will just create problems and we shouldn't.

@cboisson @woodmd - Please comment!

@cdeil

This comment has been minimized.

Show comment
Hide comment
@cdeil

cdeil Aug 30, 2018

Member

The maps spec doesn't mention how to store unit, and the example files contain no unit:
https://github.com/open-gamma-ray-astro/gamma-astro-data-formats/blob/master/source/skymaps/wcs/make_wcs_files.py
Suggest to add a mention there that the FITS standard BUNIT should be used, and add a unit to the example files.

This package seems to use BUNIT only, not read or write UNIT or SPECUNIT:
https://github.com/radio-astro-tools/spectral-cube/blob/e9c6c97fda60f460c09e8dc1c25c04578a5ad49e/spectral_cube/cube_utils.py#L431

If anyone has files that contain "UNIT" or "SPECUNIT" (e.g. from the Fermi ST or other missions) instead fo "BUNIT" that you want to use with Gammapy, please mention it here so that we know whether to support that on read.

Member

cdeil commented Aug 30, 2018

The maps spec doesn't mention how to store unit, and the example files contain no unit:
https://github.com/open-gamma-ray-astro/gamma-astro-data-formats/blob/master/source/skymaps/wcs/make_wcs_files.py
Suggest to add a mention there that the FITS standard BUNIT should be used, and add a unit to the example files.

This package seems to use BUNIT only, not read or write UNIT or SPECUNIT:
https://github.com/radio-astro-tools/spectral-cube/blob/e9c6c97fda60f460c09e8dc1c25c04578a5ad49e/spectral_cube/cube_utils.py#L431

If anyone has files that contain "UNIT" or "SPECUNIT" (e.g. from the Fermi ST or other missions) instead fo "BUNIT" that you want to use with Gammapy, please mention it here so that we know whether to support that on read.

@cdeil cdeil added the bug label Aug 30, 2018

@cdeil cdeil added this to To do in Map analysis via automation Aug 30, 2018

@cdeil

This comment has been minimized.

Show comment
Hide comment
@cdeil

cdeil Aug 31, 2018

Member

I've made the change from UNIT to BUNIT, attached in commit e6a31fb here.

The change is covered by a test for both HPX and WCS already.

This change caused one fail: in test_compute_lima_on_off_image an image is read with an invalid unit:

fits.open('../gammapy-extra/test_datasets/unbundled/hess/survey/hess_survey_snippet.fits.gz')['ON'].header
BUNIT   = 'Count   '

It fails with a ValueError with a very complex traceback:
https://gist.github.com/cdeil/83b3807db856f6e1a438355365190370
I'll complain about this in Astropy.

But the question remains: what should we do for maps with invalid units?

I think astropy.table.Table.read will always emit a warning and move on, but e.g. astropy.nddata.CCDData will raise a nice error mentioning BUNIT and how to fix it.

from astropy.nddata import CCDData
CCDData.read('../gammapy-extra/test_datasets/unbundled/hess/survey/hess_survey_snippet.fits.gz')

@adonath @woodmd @registerrier - Thoughs?

I'm +1 to emit a warning and put no unit in the map. It's not the users fault often if they get bad data, so I think we should be "lenient on read and strict on write".

Member

cdeil commented Aug 31, 2018

I've made the change from UNIT to BUNIT, attached in commit e6a31fb here.

The change is covered by a test for both HPX and WCS already.

This change caused one fail: in test_compute_lima_on_off_image an image is read with an invalid unit:

fits.open('../gammapy-extra/test_datasets/unbundled/hess/survey/hess_survey_snippet.fits.gz')['ON'].header
BUNIT   = 'Count   '

It fails with a ValueError with a very complex traceback:
https://gist.github.com/cdeil/83b3807db856f6e1a438355365190370
I'll complain about this in Astropy.

But the question remains: what should we do for maps with invalid units?

I think astropy.table.Table.read will always emit a warning and move on, but e.g. astropy.nddata.CCDData will raise a nice error mentioning BUNIT and how to fix it.

from astropy.nddata import CCDData
CCDData.read('../gammapy-extra/test_datasets/unbundled/hess/survey/hess_survey_snippet.fits.gz')

@adonath @woodmd @registerrier - Thoughs?

I'm +1 to emit a warning and put no unit in the map. It's not the users fault often if they get bad data, so I think we should be "lenient on read and strict on write".

@cdeil cdeil merged commit 2c4e535 into gammapy:master Sep 2, 2018

0 of 2 checks passed

continuous-integration/appveyor/pr Waiting for AppVeyor build to complete
Details
continuous-integration/travis-ci/pr The Travis CI build is in progress
Details

Map analysis automation moved this from To do to Done Sep 2, 2018

@cdeil cdeil changed the title from Support BUNIT FITS header key? to Change map FITS unit header key to standard "BUNIT" Sep 9, 2018

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment