In [1]:
import numpy as np
import matplotlib
import astropy.io.fits as afits
from astropy.wcs import WCS
import reproject
from astropy.visualization import ZScaleInterval
import astropy.table as at
import astropy.coordinates as coords
import astropy.units as u
from astropy.visualization.wcsaxes import WCSAxes
import astropy.visualization.wcsaxes.frame as frame
%matplotlib notebook
%pylab
In [2]:
!ls ../01_Registration/out/
In [3]:
# load the two fits images
f1 = afits.open('../01_Registration/out/wdd7.040920_0452.051_6.solved.fits')
f2 = afits.open('../01_Registration/out/wdd7.080104_0214.1025_6.solved.fits')
f1d = f1[0].data
f2d = f2[0].data
# load the two WCS solutions we computed
f1wcs = WCS(afits.Header.fromfile('../01_Registration/out/wdd7.040920_0452.051_6.wcs'))
f2wcs = WCS(afits.Header.fromfile('../01_Registration/out/wdd7.080104_0214.1025_6.wcs'))
In [4]:
f2reproj, f2footprint = reproject.reproject_exact((f2d, f2wcs), f1wcs, shape_out=f1d.shape, parallel=True)
In [5]:
zscaler = ZScaleInterval(nsamples=1000, contrast=0.25)
f1s = zscaler(f1d)
f2s = zscaler(f2reproj)
fig = figure(figsize=(10,4))
ax = fig.add_subplot(1,1,1)
ax.imshow(f1s.T, cmap='Reds')
ax.imshow(f2s.T, cmap='Blues', alpha=0.5)
tight_layout()
xlabel('x')
ylabel('y')
savefig('out/reprojected.pdf')
In [6]:
mask = (f2footprint > 0)
f1m = f1d.mean()
f2m = f2reproj[mask].mean()
In [7]:
diff = f1d/f1m - f2reproj/f2m
diffs = zscaler(diff)
In [8]:
fig2 = figure(figsize=(10,4))
ax2 = fig2.add_subplot(1,1,1, projection=f1wcs)
ax2.imshow(diffs.T, cmap='Greys', origin='lower')
# add a grey grid for what whatever we're using as the native projection
ax2.coords.grid(True, color='green', ls='solid')
ra = ax2.coords['ra']
dec = ax2.coords['dec']
ra.set_axislabel(r'$\alpha$ (J2000)')
dec.set_axislabel(r'$\delta$ (J2000)')
tight_layout()
savefig('out/kinda_difference_image.pdf')