In [12]:
import time
from datetime import datetime
import pytz
import os
import sys
import numpy as np
import matplotlib
%matplotlib ipympl
import matplotlib.pyplot as plt
import pandas as pd
from matplotlib import style
style.use('ggplot')
import multiprocessing
from multiprocessing import Pool
from pathlib import Path
import pandas as pd
import photutils
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
from mmtwfs.wfs import *
from mmtwfs.zernike import ZernikeVector
from mmtwfs.telescope import MMT
from mmtwfs.custom_exceptions import WFSAnalysisFailed
tz = pytz.timezone("America/Phoenix")
In [2]:
%load_ext line_profiler
%load_ext autoreload
%autoreload 2
In [3]:
hdr_keys = ['DATE-OBS', 'RA', 'DEC', 'AZ', 'EL', 'ROT', 'TEMP', 'FOCUS', 'CHAM_DPT', 'CHAM_T', 'OUT_T', 'WIND_W', 'WDIR_W', 'PRESSURE',
'EXPTIME', 'AIRMASS', 'OSSTEMP', 'TRANSX', 'TRANSY', 'TILTX', 'TILTY']
f5_hdr_keys = ['DATEOBS', 'UT', 'RA', 'DEC', 'AZ', 'EL', 'ROT', 'FOCUS', 'T_OUT', 'T_CHAM', 'RH_OUT', 'RH_CHAM', 'P_BARO', 'EXPTIME', 'AIRMASS',
'WIND', 'WINDDIR', 'OSSTEMP', 'TRANSX', 'TRANSY', 'TILTX', 'TILTY']
old_f9_keys = ['DATE-OBS', 'TIME-OBS', 'RA', 'DEC', 'AZ', 'EL', 'ROT', 'FOCUS', 'T_OUT', 'T_CHAM', 'RH_OUT', 'RH_CHAM', 'P_BARO', 'EXPTIME', 'AIRMASS',
'WIND', 'WINDDIR', 'OSSTEMP', 'TRANSX', 'TRANSY', 'TILTX', 'TILTY']
In [10]:
def check_image(f, wfskey=None):
hdr = {}
with fits.open(f) as hdulist:
for h in hdulist:
hdr.update(h.header)
data = hdulist[-1].data
# if wfskey is None, figure out which WFS from the header info...
if wfskey is None:
# check for MMIRS
if 'WFSNAME' in hdr:
if 'mmirs' in hdr['WFSNAME']:
wfskey = 'mmirs'
if 'mmirs' in f.name:
wfskey = 'mmirs'
# check for binospec
if 'bino' in f.name:
wfskey = 'binospec'
if 'ORIGIN' in hdr:
if 'Binospec' in hdr['ORIGIN']:
wfskey = 'binospec'
# check for new F/9
if 'f9wfs' in f.name:
wfskey = 'newf9'
if 'OBSERVER' in hdr:
if 'F/9 WFS' in hdr['OBSERVER']:
wfskey = 'newf9'
if wfskey is None and 'CAMERA' in hdr:
if 'F/9 WFS' in hdr['CAMERA']:
wfskey = 'newf9'
# check for old F/9
if 'INSTRUME' in hdr:
if 'Apogee' in hdr['INSTRUME']:
wfskey = 'oldf9'
if 'DETECTOR' in hdr:
if 'Apogee' in hdr['DETECTOR']:
wfskey = 'oldf9'
# check for F/5 (hecto)
if wfskey is None and 'SEC' in hdr: # mmirs has SEC in header as well and is caught above
if 'F5' in hdr['SEC']:
wfskey = 'f5'
if Path(f.parent / "F5").exists():
wfskey = 'f5'
if wfskey is None:
# if wfskey is still None at this point, whinge.
print(f"Can't determine WFS for {f.name}...")
if 'AIRMASS' not in hdr:
if 'SECZ' in hdr:
hdr['AIRMASS'] = hdr['SECZ']
else:
hdr['AIRMASS'] = np.nan
# we need to fix the headers in all cases to have a proper DATE-OBS entry with
# properly formatted FITS timestamp. in the meantime, this hack gets us what we need
# for analysis in pandas.
dtime = None
if 'DATEOBS' in hdr:
dateobs = hdr['DATEOBS']
if 'UT' in hdr:
ut = hdr['UT'].strip()
elif 'TIME-OBS' in hdr:
ut = hdr['TIME-OBS']
else:
ut = "07:00:00" # midnight
timestring = dateobs + " " + ut + " UTC"
if '-' in timestring:
dtime = datetime.strptime(timestring, "%Y-%m-%d %H:%M:%S %Z")
else:
dtime = datetime.strptime(timestring, "%a %b %d %Y %H:%M:%S %Z")
else:
if wfskey == "oldf9":
d = hdr['DATE-OBS']
if '/' in d:
day, month, year = d.split('/')
year = str(int(year) + 1900)
timestring = year + "-" + month + "-" + day + " " + hdr['TIME-OBS'] + " UTC"
else:
timestring = d + " " + hdr['TIME-OBS'] + " UTC"
dtime = datetime.strptime(timestring, "%Y-%m-%d %H:%M:%S %Z")
else:
if 'DATE-OBS' in hdr:
timestring = hdr['DATE-OBS'] + " UTC"
try:
dtime = datetime.strptime(timestring, "%Y-%m-%dT%H:%M:%S.%f %Z")
except:
dtime = datetime.strptime(timestring, "%Y-%m-%dT%H:%M:%S %Z")
else:
dt = datetime.fromtimestamp(f.stat().st_ctime)
local_dt = tz.localize(dt)
dtime = local_dt.astimezone(pytz.utc)
if dtime is None:
print(f"No valid timestamp in header for {f.name}...")
obstime = None
else:
obstime = dtime.isoformat().replace('+00:00', '')
hdr['WFSKEY'] = wfskey
hdr['OBS-TIME'] = obstime
return data, hdr
In [5]:
findpars = {
"oldf9": {"fwhm": 8.0, "thresh": 5.0},
"newf9": {"fwhm": 12.0, "thresh": 7.0},
"f5": {"fwhm": 9.0, "thresh": 5.0},
"mmirs": {"fwhm": 7.0, "thresh": 5.0},
"binospec": {"fwhm": 7.0, "thresh": 5.0}
}
In [6]:
#%cd /Users/tim/MMT/wfsdat/20180209
%cd /Volumes/LACIE\ SHARE/wfsdat
In [11]:
file = Path("20170107/manual_wfs_0000.fits")
data, hdr = check_image(file)
data.shape, hdr
Out[11]:
In [ ]:
def process_image(f, clobber=False):
tablefile = Path(str(f.parent / f.stem) + ".csv")
if '_ave' in str(f):
print(f"Not processing {f.name} because it's an average of multiple images")
return
if tablefile.exists() and not clobber:
print(f"{f.name} already processed...")
return
else:
print(f"Processing {f.name}...")
try:
data, hdr = check_image(f)
w = hdr['WFSKEY']
mean, median, stddev = stats.sigma_clipped_stats(data, sigma=3.0, iters=None)
data = data - median
spots, fig = wfsfind(data, fwhm=findpars[w]['fwhm'], threshold=findpars[w]['thresh'], plot=False)
except WFSAnalysisFailed as e:
try:
spots, fig = wfsfind(data, fwhm=2.*findpars[w]['fwhm'], threshold=findpars[w]['thresh'], plot=False)
except Exception as e:
print(f"Failed to find spots for {str(Path(str(f.parent / f.stem)))}: {e}")
return
except Exception as ee:
print(f"Failure in data loading or spot finding in {str(Path(str(f.parent / f.stem)))}: {ee}")
return
apsize = 15.
apers = photutils.CircularAperture(
(spots['xcentroid'], spots['ycentroid']),
r=apsize
)
masks = apers.to_mask(method="subpixel")
props = []
spot_lines = []
fit_lines = []
for m, s in zip(masks, spots):
tline = {}
subim = m.cutout(data)
try:
props_table = photutils.data_properties(subim).to_table()
except Exception as e:
print(f"Can't measure source properties for {str(Path(str(f.parent / f.stem)))}: {e}")
continue
moment_fwhm = 0.5 * (props_table['semimajor_axis_sigma'][0].value + props_table['semiminor_axis_sigma'][0].value) * stats.gaussian_sigma_to_fwhm
props.append(props_table)
spot_lines.append(s)
tline['filename'] = f.name
for k in hdr:
if 'COMMENT' not in k and k != '':
tline[k] = hdr[k]
y, x = np.mgrid[:subim.shape[0], :subim.shape[1]]
sigma = (props_table['semimajor_axis_sigma'][0].value + props_table['semiminor_axis_sigma'][0].value) / 2.
fitter = LevMarLSQFitter()
gauss_model = Gaussian2D(
amplitude=subim.max(),
x_mean=subim.shape[1]/2.,
y_mean=subim.shape[0]/2.,
x_stddev = sigma,
y_stddev = sigma
) + Polynomial2D(degree=0)
moffat_model = Moffat2D(
amplitude=subim.max(),
x_0=subim.shape[1]/2.,
y_0=subim.shape[0]/2.,
gamma=sigma
) + Polynomial2D(degree=0)
try:
gauss_fit = fitter(gauss_model, x, y, subim)
gauss_resid = subim - gauss_fit(x, y)
gauss_fwhm = 0.5 * (gauss_fit.x_stddev_0.value + gauss_fit.y_stddev_0.value) * stats.gaussian_sigma_to_fwhm
tline['gauss_x'] = gauss_fit.x_mean_0.value
tline['gauss_y'] = gauss_fit.y_mean_0.value
tline['gauss_sigx'] = gauss_fit.x_stddev_0.value
tline['gauss_sigy'] = gauss_fit.y_stddev_0.value
tline['gauss_amplitude'] = gauss_fit.amplitude_0.value
tline['gauss_theta'] = gauss_fit.theta_0.value
tline['gauss_background'] = gauss_fit.c0_0_1.value
tline['gauss_rms'] = np.nanstd(gauss_resid)
tline['gauss_fwhm'] = gauss_fwhm
except Exception as e:
print(f"Gaussian fit failed in {str(Path(str(f.parent / f.stem)))}: {e}")
tline['gauss_x'] = np.nan
tline['gauss_y'] = np.nan
tline['gauss_sigx'] = np.nan
tline['gauss_sigy'] = np.nan
tline['gauss_amplitude'] = np.nan
tline['gauss_theta'] = np.nan
tline['gauss_background'] = np.nan
tline['gauss_rms'] = np.nan
tline['gauss_fwhm'] = np.nan
try:
moffat_fit = fitter(moffat_model, x, y, subim)
moffat_resid = subim - moffat_fit(x, y)
gamma = moffat_fit.gamma_0.value
alpha = moffat_fit.alpha_0.value
moffat_fwhm = np.abs(2. * gamma * np.sqrt(2.**(1./alpha) - 1.))
tline['moffat_amplitude'] = moffat_fit.amplitude_0.value
tline['moffat_gamma'] = gamma
tline['moffat_alpha'] = alpha
tline['moffat_x'] = moffat_fit.x_0_0.value
tline['moffat_y'] = moffat_fit.y_0_0.value
tline['moffat_background'] = moffat_fit.c0_0_1.value
tline['moffat_rms'] = np.nanstd(moffat_resid)
tline['moment_fwhm'] = moment_fwhm
tline['moffat_fwhm'] = moffat_fwhm
except Exception as e:
print(f"Moffat fit failed in {str(Path(str(f.parent / f.stem)))}: {e}")
tline['moffat_amplitude'] = np.nan
tline['moffat_gamma'] = np.nan
tline['moffat_alpha'] = np.nan
tline['moffat_x'] = np.nan
tline['moffat_y'] = np.nan
tline['moffat_background'] = np.nan
tline['moffat_rms'] = np.nan
tline['moment_fwhm'] = np.nan
tline['moffat_fwhm'] = np.nan
fit_lines.append(tline)
fit_table = Table(fit_lines)
spot_table = Table(vstack(spot_lines))
prop_table = Table(vstack(props))
t = hstack([spot_table, prop_table, fit_table])
try:
if tablefile.exists():
if clobber:
t.write(tablefile, format="csv")
else:
t.write(tablefile, format="csv")
except Exception as e:
print(f"Failed to write {str(tablefile)}: {e}")
return t
In [ ]:
process_image(file, clobber=True)
In [ ]:
t.hist(column='gauss_fwhm', bins=50, range=(0,10), alpha=0.6)
t.hist(column='moffat_fwhm', bins=50, range=(0,10), alpha=0.6)
plt.show()
In [ ]:
t['gauss_background'].mean(), t['moffat_background'].mean()
In [ ]:
t.hist(column='ellipticity', bins=50, range=(0,1.0), alpha=0.6)
In [ ]:
(t['moffat_rms']/t['moffat_amplitude']).mean()
In [ ]:
%matplotlib ipympl
plt.imshow(gauss_resids[100])
plt.show()
In [ ]:
#%matplotlib ipympl
plt.imshow(moffat_resids[10])
plt.show()
In [ ]:
%lprun -f process_image process_image(file)
In [ ]:
rootdir = Path("/Volumes/LACIE SHARE/wfsdat")
dirs = sorted(list(rootdir.glob("20071*")))
In [ ]:
nproc = 8
for d in dirs:
if d.is_dir():
print(f"Working in {d.name}...")
fitsfiles = sorted(list(d.glob("*.fits")))
with Pool(processes=nproc) as pool: # my mac's i7 has 4 cores + hyperthreading so 8 virtual cores.
pool.map(process_image, fitsfiles) # plines comes out in same order as fitslines!
In [ ]:
test = rootdir / "20030321" / "end_0000.fits"
In [ ]:
str(test.parent / test.stem) + ".csv"
In [ ]:
test, fig = wfsfind("20031019/start_0001.fits", fwhm=8.)
In [ ]:
test, hdr = check_image(Path("20031019/start_0001.fits"))
mean, median, stddev = stats.sigma_clipped_stats(test, sigma=3.0, iters=None)
In [ ]:
test_table = process_image(Path("20031019/start_0001.fits"))
In [ ]:
tstfile = Path("20030517/wfs0029.fits")
In [ ]:
dt = datetime.fromtimestamp(tstfile.stat().st_ctime)
local_dt = tz.localize(dt)
utc_dt = local_dt.astimezone(pytz.utc).isoformat().replace('+00:00', '')
utc_dt
In [ ]:
test_table
In [ ]:
hdr
In [20]:
a = np.random.random(10000000)
b = np.random.random(10000000)
In [21]:
%%timeit
t = Table([a, b])
In [22]:
%%timeit
df = pd.DataFrame([a, b])
In [23]:
del df
In [ ]: