This is the Python module that has the FITS reader (and writer).
In [23]:
import pyfits
# some configuration for plotting the way I like it. Should really put this
# into my general config file.
%pylab qt
rcParams['image.interpolation'] = 'nearest'
rcParams['image.cmap'] = 'gray'
rcParams['image.aspect'] = 'auto'
going to scan5 for now.
In [17]:
fitsfiles = !ls /Users/maye/data/irtf/guidedog/scan4/fname*.fits
In [18]:
def read_scan1(fname):
data = pyfits.getdata(fname, 0)
return data[248:394, :]
In [19]:
pyfits.getheader(fitsfiles[0])['TIME_OBS']
Out[19]:
In [40]:
data = pyfits.getdata(fitsfiles[10])
In [41]:
from scipy.ndimage import median_filter
filtered = median_filter(data, 2)
In [42]:
imshow(filtered, vmax=2000)
Out[42]:
In [26]:
hist(data.flatten(), 50, log=True)
Out[26]:
In [4]:
cube = []
for fitsfile in fitsfiles:
# fish out the image number from the filename
nr = int(fitsfile.split('-')[1][:5])
cube.append(read_scan1(fitsfile))
if nr > 104:
break
# making a contiguous C-array out of the slow list container
cube = np.array(cube)
cube.shape
Out[4]:
In [ ]:
Above dimensions indicate that I have stacked 50 spectra with each 146x1024 pixels.
In [5]:
moon_master = cube.mean(axis=0)
moon_master.shape
Out[5]:
Again, the dimensions indicate the correct result that I averaged the first axis, resulting in an image with the same dimensions as each original spectral image.
Now on to some filtering.
In [6]:
from scipy.ndimage import median_filter
In [143]:
imshow(moon_master, aspect='auto')
colorbar(orientation='horizontal')
title('Moon master raw')
savefig('moon_master_original.png', dpi=200)
In [11]:
moon_master_filtered = median_filter(moon_master, 2)
In [12]:
imshow(moon_master_filtered, aspect='auto')
colorbar(orientation='horizontal')
title('Moon master median filter (size 2)')
savefig('hermite_master_median_filtered.png', dpi=200)
In [23]:
cube = []
for fitsfile in fitsfiles:
# fish out the image number from the filename
nr = int(fitsfile.split('-')[1][:5])
if nr < 107:
continue
cube.append(read_scan1(fitsfile))
# making a contiguous C-array out of the slow list container
sky_cube = np.array(cube)
sky_cube.shape
Out[23]:
In [24]:
sky_master = sky_cube.mean(axis=0)
sky_master.shape
Out[24]:
In [145]:
imshow(sky_master, aspect='auto', vmax=20000)
colorbar(orientation='horizontal')
title('sky master, raw, clipping at 20000')
savefig('sky_master_raw.png', dpi=200)
In [26]:
# do a minimal pepper/salt filter
sky_master_filtered = median_filter(sky_master, 3)
In [144]:
imshow(sky_master_filtered, aspect='auto', vmax=20000)
colorbar(orientation='horizontal')
title('sky master, median size 3 filter, clipping at 20000')
savefig('sky_master_median_filtered.png', dpi=200)
Taking an upper limit of 20000 counts to exclude hot pixels.
In [28]:
sky_master_filtered[sky_master_filtered > 20000]= np.nan
In [29]:
master = moon_master_filtered - sky_master_filtered
In [30]:
rcParams['image.aspect'] = 'auto'
In [31]:
imshow(master)
colorbar(orientation='horizontal')
title('Background("sky")-subtracted master')
savefig('master_subtracted.png',dpi=200)
Now, having the master, let's cut out regions.
In [32]:
master.shape
Out[32]:
In [33]:
region1 = master[:36, :]
region2 = master[37:73, :]
region3 = master[74:110, :]
region4 = master[110:,:]
In [34]:
for region in [region1,region2,region3,region4]:
print region.shape
In [35]:
spec1 = region1.mean(axis=0)
spec2 = region2.mean(axis=0)
spec3 = region3.mean(axis=0)
spec4 = region4.mean(axis=0)
In [36]:
for i,spec in enumerate([spec1,spec2,spec3,spec4]):
plot(spec, label=i+1)
legend(loc='best')
title('master spectra for each region')
savefig('master_spectra_in_counts.png',dpi=200)
In [37]:
rcParams['lines.linewidth'] = 1.5
In [14]:
moon_regions = []
moon_regions.append(moon_master_filtered[:36, :])
moon_regions.append(moon_master_filtered[37:73, :])
moon_regions.append(moon_master_filtered[74:110, :])
moon_regions.append(moon_master_filtered[110:,:])
Trying same with the raw master, to see if it has an effect...
In [115]:
moon_regions_raw = []
moon_regions_raw.append(moon_master[:36, :])
moon_regions_raw.append(moon_master[37:73, :])
moon_regions_raw.append(moon_master[74:110, :])
moon_regions_raw.append(moon_master[110:,:])
In [ ]:
# only if we have them, obviously
sky_regions = []
sky_regions.append(sky_master_filtered[:36, :])
sky_regions.append(sky_master_filtered[37:73, :])
sky_regions.append(sky_master_filtered[74:110, :])
sky_regions.append(sky_master_filtered[110:, :])
In [16]:
# these are the wavelength conversions from pixel number in the wavelegnh
waves = pd.read_excel('lucy_spectral_cal.xlsx', parse_cols=0)
In [123]:
plot(waves, moon_regions[3].mean(axis=0), label='salt/pepper filtered')
plot(waves, moon_regions_raw[3].mean(axis=0), label='raw')
Out[123]:
In [113]:
fig, axes = subplots(nrows=2)
for i, moon_region in enumerate(moon_regions):
axes[0].plot(waves, moon_region.mean(axis=0),
label='moon{}'.format(i+1))
axes[0].set_title('Moon regions master spectra')
axes[0].legend(loc='best')
for i, sky_region in enumerate(sky_regions):
axes[1].plot(waves, sky_region.mean(axis=0),
label='sky{}'.format(i+1))
axes[1].set_title('Sky regions master spectra')
axes[1].legend(loc='best')
axes[1].set_xlabel('Wavelengths [micron]')
axes[0].set_ylabel('Digital counts unscaled')
axes[1].set_ylabel('Digital counts unscaled')
savefig('moon_and_sky_separated_regions.png', dpi=200)
In [21]:
# copied from above, but axes2 removed,
# as we don't have sky data
for i, moon_region in enumerate(moon_regions):
plot(waves, moon_region.mean(axis=0),
label='moon{}'.format(i+1))
title('Hermite regions master spectra')
legend(loc='best')
xlabel('Wavelengths [micron]')
ylabel('Digital counts unscaled')
savefig('hermite_regions_spectra.png', dpi=200)
In [19]:
#sky_masters = []
#for sky in sky_regions:
# sky_masters.append(sky.mean(axis=0))
moon_masters = []
for moon in moon_regions:
moon_masters.append(moon.mean(axis=0))
In [125]:
fig, axes = subplots(nrows=2)
for t in [(4,3),(4,2),(4,1), (3,2), (3,1), (2,1)]:
axes[0].plot(waves, sky_masters[t[0]-1]/sky_masters[t[1]-1],
label='{}/{}'.format(*t))
axes[1].plot(waves, moon_masters[t[0]-1]/moon_masters[t[1]-1],
label='{}/{}'.format(*t))
axes[0].set_title('Sky regional master spectra ratios')
axes[1].set_title('Moon regional master spectra ratios')
for ax in axes:
ax.legend(loc='best')
axes[0].set_ylim(0.8, 1.2)
axes[1].set_ylim(0.8, 1.8)
axes[1].set_xlabel("Wavelengths [micron]")
savefig('individual_moon_and_sky_ratios.png', dpi=200)
In [20]:
# copied from above, but removed sky parts
for t in [(4,3),(4,2),(4,1), (3,2), (3,1), (2,1)]:
plot(waves, moon_masters[t[0]-1]/moon_masters[t[1]-1],
label='{}/{}'.format(*t))
title('Hermite regional master spectra ratios')
legend(loc='best')
xlabel("Wavelengths [micron]")
savefig('hermite_ratios.png', dpi=200)
In [132]:
fig, axes = subplots(nrows=1)
for i in range(4):
axes.plot(waves, moon_masters[i]/sky_masters[i], label=i)
axes.legend(loc='best')
Out[132]:
In [137]:
for i in range(0,1000, 100):
plot(moon_master[:, i:i+100].mean(axis=1), label=i)
legend(loc='best')
Out[137]:
In [ ]: