Updating FITS. headers via Sunpy/astropy

[Tom Bridgman]

I understand that SunPy FITS and AstroPy FITS do not exactly play nice together.

I've got some SDO/AIA FITS files that I want to do the equivalent of Level 1 to Level 1.5 conversion and save the resulting FITS files WITH appropriately updated HISTORY keywords (similar to the SolarSoft process).

I've tried a number of options with no success.  Often I get this mysterious message when I try to add the 'headerstring'

WARNING: SunpyMetadataWarning: The meta key comment is not valid ascii, dropping from the FITS header [sunpy.io.fits]

Yet the string seems to be perfectly fine ASCII.

Any pointers for what I'm missing?  Code snippets below of what I've tried so far.

Thanks,

Tom

----------------

Using 

from astropy.io import fits

import sunpy.map

import sunpy.io

import aiapy

from aiapy.calibrate import update_pointing

# Read FITS via astropy

hdus = fits.open(dataFile) # read extension and header

# convert to Sunpy Map

dataViewLevel1 = sunpy.map.Map(hdus[1].data, hdus[1].header)

# update pointing info

dataViewPupdate = update_pointing(dataViewLevel1) # update pointing information

headerstring = "AIApy (v%s): Update pointing"%str(aiapy.__version__)

print(headerstring)

# neither of the two lines below work

dataViewPupdate.fits_header['history'] = headerstring

#dataViewPupdate.fits_header.set('history', headerstring)

header = dataViewPupdate.fits_header

data = dataViewPupdate.data

# tried the two options below for writing the file. Either fails with error or doesn't update the header.

# sunpy.io.fits.write('data/TestOutput.fits',

# header, data,

# overwrite=True,

# hdu_type=CompImageHDU,

# checksum=True)

# optionally AstroPy FITS

fits.writeto('data/TestOutput.fits',hdus[0].data, hdus[0].header)

# append updated HDU

fits.append('data/TestOutput.fits', data,header)

[will.t...@gmail.com]

Hi Tom,

Thanks for your question. Below I've included a minimal working version of your code snippet. Note that here I am using astropy v5.0, sunpy v3.1.3 and aiapy v0.6.4.

```

import sunpy.map

import aiapy
from aiapy.calibrate import update_pointing

from astropy.io.fits import CompImageHDU

m_l1 = sunpy.map.Map(dataFile)
m_l1_pupdate = update_pointing(m_l1)
m_l1_pupdate.meta['history'] = f'AIApy (v{aiapy.__version__}): Update pointing'
m_l1_pupdate.save('data/TestOutput.fits', checksum=True, hdu_type=CompImageHDU, overwrite=True)

# Or using sunpy.io.fits
sunpy.io.fits.write('data/TestOutput.fits', m_l1_pupdate.data, m_l1_pupdate.fits_header, checksum=True, hdu_type=CompImageHDU, overwrite=True)
# Or using astropy.io.fits
CompImageHDU(m_l1_pupdate.data, m_l1_pupdate.fits_header).writetoto('data/TestOutput.fits', checksum=True, overwrite=True)

```

A few clarifying comments:

* All three ways of writing the file out to a FITS file should yield equivalent results.

* The `fits_header` property of a sunpy Map is immutable and is just a representation (in astropy.io.fits.Header format) of the `meta` attribute. This is why you're change was not updating the header. If you want to explicitly modify the underlying metadata, you should modify the keys in `map.meta`. This will then be reflected in `map.fits_header`.

* The most convenient way to write a Map to a FITS file is through the Map.save method. As I've shown above, most of the arguments that you can pass to sunpy.io.fits.write (and by extension astropy.io.fits.writeto) can also be passed to Map.save.

* Regarding your statement "I understand that SunPy FITS and AstroPy FITS do not exactly play nice together," sunpy.io.fits is a wrapper around astropy.io.fits that provides a lot of convenience when dealing with solar data. This convenience is at the cost of generality. The main purpose of sunpy.io.fits is to provide a FITS reader for 2D image data ingested through sunpy.map.Map and a writer for sunpy.map.Map.save. For a general purpose FITS reader or in cases where you want *very* fine-grained control over what gets written to your FITS file, astropy.io.fits is the recommended solution.

