Sample code of using the evq_lc_info.fits to access EVQ light curves, and reproduce some plots used in the manuscript

In [1]:
table_location='/Users/yl/Documents/R_Yue/evq_lc_info.fits'
In [2]:
import os
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from astropy.io import fits
from astropy.table import Table
In [3]:
plt.rcParams["font.family"] = "Times New Roman"
plt.rcParams["font.weight"] = 'light'
params = {'mathtext.default': 'regular' }          
plt.rcParams.update(params)
In [4]:
# read in the fits table
table=Table.read(table_location)
In [5]:
# column names of the table
table.colnames
Out[5]:
['SDSS_DR7_ID',
 'median_data_each_season',
 'std_wrt_line_fitting',
 'redshift',
 'RA',
 'DEC',
 'rms_vari_arr',
 'Bezier_curve',
 'num_of_inflection_pts',
 'category',
 'median_err_each_season',
 'median_rms_vari(excluding_0)',
 'num_of_pts_each_season',
 'complete_light_curve',
 'survey_of_each_data_point',
 'total_absolute_curvature',
 'inflection points',
 'control points']
In [6]:
# header of the fits table
fits.getheader(table_location)
Out[6]:
SIMPLE  =                    T / conforms to FITS standard                      
BITPIX  =                    8 / array data type                                
NAXIS   =                    0 / number of array dimensions                     
EXTEND  =                    T                                                  
COLUMN1 = 'SDSS_DR7_ID'                                                         
COLUMN2 = 'median_data_each_season' / array of (mjd, mag)                       
COLUMN3 = 'std_wrt_line_fitting'                                                
COLUMN4 = 'redshift'                                                            
COLUMN5 = 'RA      '                                                            
COLUMN6 = 'DEC     '                                                            
COLUMN7 = 'rms_vari_arr'       / array of (median mag, rms vari in the group)   
COLUMN8 = 'Bezier_curve'       / array (mjd, magnitude)                         
COLUMN9 = 'num_of_inflection_pts'                                               
COLUMN10= 'category'                                                            
COLUMN11= 'median_err_each_season'                                              
COLUMN12= 'median_rms_vari(excluding_0)'                                        
COLUMN13= 'num_of_pts_each_season'                                              
COLUMN14= 'complete_light_curve' / array of (mjd, mag, err)                     
COLUMN15= 'survey_of_each_data_point'                                           
COLUMN16= 'total_absolute_curvature'                                            
COLUMN17= 'infection points'   / array (mjd, magnitude)                         
COLUMN18= 'control points'     / control points used to construc the Bezier curv

plot the light curve for one of the EVQs

In [7]:
# ID number of the object in SDSS_DR7 as a string
idnum='28'

# get its index in the table
ids=np.asarray(table['SDSS_DR7_ID'],dtype=str)
index=np.where(ids==idnum)[0][0]

# access complete light curve data, in the form of (mjd, g-magnitude, magitude err)
lc=table['complete_light_curve'][index]
# strip all the nan
lc=lc[~np.isnan(lc).any(axis=1)]

# get the survey of each data point
survey=np.asarray(table['survey_of_each_data_point'][index],dtype=str)
# strip all the "nan"
survey=survey[~(survey=='nan')]

# access the smoothed representation of the light curve
bc=table['Bezier_curve'][index]
# strip all the nan
bc=bc[~np.isnan(bc).any(axis=1)]

# obtain the median point for each season
lc_mid=table['median_data_each_season'][index]
# strip all the nan
lc_mid=lc_mid[~np.isnan(lc_mid).any(axis=1)]

# obtain index of data points of different surveys, so we can plot them with different markers
sdss_index=np.where(survey=='SDSS')[0]
ps_index=np.where(survey=='PanStarrs')[0]
des_index=np.where(survey=='DES')[0]

# obtain the identified inflection points
inf=table['inflection points'][index]
# strip all the nan
inf=inf[~np.isnan(inf).any(axis=1)]

# obtain ra, dec, and z
ra=table['RA'][index]
dec=table['DEC'][index]
z=table['redshift'][index]
In [8]:
# plot the light curve
fig=plt.figure(figsize=(12,4))
plt.gca().invert_yaxis()

plt.errorbar(lc[sdss_index][:,0],lc[sdss_index][:,1],yerr=lc[sdss_index][:,2],elinewidth=1,fmt='o',ecolor='b',mfc='none', mec='k', ms=8,label='SDSS')
plt.errorbar(lc[ps_index][:,0],lc[ps_index][:,1],yerr=lc[ps_index][:,2],elinewidth=1,fmt='^',ecolor='b',mfc='none', mec='r', ms=9,label='PanStarrs')
plt.errorbar(lc[des_index][:,0],lc[des_index][:,1],yerr=lc[des_index][:,2],elinewidth=1,fmt='s',ecolor='b',mfc='none', mec='g', ms=9,label='DES')

plt.plot(lc_mid[:,0],lc_mid[:,1],'o',ms=13,color='#32CD32',label='season median',alpha=0.7)
plt.plot(bc[:,0],bc[:,1],color='b',alpha=0.8,label=r'$B\acute{e}zier$ curve')
plt.plot(inf[:,0],inf[:,1],'*',color='r', ms=17,label='inflection points')

