In [1]:
import warnings
import time
from datetime import datetime
import os
import sys
import numpy as np
import matplotlib
%matplotlib ipympl
import matplotlib.pyplot as plt
import pandas as pd
import multiprocessing
from multiprocessing import Pool
import photutils
from pathlib import Path
import astropy.units as u
from astropy import stats
from astropy.io import fits
from astropy.table import Table, hstack, vstack
from astropy.modeling.models import Gaussian2D, Polynomial2D, Moffat2D
from astropy.modeling.fitting import LevMarLSQFitter
In [2]:
%cd "/Volumes/LACIE SHARE/bino_gdr"
root = Path(".")
In [3]:
%ls
In [4]:
hdr_keys = ['DATE-OBS', 'RA', 'DEC', 'AZ', 'EL', 'ROT', 'TEMP', 'FOCUS', 'EXPTIME']
guiders = ['GDR1', 'GDR2']
darks = {}
darks['GDR1'] = fits.open(root / "GDR1" / "box1_median_dark.fits")[0].data
darks['GDR2'] = fits.open(root / "GDR2" / "box1_median_dark.fits")[0].data
# 13 um pixels binned 2x2 and 0.167 mm/arcsec plate scale
pix_scale = 2 * 0.013 / 0.167 # arcsec/pixel
apsize = int(5. / pix_scale) # 5" arcsec aperture width
print(apsize)
In [5]:
def process_image(fitsfile, retall=False):
with fits.open(fitsfile) as hdu:
hdr = hdu[0].header
data = hdu[0].data
try:
mean, median, std = stats.sigma_clipped_stats(data, sigma=4.0, iters=5)
except:
if retall:
return None, data, None, None, None
else:
return None
# in some images there's a big step in the background so we take the median along each row and subtract that vector
m = np.median(data, axis=1)
data = data - m[:, None]
daofind = photutils.DAOStarFinder(fwhm=8.0, threshold=4.*std, sharphi=0.9, exclude_border=True)
srcs = daofind(data)
nsrcs = len(srcs)
if nsrcs == 0:
print(f"No star detected in {fitsfile}")
if retall:
return None, data, None, None, None
else:
return None
srcs.sort('flux')
srcs.reverse()
apers = photutils.CircularAperture(
(srcs['xcentroid'][0], srcs['ycentroid'][0]),
r=apsize/2.
)
masks = apers.to_mask(method="subpixel")
subim = masks[0].cutout(data)
props = photutils.data_properties(subim)
tline = hstack([srcs[0], props.to_table()])
tline['filename'] = fitsfile
try:
for k in hdr_keys:
tline[k] = hdr[k]
except:
if retall:
return None, data, None, None, None
else:
return None
y, x = np.mgrid[:data.shape[0], :data.shape[1]]
sigma = (tline['semimajor_axis_sigma'][0].value + tline['semiminor_axis_sigma'][0].value) / 2.
fitter = LevMarLSQFitter()
gauss_model = Gaussian2D(
amplitude=data.max(),
x_mean=data.shape[1]/2.,
y_mean=data.shape[0]/2.,
x_stddev = sigma,
y_stddev = sigma
)
moffat_model = Moffat2D(
amplitude=data.max(),
x_0=data.shape[1]/2.,
y_0=data.shape[0]/2.,
gamma=sigma
)
gauss_fit = fitter(gauss_model, x, y, data)
moffat_fit = fitter(moffat_model, x, y, data)
if retall:
print(gauss_fit)
print(moffat_fit)
gamma = moffat_fit.gamma.value
alpha = moffat_fit.alpha.value
moffat_fwhm = pix_scale * 2. * gamma * np.sqrt(2.**(1./alpha) - 1.)
gauss_fwhm = pix_scale * 0.5 * (gauss_fit.x_stddev.value + gauss_fit.y_stddev.value) * stats.gaussian_sigma_to_fwhm
moment_fwhm = pix_scale * 0.5 * (tline['semimajor_axis_sigma'] + tline['semiminor_axis_sigma']) * stats.gaussian_sigma_to_fwhm
tline['gauss_x'] = gauss_fit.x_mean.value
tline['gauss_y'] = gauss_fit.y_mean.value
tline['gauss_sigx'] = gauss_fit.x_stddev.value
tline['gauss_sigy'] = gauss_fit.y_stddev.value
tline['gauss_amplitude'] = gauss_fit.amplitude.value
tline['gauss_theta'] = gauss_fit.theta.value
tline['moffat_amplitude'] = moffat_fit.amplitude.value
tline['moffat_gamma'] = gamma
tline['moffat_alpha'] = alpha
tline['moffat_x'] = moffat_fit.x_0.value
tline['moffat_y'] = moffat_fit.y_0.value
tline['stats_mean'] = mean
tline['stats_median'] = median
tline['stats_std'] = std
tline['moment_fwhm'] = moment_fwhm
tline['gauss_fwhm'] = gauss_fwhm
tline['moffat_fwhm'] = moffat_fwhm
resid = data - moffat_fit(x, y)
if retall:
return tline, data, resid, gauss_fit, moffat_fit
else:
del data
del hdr
del x
del y
del resid
return tline
In [6]:
fdir = root / "GDR2" / "2018.0412"
In [7]:
files = sorted(list(fdir.glob("gdr*.fits")))
In [8]:
len(files)
Out[8]:
In [9]:
t, i, resid, gf, mf = process_image("GDR2/2018.0412/gdr2_raw_box_img_2018.0412.111604.1.fits", retall=True)
In [10]:
t
Out[10]:
In [11]:
%matplotlib ipympl
plt.imshow(resid)
plt.show()
In [12]:
nproc = 8
with Pool(processes=nproc) as pool: # my mac's i7 has 4 cores + hyperthreading so 8 virtual cores.
plines = pool.map(process_image, files) # plines comes out in same order as fitslines!
plines = list(filter(None.__ne__, plines)) # trim out any None entries
df = Table(vstack(plines)).to_pandas()
In [13]:
len(df)
Out[13]:
In [14]:
%matplotlib ipympl
df.hist(column=["moffat_fwhm", 'gauss_fwhm', 'moment_fwhm'], range=(0, 4), bins=100, alpha=0.6)
df.hist(column=["moffat_alpha"], range=(0, 10), bins=100, alpha=0.6)
Out[14]:
In [15]:
df.to_csv("2018_04_12.csv")
In [59]:
df.columns
Out[59]:
In [21]:
df.hist(column=['ROT'], range=(-180, 180.0), bins=100, alpha=0.6)
Out[21]:
In [16]:
df['DATE-OBS']
Out[16]:
In [17]:
import gc
In [18]:
gc.collect()
Out[18]:
In [14]:
test = fits.open("GDR2/2018.0209/gdr2_raw_box_img_2018.0209.114833.1.fits")[0].data
In [28]:
m = np.median(test, axis=1)
c = test - m[:, None]
In [29]:
%matplotlib ipympl
plt.imshow(c)
plt.show()
In [ ]: