Repository: astropy/astropy

Issue:
New FITS compression code breaks loading Pan-STARRS images 
### Description

The code for reading compressed FITS files in 5.3 is crashing on loading the images from Pan-STARRS. It was working nicely in 5.2.

The problem manifests like that:
```
$ python -c 'from astropy.io import fits; fits.getdata("http://ps1images.stsci.edu/rings.v3.skycell/1004/016/rings.v3.skycell.1004.016.stk.g.unconv.fits")'
Traceback (most recent call last):
  File "<string>", line 1, in <module>
  File "/ltzfs/ASTRO101/anaconda3/envs/astro101/lib/python3.11/site-packages/astropy/io/fits/convenience.py", line 232, in getdata
    data = hdu.data
           ^^^^^^^^
  File "/ltzfs/ASTRO101/anaconda3/envs/astro101/lib/python3.11/site-packages/astropy/utils/decorators.py", line 837, in __get__
    val = self.fget(obj)
          ^^^^^^^^^^^^^^
  File "/ltzfs/ASTRO101/anaconda3/envs/astro101/lib/python3.11/site-packages/astropy/io/fits/hdu/compressed.py", line 1558, in data
    data = self.section[...]
           ~~~~~~~~~~~~^^^^^
  File "/ltzfs/ASTRO101/anaconda3/envs/astro101/lib/python3.11/site-packages/astropy/io/fits/hdu/compressed.py", line 2189, in __getitem__
    data = decompress_hdu_section(self.hdu, first_tile_index, last_tile_index)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/ltzfs/ASTRO101/anaconda3/envs/astro101/lib/python3.11/site-packages/astropy/io/fits/_tiled_compression/tiled_compression.py", line 442, in decompress_hdu_section
    tile_data[blank_mask] = np.nan
    ~~~~~~~~~^^^^^^^^^^^^
ValueError: cannot convert float NaN to integer
```

### Expected behavior

The images from Pan-STARRS were loading without issues in pre-5.3

### How to Reproduce

```python
from astropy.io import fits
fits.getdata("http://ps1images.stsci.edu/rings.v3.skycell/1004/016/rings.v3.skycell.1004.016.stk.g.unconv.fits")
```


### Versions

Linux-4.18.0-372.9.1.el8.x86_64-x86_64-with-glibc2.28
Python 3.11.4 (main, Jul  5 2023, 13:45:01) [GCC 11.2.0]
astropy 5.3.2
Numpy 1.25.2
pyerfa 2.0.0.3
Scipy 1.11.2
Matplotlib 3.7.2


