In [1]:
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import LinearSegmentedColormap
from astropy.io import fits
import scipy, scipy.ndimage
import os
import gpipy
import gpipy.plotutils as plotutil
colormap_ds9_b = plotutil.colormap_ds9_b
In [2]:
datapath = '/Users/mperrin/Dropbox (Personal)/Documents/my_papers/gpi/pol_and_4796/data'
# or get the files from Dropbox (GPI)/SCRATCH/Scratch/mperrin/HR4796A_paper_datafiles
In [ ]:
im_model = gpipy.read(os.path.join(datapath, 'hr4796_k1_vc2_psfsubtract_stokesdc.fits'))
im_klip = gpipy.read(os.path.join(datapath,'S20140325S0218_psf__klip99_northup.fits'))
In [3]:
plot_kwargs = {'aspect': 0.7,
'crop':180,
'what':'i',
'coords':'arcsec',
'colorbar': True,
'cmap': plotutil.colormap_ds9_b}
colorbar_kwargs = {'shrink': 0.8, 'pad':0.07, 'aspect':15}
plt.figure(figsize=(10,12))
plt.subplots_adjust(wspace=0.05)
plt.subplot(121)
ax1,cb1 = im_model.imshow(vmin=1, colorbar_kwargs=colorbar_kwargs, **plot_kwargs)
plt.subplot(122)
ax2,cb2 = im_klip.imshow(vmin=1, colorbar_kwargs=colorbar_kwargs, **plot_kwargs)
for ax in [ax1,ax2]:
pass
ax2.yaxis.set_visible(False)
ax1.set_ylabel('Offset [Arcsec]')
for cb in [cb1,cb2]:
cb.set_label('Counts')
# cb.ax.xaxis.set_major_formatter(plotutil.NicerLogFormatter())
#cb = plt.colorbar(shrink=0.5)
In [4]:
im_model = gpipy.read(os.path.join(datapath, 'hr4796_k1_vc2_psfsubtract_stokesdc.fits'))
im_klip = gpipy.read(os.path.join(datapath, 'S20140325S0218_psf__klip99_northup.fits'))
plot_kwargs2 = {'aspect': 0.8,
'crop':180,
'what':'i',
'coords':'arcsec',
'colorbar': True,
'scale': 'linear',
'smooth': 0.75,
'counts_per_sec': True,
'cmap': plotutil.colormap_ds9_b}
colorbar_kwargs = {'shrink': 0.9, 'pad':0.07, 'aspect':15}
plt.figure(figsize=(10,9))
plt.subplots_adjust(wspace=0.07, hspace=0.01, top=0.99, bottom=0.01, left=0.07, right=0.98)
plt.subplot(121)
ax1,cb1 = im_klip.imshow(vmin=-8./60, vmax=0.22*np.nanmax(im_klip.i)/60, colorbar_kwargs=colorbar_kwargs, **plot_kwargs2)
plotutil.tvcircle(25*0.0143, color='cyan', ls='--', lw=2)
plt.title('ADI + KLIP PSF subtraction', fontsize=18)
plt.subplot(122)
ax2,cb2 = im_model.imshow(vmin=-30./60, vmax = 0.04*np.nanmax(im_model.i)/60, colorbar_kwargs=colorbar_kwargs, **plot_kwargs2)
plotutil.tvcircle(50*0.0143, color='cyan', ls='--', lw=2)
plt.title('Interpolated Model PSF sub.', fontsize=18)
for ax,xc in zip([ax1,ax2],[0,0.0143]):
# offset to xc to fix an apparent off-by-one-pixel issue?
plotutil.compass((-0.65, -1.2), ax=ax, length=0.15, textscale=2, labeleast=True)
#ax.add_patch(matplotlib.patches.Circle((0,0), 0.153, alpha=1.0, facecolor='none', edgecolor='black', lw=1.5))
ax.add_patch(matplotlib.patches.Circle((xc,0), 0.153, alpha=0.5, facecolor='black',edgecolor='white', lw=1, ls='solid'))
ax.add_patch(matplotlib.patches.Circle((xc,0), 0.153, alpha=1.0, facecolor='none', edgecolor='white', lw=1, ls='dotted'))
ax2.yaxis.set_visible(False)
ax1.set_ylabel('Offset [Arcsec]')
for cb in [cb1,cb2]:
cb.locator = matplotlib.ticker.MaxNLocator(6)
cb.update_ticks()
cb.set_label('Counts second$^{-1}$ spaxel$^{-1}$ ')
# cb.ax.xaxis.set_major_formatter(plotutil.NicerLogFormatter())
#cb = plt.colorbar(shrink=0.5)
#plt.tight_layout()
plt.savefig("fig_total_intensity_subtractions_v3.pdf", transparent=True)
In [5]:
im_model = gpipy.read(os.path.join(datapath,'hr4796_k1_vc2_psfsubtract_stokesdc.fits'))
# slice 0: PSF subtracted total intensity
# slice 1: Stokes Q
# slice 2: Stokes U
# slice 3: original total intensity (no PSF sub)
plot_kwargs3 = {'aspect': 1.0,
'crop':170,
'coords':'arcsec',
'scale': 'linear',
'smooth': 0.75,
'cmap': plotutil.colormap_ds9_b}
plt.figure(figsize=(12,3))
plt.subplot(131)
im_model.imshow(what=3, vmin=100, vmax=5000, **plot_kwargs3)
plt.subplot(132)
exp=0.5
plt.imshow(im_model.data[3]**exp,vmin=100**exp, vmax=5000**exp,cmap=plotutil.colormap_ds9_b)
plt.subplot(133)
sqnorm = plotutil.SqrtNorm(vmin=100, vmax=5000)
plt.imshow(im_model.data[3], norm=sqnorm,cmap=matplotlib.cm.cubehelix)
Out[5]:
In [7]:
# Set up custom colormap for right hand panels. Tweaked for visual appeal
modified_terrain_data = (
(0.00, (0.2, 0.2, 0.6)),
(0.15, (0.0, 0.6, 1.0)),
(0.25, (0.0, 0.8, 0.4)),
(0.50, (1.0, 1.0, 0.6)),
(0.75, (1.0, 0.36, 0.33)),
(1.00, (1.0, 1.0, 1.0)))
modified_terrain = matplotlib.colors.LinearSegmentedColormap.from_list('modified_terrain', modified_terrain_data)
In [8]:
files = ['hr4796_h_vc1_psfsubtract.fits','hr4796_k1_vc1_psfsubtract.fits','hr4796_k1_vc2_psfsubtract_stokesdc.fits']
itimes = [60, 60, 60]
plot_kwargs4 = {'aspect': 1.0,
'crop':170,
'coords':'arcsec',
'scale': 'linear',
'smooth': 0.75}
# 'cmap': plotutil.colormap_ds9_b}
from copy import deepcopy
colormap_for_p = deepcopy(plotutil.colormap_ds9_b)
colormap_for_p.set_bad('white')
plt.figure(figsize=(10,8))
plt.subplots_adjust(hspace=0, wspace=0)
ncol=4
for ifile, filename in enumerate(files):
cube = gpipy.read(os.path.join(datapath, files[2])) # only the VC2 K1 as modified to look lik a regular GPI cube so gpidata can read it.
mydata = fits.getdata(os.path.join(datapath,files[ifile]))
cube.data = mydata*1.0/itimes[ifile]
ncol=4
# Total intensity, no PSB sub in left hand column
ax = plt.subplot(len(files),ncol,ncol*ifile+1)
vmax =70/1.5 if ifile > 0 else 100 #150
vmin=10 if ifile > 0 else 10
vmin = [10,12,10][ifile]
cube.imshow(what=3, vmin=vmin, vmax=vmax, cmap=plotutil.colormap_ds9_b, **plot_kwargs4)
# Total intensity, interpolated model PSB sub in 2nd left column
ax = plt.subplot(len(files),ncol,ncol*ifile+2)
cube.imshow(what=0, vmin=0, vmax=6, cmap=plotutil.colormap_ds9_b, **plot_kwargs4)
# Polarized intensity, in 3rd column
ax = plt.subplot(len(files),ncol,ncol*ifile+3)
cube.imshow(what='p', vmin=0, vmax=3.5, cmap=plotutil.colormap_ds9_b, **plot_kwargs4)
# Total intensity, interpolated model PSB sub in 2nd left column
ax = plt.subplot(len(files),ncol,ncol*ifile+4)
fraction = cube.p/cube.i
smooth=2
smoothed_q = scipy.ndimage.filters.gaussian_filter(cube[1], smooth )
smoothed_u = scipy.ndimage.filters.gaussian_filter(cube[2], smooth )
smoothed_p = np.sqrt(smoothed_q**2+smoothed_u**2)
mask = np.ones_like(cube.i)
mask[smoothed_p < 0.3] = np.nan
mask[:,200:] = np.nan
mask[:,0:50] = np.nan
mask[0:67, :] = np.nan
mask[225:, :] = np.nan
mask[smoothed_p/cube.i < 0.1] = np.nan
mask2 = np.isfinite(mask)
struct = morph.iterate_structure(morph.generate_binary_structure(2, 1), 10).astype(int)
mask2 = (morph.binary_closing(mask2,struct ))
mask3 = np.ones_like(mask2) + np.nan
mask3[mask2] = 1
cm = modified_terrain #matplotlib.cm.terrain
cm.set_under('white')
#mask[cube.i < 0.1] = np.nan
#cube.data[1][mask] = np.nan
#plt.imshow(mask)
cube.imshow(what='fraction', vmin=0, vmax=1,cmap=cm, mask=mask3,**plot_kwargs4)
for ic in range(ncol):
for ir in range(3):
ax = plt.subplot(len(files),ncol, ncol*ir+ic+1)
ax.xaxis.set_visible(False)
ax.yaxis.set_visible(False)
xc=[-0.04, -0.03, 0][ir]
yc=[0.01, 0.02, 0][ir]
if ir == 0 and ic==3:
ax.add_patch(matplotlib.patches.Circle((xc,yc), 0.2, facecolor='white',edgecolor='white', lw=1, ls='solid'))
radius = 0.153 if ir > 0 else 0.123
ax.add_patch(matplotlib.patches.Circle((xc,yc), radius, alpha=0.5, facecolor='black',edgecolor='white', lw=1, ls='solid'))
ax.add_patch(matplotlib.patches.Circle((xc,yc), radius, alpha=1.0, facecolor='none', edgecolor='black', lw=1, ls='dotted'))
ax = plt.subplot(len(files),ncol, ncol*2+ic+1)
#cb = plt.colorbar(ax.images[0], ax=ax, orientation='horizontal', shrink=0.75 )
#cb.locator = matplotlib.ticker.MaxNLocator(4)
#cb.update_ticks()
plt.savefig('gallery_revised_v2.pdf', transparent=True)
In [64]:
# Read in data
mydata = fits.open(os.path.join(datapath,'hr4796_k1_vc2_psfsubtract.fits'))
mystokes = gpipy.read(os.path.join(datapath,'hr4796_k1_vc2_psfsubtract_stokesdc.fits'))
# Normalize to units of counts/second
mydata[0].data /= mystokes.itime
mystokes.data /= mystokes.itime
def plot_image(mydata, what='total_gaspard', vmin=300, vmax=5e3,ext=0, center=None,
verbose=False, colorbar=False, smooth=0.9, scale='log',
aspect=0.8, colorbar_kwargs=None, **kwargs):
"""
Parameters
----------
ext : int
Which extension to get the datacube extension from? This is needed because Gaspard's files are output nonstandardly.
what: string
Similar.
"""
#print ext
#print mydata[ext].header
cube = mydata[ext].data
if center is None:
cenx = mydata[ext].header['PSFCENTX']
ceny = mydata[ext].header['PSFCENTY']
else:
cenx = center[1], center[0]
box = 90
if scale=='log':
norm = matplotlib.colors.LogNorm(vmin=vmin,vmax=vmax)
else:
norm = matplotlib.colors.Normalize(vmin=vmin,vmax=vmax)
if verbose:
print ("vmin: {0}\tvmax: {1}".format(vmin,vmax))
if what=='total_gaspard':
what_to_show = cube[3]
elif what=='total':
what_to_show = cube[0]
elif what=='polint':
p = np.sqrt(cube[1]**2+cube[2]**2)
what_to_show = p
elif what=='polint_smoothed':
smoothed_q = scipy.ndimage.filters.gaussian_filter(cube[1], smooth )
smoothed_u = scipy.ndimage.filters.gaussian_filter(cube[2], smooth )
p = np.sqrt(smoothed_q**2+smoothed_u**2)
what_to_show = p
elif what=='polfrac':
pdivi = np.sqrt(cube[1]**2+cube[2]**2)/cube[0]
what_to_show = pdivi
pixelscale = 0.0143
shp = what_to_show.shape
extent= np.asarray([-cenx-0.5, shp[1]-cenx-0.5, -ceny-0.5, shp[0]-ceny-0.5]) * pixelscale
plt.imshow(what_to_show, cmap=colormap_ds9_b, norm=norm, extent=extent, **kwargs)
ax = plt.gca()
ax.set_xlim( -box*aspect*pixelscale, +box*aspect*pixelscale)
ax.set_ylim( -box*pixelscale, +box*pixelscale)
#ax.set_xlim( cenx-box*aspect, cenx+box*aspect)
#ax.set_ylim( ceny-box, ceny+box)
if colorbar:
colorbar_orientation='horizontal'
if colorbar_kwargs is None:
colorbar_kwargs = {}
cb = plt.colorbar(ax.images[0], ax=ax, orientation=colorbar_orientation, drawedges=False, **colorbar_kwargs)
if scale.lower() =='log':
ticks = np.logspace(np.log10(vmin), np.log10(vmax), np.log10(vmax/vmin)+1)
if colorbar_orientation=='horizontal' and vmax==1e-1 and vmin==1e-8: ticks = [1e-8, 1e-6, 1e-4, 1e-2, 1e-1] # looks better
cb.set_ticks(ticks)
cb.set_ticklabels(ticks)
#cb.solids.set_rasterized(True) # work around bug that makes faint vertical lines on colorbar;
# see http://stackoverflow.com/questions/15003353/why-does-my-colorbar-have-lines-in-it
return (ax,cb)
else:
return ax
In [10]:
plt.figure(figsize=(10,9))
plot_image(mydata, what='polint_smoothed', colorbar=True, vmin=0, vmax=2.5, smooth=0.7, alpha=0.5, scale='linear')
mystokes.polvect(color='black', showevery=4,
pmin=0.3333, scale=2.5, lw=1,
alpha=0.8,
scalebar_location = (0.07, 0.03), scalebar_value=1, scalebar_label="P = 1 count sec$^{-1}$ spaxel$^{-1}$",scalebar_fontsize=14)
plt.xlabel('Offset [Arcsec]')
plt.ylabel('Offset [Arcsec]')
plt.title('Polarization Vectors', fontsize=18)
plt.savefig('fig_pol_vectors_VC2_v3.pdf', transparent=True)
In [11]:
plt.figure(figsize=(10,9))
plot_image(mydata, what='polint_smoothed', colorbar=True, vmin=0, vmax=3, smooth=0.85, scale='linear')
plt.savefig('fig_polint_smoothed.pdf', transparent=True)
In [12]:
plt.imshow(mystokes.p, vmin=0.3, vmax=3)
Out[12]:
In [13]:
import copy
import scipy.ndimage as ndimage
def plot_one(num=217):
sky_fn = '/Users/mperrin/GPI/data/Reduced/140325/S20140325S0219_podc.fits'
sci_fn = '/Users/mperrin/GPI/data/Reduced/140325/S20140325S0{0}_podc.fits'.format(num)
sky = fits.open(sky_fn)
sci = fits.open(sci_fn)
#sub = copy.deepcopy(sci)
subtracted = sci[1].data- sky[1].data
log = matplotlib.colors.LogNorm(vmin=20,vmax=3e3)
plt.figure(figsize=(8,8))
plt.subplots_adjust(top=1.0, bottom=0.0, left=0.0, right=1.0)
subtracted[0][~np.isfinite(subtracted[0])] = np.nanmin(subtracted[0])
rotated = ndimage.rotate(subtracted[0], 24.5)
#plt.imshow(subtracted[0], cmap=colormap_ds9_b, norm=log)
plt.imshow(rotated, cmap=colormap_ds9_b, norm=log)
ax = plt.gca()
ax.set_xlim(100,285)
ax.set_ylim(95,280)
ax.xaxis.set_visible(False)
ax.yaxis.set_visible(False)
plt.savefig('single_hr4796a_im{0}.png'.format(num))
#plot_image(sub) #, what='polint',vmin=10,vmax=5e2)
#for i in range(200,206)+range(215,218):
for i in [200]:
plot_one(i)
#plot_one(217)
# to make a simple movie:
#convert -loop 5 -delay 30 single_hr4796a_im2*png raw_hr4796.gif
# PA from header read out using GETROT in IDL: -69.8
#ax = plt.gca()
#plotutil.compass((0,0), ax=ax, length=0.15, textscale=2, color='black')
In [65]:
plt.figure(figsize=(14,8))
plt.subplots_adjust(top=1.0, bottom=0.0, left=0.0, right=1.0)
plt.subplots_adjust(wspace=0)
mydata = fits.open(os.path.join(datapath,'hr4796_k1_vc2_psfsubtract.fits'))
ax2 = plt.subplot(121)
plot_image(mydata) #, vmin=10,vmax=5e4)
ax3 = plt.subplot(122)
plot_image(mydata, what='polint',vmin=10,vmax=5e2)
for ax in [ax2,ax3]:
ax.xaxis.set_visible(False)
ax.yaxis.set_visible(False)
plt.savefig('hr4796a_k1_vc2_display.png', transparent=True)
In [14]:
plt.figure(figsize=(10,9))
scale='linear'
if scale == 'log':
vmin=5./60
vmax=350./60
else:
vmin=0
vmax=2.5
plot_image(mydata, what='polint_smoothed', scale=scale, colorbar=True, vmin=vmin,vmax=vmax, smooth=0.85)
plt.xlabel('Offset [Arcsec]')
plt.ylabel('Offset [Arcsec]')
#plt.title('Polarization Vectors', fontsize=18)
plt.savefig('fig_polint_smoothed.pdf', transparent=True)
In [15]:
plt.figure(figsize=(10,9))
plt.subplots_adjust(wspace=0.07, hspace=0.01, top=0.99, bottom=0.01, left=0.07, right=0.98)
plt.subplot(121)
ax1,cb1 = im_klip.imshow(vmin=-8./60, vmax=0.22*np.nanmax(im_klip.i)/60,
colorbar_kwargs=colorbar_kwargs, **plot_kwargs2)
plotutil.tvcircle(25*0.0143, color='cyan', ls='--', lw=2)
plt.text(-0.95, 1.15, 'Total Intensity (ADI + KLIP)', fontsize=16, color='white')
plt.subplot(122)
ax2, cb2 = plot_image(mydata, what='polint_smoothed', colorbar=True, scale='linear',
vmin=0, vmax=2.5, smooth=0.85,
colorbar_kwargs=colorbar_kwargs)
plt.text(-0.95, 1.15, 'Polarized Intensity', fontsize=16, color='white')
#plt.xlabel('Offset [Arcsec]')
#plt.ylabel('Offset [Arcsec]')
ax2.yaxis.set_visible(False)
ax1.set_ylabel('Offset [Arcsec]')
for cb in [cb1,cb2]:
cb.locator = matplotlib.ticker.MaxNLocator(6)
cb.update_ticks()
cb1.set_label('Total Intensity [counts s$^{-1}$ pixel$^{-1}$]')
cb2.set_label('Polarized Intensity [counts s$^{-1}$ pixel$^{-1}$]')
for ax in [ax1, ax2]:
ax.add_patch(matplotlib.patches.Circle((0,0), 0.153, alpha=0.5, facecolor='black', edgecolor='white', linestyle='dashed', lw=2))
plt.savefig('fig_combined_I_P_smoothed.pdf', transparent=True)
In [15]:
In [ ]: