The FITS standard states:
3.3.2. Primary data array ... The individual data values shall be stored in big-endian byte order such that the byte containing the most-significant bits of the value appears first in the FITS file, followed by the remaining bytes, if any, in decreasing order of significance.
...
5. Data representation ... Bytes are in big-endian order of decreasing significance.
7.3.3.1. Main data table ...
16-Bit integer... The most-significant byte shall be first (big-endian byte order). ...
Single-precision floating point ... in big-endian byte order ...
Problem: Pillow does not seem to handle this big endian byte order. That's maybe OK for 8bit images. But #7894 already added a 16bit test image Tests/images/m13.fits, and loading and viewing that file demonstrates the problem.
The problem is even worse for floating point images. Example files can be downloaded here.
A workaround is using numpy's byteswap() function after loading the image.
With pillow (wrong byte order):
import PIL.Image
from matplotlib import pyplot as plt
filename = "502nmos.fits"
arr = numpy.asarray(PIL.Image.open(filename))
plt.imshow(arr, cmap='gray')
plt.show()
With byteswap (corrected byte order):
import PIL.Image
from matplotlib import pyplot as plt
filename = "502nmos.fits"
arr = numpy.asarray(PIL.Image.open(filename))
arr = arr.byteswap() # workaround to correct Big Endian byte order
plt.imshow(arr, cmap='gray')
plt.show()
With byteswap (corrected byte order) and normalized (enhanced contrast):
import PIL.Image
from matplotlib import pyplot as plt
import numpy as np
filename = "502nmos.fits"
arr = numpy.asarray(PIL.Image.open(filename))
arr = arr.byteswap() # workaround to correct Big Endian byte order
arr = np.clip(arr, np.percentile(arr, 0.4), np.percentile(arr, 99.6)) # normalize
plt.imshow(arr, cmap="gray")
plt.show()

The FITS standard states:
Problem: Pillow does not seem to handle this big endian byte order. That's maybe OK for 8bit images. But #7894 already added a 16bit test image Tests/images/m13.fits, and loading and viewing that file demonstrates the problem.
The problem is even worse for floating point images. Example files can be downloaded here.
A workaround is using numpy's
byteswap()function after loading the image.With pillow (wrong byte order):
With byteswap (corrected byte order):
With byteswap (corrected byte order) and normalized (enhanced contrast):