Discussion:
github-actions[bot]: Thank you for your contribution to Astropy! 🌌 This checklist is meant to remind the package maintainers who will review this pull request of some common things to look for.

  - [ ] Do the proposed changes actually accomplish desired goals?
  - [ ] Do the proposed changes follow the [Astropy coding guidelines](https://docs.astropy.org/en/latest/development/codeguide.html)?
  - [ ] Are tests added/updated as required? If so, do they follow the [Astropy testing guidelines](https://docs.astropy.org/en/latest/development/testguide.html)?
  - [ ] Are docs added/updated as required? If so, do they follow the [Astropy documentation guidelines](https://docs.astropy.org/en/latest/development/docguide.html#astropy-documentation-rules-and-guidelines)?
  - [ ] Is rebase and/or squash necessary? If so, please provide the author with appropriate instructions. Also see instructions for [rebase](https://docs.astropy.org/en/latest/development/workflow/development_workflow.html#rebase-if-necessary) and [squash](https://docs.astropy.org/en/latest/development/workflow/development_workflow.html#squash-if-necessary).
  - [ ] Did the CI pass? If no, are the failures related? If you need to run daily and weekly cron jobs as part of the PR, please apply the "Extra CI" label. Codestyle issues can be fixed by the [bot](https://docs.astropy.org/en/latest/development/workflow/development_workflow.html#pre-commit).
  - [ ] Is a change log needed? If yes, did the change log check pass? If no, add the "no-changelog-entry-needed" label. If this is a manual backport, use the "skip-changelog-checks" label unless special changelog handling is necessary.
  - [ ] Is this a big PR that makes a "What's new?" entry worthwhile and if so, is (1) a "what's new" entry included in this PR and (2) the "whatsnew-needed" label applied?
  - [ ] At the time of adding the milestone, if the milestone set requires a backport to release branch(es), apply the appropriate "backport-X.Y.x" label(s) *before* merge.
neutrinoceros: This would supersede #15756 and fix #15236
@rlwastro please make sure formatting is green (you may comment `pre-commit.ci autofix` here or fix it locally and push), and please include a changelog entry (see `docs/changes/README.rst`)
rlwastro: > This would supersede #15756 and fix #15236 @rlwastro please make sure formatting is green (you may comment `pre-commit.ci autofix` here or fix it locally and push), and please include a changelog entry (see `docs/changes/README.rst`)

I fixed the formatting and added a changelog entry.  Let me know if you see other changes that are required.

I'm working on creating a small compressed image that can be used for testing.
rlwastro: @neutrinoceros I'm testing this fix on various images and am beginning to wonder whether I have the correct fix for this issue.

Currently if we read a FITS file with an _uncompressed_ integer image that has a `BLANK` keyword in the header, the image gets converted from an integer datatype to a `float32`, and the blanked pixels are set to `np.nan`.  I could argue with that behavior (perhaps a masked integer array would be better than a `float32` array).  But that is the way it acts now, and as far as I know people are OK with it.

With my version of the code, a FITS file with a _compressed_ integer image that has a `BLANK` keyword in the header does _not_ get converted to a float.  Instead it is kept as an integer array, with the blank values retained as the special integer value.

I think for consistency with the uncompressed image behavior, this should instead convert the array to `float32` and replace the blanks with `np.nan`.  If you agree, I can withdraw this pull request and submit a new one.  But let me know if you think it needs more discussion before doing that.

This was a bit tricky to discover because the Pan-STARRS images in the original bug report also have `BZERO` and `BSCALE` keywords, and it turns out they do in the end get converted to `float32` with `np.nan` for blanks.  So the change I'm proposing will in the end give the same result for the PS images.  But I'm pretty sure the new approach would be better for other images.

I do think this is important to fix soon!  Currently it simply dies with an exception when one of these images is opened.  That is clearly bad.
neutrinoceros: > I think for consistency with the uncompressed image behavior, this should instead convert the array to float32 and replace the blanks with np.nan. If you agree, I can withdraw this pull request and submit a new one. But let me know if you think it needs more discussion before doing that.

I can see both ways as satisfying solutions, albeit non ideal. I think @saimn and/or @astrofrog should weigh in here.
rlwastro: @neutrinoceros I figured out the right place to make this change.  The code that applies the `BSCALE/BZERO` scaling and replaces `BLANK` values with `NaN` (in `compressed.py:CompImageHDU._scale_data()`) was skipped if the `BSCALE` and `BZERO` values were omitted.  It is now executed if `BLANK` is specified (even in the absence of `BSCALE/BZERO`).

I updated the description to reflect these additional changes.  This version should have the correct behavior (as similar to uncompressed image HDUs as possible).

There is room for improvement to the handling of errors in compressed image headers.  I have not tried to make these changes because the code in `header.py` is pretty scary.  My changes do not make these problems any worse at least!
* For integer images, if the `ZBLANK` header keyword is specified but there is no `BLANK` keyword, `ZBLANK` should be copied to `BLANK` (with a warning).
* For integer images, if neither `ZBLANK` nor `BLANK` appears in the header, but there is a `ZBLANK` column in the table, then `BLANK` should be set to a default value (with a warning).
rlwastro: @neutrinoceros I made some further improvements in the code that cleaned up some remaining issues with unusual file formats.  I believe this patch is complete now.

I added new tests in the `test_tiled_compression.py` module that I believe should test all the added code.  I added some additional (small) test files as well to support those tests. There are still some warnings about code coverage.  I'm not sure what the problem is there.  I can't tell whether the new tests have actually been run yet or not.

I believe this is very close to complete, so if you or others can help get it through the final stages, that would be very helpful.
rlwastro: Thanks @pllim !  I fixed some filename typos in my test code, and now the patch coverage is 100% complete.  All checks now pass.
pllim: Alas, we just updated branch protection rules due to matplotlib bump. Would you mind rebasing? Thanks!
rlwastro: @pllim Ah crap, I think I screwed that up.  This is very much out of my comfort zone in git!  I'll see if I can fix it.  If not I may have to start over.
pllim: Maybe try this?

```
git fetch upstream main
git rebase upstream/main
```

And see if the history looks better before you force push?
rlwastro: Whew, thanks!  That did the trick.  I think things look like they will be OK now.  I'm just waiting to confirm the checks are all still OK.
rlwastro: > First of, I want to point out this seems like great work to me, though again I'm no authority when it comes to `io.fits`, but your work is greatly appreciated.
> 
> Since you added tests (which is great) alongside some binary files, it would be important to also include the code that was used to generate them if possible. See #16371 for context.

@neutrinoceros @pllim 
I can do that, and I agree it seems like a good idea.  But I don't know where to put the script.  I tried putting it in the `data/` directory, but the pre-commit check complains that it needs an `__init__.py` file.  I'm pretty sure we don't want that.  I could put it in `tests/`, but then pytest might try to run it?  I poked around in the astropy directories and did not see any obvious place to use for this sort of thing.

I have the code ready to check in but need some guidance before I can complete the process.
rlwastro: @pllim @neutrinoceros 

Maybe I should try to be clearer about my question: if I have a script `build_my_test_data.py` that creates the files used for testing, should I just stick it in the `tests/` directory?

I do not want it to run as part of testing, but I assume it should be somewhere in the vicinity of the `test` directory since that is what it is for.
dhomeier: Thanks for your work on this! Since saving the code to generate binary data has only been declared a desired policy in #16371, but no decision been made how to put that policy into effect, probably needs further discussion where to put those scripts.
However looking at the test from the [previous attempt to fix this this](https://github.com/astropy/astropy/pull/15756/files#diff-3d706f4a6b686b1787aeb0b2332764235716d20699794cb78f2b0709bd059582) I am wondering if _some_ of the test files could not be created directly in and by the test? As I understand, some of your examples are technically invalid FITS, so probably will not be possible to create (easily) with `io.fits`, but the remainder should.
Still trying to bring this into context with the definition of `BLANK` in the standard; probably will still take someone more familiar with its usage and especially the `ZBLANK` and table columns cases.
Last but not least note that much of this code is scheduled for a significant refactor in #15474, so it might better be hold until that PR is merged.
neutrinoceros: Agreed with @dhomeier that generating data as part of the test seems preferable if it's quick enough.
rlwastro: @dhomeier @neutrinoceros 
The files are generated by running fpack (plus some modifications to the fpack file using astropy.io.fits).  My opinion is that it is much better to rely on a file generated externally for the test than on generating a file on-the-fly.  Changes in the code could result in an error during compression that is compensated by an error in uncompression.  It seems dangerous to assume that cannot happen in this kind of code.

I could modify the tests so that there is only one additional test file.  I would be very reluctant to try to get to zero test files though.

I am very familiar with the FITS compression standard since [I wrote it](https://ui.adsabs.harvard.edu/abs/2012arXiv1201.1336W/abstract)!  So I am confident that my interpretation of the standard and the approach I've taken is correct.  Of course that doesn't mean it is not a good idea to find someone else to review it.

If you think the fix for #15474 will happen soon, I can see the benefit of waiting for this one.  I can't tell from the discussion there whether that major change is really imminent or not.  It is not obvious whether the discussion in that pull request is converging.  It's been 3 months since the request got switched back to a draft.
neutrinoceros: > If you think the fix for https://github.com/astropy/astropy/pull/15474 will happen soon

@astrofrog and @saimn would be in a better position to answer that question
rlwastro: After further thought, I decided that it would be OK to replace the test files with files that are completely generated within the `test_tiles_compression.py` script.  That makes the changes to the test script much longer, but now the additional FITS files are not part of the pull request.

Incidentally, in response to the comment by @dhomeier : none of the compressed files that are tested are actually invalid, as far as I can see.  The FITS spec is reluctant to say _must_ instead of _should_ for any of its requirements, so these all fit in the category of formats that are not the preferred approach but nonetheless are acceptable.  I think that makes them worth testing. It is quite possible that `astropy.io.fits` is more robust to these variations than `cfitsio` and `fpack/funpack` are.
neutrinoceros: @rlwastro we'll need to rebase this branch as some point to fix docs builds. Are you comfortable with the process ?
rlwastro: > @rlwastro we'll need to rebase this branch as some point to fix docs builds. Are you comfortable with the process ?

@neutrinoceros I think I can probably figure it out.  I struggled last time but pllim helped.
dhomeier: > Incidentally, in response to the comment by @dhomeier : none of the compressed files that are tested are actually invalid, as far as I can see. The FITS spec is reluctant to say _must_ instead of _should_ for any of its requirements, so these all fit in the category of formats that are not the preferred approach but nonetheless are acceptable. I think that makes them worth testing. It is quite possible that `astropy.io.fits` is more robust to these variations than `cfitsio` and `fpack/funpack` are.

If it can create them, all the better, as it is by default strict on writing – but something that just _should_ not be in there should indeed at most raise a warning. If you feel that something needs to be tested independently of the `io.fits` writer, it is definitely acceptable to add new files in `tests/data`, but I don't see any obvious case here, so the default should be to generate those in the test run (there are a handful of files stored there, but some 100 created by the test suite).
dhomeier: > I could try to replicate the weird discrepancy between the `uint` and `int` behavior, but that seems like a bad idea because (1) it appears to be wrong and (2) there are big changes coming anyway from changing the `CompImageHDU` to inherit from `ImageHDU` instead of `BinTableHDU`. It does not seem worthwhile to try to do any behavior cleanup of this now.

I tend to agree with (1) now, and quite possibly this combination of `uint` and `BLANK` simply has not been used by anyone since I have seen no questions about the weird behaviour. And the doc clearly states that for `ignore_blank=False` blank values _will_ be converted to NaN, so not doing so in this case can be considered a bug as well.
It seems that #15474 is relegating most or all of this scaling to the functions in [image.py](https://github.com/astropy/astropy/pull/15474/files#diff-1daacadfe9276932278eb86f8b21a7642fb2560cf5408d3e64e2f5ee5dff48c1), this should automatically ensure more consistent behaviour between uncompressed and compressed data. But this will only come in with a major release due to the API changes, so I agree now to bring this bug fix in as soon as possible, and then probably fix the corresponding `image.py` code so everything will still be consistent afterwards.
pllim: Hi. Just catching up on notifications and saw https://github.com/astropy/astropy/pull/16550#issuecomment-2174532227 but not sure if it has been resolved or not.
dhomeier: > Hi. Just catching up on notifications and saw [#16550 (comment)](https://github.com/astropy/astropy/pull/16550#issuecomment-2174532227) but not sure if it has been resolved or not.

There are no new test files added anymore, so no need for external scripts at this point. I want to take a closer look at the current test coverage, but after rebase this looks about ready to go for me.

pllim: Regular CI is all green. But if maintainers also want to run optional CI, please apply `Extra CI` label too. Otherwise, I think this is ready for final review as long as we agree to use the "Squash and Merge" button. Thanks!
rlwastro: > Regular CI is all green. But if maintainers also want to run optional CI, please apply `Extra CI` label too. Otherwise, I think this is ready for final review as long as we agree to use the "Squash and Merge" button. Thanks!

"Squash and Merge" is fine with me.  Thanks!
astrofrog: I would like to make a suggestion - this is clearly an important bug to fix, but if we wait until after #15474 then any bug fix merged after to the compressed FITS code won't be backportable anyway. Therefore one option would be to change this PR to be against the v6.1.x branch directly, and then once #15474 is merged we can cherry pick the regression test and re-do a new fix for the bug in v7.0.x?
rlwastro: @astrofrog @saimn @dhomeier 

> I would like to make a suggestion - this is clearly an important bug to fix, but if we wait until after #15474 then any bug fix merged after to the compressed FITS code won't be backportable anyway. Therefore one option would be to change this PR to be against the v6.1.x branch directly, and then once #15474 is merged we can cherry pick the regression test and re-do a new fix for the bug in v7.0.x?

I would be willing to give that a try if it seems like an acceptable solution.  (This is all pushing the limits of my comfort with git, but it's good to learn new things.)  I won't try it unless there is agreement that this would be a good approach.
astrofrog: @saimn - what do you think?
astrofrog: An alternative is that we do merge it into main but I would add a revert as my first commit in #15474 and re-implement a fix at the end of the git history in that PR as actually rebasing would be very tricky.
neutrinoceros: FWIW I was about to make the exact same suggestion Tom did in https://github.com/astropy/astropy/pull/16550#issuecomment-2189916395 
saimn: Having the bug fixed would be great so yes making the PR against the v6.1.x branch seems a good idea.
neutrinoceros: @rlwastro here are some instructions for this operation:
```
git remote add upstream https://github.com/astropy/astropy.git
git fetch upstream v6.1.x
git checkout rlwastro-compressionFix
git rebase -i v6.1.x
# ... rebasing ...
git push --force
```

Then on this page, click the "Edit" button next to the PR title: you should be able to change the target branch from main to v6.1.x (the dropdown should be next to "open")

Patch:
diff --git a/astropy/io/fits/hdu/compressed/_tiled_compression.py b/astropy/io/fits/hdu/compressed/_tiled_compression.py
index cc6909aafb..4901b936e8 100644
--- a/astropy/io/fits/hdu/compressed/_tiled_compression.py
+++ b/astropy/io/fits/hdu/compressed/_tiled_compression.py
@@ -359,7 +359,9 @@ def decompress_image_data_section(
     else:
         zzero_column = None
 
-    zblank_header = compressed_header.get("ZBLANK", None)
+    zblank_header = compressed_header.get(
+        "ZBLANK", compressed_header.get("BLANK", None)
+    )
 
     gzip_compressed_data_column = None
     gzip_compressed_data_dtype = None
@@ -465,7 +467,7 @@ def decompress_image_data_section(
             if zblank is not None:
                 if not tile_data.flags.writeable:
                     tile_data = tile_data.copy()
-                tile_data[blank_mask] = np.nan
+                tile_data[blank_mask] = zblank_header if zbitpix > 0 else np.nan
 
         image_data[tile_slices] = tile_data
 
diff --git a/astropy/io/fits/hdu/compressed/compressed.py b/astropy/io/fits/hdu/compressed/compressed.py
index 470e35408b..9df53d1442 100644
--- a/astropy/io/fits/hdu/compressed/compressed.py
+++ b/astropy/io/fits/hdu/compressed/compressed.py
@@ -3,6 +3,7 @@
 import ctypes
 import math
 import time
+import warnings
 from contextlib import suppress
 
 import numpy as np
@@ -21,6 +22,7 @@ from astropy.io.fits.util import (
 )
 from astropy.utils import lazyproperty
 from astropy.utils.decorators import deprecated_renamed_argument
+from astropy.utils.exceptions import AstropyUserWarning
 
 from .header import (
     CompImageHeader,
@@ -362,6 +364,28 @@ class CompImageHDU(BinTableHDU):
         self._orig_bscale = self._bscale
         self._orig_bitpix = self._bitpix
 
+        if (
+            self._bitpix > 0
+            and "BLANK" not in self._header
+            and "ZBLANK" not in self._header
+        ):
+            # check for column named "ZBLANK"
+            for i in range(1, self._header["TFIELDS"] + 1):
+                if self._header[f"TTYPE{i}"] == "ZBLANK":
+                    # required BLANK keyword is missing
+                    # use most negative value as default
+                    self._header["BLANK"] = -(1 << (self._bitpix - 1))
+                    warnings.warn(
+                        f"Setting default value {self._header['BLANK']} for missing BLANK keyword in compressed extension",
+                        AstropyUserWarning,
+                    )
+                    break
+
+        if self._bitpix > 0:
+            self._blank = self._header.get("BLANK", self._header.get("ZBLANK"))
+        else:
+            self._blank = None
+
     def _remove_unnecessary_default_extnames(self, header):
         """Remove default EXTNAME values if they are unnecessary.
 
@@ -537,15 +561,15 @@ class CompImageHDU(BinTableHDU):
             self.name = name
 
     def _scale_data(self, data):
-        if self._orig_bzero != 0 or self._orig_bscale != 1:
-            new_dtype = self._dtype_for_bitpix()
-            data = np.array(data, dtype=new_dtype)
-
-            if "BLANK" in self._header:
-                blanks = data == np.array(self._header["BLANK"], dtype="int32")
+        if self._orig_bzero != 0 or self._orig_bscale != 1 or self._blank is not None:
+            if self._blank is not None:
+                blanks = data == np.array(self._blank, dtype=data.dtype)
             else:
                 blanks = None
 
+            new_dtype = self._dtype_for_bitpix()
+            data = np.array(data, dtype=new_dtype)
+
             if self._orig_bscale != 1:
                 np.multiply(data, self._orig_bscale, data)
             if self._orig_bzero != 0:
@@ -556,7 +580,10 @@ class CompImageHDU(BinTableHDU):
                 np.add(data, self._orig_bzero, out=data, casting="unsafe")
 
             if blanks is not None:
-                data = np.where(blanks, np.nan, data)
+                # use float32 version of nan to reduce data size for uint conversion
+                # result will still be float64 for larger uint sizes (e.g., uint32)
+                # for float types np.where retains the type of the data array
+                data = np.where(blanks, np.float32(np.nan), data)
 
         return data
 
diff --git a/astropy/io/fits/hdu/compressed/header.py b/astropy/io/fits/hdu/compressed/header.py
index 62cf3792d7..cbc811d67e 100644
--- a/astropy/io/fits/hdu/compressed/header.py
+++ b/astropy/io/fits/hdu/compressed/header.py
@@ -467,6 +467,14 @@ def _bintable_header_to_image_header(bintable_header):
         del image_header["PCOUNT"]
         del image_header["GCOUNT"]
 
+    # Fill in BLANK keyword if necessary
+    if (
+        image_header["BITPIX"] > 0
+        and "BLANK" not in image_header
+        and "ZBLANK" in bintable_header
+    ):
+        image_header["BLANK"] = bintable_header["ZBLANK"]
+
     # Look to see if there are any blank cards in the table
     # header.  If there are, there should be the same number
     # of blank cards in the image header.  Add blank cards to
diff --git a/astropy/io/fits/hdu/hdulist.py b/astropy/io/fits/hdu/hdulist.py
index ba9c83c00f..92f53915ed 100644
--- a/astropy/io/fits/hdu/hdulist.py
+++ b/astropy/io/fits/hdu/hdulist.py
@@ -1507,7 +1507,7 @@ class HDUList(list, _Verify):
                     del hdu.data
                 hdu._file = ffo
 
-            if sys.platform.startswith("win"):
+            if sys.platform.startswith("win") and NUMPY_LT_2_0:
                 # On Windows, all the original data mmaps were closed above.
                 # However, it's possible that the user still has references to
                 # the old data which would no longer work (possibly even cause
@@ -1518,27 +1518,14 @@ class HDUList(list, _Verify):
                 # lead to odd behavior in practice.  Better to just not keep
                 # references to data from files that had to be resized upon
                 # flushing (on Windows--again, this is no problem on Linux).
+                # Note that this hack is only possible on numpy 1.x:
+                # in 2.x, we cannot write directly to the data attribute
                 for idx, mmap, arr in mmaps:
-                    if mmap is None:
-                        continue
-                    if NUMPY_LT_2_0:
-                        # Note that this hack is only possible on numpy 1.x:
-                        # in 2.x, we cannot write directly to the data attribute
+                    if mmap is not None:
                         # https://github.com/numpy/numpy/issues/8628
                         with warnings.catch_warnings():
                             warnings.simplefilter("ignore", category=DeprecationWarning)
                             arr.data = self[idx].data.data
-                    elif sys.getrefcount(arr) > 2:
-                        # 2 is the minimum number of references to this object
-                        # counting `arr` and a reference as an argument to getrefcount(),
-                        # see  https://docs.python.org/3/library/sys.html#sys.getrefcount
-                        warnings.warn(
-                            "Memory map object was closed but appears to still "
-                            "be referenced. Further access will result in undefined "
-                            "behavior (possibly including segmentation faults).",
-                            category=UserWarning,
-                            stacklevel=2,
-                        )
                 del mmaps  # Just to be sure
 
         else:
diff --git a/docs/changes/io.fits/16550.bugfix.rst b/docs/changes/io.fits/16550.bugfix.rst
new file mode 100644
index 0000000000..cff2187839
--- /dev/null
+++ b/docs/changes/io.fits/16550.bugfix.rst
@@ -0,0 +1,1 @@
+Fix a spurious exception when reading integer compressed images with blanks.
diff --git a/docs/changes/io.fits/16581.bugfix.rst b/docs/changes/io.fits/16581.bugfix.rst
index cbfa55604e..41edd54473 100644
--- a/docs/changes/io.fits/16581.bugfix.rst
+++ b/docs/changes/io.fits/16581.bugfix.rst
@@ -1,5 +1,5 @@
 Fixed a crash that occurred for files opened via
 ``fits.open(..., mode='update')``, on Windows, and with numpy 2.0 installed.
-A warning is now emitted in cases most likely to escalate into
-undefined behavior (e.g., segfaults), i.e., when a closed memory map object is
-still referenced by external code. Please report any regression found.
+It is possible, though unlikely, that this patch reveals memory issues
+for downstream code run under these exact conditions. Please report any
+regression found.


Test patch:
diff --git a/astropy/io/fits/hdu/compressed/tests/test_tiled_compression.py b/astropy/io/fits/hdu/compressed/tests/test_tiled_compression.py
index 9d79ec3a35..4985ca1df0 100644
--- a/astropy/io/fits/hdu/compressed/tests/test_tiled_compression.py
+++ b/astropy/io/fits/hdu/compressed/tests/test_tiled_compression.py
@@ -8,6 +8,7 @@ from numpy.testing import assert_allclose, assert_equal
 from astropy.io import fits
 from astropy.io.fits.hdu.compressed._codecs import PLIO1
 from astropy.io.fits.hdu.compressed._compression import CfitsioException
+from astropy.utils.exceptions import AstropyUserWarning
 
 from .conftest import fitsio_param_to_astropy_param
 
@@ -65,6 +66,97 @@ def test_zblank_support(canonical_data_base_path, tmp_path):
         assert_equal(np.round(hdul[1].data), reference)
 
 
+def test_integer_blank_support(canonical_data_base_path, tmp_path):
+    # This uses a test 64x64 section of the m13 image to which three NaN values are added.
+    # It is converted to an integer image with blank values and compressed.
+    # Then 3 variations of the BLANK/ZBLANK specification are created by
+    # editing the header and compressed data table.
+    # When these are read in, the result is a float image that should be the same
+    # as the original (including NaNs for blanks).
+
+    with fits.open(canonical_data_base_path / "m13.fits") as hdul:
+        # pick a 64x64 pixel subset
+        data = hdul[0].data[116:180, 116:180]
+        header = hdul[0].header.copy()
+        del header["CHECKSUM"]
+        del header["DATASUM"]
+
+    # float version of image with NaNs as blanks
+    reference = data.astype(np.float32)
+    reference[[1, 2, 3], [1, 2, 3]] = np.nan
+
+    # set blanks in the int16 image and BLANK keyword in header
+    # choose a blank value that differs from the default
+    blank = -16384
+    data[np.isnan(reference)] = blank
+    header["BLANK"] = blank
+
+    # create the compressed file
+    cfile1 = tmp_path / "compressed_with_blank.fits"
+    hdu = fits.CompImageHDU(
+        data=data, header=header, compression_type="RICE_1", tile_shape=(16, 16)
+    )
+    hdu.writeto(cfile1, overwrite=True, checksum=False)
+
+    # replace BLANK keyword in header with ZBLANK keyword
+    cfile2 = tmp_path / "compressed_with_zblank.fits"
+    with fits.open(cfile1, disable_image_compression=True) as hdul:
+        assert "ZBLANK" not in hdul[1].header
+        hdul[1].header["ZBLANK"] = hdul[1].header["BLANK"]
+        del hdul[1].header["BLANK"]
+        hdul.writeto(cfile2, overwrite=True, checksum=False)
+
+    # replace ZBLANK in header with with ZBLANK in table column
+    # This creates a file structure that is unlikely to be encountered in the wild
+    # but that is apparently allowed by the FITS standard.
+    # Two versions are created, one without the BLANK keyword and one with it.
+    cfile3 = tmp_path / "compressed_with_zblank_column.fits"
+    cfile4 = tmp_path / "compressed_with_zblank_column_and_blank.fits"
+    with fits.open(cfile2, disable_image_compression=True) as hdul:
+        phdu = hdul[0]
+        thdu = hdul[1]
+        orig_table = hdul[1].data
+        orig_cols = orig_table.columns
+        zblank = thdu.header["ZBLANK"]
+        new_cols = fits.ColDefs(
+            [
+                fits.Column(
+                    name="COMPRESSED_DATA",
+                    format="1PB()",
+                    array=thdu.data.field("COMPRESSED_DATA"),
+                ),
+                fits.Column(
+                    name="ZBLANK",
+                    format="I",
+                    array=np.zeros(len(orig_table), dtype=np.int32) + zblank,
+                ),
+            ]
+        )
+        new_thdu = fits.BinTableHDU.from_columns(new_cols, header=thdu.header)
+        del new_thdu.header["ZBLANK"]
+        new_hdul = fits.HDUList([phdu, new_thdu])
+        new_hdul.writeto(cfile3, overwrite=True, checksum=False)
+        new_thdu.header["BLANK"] = blank
+        new_hdul = fits.HDUList([phdu, new_thdu])
+        new_hdul.writeto(cfile4, overwrite=True, checksum=False)
+
+    # now test the 4 files to confirm they all uncompress correctly
+    for filename in (cfile1, cfile2, cfile4):
+        with fits.open(canonical_data_base_path / filename) as hdul:
+            assert_equal(hdul[1].data, reference)
+            # ensure that the uncompressed header is created with BLANK keyword
+            assert hdul[1].header.get("BLANK") == blank
+
+    # this one generates an expected warning
+    with pytest.warns(
+        AstropyUserWarning,
+        match="Setting default value -32768 for missing BLANK keyword in compressed extension",
+    ):
+        with fits.open(cfile3) as hdul:
+            assert_equal(hdul[1].data, reference)
+            assert hdul[1].header.get("BLANK") == -32768
+
+
 @pytest.mark.parametrize(
     ("shape", "tile_shape"),
     (