* I am unsure as to what is causing the SunpyMetadataWarning. I was not able to reproduce that behavior.

Best,

Will

[Tom Bridgman]

Thanks.

Is there anything about the meta method/attribute connection to updating headers in the sunpy docs? 

What I find appears to be for a different use.

Tom

[will.t...@gmail.com]

Hi Tom,

It appears we don't have much in the way of comprehensive documentation regarding manipulating metadata. There are some cursory comments about modifying keys in `.meta` here: https://docs.sunpy.org/en/stable/code_ref/map.html#fixing-map-metadata.

This omission is partially intentional. The majority of the attributes on map are derived from the underlying metadata as stored in `.meta`. These include those related to the WCS (e.g. `.reference_coordinate`), the coordinate frame (e.g. `observer_coordinate`), or other properties such as the observatory or processing level. By modifying the underlying  metadata, there is a risk of (incorrectly) altering these higher-level properties. Additionally, there is also no guarantee that the modified metadata is now serializable back to FITS.

Of course, there are cases, such as the one you've described here, where modifying a key in the metadata is innocuous and also helpful. The issue with the way our metadata is currently structured is that everything is in one dictionary (`.meta`) and the distinction between keys that shouldn't be (or at least not without great care) manually manipulated (e.g. CRPIX1) and keys that are easily modified (e.g. HISTORY, COMMENTS) is not clear. How exactly to alter out metadata model to make it "safer" but also more flexible (e.g. providing some documented entrypoint for appending to the history) has been an ongoing discussion for a while now. If you have any thoughts or suggestions on this topic, any feedback is greatly appreciated, either here or as an issue on GitHub.

Best,

Will

[Tom Bridgman]

Finally got to test this fix in my code.  Because I was adding *multiple* HISTORY keywords to indicate the different correction steps, only the last one survived to be written to the file.

I tried several variations using the different file writers with still no luck.

This is clearly an issue with maintaining the FITS standard in Python.

[will.t...@gmail.com]

Hi Tom,

The problem you're describing is actually due to the way sunpy handles metadata attached to a map rather than some deeper incompatibility with Python/astropy and the FITS standard.

The `.meta` attribute is a (subclass of a) dictionary which, by definition, does not permit more than one key of the same name. As you point out, the FITS standard does explicitly allow for any number of occurrences of HISTORY and COMMENT. Astropy does actually provide a special case for these keywords in their astropy.io.fits.Header object through the add_history and add_comment methods. Below, I've included a modified version of my original script that uses this to append two HISTORY entries and then save the file out in the same way.

```

import sunpy.map

import aiapy
from aiapy.calibrate import update_pointing, register

from astropy.io.fits import CompImageHDU

m_l1 = sunpy.map.Map(dataFile)
m_l1_pupdate = update_pointing(m_l1)

m_l15 = register(m_l1_pupdate)
header = m_l15.fits_header
header.add_history(f'AIApy (v{aiapy.__version__}): Update pointing')
header.add_history(f'AIApy (v{aiapy.__version__}): Register image')
CompImageHDU(m_l15.data, header).writeto('data/TestOutput.fits', checksum=True, overwrite=True)

```

The resulting file now includes both history entries. Note that if you read this file directly into a sunpy Map, the two keys will be combined into one history entry in the `.meta` attribute. I hope that is helpful!

Additionally, there are a number of open issues on both sunpy and aiapy related to this problem that may be of interest to you and are related to the larger problem of a more advanced metadata object that I mentioned in my previous message:

* https://github.com/sunpy/sunpy/issues/997

* https://github.com/sunpy/sunpy/issues/3876

* https://gitlab.com/LMSAL_HUB/aia_hub/aiapy/-/issues/56

Best,

Will

[Tom Bridgman]

Ah!  THAT worked.  Thanks.

It seems that you cannot update headers when they are still within a Map object, yet you get no error message when you try.  I'm wondering if this is another post Python 1.5 feature that I never really grokked.

Good to see that others are concerned about preserving some modification history documentation.

One day I'll get back to activating the new security features on my GitHub account and maybe add some of these types of examples to the docs.

Tom