astropy: Visualization: Image displays in duplicate?

Description

In this image, the central grayscale image shows up twice, once where it should be and once wrapping at the corners:

Expected behavior

The original image looks like this:

so it should appear that way (rotated) on the Galactic frame

How to Reproduce

from astropy.io import fits
from astropy.utils.data import download_file
import matplotlib.pyplot as plt
import os
from astropy.wcs import WCS
from astropy.visualization import simple_norm

fn_meerkat = 'MeerKAT_Galactic_Centre_1284MHz-StokesI.fits'
if not os.path.exists(fn_meerkat):
    meerkat = download_file(f"https://archive-gw-1.kat.ac.za/public/repository/10.48479/fyst-hj47/data/{fn_meerkat}")
    shutil.move(meerkat, fn_meerkat)
meerkat = fits.open(fn_meerkat)

target_header_cmz = fits.Header.fromstring("""
NAXIS   =                    2
NAXIS1  =                 1920
NAXIS2  =                  960
CTYPE1  = 'GLON-MOL'
CRPIX1  =                960.5
CRVAL1  =                  0.0
CDELT1  =                -0.05
CUNIT1  = 'deg     '
CTYPE2  = 'GLAT-MOL'
CRPIX2  =                480.5
CRVAL2  =                  0.0
CDELT2  =                 0.05
CUNIT2  = 'deg     '
COORDSYS= 'icrs    '
""", sep='\n')

plt.figure(figsize=(10,5), dpi=200)
ax = plt.subplot(1,1,1, projection=WCS(target_header_cmz),)
                 #frame_class=EllipticalFrame)
ax.set_xlim(840, 1080)
ax.set_ylim(360, 600)

ax.imshow(meerkat[0].data.squeeze(),
          norm=simple_norm(meerkat[0].data, min_percent=0.1, max_percent=99.9, stretch='log'),
          transform=ax.get_transform(WCS(meerkat[0].header).celestial),
          cmap='gray')

Versions

Linux-3.10.0-1160.80.1.el7.x86_64-x86_64-with-glibc2.17
Python 3.9.13 | packaged by conda-forge | (main, May 27 2022, 16:56:21) 
[GCC 10.3.0]
astropy 5.1
Numpy 1.23.1
pyerfa 2.0.0.1
Scipy 1.9.0
Matplotlib 3.5.2

About this issue

Original URL
State: open
Created a year ago
Comments: 16 (16 by maintainers)

Most upvoted comments

No, don’t close this issue yet - if this usage isn’t supported, we need to clarify that in the documentation and/or add an error or warning message.

The specific behavior in this bug report is certainly incorrect and unexpected and points to something weird happening under the hood, too. Maybe the solution is to warn users away from using imshow, but some change is needed.

keflavich on Feb 21, 2023

Here’s a more minimal example than using the ~850 MB 100-megapixel data array above. The middle panel shows how imshow() makes up stuff beyond the bounds of the data due to the rotation in the transform. The right panel shows how pcolormesh() can be used instead.

import matplotlib.pyplot as plt
import numpy as np

import astropy.units as u
from astropy.visualization.wcsaxes import Quadrangle
from astropy.wcs import WCS

# Create a simple gradient array
data = np.broadcast_to(np.arange(10), (5, 10))

# Source WCS for data array
wcs = WCS({
    'naxis': 2,
    'naxis1': 10,
    'naxis2': 5,
    'ctype1': 'RA---TAN',
    'ctype2': 'DEC--TAN',
    'crval1': 0,
    'crval2': 0,
    'crpix1': 5.5,
    'crpix2': 3,
    'cdelt1': 1,
    'cdelt2': 1,
})

# Target WCS that is rotated by 30 deg
wcs_30 = WCS({
    'naxis': 2,
    'naxis1': 20,
    'naxis2': 20,
    'ctype1': 'RA---TAN',
    'ctype2': 'DEC--TAN',
    'crval1': 0,
    'crval2': 0,
    'crpix1': 10.5,
    'crpix2': 10.5,
    'cdelt1': 1,
    'cdelt2': 1,
    'crota2': 30,
})

fig = plt.figure(figsize=(12, 3))

ax1 = fig.add_subplot(131, projection=wcs)
ax1.imshow(data, origin='lower')
ax1.set_title('On non-rotated WCS')

ax2 = fig.add_subplot(132, projection=wcs_30)
ax2.imshow(data, origin='lower', transform=ax2.get_transform(wcs))
ax2.add_patch(Quadrangle(
    (-5*u.deg, -2.5*u.deg), 10*u.deg, 5*u.deg,
    edgecolor='red', facecolor='none',
    transform=ax2.get_transform('icrs')))
ax2.set_xlim(-0.5, 19.5)
ax2.set_ylim(-0.5, 19.5)
ax2.set_title('Using imshow')

ax3 = fig.add_subplot(133, projection=wcs_30)
# Need to provide pixel edges to pcolormesh()
ax3.pcolormesh(np.arange(data.shape[1] + 1) - 0.5,
               np.arange(data.shape[0] + 1) - 0.5,
               data,
               transform=ax3.get_transform(wcs))
ax3.grid(False)
ax3.add_patch(Quadrangle(
    (-5*u.deg, -2.5*u.deg), 10*u.deg, 5*u.deg, 
    edgecolor='red', facecolor='none',
    transform=ax3.get_transform('icrs')))
ax3.set_xlim(-0.5, 19.5)
ax3.set_ylim(-0.5, 19.5)
ax3.set_aspect('equal')
ax3.set_title('Using pcolormesh')

plt.show()

example

ayshih on Feb 24, 2023