plt.title('z={},  SDSS_DR7_ID={},  RA={:.5f},  DEC= {:.5f},  $\Delta$g={:.2f},  g$_{{median}}$={:.2f}'.format(z,idnum,ra,dec,np.max(lc[:,1])-np.min(lc[:,1]),np.median(lc[:,1])),
          y=0.999,fontsize=16.6)
plt.legend(prop={'size': 15.5},bbox_to_anchor=(1., 0.49),loc='center left')
plt.tick_params(axis='both', which='major', labelsize=17)
plt.xlabel('MJD',fontsize=16)
plt.ylabel('g-band magnitude',fontsize=18)
Out[8]:
Text(0, 0.5, 'g-band magnitude')

plot rms magnitude vs redshift for all EVQs with categories

In [9]:
# obtain index for different category, so we can plot them with different markers
category=np.asarray(table['category'],dtype=str)
mono_index=np.where(category=='monotonic')[0]
com_index=np.where(category=='complex')[0]
pv_index=np.where(category=='peaked')[0]

# obtain redshift
redshift=np.asarray(table['redshift'],dtype=float)

# obtain rms magnitude
rms_mag=np.asarray(table['median_rms_vari(excluding_0)'],dtype=float)
In [10]:
fig = plt.figure(figsize=(8, 8))
grid = plt.GridSpec(4, 4, hspace=0.3, wspace=0.55)
main_ax = fig.add_subplot(grid[:-1, 1:])

y_hist = fig.add_subplot(grid[:-1, 0], sharey=main_ax)
x_hist = fig.add_subplot(grid[-1, 1:],sharex=main_ax)
fig.text(0.6, 0.06, 'redshift', ha='center', va='center',fontsize=20.7)
fig.text(0.01, 0.6, '$\sqrt{g^2_{std}-g^2_{error}}}$', ha='center', va='center', rotation='vertical',fontsize=20.7)

# make convenient names for quantities to be plotted on the x and y axes
xpv=redshift[pv_index]
ypv=rms_mag[pv_index]
xcom=redshift[com_index]
ycom=rms_mag[com_index]
xm=redshift[mono_index]
ym=rms_mag[mono_index]

# scatter points on the main axes
main_ax.plot(xpv, ypv,'o',color='#A9A9A9', ms=7,label='peak/valley')
main_ax.plot(xcom, ycom,'^',color='#1E90FF', ms=8,alpha=0.7,label='complex')
main_ax.plot(xm, ym,'s',color='#FF4500', ms=7,alpha=0.7,label='monotonic')
main_ax.legend(prop={'size': 16},loc=0)
main_ax.tick_params(axis='both', which='major', labelsize=17)

# calculate weights for normalizing the histogram
wpvx = np.ones_like(xpv)/float(len(xpv))
wmx = np.ones_like(xm)/float(len(xm))
wcomx = np.ones_like(xcom)/float(len(xcom))

# histogram on the attached axes
bins1=np.linspace(min(np.r_[xm,xpv,xcom]),max(np.r_[xm,xpv,xcom]),30)
x_hist.hist(xpv, bins1, histtype='step',facecolor='None',linewidth=1.4,weights=wpvx,
            orientation='vertical', edgecolor='#808080')
x_hist.hist(xcom, bins1, histtype='step',facecolor='None',linewidth=1.4,weights=wcomx,
            orientation='vertical', edgecolor='#1E90FF',alpha=0.9)
x_hist.hist(xm, bins1, histtype='step',facecolor='None',linewidth=1.4,weights=wmx,
            orientation='vertical', edgecolor='#FF4500',alpha=0.9)
x_hist.tick_params(axis='both', which='major', labelsize=17)

# calculate weights for normalizing the histogram
# get rid of nan when necessary
wpvy = np.ones_like(ypv[np.where(~np.isnan(ypv))])/float(len(ypv[np.where(~np.isnan(ypv))]))
wmy = np.ones_like(ym[np.where(~np.isnan(ym))])/float(len(ym[np.where(~np.isnan(ym))]))
wcomy = np.ones_like(ycom)/float(len(ycom))

bins2=np.linspace(np.nanmin(np.array(np.r_[ym,ypv,ycom])),np.nanmax(np.array(np.r_[ym,ypv,ycom])),30)
y_hist.hist(ypv[np.where(~np.isnan(ypv))], bins2, histtype='step',linewidth=1.4,
            orientation='horizontal', edgecolor='#808080',weights=wpvy)
y_hist.hist(ycom, bins2, histtype='step',linewidth=1.4,weights=wcomy,
            orientation='horizontal', edgecolor='#1E90FF',alpha=0.9)
y_hist.hist(ym[np.where(~np.isnan(ym))], bins2, histtype='step',linewidth=1.4,
            orientation='horizontal', edgecolor='#FF4500',alpha=0.9,weights=wmy)
y_hist.invert_xaxis()
y_hist.tick_params(axis='both', which='major', labelsize=17)