leejjoon / pywcsgrid2

astromonical fits image for python matplotlib
http://leejjoon.github.com/pywcsgrid2/
MIT License
25 stars 13 forks source link

How to display an image in equatorial frame (north pointing up)? #14

Open will-henney opened 11 years ago

will-henney commented 11 years ago

Hi,

I have managed to successfully use pywcsgrid2 for displaying a FITS image in the pixel reference frame (see below), but I would like to display it in the standard equatorial frame, where the negative x-axis is RA and the positive y-axis is Dec.

Is it possible to do this easily with pywcsgrid2? From looking at the docs, the only way I could see would be to make a fake FITS header that has the axes rotated, and then use that in a similar way to what is described here. But my first attempt at this failed, so I thought I would first ask if there might be an easier way.

Cheers

Will

pywcsgrid2-test

will-henney commented 11 years ago

Just to say that I put in a bit more effort and did get it to work using parasitic axes. I give an example below. A few details:

  1. It is necessary to explicitly set the plot window with set_xlim, set_ylim.
  2. It doesn't work with saving to a JPG file, but PDF and PNG are fine.

I'd still be interested in any way this could be made easier.

import numpy as np
import astropy.wcs as wcs
from astropy.io import fits
import pywcsgrid2
import matplotlib.pyplot as plt

hdu = fits.open("j8oc01010_drz/LL1-extract.fits")["SCI"]

# Data to plot is in rotated frame (2)
h2 = hdu.header
w2 = wcs.WCS(h2)

# To make the rotation easier, first shift the WCS reference pixel to
# the center of the grid
nx, ny = h2["naxis1"], h2["naxis2"]
i0, j0 = (float(nx) + 1.0)/2, (float(ny) + 1.0)/2
[ra0, dec0], = w2.wcs_pix2world([[i0, j0]], 1)
h2.update(crpix1=i0, crpix2=j0, crval1=ra0, crval2=dec0)

# First plot the image in the original (rotated) frame
ax2 = pywcsgrid2.axes(wcs=h2)
arcsec = 1./3600
ax2.imshow(hdu.data, origin='lower', 
           vmin=6, vmax=16, cmap=plt.cm.gray_r)
ax2["fk5"].plot([ra0, ra0], 
                [dec0 + 10*arcsec, dec0 + 15*arcsec], "r")
ax2.grid()
ax2.figure.savefig("pywcsgrid2-rotate-test-pixframe.pdf", dpi=150)

# Now construct a header that is in the equatorial frame (1)
h1 = h2.copy()
h1.update(
    cd1_1=-np.hypot(h2["CD1_1"], h2["CD1_2"]), 
    cd1_2=0.0, 
    cd2_1=0.0, 
    cd2_2=np.hypot(h2["CD2_1"], h2["CD2_2"]), 
    orientat=0.0)
# Finally plot the image in the new frame
ax = pywcsgrid2.axes(wcs=h1, aspect=1)
ax.set_xlim(-nx/4, 5*nx/4)
ax.set_ylim(-ny/4, 5*ny/4)
ax["fk5"].plot([ra0, ra0], 
               [dec0 + 10*arcsec, dec0 + 15*arcsec], "r")
ax[h2].imshow_affine(hdu.data, origin='lower', 
                     vmin=6, vmax=16, cmap=plt.cm.gray_r)
ax.grid()
ax.figure.savefig("pywcsgrid2-rotate-test-eqframe.pdf", dpi=150)

pywcsgrid2-rotate-test-eqframe

leejjoon commented 11 years ago

Good to hear that you figured out how to do this. And the way you do is in fact what I would recommend.

Issue 1 - Yes, inshow_affine commands which acts on parasite axes (e.g., ax[h2]) does NOT attempt to adjust the axes limits of its host axes (e.g. ax). This is just to make things simple. I just opened a new issue for this as a wishlist (#15).

Issue 2 - Not sure why that happens. JPG should use the same backend as the PNG one. Can you describe more about the problem? Saving JPG works for other script that do not use imshow_affine?