Fix sunpy.io to update header info after scaling applied when data is accessed

changelog/7141.bugfix.rst

@@ -0,0 +1 @@
+Fixed a bug where fits headers did not match the current state of the corresponding astropy.io.fits headers.

sunpy/data/test/__init__.py

@@ -112,7 +112,9 @@ def get_dummy_map_from_header(filename):
     """
     filepath = get_test_filepath(filename)
     header = _fits.format_comments_and_history(astropy.io.fits.Header.fromtextfile(filepath))
-    data = np.random.rand(header['NAXIS2'], header['NAXIS1'])
+    data = np.random.randint(low=0, high=100, size = (header['NAXIS2'], header['NAXIS1']))
+    if "BSCALE" in header:
+        data = data * header["BSCALE"] + header["BZERO"]
     if 'BITPIX' in header:
         data = data.astype(astropy.io.fits.BITPIX2DTYPE[header['BITPIX']])
     # NOTE: by reading straight from the data header pair, we are skipping

sunpy/data/test/tests/test_funcs.py

@@ -0,0 +1,17 @@
+from pathlib import Path
+
+import numpy as np
+
+from sunpy.data.test import get_dummy_map_from_header, get_test_filepath
+
+TEST_IRIS_SJI_HEADER = get_test_filepath('iris_l2_20130801_074720_4040000014_SJI_1400_t000.header')
+
+def test_get_dummy_map_from_header_b_scale_zero_applied():
+    # We use the IRIS one since it contains a BZERO and BSCALE in the test header
+    iris_map = get_dummy_map_from_header(Path(TEST_IRIS_SJI_HEADER))
+    assert iris_map.meta["bscale"] == 0.25
+    assert iris_map.meta["bzero"] == 7992
+    # Data from the test header is scaled by BSCALE and BZERO
+    # but the max value it will return is 100 and converted to int
+    # So its 24 instead of 25 due to the rounding
+    assert np.max(iris_map.data) == 8016

sunpy/io/_fits.py

@@ -44,14 +44,17 @@
 
 def read(filepath, hdus=None, memmap=None, **kwargs):
     """
-    Read a fits file.
+    Read a FITS file.
 
     Parameters
     ----------
     filepath : `str`
         The fits file to be read.
     hdus : `int` or iterable
         The HDU indexes to read from the file.
+    memmap : `bool`, optional
+        Set to use memory mapping.
+        Default is `None`.
     **kwargs : `dict`, optional
         Passed to `astropy.io.fits.open`.
 
@@ -74,17 +77,16 @@ def read(filepath, hdus=None, memmap=None, **kwargs):
                 hdulist = hdulist[hdus]
             elif isinstance(hdus, collections.abc.Iterable):
                 hdulist = [hdulist[i] for i in hdus]
-
         hdulist = fits.hdu.HDUList(hdulist)
         for h in hdulist:
             h.verify('silentfix+warn')
-
-        headers = get_header(hdulist)
         pairs = []
-
-        for i, (hdu, header) in enumerate(zip(hdulist, headers)):
+        for i, hdu in enumerate(hdulist):
             try:
-                pairs.append(HDPair(hdu.data, header))
+                # Access data to force any BSCALE/BZERO scaling to be applied
+                # which will update the keywords in the header.
+                data = hdu.data
+                pairs.append(HDPair(data, get_header(hdu)[-1]))
             except (KeyError, ValueError) as e:
                 message = f"Error when reading HDU {i}. Skipping.\n"
                 for line in traceback.format_tb(sys.exc_info()[2]):
@@ -98,12 +100,12 @@ def read(filepath, hdus=None, memmap=None, **kwargs):
 
 def get_header(afile):
     """
-    Read a fits file and return just the headers for all HDU's.
+    Gets the header from a FITS file or HDUList or individual HDU.
 
     Parameters
     ----------
-    afile : `str` or `astropy.io.fits.HDUList`
-        The file to be read, or HDUList to process.
+    afile : `str` or `astropy.io.fits.HDUList` or `astropy.io.fits.PrimaryHDU`
+        The file to be read or HDUList or PrimaryHDU to process.
 
     Returns
     -------
@@ -113,15 +115,15 @@ def get_header(afile):
     if isinstance(afile, fits.HDUList):
         hdulist = afile
         close = False
+    elif isinstance(afile, (fits.PrimaryHDU | fits.hdu.base.ExtensionHDU)):
+        hdulist = [afile]
+        close = False
     else:
         hdulist = fits.open(afile, ignore_blank=True)
         hdulist.verify('silentfix')
         close = True
-
     try:
-        headers = []
-        for hdu in hdulist:
-            headers.append(format_comments_and_history(hdu.header))
+        headers = [format_comments_and_history(hdu.header) for hdu in hdulist]
     finally:
         if close:
             hdulist.close()
@@ -156,11 +158,11 @@ def format_comments_and_history(input_header):
     header['COMMENT'] = comment
     header['HISTORY'] = history
 
-    # Strip out KEYCOMMENTS to a dict, the hard way
-    keydict = {}
-    for card in input_header.cards:
-        if card.comment != '':
-            keydict.update({card.keyword: card.comment})
+    keydict = {
+        card.keyword: card.comment
+        for card in input_header.cards
+        if card.comment != ''
+    }
     header['KEYCOMMENTS'] = keydict
     waveunit = extract_waveunit(header)
     if waveunit is not None:

sunpy/io/tests/test_fits.py

@@ -17,8 +17,7 @@
 TEST_EIT_HEADER = get_test_filepath('EIT_header/efz20040301.000010_s.header')
 TEST_SWAP_HEADER = get_test_filepath('SWAP/resampled1_swap.header')
 TEST_GONG_HEADER = get_test_filepath('gong_synoptic.header')
-# Some of the tests images contain an invalid BLANK keyword
-# ignore the warning raised by this
+# Some of the tests images contain an invalid BLANK keyword ignore the warning raised by this
 pytestmark = pytest.mark.filterwarnings("ignore:Invalid 'BLANK' keyword in header")
 
 
@@ -133,8 +132,6 @@ def test_write_with_metadict_header_astropy(tmpdir):
 
 # Various warnings are thrown in this test, but we just want to check that the code
 # works without exceptions
-
-
 @pytest.mark.filterwarnings('ignore')
 def test_fitsheader():
     """Test that all test data can be converted back to a FITS header."""
@@ -180,6 +177,7 @@ def test_warn_longkey():
 def test_read_memmap():
     data, _ = _fits.read(TEST_AIA_IMAGE, memmap=True)[0]
     assert data.base is not None
+    # Simple check to see if the base is a memory map
     assert isinstance(data.base, mmap.mmap)
 
     data, _ = _fits.read(TEST_AIA_IMAGE, memmap=False)[0]