In [1]:
import os
import numpy as np
from astropy import units as u
from astropy.coordinates import SkyCoord
from astropy import coordinates
from astropy.io import fits
import matplotlib.pyplot as plt
import aplpy
%matplotlib inline
In [2]:
imagename = "./images/uid___A002_X969646_X1aa8.ms.split.cal-CALIBRATE_BANDPASS-J0501-0159.ms.self3.substracted.cont.image.fits"
In [3]:
hdu_list = fits.open(imagename)
In [4]:
hdu_list.info()
In [5]:
image_data = hdu_list[0].data
print(type(image_data))
print(image_data.shape)
img = image_data[0][0]
In [6]:
image_data.max(), image_data.min()
Out[6]:
In [7]:
fig = aplpy.FITSFigure(imagename)
fig.add_beam()
fig.beam.set_color('red')
fig.add_scalebar(10 * u.arcsecond)
fig.scalebar.set_label('10 arcsec')
fig.scalebar.set_color("white")
fig.show_colorscale(vmin=0, vmax=0.0025)
In [8]:
from wavelet import wt
wt = wt()
In [9]:
result = wt.atrous(img, 6)
In [10]:
plt.figure(figsize=(16,6))
plt.subplot(241)
plt.imshow(img)
plt.xticks([]), plt.yticks([])
for i,lvl in enumerate(result):
plt.subplot(24*10+i+2)
plt.imshow(lvl)
plt.xticks([]), plt.yticks([])
plt.show()
In [11]:
filteredplane = wt.filtering(result, threshold=5)#, mask=(192, 320, 192, 320))
#filteredplane = wt.filtering(result, waveletNoise=True, imageNoise = 8.4E-05) #waveletNoise=True,
In [12]:
plt.figure(figsize=(16,6))
plt.subplot(241)
plt.imshow(img)
plt.xticks([]), plt.yticks([])
for i,lvl in enumerate(filteredplane):
plt.subplot(24*10+i+2)
plt.imshow(lvl)
plt.xticks([]), plt.yticks([])
plt.show()
In [13]:
filtered = wt.restore(filteredplane, 0, 6)
In [14]:
plt.figure(figsize=(16,12))
plt.subplot(121)
plt.imshow(img)
plt.title("original image")
plt.xticks([]), plt.yticks([])
plt.subplot(122)
plt.title("filtered")
plt.imshow(filtered)
plt.xticks([]), plt.yticks([])
plt.show()
In [15]:
hdu_list[0].data[0][0] = filtered
In [16]:
outputname = 'filtered.fits'
if os.path.exists(outputname):
os.remove(outputname)
hdu_list.writeto(outputname)
In [17]:
rms = 8.4e-05
c = SkyCoord('05h01m12.8s -01d59m14s', unit=(u.hourangle, u.deg), frame='icrs')
center_x, center_y = [c.ra.value, c.dec.value]
PB = 58./3600.0
multp = np.array([4, 10, 50, 100, 500, 1000])
lvl = rms*multp
image = aplpy.FITSFigure(outputname)
#image.tick_labels.set_font(size='small')
#image.show_contour(colors='white', levels=lvl, alpha=0.75)
image.hide_xaxis_label()
image.hide_yaxis_label()
image.hide_tick_labels()
image.add_beam()
image.beam.set_color('gray')
image.add_scalebar(1 * u.arcsecond)
image.scalebar.set_label('1 arcsec')
image.scalebar.set_color("white")
# marker_size = 100
# image.show_markers(ra, dec, edgecolor='red', facecolor='none', marker='o', s=marker_size, alpha=0.4)
# primary beam circle
image.show_circles(center_x, center_y, PB/2.0, edgecolor='gray')
image.show_colorscale(vmin=0, vmax=0.0025)
In [18]:
image.savefig("filtered.png", dpi=300)
In [19]:
hdu_list = fits.open("filtered.fits")
image_data = hdu_list[0].data
image_data.max()
Out[19]:
In [20]:
image_data.max()
Out[20]:
In [21]:
image_data.min()
Out[21]:
In [ ]: