Subarray determination

As we took the data in a sub-array configuration, we first need to determine where in the FITS data the subarray is positioned.



In [1]:

    
# 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'] = 'hot'









    



Populating the interactive namespace from numpy and matplotlib

This is the Python module that has the FITS reader (and writer).



In [2]:

    
import pyfits

Reading in filenames that are fits

For the run scan5 we forgot to rename the filenames, so the data files are still called limb-xxx (from the run scan4) even so the data is from the Hermite crater, taken with the prism (in contrast to the grate we used for scan3 on day 1.)

going to scan5 for now.



In [3]:

    
fitsfiles = !ls /Users/maye/data/irtf/bigdog/scan5/limb*.fits



In [4]:

    
# plot 1 test image
data = pyfits.getdata(fitsfiles[0])
imshow(data)
plt.axis('off')









    Out[4]:





(-0.5, 1023.5, 649.5, -0.5)

The next plot determines the best top row for where the data starts. As soon as the rows become indistinguishable no refraction from the slit should be visible anymore.



In [10]:

    
for line in range(245,255):
    plot(data[line, :], label=line)
legend(loc='best')









    Out[10]:





<matplotlib.legend.Legend at 0x10dc26b10>

We believe row 247 (the red line) still shows too many differences so we should include only above it. Now the same for the higher row number:



In [11]:

    
for line in range(390, 400):
    plot(data[line, :], label=line)
legend(loc='best')









    Out[11]:





<matplotlib.legend.Legend at 0x10dc26ad0>

It looks like up to row 394 (0-indexing) the data is reliable. Creating an appropriate subframe reader:



In [6]:

    
def read_scan1(fname):
    data = pyfits.getdata(fname, 0)
    return data[248:394, :]



In [7]:

    
data = read_scan1(fitsfiles[0])



In [8]:

    
data.shape









    Out[8]:





(146, 1024)

Moon master

Above shape is the shape of one spectral image. Now reading in all data into a 3D cube. The on_moon data goes until image number 104.



In [9]:

    
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[9]:





(50, 146, 1024)

Above dimensions indicate that I have stacked 50 spectra with each 146x1024 pixels.



In [10]:

    
moon_master = cube.mean(axis=0)
moon_master.shape









    Out[10]:





(146, 1024)

Again, the dimensions indicate the correct result that I averaged the first axis, and therefore have a 'moon_master' image:



In [13]:

    
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 [14]:

    
moon_master_filtered = median_filter(moon_master, 2)



In [16]:

    
imshow(moon_master_filtered, aspect='auto')
colorbar(orientation='horizontal')
title('Hermite master median filter (size 2)')
savefig('hermite_master_median_filtered.png', dpi=200)

Sky master

Starting to read from 107 (106 might be okay as well).



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]:





(33, 146, 1024)



In [24]:

    
sky_master = sky_cube.mean(axis=0)
sky_master.shape









    Out[24]:





(146, 1024)



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

Regions

For comparison with Raquel's analysis, I also will split the spatial region into 4 bands, each with 36 pixels width. First, let's substract the sky.



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]:





(146, 1024)



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









    



(36, 1024)
(36, 1024)
(36, 1024)
(36, 1024)



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

Regions before subtracting.



In [17]:

    
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 [47]:

    
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 [18]:

    
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]:





[<matplotlib.lines.Line2D at 0x125412090>]



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 [20]:

    
moon_regions[0].shape









    Out[20]:





(36, 1024)



In [21]:

    
spectrum = moon_regions[0].mean(axis=0)



In [24]:

    
spectrum.tofile?



In [ ]:

    
[(i, j) for i in range(4) and j in range(3)]



In [67]:

    
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 [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]:





<matplotlib.legend.Legend at 0x11f7adf90>



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]:





<matplotlib.legend.Legend at 0x126a1b950>



In [ ]: