In [133]:
import numpy as np
import matplotlib.pyplot as plt
%matplotlib notebook
from astropy.table import Table
from cesium import featurize
import sqlite3
from astropy.stats import sigma_clip as sc
from matplotlib.ticker import MultipleLocator
# plt.ioff()
In [24]:
filename = input("Enter .tbl filename: ")
dataTable = Table.read(filename, format='ipac')
times = dataTable["obsmjd"]
values = dataTable["mag_autocorr"]
errors = dataTable["magerr_auto"]
oids = dataTable["oid"]
fids = dataTable["fid"]
printFrequency(times, values, errors)
frequency = float(input("Enter the frequency: "))
period = 1/frequency
In [8]:
def printFrequency(times, values, errors):
""" This function obtains the frequency using cesium, which implies the period.
Arguments:
dataTable (array-like): this is the .tbl data file downloaded from Caltech
IRSA website
"""
# freq1_amplitude1: Get the amplitude of the jth harmonic of the ith frequency from a fitted Lomb-Scargle model.
features_to_use = ["freq1_freq", "amplitude", "freq1_amplitude1"]
fset_cesium = featurize.featurize_time_series(times, values, errors, features_to_use=features_to_use)
print(fset_cesium)
In [3]:
def zoo_plot(times, values, errors, indexesTempOID, tempOID, fmtToBeUsed, period):
fig, (ax1, ax2) = plt.subplots(1, 2, figsize = (10, 5))
# raw light curve
ax1.errorbar(times[indexesTempOID], values[indexesTempOID], yerr=errors[indexesTempOID], fmt=fmtToBeUsed,
markersize=5.75, mec="k", mew=0.59)
ax1.set_ylim(ax1.get_ylim()[::-1])
ax1.set_title("LC for oid={} (raw data)".format(tempOID))
ax1.set_xlabel("time")
ax1.set_ylabel("brightness (mag)")
ax1.yaxis.grid(color='gray', linestyle='dashed')
ax1.xaxis.grid(color='gray', linestyle='dashed')
# phase-folded light curve
tempLength = len(indexesTempOID[0])
timesTemp = np.empty(tempLength)
timesPhaseFolded = np.empty(tempLength)
timesTemp = times[indexesTempOID]
for i in range(tempLength):
timesPhaseFolded[i] = (timesTemp[i] % period) / period
timesPhaseFoldedPlusOne = np.empty(tempLength)
for j in range(tempLength):
timesPhaseFoldedPlusOne[j] = timesPhaseFolded[j]+1
ax2.errorbar(timesPhaseFolded, values[indexesTempOID], yerr=errors[indexesTempOID], fmt=fmtToBeUsed,
markersize=5.75, mec="k", mew=0.59)
ax2.errorbar(timesPhaseFoldedPlusOne, values[indexesTempOID], yerr=errors[indexesTempOID], fmt=fmtToBeUsed,
markersize=5.75, mec="k", mew=0.59)
ax2.set_ylim(ax2.get_ylim()[::-1]) # reversing the y-axis as desired
ax2.set_title("LC for oid={} (phase-folded)".format(tempOID))
ax2.set_xlabel("phase")
ax2.yaxis.grid(color='gray', linestyle='dashed')
ax2.xaxis.grid(color='gray', linestyle='dashed')
fig.tight_layout()
In [27]:
def IPAC_Table_Parser(times, values, errors, oids, fids, period): # took out period as argument
""" This Zooniverse light curve plotting
function produces phase-folded light
curves and non-phase-folded light curves
(i.e. the raw data) with error bars for each
unique object ID from data (.tbl) file obtained
from Caltech PTF database.
Arguments:
times: array of time values
values: array of magnitude values
errors: array of errors on the magnitudes
oids: array of object ID's
fids: array of filter ID's
period (float): the period of the light
curve for phase folding (obtained using cesium)
"""
length = len(oids)
oidsArray = np.empty(length)
for i in range(length):
oidsArray[i] = oids[i]
oidsArraySorted = np.empty(length)
oidsArraySorted = np.sort(oidsArray)
fmtToBeUsed = 'ro'
index = 0
swapVar = 0
tempOID = 0
newOID = int(oidsArraySorted[index])
while newOID != tempOID: # while loop because number of unique OID's is not fixed
swapVar = newOID
tempOID = swapVar
indexesTempOID = 0
indexesTempOID = np.where(dataTable[3][:] == tempOID)
fmtToBeUsedNumber = fids[indexesTempOID[0][0]]
# Filter identifier (1 = g; 2 = R)
if fmtToBeUsedNumber == 2:
fmtToBeUsed = 'ro'
else:
fmtToBeUsed = 'go'
# THIS IS WHERE THE ERROR IS! TRYING TO APPLY CESIUM TO EACH OID SEPARATELY
#print("times[indexesTempOID]:", times[indexesTempOID])
#print("values[indexesTempOID]:", values[indexesTempOID])
#print("errors[indexesTempOID]:", errors[indexesTempOID])
zoo_plot(times, values, errors, indexesTempOID, tempOID, fmtToBeUsed, period)
index += 1
newOID = int(oidsArraySorted[index])
# at this point, newOID will equal tempOID, unless the there is only one oid of that variety
while (newOID == tempOID and index < (length-1)):
index += 1
newOID = int(oidsArraySorted[index])
In [28]:
IPAC_Table_Parser(times, values, errors, oids, fids, period)
In [97]:
def zoo_plot2(t, m, m_unc, period = 1.0, ptf_filt = "R"):
"""Plots PTF data for Zooniverse"""
color_dict = {"g": "MediumAquaMarine", "R": "FireBrick"}
#fig, (ax1, ax2) = plt.subplots(2, 1, figsize = (8,8), gridspec_kw = {'height_ratios':[1, 3]}) this works
fig, (ax1, ax2) = plt.subplots(1, 2, figsize = (10, 5))
# raw light curve
ax1.errorbar(t, m, m_unc,
fmt = "o", mec="k", mew=0.5,
color = color_dict[ptf_filt])
ax1.set_ylim(ax1.get_ylim()[::-1])
ax1.yaxis.grid(color='gray', linestyle='dashed')
ax1.xaxis.grid(color='gray', linestyle='dashed')
ax1.set_title("Light Curve (raw data)")
ax1.set_xlabel("Time (d)")
ax1.set_ylabel("Brightness (mag)")
# phase folded light curve
ax2.errorbar(t / period % 1, m, m_unc,
fmt = "o", mec="k", mew=0.5,
color = color_dict[ptf_filt])
ax2.errorbar((t / period % 1)+1, m, m_unc,
fmt = "o", mec="k", mew=0.5,
color = color_dict[ptf_filt])
ax2.set_ylim(ax2.get_ylim()[::-1])
ax2.yaxis.grid(color='gray', linestyle='dashed')
ax2.xaxis.grid(color='gray', linestyle='dashed')
ax2.set_title("Light Curve for Object (phase-folded)")
ax2.set_xlabel("Phase")
#ax2.set_ylabel("brightness (mag)")
fig.tight_layout()
return fig
In [6]:
conn = sqlite3.connect('/home/nick/Desktop/NUREU17/LSST/VariableStarClassification/features.db')
cur = conn.cursor()
In [106]:
cur.execute("""select oid, freq1_freq from sigfeats""")
raw = cur.fetchall()
freqs = list(zip(*raw))
for i in range(len(freqs[0])):
ipac_lc = Table.read('/home/nick/Desktop/NUREU17/LSST/VariableStarClassification/scripts/ptf_query/byOid/curves_oid{:_>17}.tbl'.format(freqs[0][i]), format = 'ipac')
color = ipac_lc['fid'][0]
mask = ~(sc(ipac_lc['mag_autocorr'])).mask
if color == 1:
fig = zoo_plot2(ipac_lc['obsmjd'][mask], ipac_lc['mag_autocorr'][mask], ipac_lc['magerr_auto'][mask], period = 1/freqs[1][i], ptf_filt = "g")
elif color == 2:
fig = zoo_plot2(ipac_lc['obsmjd'][mask], ipac_lc['mag_autocorr'][mask], ipac_lc['magerr_auto'][mask], period = 1/freqs[1][i], ptf_filt = "R")
plt.savefig('/home/nick/Desktop/NUREU17/LSST/VariableStarClassification/scripts/ptf_query/New_Data_Files/Figs/v01/figs_v01_{:_>17}.png'.format(freqs[0][i]))
plt.close(fig)
In [61]:
i = 0
print(freqs[1][i])
In [33]:
zoo_plot2(ipac_lc["obsmjd"], ipac_lc["mag_autocorr"], ipac_lc["magerr_auto"], period = 1.09296246931)
In [9]:
def zoo_plot3(dataTable, period):
""" This is a fully general light curve plotting function that works with the dataTable
downloaded from the Caltech IRSA website to produce phase-folded light curves and non-phase-folded
light curves (i.e. the raw data) with error bars for each observation.
Caltech IRSA website: http://irsa.ipac.caltech.edu/cgi-bin/Gator/nph-scan?projshort=PTF
Arguments:
dataTable (array-like) : data table from Caltech IRSA website (.tbl)
period (float) : the period of the light curve for phase folding (obtained with cesium)
"""
times = dataTable["obsmjd"]
values = dataTable["mag_autocorr"]
errors = dataTable["magerr_auto"]
oids = dataTable["oid"]
fids = dataTable["fid"]
length = len(oids)
oidsArray = np.empty(length)
for i in range(length):
oidsArray[i] = oids[i]
oidsArraySorted = np.empty(length)
oidsArraySorted = np.sort(oidsArray)
fmtToBeUsed = 'ro'
index = 0
swapVar = 0
tempOID = 0
newOID = int(oidsArraySorted[index])
while newOID != tempOID: # while loop because number of unique OID's is not fixed
swapVar = newOID
tempOID = swapVar
indexesTempOID = np.where(dataTable[3][:] == tempOID)
fmtToBeUsedNumber = fids[indexesTempOID[0][0]]
# Filter identifier (1 = g; 2 = R)
if fmtToBeUsedNumber == 2:
fmtToBeUsed = 'ro'
else:
fmtToBeUsed = 'go'
# returns figure object and axis object
fig, (ax1, ax2) = plt.subplots(2, 1, figsize = (8,8), gridspec_kw = {'height_ratios':[1, 3]})
# raw light curve
ax1.errorbar(times[indexesTempOID], values[indexesTempOID], yerr=errors[indexesTempOID], fmt=fmtToBeUsed,
markersize=3, mec="k", mew=0.59)
ax1.set_ylim(ax1.get_ylim()[::-1])
ax1.set_title("Light Curve for Object {} (raw data)".format(tempOID))
ax1.set_xlabel("time")
ax1.set_ylabel("brightness (mag)")
ax1.yaxis.grid(color='gray', linestyle='dashed')
ax1.xaxis.grid(color='gray', linestyle='dashed')
# phase-folded light curve
tempLength = len(indexesTempOID[0])
timesTemp = np.empty(tempLength)
timesPhaseFolded = np.empty(tempLength)
timesTemp = times[indexesTempOID]
for i in range(tempLength):
timesPhaseFolded[i] = (timesTemp[i] % period) / period
timesPhaseFoldedPlusOne = np.empty(tempLength)
for j in range(tempLength):
timesPhaseFoldedPlusOne[j] = timesPhaseFolded[j]+1
ax2.errorbar(timesPhaseFolded, values[indexesTempOID], yerr=errors[indexesTempOID], fmt=fmtToBeUsed,
markersize=3, mec="k", mew=0.59)
ax2.set_ylim(ax2.get_ylim()[::-1]) # reversing the y-axis as desired
ax2.set_title("Light Curve for Object {} (phase-folded)".format(tempOID))
ax2.set_xlabel("phase")
ax2.set_ylabel("brightness (mag)")
ax2.yaxis.grid(color='gray', linestyle='dashed')
ax2.xaxis.grid(color='gray', linestyle='dashed')
fig.tight_layout()
plt.show()
index += 1
newOID = int(oidsArraySorted[index])
# at this point, newOID will equal tempOID, unless the there is only one oid of that variety
while (newOID == tempOID and index < (length-1)):
index += 1
newOID = int(oidsArraySorted[index])
In [11]:
filename = input("Enter .tbl filename: ")
dataTable = Table.read(filename, format='ipac')
printFrequency(dataTable)
frequency = float(input("Enter the frequency: "))
period = 1/frequency
zoo_plot3(dataTable, period)
print("'period':", period)
In [2]:
filename = input("Enter .tbl filename: ")
In [26]:
ipac_lc = Table.read(filename, format='ipac')
hjd = np.array(ipac_lc["obsmjd"])
mag = np.array(ipac_lc["mag_autocorr"])
mag_unc = np.array(ipac_lc["magerr_auto"])
clipped_obs = sc(mag).mask
cs_feats = featurize.featurize_time_series(hjd[~clipped_obs],
mag[~clipped_obs],
mag_unc[~clipped_obs],
features_to_use=["freq1_freq",
"amplitude",
"freq1_amplitude1"])
period = 1/float(cs_feats["freq1_freq"][0])
In [44]:
def strateva(mag, a = 0.027665,
b = 1.285193044293278e-09,
c = 1.0999892234090642e-19):
return a + b*10**(0.4*mag) + c*10**(0.8*mag)
In [92]:
def lc_breaks(time, break_width=30):
breaks = np.where(np.diff(time[np.argsort(time)]) > break_width)
return breaks
In [174]:
def aam_plot(time, mag, mag_unc,
period=1, ptf_filt=2,
):
"""Plot a light curve suitable for Zooniverse classification
"""
color_dict = {1: "LightSeaGreen",
2: "Crimson"}
model_unc = strateva(mag)
mean_mag = np.mean(mag)
obs_breaks = lc_breaks(hjd[~clipped_obs])[0] + 1
seasons = len(obs_breaks) + 1
season_dict = {}
total_obs_duration = 0
for season_num in range(seasons):
if season_num == 0:
season = time[np.argsort(time)][0:obs_breaks[season_num]]
elif season_num == seasons - 1:
season = time[np.argsort(time)][obs_breaks[season_num-1]:]
else:
season = time[np.argsort(time)][obs_breaks[season_num-1]:obs_breaks[season_num]]
season_dict["start{:d}".format(season_num)] = min(season)
season_dict["end{:d}".format(season_num)] = max(season)
season_dict["length{:d}".format(season_num)] = np.round(np.ptp(season))
total_obs_duration += np.round(max(season) - min(season))
###### FIGURE #######
fig = plt.figure()
gs = GridSpec(4, int(total_obs_duration))
amplitude = np.ptp(mean_mag - mag)
ymin = min(mean_mag - mag) - 0.05*amplitude
ymax = max(mean_mag - mag) + 0.05*amplitude
ymajor = np.round(4*np.ptp(mean_mag - mag))/16
yminor = np.round(4*np.ptp(mean_mag - mag))/80
for season_num in range(seasons):
obs_length = season_dict["length{:d}".format(season_num)]
obs_start = season_dict["start{:d}".format(season_num)]
obs_end = season_dict["end{:d}".format(season_num)]
if season_num == 0:
ax_start = 0
ax_end = int(obs_length)
else:
ax_start = ax_end
ax_end = ax_start + int(obs_length)
ax = plt.subplot(gs[0, ax_start:ax_end])
obs_this_season = np.logical_and(time >= obs_start,
time <= obs_end)
ax.errorbar(time[obs_this_season] - 54900.5,
mean_mag - mag[obs_this_season],
yerr=model_unc[obs_this_season],
fmt="o", ms=5, mfc=color_dict[ptf_filt],
ecolor=color_dict[ptf_filt], elinewidth=0.5,
mec="0.2", mew=0.2)
ax.xaxis.set_major_locator(MultipleLocator(25))
ax.xaxis.set_minor_locator(MultipleLocator(5))
ax.set_ylim(ymin,ymax)
ax.yaxis.set_major_locator(MultipleLocator(ymajor))
ax.yaxis.set_minor_locator(MultipleLocator(yminor))
if season_num != 0:
ax.set_yticklabels([])
ax_phase = plt.subplot(gs[1:, :])
for repeat in [-1, 0, 1]:
ax_phase.errorbar((time/period) % 1 + repeat,
mean_mag - mag,
yerr=model_unc,
fmt="o", ms=5, mfc=color_dict[ptf_filt],
ecolor=color_dict[ptf_filt], elinewidth=0.5,
mec="0.2", mew=0.2)
ax_phase.set_xlabel("Phase")
ax_phase.set_ylabel("brightness (mag)")
ax_phase.axvline(x=0, ls='--', color='0.8', lw=1)
ax_phase.axvline(x=1, ls='--', color='0.8', lw=1)
ax_phase.set_xlim(-0.2,1.2)
ax_phase.set_ylim(ymin,ymax)
ax_phase.yaxis.set_major_locator(MultipleLocator(ymajor))
ax_phase.yaxis.set_minor_locator(MultipleLocator(yminor))
fig.tight_layout()
fig.show()
In [175]:
aam_plot(hjd[~clipped_obs], mag[~clipped_obs], mag_unc[~clipped_obs],
period=period)
In [176]:
filename = input("Enter .tbl filename: ")
ipac_lc = Table.read(filename, format='ipac')
hjd = np.array(ipac_lc["obsmjd"])
mag = np.array(ipac_lc["mag_autocorr"])
mag_unc = np.array(ipac_lc["magerr_auto"])
clipped_obs = sc(mag).mask
cs_feats = featurize.featurize_time_series(hjd[~clipped_obs],
mag[~clipped_obs],
mag_unc[~clipped_obs],
features_to_use=["freq1_freq",
"amplitude",
"freq1_amplitude1"])
period = 1/float(cs_feats["freq1_freq"][0])
In [177]:
aam_plot(hjd[~clipped_obs], mag[~clipped_obs], mag_unc[~clipped_obs],
period=period)
In [178]:
period
Out[178]:
In [ ]: