WISE_ZTF_download_plot
Published:
Import packages
from astroquery.ipac.ned import Ned
import ztfquery
from ztfquery import lightcurve
import numpy as np
import pandas as pd
import os,sys
from astropy.table import Table
from astropy import config as _config
from astropy import units as u
from astropy.coordinates import SkyCoord
from astroquery.irsa import Irsa
from scipy.stats import spearmanr#
from scipy.stats.stats import pearsonr
from astropy.table import Table
from astropy.cosmology import FlatLambdaCDM,Planck13,Planck15,z_at_value
from astropy import units as u
from astropy.cosmology import LambdaCDM
cosmo = LambdaCDM(H0=70, Om0=0.27, Ode0=0.73)
from astropy.time import Time
from astropy.time import TimeYearDayTime
from datetime import datetime
import time
from time import strftime,strptime
import calendar
from dateutil.parser import parse
#from adjustText import adjust_text
import matplotlib as mpl
import matplotlib.pyplot as plt
from pylab import cm
from collections import OrderedDict
#from adjustText import adjust_text
from astropy.table import Table
from astropy.io import fits
from astropy import units as u
from astropy import constants
from astropy.coordinates import SkyCoord
from astropy.visualization import hist
from astroML.datasets import fetch_imaging_sample, fetch_sdss_S82standards
from astroML.crossmatch import crossmatch_angular
%matplotlib inline
%config InlineBackend.figure_format='svg'
# Edit the font, font size, and axes width
mpl.rcParams['font.family'] = 'Times New Roman' #'Avenir'
plt.rcParams['font.size'] = 18
plt.rcParams['axes.linewidth'] = 1.5
#from adjustText import adjust_text
import matplotlib as mpl
def set_ax_tick(ax):
ax.xaxis.set_tick_params(which='major', size=10, width=2, direction='in', top='on',)
ax.xaxis.set_tick_params(which='minor', size=5, width=2, direction='in', top='on')
ax.yaxis.set_tick_params(which='major', size=10, width=2, direction='in', right='on')
ax.yaxis.set_tick_params(which='minor', size=5, width=2, direction='in', right='on')
def set_ax_locator(ax,xma=1,xmi=0.2,yma=1,ymi=0.2):
ax.xaxis.set_major_locator(mpl.ticker.MultipleLocator(xma))
ax.xaxis.set_minor_locator(mpl.ticker.MultipleLocator(xmi))
ax.yaxis.set_major_locator(mpl.ticker.MultipleLocator(yma))
ax.yaxis.set_minor_locator(mpl.ticker.MultipleLocator(ymi))
def plot_secax(ax,mi_interval=365,ma_interval=365*2,rotation=30,):
secax1 = ax.secondary_xaxis('top', functions=(mjd2numsdate,datenums2mjd))
secax1.xaxis.set_major_formatter(mdates.DateFormatter("%Y-%m"))
secax1.xaxis.set_major_locator(mdates.DayLocator(interval=ma_interval))
secax1.xaxis.set_minor_locator(mdates.DayLocator(interval=mi_interval))
secax1.xaxis.set_tick_params(which='major', size=10, width=2, direction='out')
secax1.xaxis.set_tick_params(which='minor', size=5, width=2, direction='out')
secax1.xaxis.set_tick_params(which='major', size=10, width=2, direction='out', rotation=rotation,)
def set_ax_legend(ax,bbox_to_anchor=(0.01, 0.99)):
ax.xaxis.set_tick_params(which='major', size=10, width=2, rotation=0,)
handles, labels = ax.get_legend_handles_labels()
# remove the errorbars
#hdl = [h[0] for h in handles]
hdl = handles
labels_dict=dict(zip(labels, hdl)) #key,values
by_label=OrderedDict(sorted(labels_dict.items(),key=lambda t:t[0]))
#by_label = OrderedDict(zip(labels, handles))
ax.legend(by_label.values(), by_label.keys(), bbox_to_anchor=bbox_to_anchor,
loc=2, numpoints=1,ncol=1,fontsize=11.)
def set_ax_legend_sequence(ax,bbox_to_anchor=(0.01, 0.99)):
ax.xaxis.set_tick_params(which='major', size=10, width=2, rotation=0,)
handles, labels = ax.get_legend_handles_labels()
# remove the errorbars
#hdl = [h[0] for h in handles]
hdl = handles
labels_dict=dict(zip(labels, hdl)) #key,values
#by_label=OrderedDict(sorted(labels_dict.items(),key=lambda t:t[0]))
by_label = OrderedDict(zip(labels, handles))
ax.legend(by_label.values(), by_label.keys(), bbox_to_anchor=bbox_to_anchor,
loc=2, numpoints=1,ncol=1,fontsize=11.)
def set_mag_ylim(ax):
bottom, top = ax.set_ylim()
if bottom< top:
ax.set_ylim(top,bottom)
def set_mag_xlim(ax):
bottom, top = ax.set_xlim()
if bottom< top:
ax.set_xlim(top,bottom)
def drop_index(data):
data=data.reset_index(drop=True)
return data
#import pandas as pd
#import numpy as np
#import os
def get_obsids(path):
dirname=os.listdir(path)
obsids=[]
for i in dirname:
if i.isdigit():
obsids.append(i)
obsids.sort()
return obsids
**Time_convert **
from astropy.time import Time
from astropy.io import fits
import time
import datetime
from matplotlib.pyplot import MultipleLocator
import matplotlib.dates as mdates
def datetime2mjd(x):
mjd_ref=59000
mjd_minus_mdates_num=mdates.date2num(convert_xaxis_time(mjd_ref))-mjd_ref
x=mdates.date2num(x)
y = x - mjd_minus_mdates_num
return y
def mjd2datetime(x):
mjd_ref=59000
mjd_minus_mdates_num=mdates.date2num(convert_xaxis_time(mjd_ref))-mjd_ref
y= x + mjd_minus_mdates_num
y= mdates.num2date(y)
return y
def datenums2mjd(x):
#x=mdates.date2num(x)
mjd_ref=59000
mjd_minus_mdates_num=mdates.date2num(convert_xaxis_time(mjd_ref))-mjd_ref
y = x - mjd_minus_mdates_num
return y
def mjd2numsdate(x):
mjd_ref=59000
mjd_minus_mdates_num=mdates.date2num(convert_xaxis_time(mjd_ref))-mjd_ref
y= x + mjd_minus_mdates_num
#y= mdates.num2date(y)
return y
def convert_xaxis_mjd(time):
return Time(time).mjd
def convert_xaxis_time(mjd):
return Time(mjd,format='mjd').to_datetime()
def date2yday(x):
"""
x is in matplotlib datenums, so they are floats.
"""
y = x - mdates.date2num(datetime(2018, 1, 1))
return y
def yday2date(x):
"""
return a matplotlib datenum (x is days since start of year of 2018)
"""
y = x + mdates.date2num(datetime(2018, 1, 1))
return y
def convert_partial_year(numbers):
datetimes=[]
for number in numbers:
year = int(number)
d = timedelta(days=(number - year)*(365 + is_leap(year)))
day_one = datetime(year,1,1)
date = d + day_one
datetimes.append(date)
return datetimes
def is_leap(year):
if not year%4 and year%100 or not year%400:
return True
return False
def convert_mjd(times):
timesmjd=[]
for i in times:
timesmjd.append(Time(i).mjd)
return timesmjd
def convert_date(times):
timesdate=[]
for i in times:
timesdate.append(Time(i,format='mjd').datetime)
return timesdate
def convert_date_single(time):
timedate=Time(time,format='mjd').datetime
return timedate
**def get_mag_weighted_mean_mag **
import numpy as np
def get_mag_weighted_mean_mag(mags,e_mags):
N=len(mags)
mag_weighted_mean_mag =0
w_n=np.sum(1.0/e_mags**2)
for mag,e_mag in zip(mags,e_mags):
w_i=1.0/e_mag**2
w=w_i/w_n
mag_weighted_mean_mag=mag_weighted_mean_mag+w*mag
return mag_weighted_mean_mag
def get_mag_weighted_mean_e_mag(mags,e_mags):
N=len(mags)
w_n=np.sum(1.0/e_mags**2)
mag_weighted_mean_mag = get_mag_weighted_mean_mag(mags,e_mags)
mag_weighted_mean_e_mag =0
for mag,e_mag in zip(mags,e_mags):
w_i=1.0/e_mag**2
w=w_i/w_n
mag_weighted_mean_e_mag=mag_weighted_mean_e_mag+w*(mag-mag_weighted_mean_mag)**2
mag_weighted_mean_e_mag_=(mag_weighted_mean_e_mag/(N-1))**0.5
return mag_weighted_mean_e_mag_
**def median_lc **
def median_lc(data,timelabel,maglabel,e_maglabel,interval=180,sigma_outlier=100):
#data=data.sort_values(timelabel)
data=drop_index(data)
mean_mjds=[]
mean_mags=[]
mean_mags_e=[]
if len(data)>2:
mjd_t0=data[timelabel][0]-10
mjd_tend=data[timelabel][len(data)-1]+10
N=int((mjd_tend-mjd_t0)//interval+1)
#print(N)
time_interval=np.linspace(mjd_t0,interval*N+mjd_t0,N+1)
#print(time_interval)
for index,i in enumerate(time_interval):
if index < N:
data_select=data[data[timelabel]>=time_interval[index]]
data_select=data_select[data_select[timelabel]<time_interval[index+1]]
data_select=drop_index(data_select)
if len(data_select)>2:
#mean_mjd=np.mean(data_select[timelabel])
mags=data_select[maglabel]
e_mags=data_select[e_maglabel]
mean_mag=get_mag_weighted_mean_mag(mags,e_mags)
mag_weighted_mean_e_mag=get_mag_weighted_mean_e_mag(mags,e_mags)
#mean_mag_e=np.mean(data_select[e_maglabel])
data_select=data_select[abs(data_select[maglabel]-mean_mag)<sigma_outlier*mag_weighted_mean_e_mag]
data_select=drop_index(data_select)
#print(len(data_select))
mags=data_select[maglabel]
e_mags=data_select[e_maglabel]
if len(data_select)>2:
mean_mag=get_mag_weighted_mean_mag(mags,e_mags)
mean_mag_e=get_mag_weighted_mean_e_mag(mags,e_mags)
mean_mjd=np.mean(data_select[timelabel])
mean_mjds.append(mean_mjd)
mean_mags.append(mean_mag)
mean_mags_e.append(mean_mag_e)
mean_mjds,mean_mags,mean_mags_e=(np.array(mean_mjds),np.array(mean_mags),np.array(mean_mags_e))
return mean_mjds,mean_mags,mean_mags_e
**def get_outlier_remove_data **
def get_outlier_remove_data(data,timelabel,maglabel,e_maglabel,interval=180,sigma_outlier=10):
data_remove=pd.DataFrame([])
#data_outlier=pd.DataFrame([])
#data=data.sort_values(timelabel)
data=drop_index(data)
mean_mjds=[]
mean_mags=[]
mean_mags_e=[]
if len(data)>2:
mjd_t0=data[timelabel][0]-20
mjd_tend=data[timelabel][len(data)-1]+20
N=int((mjd_tend-mjd_t0)//interval+1)
#print(N)
time_interval=np.linspace(mjd_t0,interval*N+mjd_t0,N+1)
#print(time_interval)
for index,i in enumerate(time_interval):
if index < N:
data_select=data[data[timelabel]>=time_interval[index]]
data_select=data_select[data_select[timelabel]<time_interval[index+1]]
data_select=drop_index(data_select)
if len(data_select)>2:
#mean_mjd=np.mean(data_select[timelabel])
mags=data_select[maglabel]
e_mags=data_select[e_maglabel]
mean_mag=get_mag_weighted_mean_mag(mags,e_mags)
mag_weighted_mean_e_mag=get_mag_weighted_mean_e_mag(mags,e_mags)
#mean_mag_e=np.mean(data_select[e_maglabel])
data_remove_outlier=data_select[abs(data_select[maglabel]-mean_mag)<=sigma_outlier*mag_weighted_mean_e_mag]
data_remove_outlier=drop_index(data_remove_outlier)
data_remove=data_remove.append(data_remove_outlier,ignore_index=True)
'''
data_out=data_select[abs(data_select[maglabel]-mean_mag)>sigma_outlier*mag_weighted_mean_e_mag]
data_out=drop_index(data_out)
data_outlier.append(data_out,ignore_index=True)
'''
#data_remove=data_remove.sort_values(timelabel)
data_remove=drop_index(data_remove)
return data_remove#,data_outlier
def get_outlier_data(data,timelabel,maglabel,e_maglabel,interval=180,sigma_outlier=10):
data_remove=pd.DataFrame([])
data_outlier=pd.DataFrame([])
#data=data.sort_values(timelabel)
data=drop_index(data)
mean_mjds=[]
mean_mags=[]
mean_mags_e=[]
if len(data)>2:
mjd_t0=data[timelabel][0]-20
mjd_tend=data[timelabel][len(data)-1]+20
N=int((mjd_tend-mjd_t0)//interval+1)
#print(N)
time_interval=np.linspace(mjd_t0,interval*N+mjd_t0,N+1)
#print(time_interval)
for index,i in enumerate(time_interval):
if index < N:
data_select=data[data[timelabel]>=time_interval[index]]
data_select=data_select[data_select[timelabel]<time_interval[index+1]]
data_select=drop_index(data_select)
if len(data_select)>2:
#mean_mjd=np.mean(data_select[timelabel])
mags=data_select[maglabel]
e_mags=data_select[e_maglabel]
mean_mag=get_mag_weighted_mean_mag(mags,e_mags)
mag_weighted_mean_e_mag=get_mag_weighted_mean_e_mag(mags,e_mags)
#mean_mag_e=np.mean(data_select[e_maglabel])
'''data_remove_outlier=data_select[abs(data_select[maglabel]-mean_mag)<=sigma_outlier*mag_weighted_mean_e_mag]
data_remove_outlier=drop_index(data_remove_outlier)
data_remove=data_remove.append(data_remove_outlier,ignore_index=True) '''
data_out=data_select[abs(data_select[maglabel]-mean_mag)>sigma_outlier*mag_weighted_mean_e_mag]
data_out=drop_index(data_out)
data_outlier=data_outlier.append(data_out,ignore_index=True)
#data_outlier=data_outlier.sort_values(timelabel)
#data_outlier=drop_index(data_outlier)
return data_outlier#,data_outlier
def get_remove_neo_mep(data_test_neo,data_test_mep):
#neo
if len(data_test_neo)>0:
timelabel,maglabel,e_maglabel=('mjd','w1mpro','w1sigmpro')
data_test_neo_remove=get_outlier_remove_data(data_test_neo,timelabel,maglabel,e_maglabel,interval=180,sigma_outlier=10)
timelabel,maglabel,e_maglabel=('mjd','w2mpro','w2sigmpro')
data_test_neo_remove=get_outlier_remove_data(data_test_neo_remove,timelabel,maglabel,e_maglabel,interval=180,sigma_outlier=10)
else:
data_test_neo_remove=[]
#mep
if len(data_test_mep)>0:
timelabel,maglabel,e_maglabel=('mjd','w1mpro_ep','w1sigmpro_ep')
data_test_mep_remove=get_outlier_remove_data(data_test_mep,timelabel,maglabel,e_maglabel,interval=180,sigma_outlier=10)
timelabel,maglabel,e_maglabel=('mjd','w2mpro_ep','w2sigmpro_ep')
data_test_mep_remove=get_outlier_remove_data(data_test_mep_remove,timelabel,maglabel,e_maglabel,interval=180,sigma_outlier=10)
else:
data_test_mep_remove=[]
return data_test_neo_remove,data_test_mep_remove
def get_outlier_neo_mep(data_test_neo,data_test_mep):
#neo
if len(data_test_neo)>0:
timelabel,maglabel,e_maglabel=('mjd','w1mpro','w1sigmpro')
data_test_neo_outiler=get_outlier_data(data_test_neo,timelabel,maglabel,e_maglabel,interval=180,sigma_outlier=10)
timelabel,maglabel,e_maglabel=('mjd','w2mpro','w2sigmpro')
data_test_neo_outiler=get_outlier_data(data_test_neo_outiler,timelabel,maglabel,e_maglabel,interval=180,sigma_outlier=10)
else:
data_test_neo_outiler=[]
#mep
if len(data_test_mep)>0:
timelabel,maglabel,e_maglabel=('mjd','w1mpro_ep','w1sigmpro_ep')
data_test_mep_outiler=get_outlier_data(data_test_mep,timelabel,maglabel,e_maglabel,interval=180,sigma_outlier=10)
timelabel,maglabel,e_maglabel=('mjd','w2mpro_ep','w2sigmpro_ep')
data_test_mep_outiler=get_outlier_data(data_test_mep_outiler,timelabel,maglabel,e_maglabel,interval=180,sigma_outlier=10)
else:
data_test_mep_outiler=[]
return data_test_neo_outiler,data_test_mep_outiler
def plot_outlier_median_lc(ax,data,timelabel,maglabel,e_maglabel,interval=180,sigma_outlier=10):
#data=data.sort_values(timelabel)
data=drop_index(data)
mean_mjds=[]
mean_mags=[]
mean_mags_e=[]
if len(data)>2:
mjd_t0=data[timelabel][0]-20
mjd_tend=data[timelabel][len(data)-1]+20
N=int((mjd_tend-mjd_t0)//interval+1)
#print(N)
time_interval=np.linspace(mjd_t0,interval*N+mjd_t0,N+1)
#print(time_interval)
for index,i in enumerate(time_interval):
if index < N:
data_select=data[data[timelabel]>=time_interval[index]]
data_select=data_select[data_select[timelabel]<time_interval[index+1]]
data_select=drop_index(data_select)
if len(data_select)>2:
#mean_mjd=np.mean(data_select[timelabel])
mags=data_select[maglabel]
e_mags=data_select[e_maglabel]
mean_mag=get_mag_weighted_mean_mag(mags,e_mags)
mag_weighted_mean_e_mag=get_mag_weighted_mean_e_mag(mags,e_mags)
#mean_mag_e=np.mean(data_select[e_maglabel])
data_select_outlier=data_select[abs(data_select[maglabel]-mean_mag)>=sigma_outlier*mag_weighted_mean_e_mag]
data_select_outlier=drop_index(data_select_outlier)
if len(data_select_outlier)>0:
ax.scatter(data_select_outlier[timelabel],data_select_outlier[maglabel],
color='red',zorder=3,s=6)
**def get_wise_data **
def get_wise_data(name_test,datawise_dir=CLAGN_wisedata_dir):
#print(name_test)
search_name=name_test.replace(' ','')
save_name=search_name.strip()
if os.path.exists(os.path.join(datawise_dir,'wise_%s_neo.csv'%save_name)):
data_test_neo=pd.read_csv(os.path.join(datawise_dir,'wise_%s_neo.csv'%save_name))
data_test_neo=data_test_neo[data_test_neo['w1mpro'].notnull()]
data_test_neo=data_test_neo[data_test_neo['w2mpro'].notnull()]
data_test_neo=data_test_neo[data_test_neo['w1sigmpro'].notnull()]
data_test_neo=data_test_neo[data_test_neo['w2sigmpro'].notnull()]
data_test_neo['W1-W2']=data_test_neo['w1mpro']-data_test_neo['w2mpro']
#data_test_neo=data_test_neo[data_test_neo['w1mpro']<15]
#data_test_neo=data_test_neo[data_test_neo['w2mpro']<13]
#data_test_neo=data_test_neo[data_test_neo['w1rchi2']<2]
#data_test_neo=data_test_neo[data_test_neo['w2rchi2']<2]
#data_test_neo=data_test_neo[data_test_neo['cc_flags']=="b'0000'"]
#data_test_neo=data_test_neo[data_test_neo['nb']<3]
#data_test_neo=data_test_neo[data_test_neo['na']==0]
#data_test_neo=data_test_neo[data_test_neo['w1snr']>10]
#data_test_neo=data_test_neo[data_test_neo['w2snr']>10]
data_test_neo=data_test_neo[data_test_neo['qual_frame']>0]
#print(len(data_test_neo))
data_test_neo=(data_test_neo.loc[:, ['mjd', 'w1mpro','w2mpro',
'w1sigmpro','w2sigmpro','W1-W2',
]].sort_values('mjd'))
data_test_neo=data_test_neo.reset_index(drop=True)
#print(data_test_neo.columns)
else:
data_test_neo=[]
if os.path.exists(os.path.join(datawise_dir,'wise_%s_mep.csv'%save_name)):
data_test_mep=pd.read_csv(os.path.join(datawise_dir,'wise_%s_mep.csv'%save_name))
data_test_mep=data_test_mep[data_test_mep['w1mpro_ep'].notnull()]
data_test_mep=data_test_mep[data_test_mep['w2mpro_ep'].notnull()]
data_test_mep=data_test_mep[data_test_mep['w1sigmpro_ep'].notnull()]
data_test_mep=data_test_mep[data_test_mep['w2sigmpro_ep'].notnull()]
data_test_mep['W1-W2']=data_test_mep['w1mpro_ep']-data_test_mep['w2mpro_ep']
#data_test_mep=data_test_mep[data_test_mep['qual_frame']>0]
data_test_mep=data_test_mep[data_test_mep['qi_fact']>0]
data_test_mep=data_test_mep[data_test_mep['saa_sep']>0]
#data_test_mep=data_test_mep[data_test_mep['cc_flags']=="b'0000'"]
#data_test_mep=data_test_mep[data_test_mep['nb']<3]
#data_test_mep=data_test_mep[data_test_mep['na']==0]
#print(len(data_test_mep))
data_test_mep=(data_test_mep.loc[:, ['mjd', 'w1mpro_ep','w2mpro_ep',
'w1sigmpro_ep','w2sigmpro_ep',
'W1-W2',
]].sort_values('mjd'))
data_test_mep=data_test_mep.reset_index(drop=True)
#print(data_test_mep.columns)
#if not os.path.exists(save_lc_img_path):
else:
data_test_mep=[]
return data_test_neo,data_test_mep
def get_wise_radec_data(radeg,decdeg,datawise_dir=CLAGN_wisedata_dir):
#print(name_test)
save_name='RA%.5f_DEC%.5f'%(radeg,decdeg)
#name=save_name
data_test_neo=pd.read_csv(os.path.join(datawise_dir,'wise_%s_neo.csv'%save_name))
data_test_neo=data_test_neo[data_test_neo['w1mpro'].notnull()]
data_test_neo=data_test_neo[data_test_neo['w2mpro'].notnull()]
data_test_neo=data_test_neo[data_test_neo['w1sigmpro'].notnull()]
data_test_neo=data_test_neo[data_test_neo['w2sigmpro'].notnull()]
data_test_neo['W1-W2']=data_test_neo['w1mpro']-data_test_neo['w2mpro']
#data_test_neo=data_test_neo[data_test_neo['w1mpro']<15]
#data_test_neo=data_test_neo[data_test_neo['w2mpro']<13]
#data_test_neo=data_test_neo[data_test_neo['w1rchi2']<2]
#data_test_neo=data_test_neo[data_test_neo['w2rchi2']<2]
#data_test_neo=data_test_neo[data_test_neo['cc_flags']=="b'0000'"]
#data_test_neo=data_test_neo[data_test_neo['cc_flags']=='0000']
#data_test_neo=data_test_neo[data_test_neo['nb']<3]
#data_test_neo=data_test_neo[data_test_neo['na']==0]
data_test_neo=data_test_neo[data_test_neo['w1snr']>3]
data_test_neo=data_test_neo[data_test_neo['w2snr']>3]
#data_test_neo=data_test_neo[data_test_neo['qual_frame']>0]
#print(len(data_test_neo))
data_test_neo=(data_test_neo.loc[:, ['mjd', 'w1mpro','w2mpro',
'w1sigmpro','w2sigmpro','W1-W2',
]].sort_values('mjd'))
data_test_neo=data_test_neo.reset_index(drop=True)
#print(data_test_neo.columns)
data_test_mep=pd.read_csv(os.path.join(datawise_dir,'wise_%s_mep.csv'%save_name))
data_test_mep=data_test_mep[data_test_mep['w1mpro_ep'].notnull()]
data_test_mep=data_test_mep[data_test_mep['w2mpro_ep'].notnull()]
data_test_mep=data_test_mep[data_test_mep['w1sigmpro_ep'].notnull()]
data_test_mep=data_test_mep[data_test_mep['w2sigmpro_ep'].notnull()]
data_test_mep['W1-W2']=data_test_mep['w1mpro_ep']-data_test_mep['w2mpro_ep']
#data_test_mep=data_test_mep[data_test_mep['qual_frame']>5]
data_test_mep=data_test_mep[data_test_mep['qi_fact']>0]
data_test_mep=data_test_mep[data_test_mep['saa_sep']>0]
#data_test_mep=data_test_mep[data_test_mep['cc_flags']=="b'0000'"]
#data_test_mep=data_test_mep[data_test_mep['cc_flags']==0]
#data_test_mep=data_test_mep[data_test_mep['nb']<3]
#data_test_mep=data_test_mep[data_test_mep['na']==0]
#print(len(data_test_mep))
data_test_mep=(data_test_mep.loc[:, ['mjd', 'w1mpro_ep','w2mpro_ep',
'w1sigmpro_ep','w2sigmpro_ep',
'W1-W2',
]].sort_values('mjd'))
data_test_mep=data_test_mep.reset_index(drop=True)
#print(data_test_mep.columns)
#if not os.path.exists(save_lc_img_path):
return data_test_neo,data_test_mep
def get_M1minM2_data(data_test_neo,data_test_mep,neo_maglabel_1='w1mpro',neo_maglabel_2='w2mpro',\
mep_maglabel_1='w1mpro_ep',mep_maglabel_2='w2mpro_ep'):
data_test_neo['delta_M']=data_test_neo[neo_maglabel_1]-data_test_neo[neo_maglabel_2]
data_test_mep['delta_M']=data_test_mep[neo_maglabel_1]-data_test_mep[neo_maglabel_2]
return data_test_neo['mjd','delta_M'],data_test_mep['mjd','delta_M']
epsilon_s_w1=0.024
epsilon_s_w2=0.028
def get_intrinsic_var(data,mag_label,magerr_label,epsilon_s):
mags=data[mag_label]
e_mags=data[magerr_label]
mag_weighted_mean_mag = get_mag_weighted_mean_mag(mags,e_mags)
mag_weighted_mean_e_mag_= get_mag_weighted_mean_e_mag(mags,e_mags)
mag_mean=mag_weighted_mean_mag
length_N=len(data[mag_label])
epsilon_square=np.average(data[magerr_label]**2)+epsilon_s**2
Sigma_square=np.sum((data[mag_label]-mag_mean)**2)/(length_N-1)
if Sigma_square>epsilon_square:
sigma_m=np.sqrt(Sigma_square-epsilon_square)
else:
sigma_m=0
if length_N<2:
sigma_m=-1#print(length_N)
return sigma_m
def get_intrinsic_var_rebin(data_rebin_neo,data_rebin_mep,epsilon_s):
mags=data[mag_label]
e_mags=data[magerr_label]
mag_weighted_mean_mag = get_mag_weighted_mean_mag(mags,e_mags)
mag_weighted_mean_e_mag_= get_mag_weighted_mean_e_mag(mags,e_mags)
mag_mean=mag_weighted_mean_mag
length_N=len(mags)
epsilon_square=np.average(data[magerr_label]**2)+epsilon_s**2
Sigma_square=np.sum((data[mag_label]-mag_mean)**2)/(length_N-1)
if Sigma_square>epsilon_square:
sigma_m=np.sqrt(Sigma_square-epsilon_square)
else:
sigma_m=0
if length_N<2:
sigma_m=-1#print(length_N)
return sigma_m
def get_minmax_mag_name_one(data_test_neo,data_test_mep):
maxmag=0
minmag=20
maxmag_var=0
if len(data_test_neo)>20:
mean_mjds_neo,mean_mags_neo,mean_mags_e_neo=median_lc(data_test_neo,'mjd','w1mpro','w1sigmpro',interval=180,)
maxmag=max(max(mean_mags_neo),maxmag)
minmag=min(min(mean_mags_neo),minmag)
if len(data_test_mep)>20:
mean_mjds_mep,mean_mags_mep,mean_mags_e_mep=median_lc(data_test_mep,'mjd','w1mpro_ep','w1sigmpro_ep',interval=180,)
maxmag=max(max(mean_mags_mep),maxmag)
minmag=min(min(mean_mags_mep),minmag)
maxmag_var=maxmag-minmag
return maxmag_var
def W1mag_to_Lum(w1,d_cm):
Lum=(-48.6-w1-2.699)/2.5+np.log10(4*np.pi*d_cm**2)+np.log10((constants.c/(3.4*u.um)).to(u.Hz).value)
return Lum #logLum
def W2mag_to_Lum(w2,d_cm):
Lum=(-48.6-w2-3.339)/2.5+np.log10(4*np.pi*d_cm**2)+np.log10((constants.c/(4.6*u.um)).to(u.Hz).value)
return Lum #logLum
def W3mag_to_Lum(w3,d_cm):
Lum=(-48.6-w3-5.174)/2.5+np.log10(4*np.pi*d_cm**2)+np.log10((constants.c/(12*u.um)).to(u.Hz).value)
return Lum #logLum
def W4mag_to_Lum(w4,d_cm):
Lum=(-48.6-w4-6.620)/2.5+np.log10(4*np.pi*d_cm**2)+np.log10((constants.c/(22*u.um)).to(u.Hz).value)
return Lum #logLum
def get_rebinwise_data(data_test_neo,data_test_mep,sigma_outlier=100):
#print(name_test)
#search_name=name_test.replace(' ','')
#save_name=search_name.strip()
data_rebin_neo=pd.DataFrame([])
data_rebin_mep=pd.DataFrame([])
timelabel,maglabel,e_maglabel=('mjd','w1mpro','w1sigmpro')
band='W1'
mean_mjds,mean_mags,mean_mags_e=median_lc(data_test_neo,timelabel,maglabel,e_maglabel,interval=180,sigma_outlier=sigma_outlier)
data_rebin_neo['t_W1']=mean_mjds
data_rebin_neo['W1']=mean_mags
data_rebin_neo['e_W1']=mean_mags_e
timelabel,maglabel,e_maglabel=('mjd','w2mpro','w2sigmpro')
band='W2'
mean_mjds,mean_mags,mean_mags_e=median_lc(data_test_neo,timelabel,maglabel,e_maglabel,interval=180,sigma_outlier=sigma_outlier)
data_rebin_neo['t_W2']=mean_mjds
data_rebin_neo['W2']=mean_mags
data_rebin_neo['e_W2']=mean_mags_e
#############################################################
timelabel,maglabel,e_maglabel=('mjd','w1mpro_ep','w1sigmpro_ep')
band='W1'
mean_mjds,mean_mags,mean_mags_e=median_lc(data_test_mep,timelabel,maglabel,e_maglabel,interval=180,sigma_outlier=sigma_outlier)
data_rebin_mep['t_W1']=mean_mjds
data_rebin_mep['W1']=mean_mags
data_rebin_mep['e_W1']=mean_mags_e
timelabel,maglabel,e_maglabel=('mjd','w2mpro_ep','w2sigmpro_ep')
band='W2'
mean_mjds,mean_mags,mean_mags_e=median_lc(data_test_mep,timelabel,maglabel,e_maglabel,interval=180,sigma_outlier=sigma_outlier)
data_rebin_mep['t_W2']=mean_mjds
data_rebin_mep['W2']=mean_mags
data_rebin_mep['e_W2']=mean_mags_e
return data_rebin_neo,data_rebin_mep
def plot_mag_w1(ax,data_test_neo,data_test_mep):
if len(data_test_neo)>0:
ax.errorbar(data_test_neo['mjd'],data_test_neo['w1mpro'],data_test_neo['w1sigmpro'],
fillstyle="none",
ls='',marker='o',color='red',label='W1',zorder=0)
if len(data_test_mep)>0:
ax.errorbar(data_test_mep['mjd'],data_test_mep['w1mpro_ep'],data_test_mep['w1sigmpro_ep'],
fillstyle="none",
ls='',marker='o',color='red',label='W1',zorder=0)
#ax.scatter(data_test_neo_w1_outlier['mjd'],data_test_neo_w1_outlier['w1mpro'],marker='.',color='red',zorder=1)
#ax.scatter(data_test_mep_w1_outlier['mjd'],data_test_mep_w1_outlier['w1mpro_ep'],marker='.',color='red',zorder=1)
def plot_mag_w2(ax,data_test_neo,data_test_mep):
if len(data_test_neo)>0:
ax.errorbar(data_test_neo['mjd'],data_test_neo['w2mpro'],data_test_neo['w2sigmpro'],
fillstyle="none",
ls='',marker='o',color='blue',label='W2',zorder=0)
if len(data_test_mep)>0:
ax.errorbar(data_test_mep['mjd'],data_test_mep['w2mpro_ep'],data_test_mep['w2sigmpro_ep'],
fillstyle="none",
ls='',marker='o',color='blue',label='W2',zorder=0)
#ax.scatter(data_test_neo_w2_outlier['mjd'],data_test_neo_w2_outlier['w2mpro'],marker='.',color='red',zorder=1)
#ax.scatter(data_test_mep_w2_outlier['mjd'],data_test_mep_w2_outlier['w2mpro_ep'],marker='.',color='red',zorder=1)
def plot_rebinmag_w1(ax,data_rebin_neo,data_rebin_mep):
if len(data_rebin_neo)>0:
ax.errorbar(data_rebin_neo['t_W1'],data_rebin_neo['W1'],data_rebin_neo['e_W1'],
ls='',marker='^',color='red',label='W1',zorder=2)
if len(data_rebin_mep)>0:
ax.errorbar(data_rebin_mep['t_W1'],data_rebin_mep['W1'],data_rebin_mep['e_W1'],
ls='',marker='^',color='red',label='W1',zorder=2)
def plot_rebinmag_w2(ax,data_rebin_neo,data_rebin_mep):
if len(data_rebin_neo)>0:
ax.errorbar(data_rebin_neo['t_W2'],data_rebin_neo['W2'],data_rebin_neo['e_W2'],
ls='',marker='+',color='blue',label='W2',zorder=2)
if len(data_rebin_mep)>0:
ax.errorbar(data_rebin_mep['t_W2'],data_rebin_mep['W2'],data_rebin_mep['e_W2'],
ls='',marker='+',color='blue',label='W2',zorder=2)
def plot_mag_w1_filled(ax,data_test_neo,data_test_mep):
if len(data_test_neo)>0:
ax.errorbar(data_test_neo['mjd'],data_test_neo['w1mpro'],data_test_neo['w1sigmpro'],
markersize=3,
ls='',marker='.',color='red',label='W1',zorder=0)
if len(data_test_mep)>0:
ax.errorbar(data_test_mep['mjd'],data_test_mep['w1mpro_ep'],data_test_mep['w1sigmpro_ep'],
markersize=3,
ls='',marker='.',color='red',label='W1',zorder=0)
#ax.scatter(data_test_neo_w1_outlier['mjd'],data_test_neo_w1_outlier['w1mpro'],marker='.',color='red',zorder=1)
#ax.scatter(data_test_mep_w1_outlier['mjd'],data_test_mep_w1_outlier['w1mpro_ep'],marker='.',color='red',zorder=1)
def plot_mag_w2_filled(ax,data_test_neo,data_test_mep):
if len(data_test_neo)>0:
ax.errorbar(data_test_neo['mjd'],data_test_neo['w2mpro'],data_test_neo['w2sigmpro'],
markersize=3,
ls='',marker='.',color='blue',label='W2',zorder=0)
if len(data_test_mep)>0:
ax.errorbar(data_test_mep['mjd'],data_test_mep['w2mpro_ep'],data_test_mep['w2sigmpro_ep'],
markersize=3,
ls='',marker='.',color='blue',label='W2',zorder=0)
def plot_mag_w1_outlier(ax,data_test_neo,data_test_mep):
if len(data_test_neo)>0:
ax.errorbar(data_test_neo['mjd'],data_test_neo['w1mpro'],data_test_neo['w1sigmpro'],
markersize=4,
ls='',marker='*',color='red',label='W1',zorder=5)
if len(data_test_mep)>0:
ax.errorbar(data_test_mep['mjd'],data_test_mep['w1mpro_ep'],data_test_mep['w1sigmpro_ep'],
markersize=4,
ls='',marker='*',color='red',label='W1',zorder=5)
#ax.scatter(data_test_neo_w1_outlier['mjd'],data_test_neo_w1_outlier['w1mpro'],marker='.',color='red',zorder=1)
#ax.scatter(data_test_mep_w1_outlier['mjd'],data_test_mep_w1_outlier['w1mpro_ep'],marker='.',color='red',zorder=1)
def plot_mag_w2_outlier(ax,data_test_neo,data_test_mep):
if len(data_test_neo)>0:
ax.errorbar(data_test_neo['mjd'],data_test_neo['w2mpro'],data_test_neo['w2sigmpro'],
markersize=4,
ls='',marker='*',color='blue',label='W2',zorder=5)
if len(data_test_mep)>0:
ax.errorbar(data_test_mep['mjd'],data_test_mep['w2mpro_ep'],data_test_mep['w2sigmpro_ep'],
markersize=4,
ls='',marker='*',color='blue',label='W2',zorder=5)
**def wise_plot_lc **
def wise_plot_lc(name_test,save_lc_img_dir=save_lc_img_dir,datawise_dir='/Users/lyubing/blog/brettlv.github.io/pythoncode/changinglookAGN/CLAGN_wise_data/'):
data_test_neo,data_test_mep=get_wise_data(name_test,datawise_dir)
#data_rebin_neo,data_rebin_mep=get_rebinwise_data(data_test_neo,data_test_mep,sigma_outlier=10)
search_name=name_test.replace(' ','')
save_name=search_name.strip()
save_lc_img_path=os.path.join(save_lc_img_dir,'%s_wise.png'%save_name)
fig = plt.figure(figsize=(8,4))
fig.subplots_adjust(hspace=0.0, wspace = 0.0)
ax = fig.add_subplot(111)
#plot_mag_w1(ax,data_test_neo,data_test_mep)
#plot_mag_w2(ax,data_test_neo,data_test_mep)
plot_mag_w1_filled(ax,data_test_neo,data_test_mep)
plot_mag_w2_filled(ax,data_test_neo,data_test_mep)
#plot_rebinmag_w1(ax,data_rebin_neo,data_rebin_mep)
#plot_rebinmag_w2(ax,data_rebin_neo,data_rebin_mep)
'''
data_test_neo_remove,data_test_mep_remove=get_remove_neo_mep(data_test_neo,data_test_mep)
plot_mag_w1_filled(ax,data_test_neo_remove,data_test_mep_remove)
plot_mag_w2_filled(ax,data_test_neo_remove,data_test_mep_remove)
data_test_neo_outlier,data_test_mep_outlier=get_outlier_neo_mep(data_test_neo,data_test_mep)
plot_mag_w1_outlier(ax,data_test_neo_outlier,data_test_mep_outlier)
plot_mag_w2_outlier(ax,data_test_neo_outlier,data_test_mep_outlier)
'''
#ax.scatter(data_test_neo['mjd'],data_test_neo['w1mpro'],color='red',label='W1')
#ax.scatter(data_test_mep['mjd'],data_test_mep['w1mpro_ep'],color='red',label='W1')
#ax.scatter(data_test_neo['mjd'],data_test_neo['w2mpro'],color='blue',label='W2')
#ax.scatter(data_test_mep['mjd'],data_test_mep['w2mpro_ep'],color='blue',label='W2')
plot_secax(ax,mi_interval=365,ma_interval=365*2)
handles, labels = ax.get_legend_handles_labels()
by_label = OrderedDict(zip(labels, handles))
ax.legend(by_label.values(), by_label.keys(),bbox_to_anchor=(0.8, 0.9),ncol=1,
loc=2,fontsize=10)
bottom, top = ax.set_ylim()
if bottom< top:
ax.set_ylim(top,bottom)
plt.xlabel('mjd')
plt.ylabel(r'$mag$')
ax.text(0.5, 0.9,name_test , horizontalalignment = 'center', verticalalignment = 'center',transform = ax.transAxes, fontsize = 10)
if not os.path.exists(save_lc_img_path):
plt.savefig(save_lc_img_path,dpi=400, transparent=False, bbox_inches='tight')
plt.close()
def wise_plotradec_lc(radeg,decdeg,datawise_dir,plotwise_dir):
data_test_neo,data_test_mep=get_wise_radec_data(radeg,decdeg,datawise_dir)
#data_rebin_neo,data_rebin_mep=get_rebinwise_data(data_test_neo,data_test_mep,sigma_outlier=10)
save_name='RA%.5f_DEC%.5f'%(radeg,decdeg)
name_test=save_name
save_lc_img_path=os.path.join(plotwise_dir,'%s_wise.png'%save_name)
fig = plt.figure(figsize=(8,4))
fig.subplots_adjust(hspace=0.0, wspace = 0.0)
ax = fig.add_subplot(111)
#plot_mag_w1(ax,data_test_neo,data_test_mep)
#plot_mag_w2(ax,data_test_neo,data_test_mep)
data_test_neo_remove,data_test_mep_remove=get_remove_neo_mep(data_test_neo,data_test_mep)
plot_mag_w1_filled(ax,data_test_neo_remove,data_test_mep_remove)
plot_mag_w2_filled(ax,data_test_neo_remove,data_test_mep_remove)
#plot_rebinmag_w1(ax,data_rebin_neo,data_rebin_mep)
#plot_rebinmag_w2(ax,data_rebin_neo,data_rebin_mep)
data_test_neo_outlier,data_test_mep_outlier=get_outlier_neo_mep(data_test_neo,data_test_mep)
plot_mag_w1_outlier(ax,data_test_neo_outlier,data_test_mep_outlier)
plot_mag_w2_outlier(ax,data_test_neo_outlier,data_test_mep_outlier)
#ax.scatter(data_test_neo['mjd'],data_test_neo['w1mpro'],color='red',label='W1')
#ax.scatter(data_test_mep['mjd'],data_test_mep['w1mpro_ep'],color='red',label='W1')
#ax.scatter(data_test_neo['mjd'],data_test_neo['w2mpro'],color='blue',label='W2')
#ax.scatter(data_test_mep['mjd'],data_test_mep['w2mpro_ep'],color='blue',label='W2')
plot_secax(ax,mi_interval=365,ma_interval=365*2)
handles, labels = ax.get_legend_handles_labels()
by_label = OrderedDict(zip(labels, handles))
ax.legend(by_label.values(), by_label.keys(),bbox_to_anchor=(0.8, 0.9),ncol=1,
loc=2,fontsize=10)
bottom, top = ax.set_ylim()
if bottom< top:
ax.set_ylim(top,bottom)
plt.xlabel('mjd')
plt.ylabel(r'$mag$')
ax.text(0.5, 0.9,name_test , horizontalalignment = 'center', verticalalignment = 'center',transform = ax.transAxes, fontsize = 10)
plt.savefig(save_lc_img_path,dpi=400, transparent=False, bbox_inches='tight')
plt.close()
**def wise_rebin_lc **
def wise_rebin_lc(name_test,save_lc_img_dir=save_lc_img_dir,datawise_dir='/Users/lyubing/blog/brettlv.github.io/pythoncode/changinglookAGN/CLAGN_wise_data/'):
data_test_neo,data_test_mep=get_wise_data(name_test,datawise_dir)
if len(data_test_neo)>20:
data_test_neo,data_test_mep=get_remove_neo_mep(data_test_neo,data_test_mep)
data_rebin_neo,data_rebin_mep=get_rebinwise_data(data_test_neo,data_test_mep,sigma_outlier=100)
#data_rebin_neo,data_rebin_mep=get_rebinwise_data(data_test_neo,data_test_mep,sigma_outlier=10)
search_name=name_test.replace(' ','')
save_name=search_name.strip()
save_lc_img_path=os.path.join(save_lc_img_dir,'%s_rebin.png'%save_name)
fig = plt.figure(figsize=(8,4))
fig.subplots_adjust(hspace=0.0, wspace = 0.0)
ax = fig.add_subplot(111)
#plot_mag_w1(ax,data_test_neo,data_test_mep)
#plot_mag_w2(ax,data_test_neo,data_test_mep)
plot_rebinmag_w1(ax,data_rebin_neo,data_rebin_mep)
plot_rebinmag_w2(ax,data_rebin_neo,data_rebin_mep)
#ax.scatter(data_test_neo['mjd'],data_test_neo['w1mpro'],color='red',label='W1')
#ax.scatter(data_test_mep['mjd'],data_test_mep['w1mpro_ep'],color='red',label='W1')
plot_secax(ax,mi_interval=365,ma_interval=365*2)
#ax.scatter(data_test_neo['mjd'],data_test_neo['w2mpro'],color='blue',label='W2')
#ax.scatter(data_test_mep['mjd'],data_test_mep['w2mpro_ep'],color='blue',label='W2')
handles, labels = ax.get_legend_handles_labels()
by_label = OrderedDict(zip(labels, handles))
ax.legend(by_label.values(), by_label.keys(),bbox_to_anchor=(0.2, 0.99),ncol=1,
loc=2,fontsize=10)
bottom, top = ax.set_ylim()
if bottom< top:
ax.set_ylim(top,bottom)
plt.xlabel('mjd')
plt.ylabel(r'$mag$')
ax.text(0.5, 0.9,name_test , horizontalalignment = 'center', verticalalignment = 'center',transform = ax.transAxes, fontsize = 10)
if not os.path.exists(save_lc_img_path):
plt.savefig(save_lc_img_path,dpi=400, transparent=False, bbox_inches='tight')
plt.close()
def wise_rebin_lcandcolor(name_test,save_lc_img_dir=save_lc_img_dir,datawise_dir='/Users/lyubing/blog/brettlv.github.io/pythoncode/changinglookAGN/CLAGN_wise_data/'):
data_test_neo,data_test_mep=get_wise_data(name_test,datawise_dir)
if len(data_test_neo)>20:
data_test_neo,data_test_mep=get_remove_neo_mep(data_test_neo,data_test_mep)
data_rebin_neo,data_rebin_mep=get_rebinwise_data(data_test_neo,data_test_mep,sigma_outlier=100)
#data_rebin_neo,data_rebin_mep=get_rebinwise_data(data_test_neo,data_test_mep,sigma_outlier=10)
search_name=name_test.replace(' ','')
save_name=search_name.strip()
save_lc_img_path=os.path.join(save_lc_img_dir,'%s_rebin_lc_color.png'%save_name)
fig = plt.figure(figsize=(8,8))
fig.subplots_adjust(hspace=0.0, wspace = 0.0)
ax = fig.add_subplot(211)
#plot_mag_w1(ax,data_test_neo,data_test_mep)
#plot_mag_w2(ax,data_test_neo,data_test_mep)
plot_rebinmag_w1(ax,data_rebin_neo,data_rebin_mep)
plot_rebinmag_w2(ax,data_rebin_neo,data_rebin_mep)
#ax.scatter(data_test_neo['mjd'],data_test_neo['w1mpro'],color='red',label='W1')
#ax.scatter(data_test_mep['mjd'],data_test_mep['w1mpro_ep'],color='red',label='W1')
plot_secax(ax,mi_interval=365,ma_interval=365*2)
#ax.scatter(data_test_neo['mjd'],data_test_neo['w2mpro'],color='blue',label='W2')
#ax.scatter(data_test_mep['mjd'],data_test_mep['w2mpro_ep'],color='blue',label='W2')
handles, labels = ax.get_legend_handles_labels()
by_label = OrderedDict(zip(labels, handles))
ax.legend(by_label.values(), by_label.keys(),bbox_to_anchor=(0.2, 0.99),ncol=1,
loc=2,fontsize=10)
bottom, top = ax.set_ylim()
if bottom< top:
ax.set_ylim(top,bottom)
ax.set_ylabel(r'$mag$')
ax.text(0.5, 0.9,name_test , horizontalalignment = 'center', verticalalignment = 'center',transform = ax.transAxes, fontsize = 10)
ax.xaxis.set_tick_params(which='major', size=10, width=2, direction='in', top='on',labelsize=0)
ax1 = fig.add_subplot(212,sharex=ax)
ax1.set_xlabel('mjd')
ax1.errorbar(data_rebin_neo['t_W1'],data_rebin_neo['W1']-data_rebin_neo['W2'],ls='',marker='.',color='black')
ax1.errorbar(data_rebin_mep['t_W1'],data_rebin_mep['W1']-data_rebin_mep['W2'],ls='',marker='.',color='black')
ax1.set_ylabel(r'$W1-W2$')
set_ax_tick(ax1)
if not os.path.exists(save_lc_img_path):
plt.savefig(save_lc_img_path,dpi=400, transparent=False, bbox_inches='tight')
plt.close()
def wise_radec_rebin_lc(radeg,decdeg,datawise_dir,plotwise_dir):
data_test_neo,data_test_mep=get_wise_radec_data(radeg,decdeg,datawise_dir)
if len(data_test_neo)>20:
data_test_neo,data_test_mep=get_remove_neo_mep(data_test_neo,data_test_mep)
data_rebin_neo,data_rebin_mep=get_rebinwise_data(data_test_neo,data_test_mep,sigma_outlier=100)
#data_rebin_neo,data_rebin_mep=get_rebinwise_data(data_test_neo,data_test_mep,sigma_outlier=10)
save_name='RA%.5f_DEC%.5f'%(radeg,decdeg)
name_test=save_name
save_lc_img_path=os.path.join(plotwise_dir,'%s_rebin.png'%save_name)
fig = plt.figure(figsize=(8,4))
fig.subplots_adjust(hspace=0.0, wspace = 0.0)
ax = fig.add_subplot(111)
#plot_mag_w1(ax,data_test_neo,data_test_mep)
#plot_mag_w2(ax,data_test_neo,data_test_mep)
plot_rebinmag_w1(ax,data_rebin_neo,data_rebin_mep)
plot_rebinmag_w2(ax,data_rebin_neo,data_rebin_mep)
#ax.scatter(data_test_neo['mjd'],data_test_neo['w1mpro'],color='red',label='W1')
#ax.scatter(data_test_mep['mjd'],data_test_mep['w1mpro_ep'],color='red',label='W1')
plot_secax(ax,mi_interval=365,ma_interval=365*2)
#ax.scatter(data_test_neo['mjd'],data_test_neo['w2mpro'],color='blue',label='W2')
#ax.scatter(data_test_mep['mjd'],data_test_mep['w2mpro_ep'],color='blue',label='W2')
handles, labels = ax.get_legend_handles_labels()
by_label = OrderedDict(zip(labels, handles))
ax.legend(by_label.values(), by_label.keys(),bbox_to_anchor=(0.2, 0.99),ncol=1,
loc=2,fontsize=10)
bottom, top = ax.set_ylim()
if bottom< top:
ax.set_ylim(top,bottom)
plt.xlabel('mjd')
plt.ylabel(r'$mag$')
ax.text(0.5, 0.9,name_test , horizontalalignment = 'center', verticalalignment = 'center',transform = ax.transAxes, fontsize = 10)
if not os.path.exists(save_lc_img_path):
plt.savefig(save_lc_img_path,dpi=400, transparent=False, bbox_inches='tight')
plt.close()
**def ztfdata_download **
def ztfdata_downloadradec(ra,dec,pathztf):
data = lightcurve.LCQuery.download_data(circle=[ra,dec,2.0/3600],BAD_CATFLAGS_MASK=32768)
save_name= 'ra%.5f_dec%.5f'%(ra,dec)
data.to_csv('%s/%s.csv'%(pathztf,save_name),index=False)
#return data
def ztfdata_download_name(name,ra,dec,pathztf):
data = lightcurve.LCQuery.download_data(circle=[ra,dec,2.0/3600],BAD_CATFLAGS_MASK=65535)#32768
#save_name= 'ra%.5f_dec%.5f'%(ra,dec)
save_name= name.strip().replace(' ','')
data.to_csv('%s/%s.csv'%(pathztf,save_name),index=False)
#return data
**def download_wise **
#Irsa.ROW_LIMIT = 1000 # value of new row limit here.
#Irsa.TIMEOUT = 120
#radius=2*u.arcsec default 10*u.arcsec
def download_wise(name,wisedata_path):
#print(i)
search_name=name.replace(' ','')
save_name=name.strip().replace(' ','')
if not os.path.exists('%s/wise_%s_mep.csv'%(wisedata_path,save_name)):
print(search_name)
try:
Irsa.ROW_LIMIT = 1000
Irsa.TIMEOUT = 120
table_wise_mep=Irsa.query_region(search_name,catalog='allwise_p3as_mep',spatial='Cone',radius=10*u.arcsec,)
table_wise_mep.write('%s/wise_%s_mep.csv'%(wisedata_path,save_name), format='csv')
time.sleep(10+np.random.randint(5,10))
except Exception:
print(name,'mep')
#continue
#raise
if not os.path.exists('%s/wise_%s_neo.csv'%(wisedata_path,save_name)):
try:
Irsa.ROW_LIMIT = 1000
Irsa.TIMEOUT = 120
table_wise_neo=Irsa.query_region(search_name,catalog='neowiser_p1bs_psd',spatial='Cone',radius=10*u.arcsec,)
table_wise_neo.write('%s/wise_%s_neo.csv'%(wisedata_path,save_name), format='csv')
time.sleep(10+np.random.randint(5,10))
except Exception:
#continue
print(name,'neo')
#raise
#Irsa.ROW_LIMIT = 1000 # value of new row limit here.
#Irsa.TIMEOUT = 120
#radius=2*u.arcsec default 10*u.arcsec
def download_wise_withradec(radeg,decdeg,wisedata_path):
#print(i)
#search_name=name.replace(' ','')
#save_name=name.strip().replace(' ','')
save_name='RA%.5f_DEC%.5f'%(radeg,decdeg)
name=save_name
if not os.path.exists('%s/wise_%s_mep.csv'%(wisedata_path,save_name)):
#print(search_name)
try:
Irsa.ROW_LIMIT = 1000
Irsa.TIMEOUT = 120
table_wise_mep=Irsa.query_region(SkyCoord(radeg,
decdeg, unit=(u.deg,u.deg),frame='icrs'),
catalog='allwise_p3as_mep',
spatial='Cone',radius=2*u.arcsec,)
table_wise_mep.write('%s/wise_%s_mep.csv'%(wisedata_path,save_name), format='csv')
time.sleep(10+np.random.randint(5,10))
except Exception as e:
print(name,'mep')
print(e)
#continue
#raise
if not os.path.exists('%s/wise_%s_neo.csv'%(wisedata_path,save_name)):
try:
Irsa.ROW_LIMIT = 1000
Irsa.TIMEOUT = 120
table_wise_neo=Irsa.query_region(SkyCoord(radeg,
decdeg, unit=(u.deg,u.deg),frame='icrs'),
catalog='neowiser_p1bs_psd',spatial='Cone',radius=2*u.arcsec,)
table_wise_neo.write('%s/wise_%s_neo.csv'%(wisedata_path,save_name), format='csv')
time.sleep(10+np.random.randint(5,10))
except Exception:
#continue
print(name,'neo')
#raise
**ZTF_mag2flux **
def ZTF_mag2flux(mag):
flux=10**(-0.4*(mag-8.9)) *1e-23 # erg s-1 cm-2 Hz-1
return flux
def ZTF_mag2logflux(mag):
flux=10**(-0.4*(mag-8.9)) *1e-23 # erg s-1 cm-2 Hz-1
return np.log10(flux)
**def ztflc_show_in_jupyter **
def ztflc_show_in_jupyter(data,name,figsize=[7,4],showtoday=False, show_upperlimits=False, formattime=False,
marker=".", mec="0.7", ms=6, ecolor="0.7", ls="", zorder=4, **kwargs):
""" kwargs goes to matplotlib's errorbar() """
FILTER_COLORS = ["C2","C3","C1"]
FILTER_CODE = ["zg","zr","zi"]
marker_CODE = [".","^","*"]
save_name=name.strip().replace(' ','_')
from astropy import time
import matplotlib.pyplot as mpl
# ----------- #
# Global #
# ----------- #
lc_dataframe = data
prop = dict(marker=marker, mec=mec, ms=ms, ecolor=ecolor, ls=ls, zorder=zorder)
# ----------- #
#
# ----------- #
fig = mpl.figure(figsize=figsize)
ax = fig.add_axes([0.1,0.15,0.8,0.75])
for filter_ in np.unique(lc_dataframe["filtercode"]):
d = lc_dataframe[lc_dataframe["filtercode"]==filter_]
if filter_ in FILTER_CODE:
prop["color"] = np.asarray(FILTER_COLORS)[np.where(np.asarray(FILTER_CODE)==filter_)][0]
else:
prop["color"] = "0.7"
if formattime:
dates = time.Time(np.asarray(d["mjd"].values, dtype="float"), format="mjd").datetime
else:
dates = d["mjd"]
ax.errorbar(dates, d["mag"],
yerr=d["magerr"],label=filter_,
**{**prop,**kwargs})
ax.invert_yaxis()
if show_upperlimits:
for filter_ in np.unique(lc_dataframe["filtercode"]):
d = lc_dataframe[lc_dataframe["filtercode"]==filter_]
if filter_ in FILTER_CODE:
color = np.asarray(FILTER_COLORS)[np.where(np.asarray(FILTER_CODE)==filter_)][0]
else:
color = "0.7"
if formattime:
dates = time.Time(np.asarray(d["mjd"].values, dtype="float"), format="mjd").datetime
else:
dates = d["mjd"]
ax.errorbar(dates, d["limitmag"],
yerr=0.1, lolims=True,marker="None", ls="None",
color=color, alpha=0.1,label=filter_
)
if showtoday:
today_color = "0.7"
import datetime
today = time.Time(datetime.date.today().isoformat(),format="iso").mjd
ax.axvline(today, ls="--", color=today_color, zorder=1, lw=1)
ax.text(today, ax.get_ylim()[0]-0.05, "Today", va="bottom", ha="right", rotation=90, color=today_color)
maxmjd= max(lc_dataframe['mjd'])
minmjd= min(lc_dataframe['mjd'])
delta_mjd=maxmjd-minmjd
if delta_mjd>750:
ma_intervals= 365
min_intervals= 90
if delta_mjd<750:
ma_intervals= delta_mjd//3
mi_intervals= delta_mjd//6
ax.set_xlabel("MJD")
plot_secax(ax,rotation=30,)
# - Labels
ax.set_ylabel("mag")
#ax.set_title('%s_%s'%(name,band))
ax.set_title('%s'%(name))
set_ax_legend(ax,bbox_to_anchor=(1.0, 0.99))
return fig,ax
**def plot_wise_ztf_lc **
def plot_wise_ztf_lc(name,ra,dec,ztfdir,wisedir,plotdir,marker=".", mec="0.7", ms=6, ecolor="0.7", ls="", zorder=4,**kwargs):
#filename_lamost='%d.fits'%obsid
prop = dict(color='green',marker=marker, mec='green', ms=ms, ecolor='green', ls=ls, zorder=zorder)
prop1 = dict(color='red',marker='.', mec='red', ms=6, ecolor='red', ls=ls, zorder=zorder)
prop2 = dict(color='blue',marker='.', mec='blue', ms=6, ecolor='blue', ls=ls, zorder=zorder)
ztfsave_name= name.strip().replace(' ','')
ztf_file_path='%s/%s.csv'%(ztfdir,ztfsave_name)
data=pd.read_csv(ztf_file_path)
data_test_neo,data_test_mep=get_wise_data(name,datawise_dir=wisedir)
data_rebin_neo,data_rebin_mep=get_rebinwise_data(data_test_neo,data_test_mep,sigma_outlier=100)
fig=plt.figure()
fig = plt.figure(figsize=(12,6))
fig.subplots_adjust(hspace=0.0, wspace = 0.2)
ax1 = fig.add_subplot(211)
lc_dataframe = data[data['filtercode']=='zg']
lc_dataframe=drop_index(lc_dataframe)
lc_dataframe_zr = data[data['filtercode']=='zr']
lc_dataframe_zr =drop_index(lc_dataframe_zr)
lc_dataframe_zi = data[data['filtercode']=='zi']
lc_dataframe_zi = drop_index(lc_dataframe_zi)
maxmjd= max(lc_dataframe['mjd'])
minmjd= min(lc_dataframe['mjd'])
delta_mjd=maxmjd-minmjd
if delta_mjd>750:
max_intervals= 365
min_intervals= 30
if delta_mjd<750:
max_intervals= delta_mjd//3
min_intervals= delta_mjd//6
ax1.errorbar(lc_dataframe['mjd'], lc_dataframe["mag"], yerr=lc_dataframe["magerr"],label='zg',
**{**prop,**kwargs})
ax1.errorbar(lc_dataframe_zr['mjd'], lc_dataframe_zr["mag"], yerr=lc_dataframe_zr["magerr"],label='zr',
**{**prop1,**kwargs})
ax1.errorbar(lc_dataframe_zi['mjd'], lc_dataframe_zi["mag"], yerr=lc_dataframe_zi["magerr"],label='zi',
**{**prop2,**kwargs})
ax1.set_xlabel("MJD")
#plot_secax(ax1,mi_interval=min_intervals,ma_interval=max_intervals,rotation=30,)
plot_secax(ax1,mi_interval=min_intervals,ma_interval=max_intervals,rotation=30,)
#ax1.axvline(mjd,color='black',ls=':',linewidth=1)
handles, labels = ax1.get_legend_handles_labels()
by_label = OrderedDict(zip(labels, handles))
ax1.legend(by_label.values(), by_label.keys(),bbox_to_anchor=(0.2, 0.99),ncol=1,
loc=2,fontsize=10)
#ax1.set_xlim(58100,59999)
ax2 = fig.add_subplot(212,sharex=ax1)
plot_rebinmag_w1(ax2,data_rebin_neo,data_rebin_mep)
plot_rebinmag_w2(ax2,data_rebin_neo,data_rebin_mep)
#ax.scatter(data_test_neo['mjd'],data_test_neo['w1mpro'],color='red',label='W1')
#ax.scatter(data_test_mep['mjd'],data_test_mep['w1mpro_ep'],color='red',label='W1')
#plot_secax(ax2,mi_interval=365,ma_interval=365*2)
#ax.scatter(data_test_neo['mjd'],data_test_neo['w2mpro'],color='blue',label='W2')
#ax.scatter(data_test_mep['mjd'],data_test_mep['w2mpro_ep'],color='blue',label='W2')
handles, labels = ax2.get_legend_handles_labels()
by_label = OrderedDict(zip(labels, handles))
ax2.legend(by_label.values(), by_label.keys(),bbox_to_anchor=(0.2, 0.99),ncol=1,
loc=2,fontsize=10)
ax1.invert_yaxis()
ax2.invert_yaxis()
ax2.set_xlabel('mjd')
ax1.set_ylabel(r'$mag$')
ax2.set_ylabel(r'$mag$')
ax1.set_title('%s'%(ztfsave_name))
save_lc_img_path=os.path.join(plotdir,'%s.png'%(ztfsave_name))
#if not os.path.exists(save_lc_img_path):
#plt.savefig(save_lc_img_path,dpi=400, transparent=False, bbox_inches='tight')
#plt.close()
return fig,ax1,ax2
def plot_wise_ztf_lc_timerange(name,ra,dec,t0,t1,ztfdir,wisedir,plotdir,marker=".", mec="0.7", ms=6, ecolor="0.7", ls="", zorder=4,**kwargs):
#filename_lamost='%d.fits'%obsid
prop = dict(color='green',marker=marker, mec='green', ms=ms, ecolor='green', ls=ls, zorder=zorder)
prop1 = dict(color='red',marker='.', mec='red', ms=6, ecolor='red', ls=ls, zorder=zorder)
prop2 = dict(color='blue',marker='.', mec='blue', ms=6, ecolor='blue', ls=ls, zorder=zorder)
ztfsave_name= name.strip().replace(' ','')
ztf_file_path='%s/%s.csv'%(ztfdir,ztfsave_name)
data=pd.read_csv(ztf_file_path)
data_test_neo,data_test_mep=get_wise_data(name,datawise_dir=wisedir)
data_rebin_neo,data_rebin_mep=get_rebinwise_data(data_test_neo,data_test_mep,sigma_outlier=100)
fig=plt.figure()
fig = plt.figure(figsize=(12,6))
fig.subplots_adjust(hspace=0.0, wspace = 0.2)
ax1 = fig.add_subplot(211)
lc_dataframe = data[data['filtercode']=='zg']
lc_dataframe=drop_index(lc_dataframe)
lc_dataframe_zr = data[data['filtercode']=='zr']
lc_dataframe_zr =drop_index(lc_dataframe_zr)
lc_dataframe_zi = data[data['filtercode']=='zi']
lc_dataframe_zi = drop_index(lc_dataframe_zi)
ax1.set_xlim(t0,t1)
#maxmjd= max(lc_dataframe['mjd'])
#minmjd= min(lc_dataframe['mjd'])
maxmjd=t1
minmjd=t0
delta_mjd=maxmjd-minmjd
if delta_mjd>750:
max_intervals= 365
min_intervals= 30
if delta_mjd<750:
max_intervals= delta_mjd//3
min_intervals= delta_mjd//6
ax1.errorbar(lc_dataframe['mjd'], lc_dataframe["mag"], yerr=lc_dataframe["magerr"],label='zg',
**{**prop,**kwargs})
ax1.errorbar(lc_dataframe_zr['mjd'], lc_dataframe_zr["mag"], yerr=lc_dataframe_zr["magerr"],label='zr',
**{**prop1,**kwargs})
ax1.errorbar(lc_dataframe_zi['mjd'], lc_dataframe_zi["mag"], yerr=lc_dataframe_zi["magerr"],label='zi',
**{**prop2,**kwargs})
ax1.set_xlabel("MJD")
#plot_secax(ax1,mi_interval=min_intervals,ma_interval=max_intervals,rotation=30,)
plot_secax(ax1,mi_interval=min_intervals,ma_interval=max_intervals,rotation=30,)
#ax1.axvline(mjd,color='black',ls=':',linewidth=1)
handles, labels = ax1.get_legend_handles_labels()
by_label = OrderedDict(zip(labels, handles))
ax1.legend(by_label.values(), by_label.keys(),bbox_to_anchor=(0.2, 0.99),ncol=1,
loc=2,fontsize=10)
#ax1.set_xlim(58100,59999)
ax2 = fig.add_subplot(212,sharex=ax1)
plot_rebinmag_w1(ax2,data_rebin_neo,data_rebin_mep)
plot_rebinmag_w2(ax2,data_rebin_neo,data_rebin_mep)
#ax.scatter(data_test_neo['mjd'],data_test_neo['w1mpro'],color='red',label='W1')
#ax.scatter(data_test_mep['mjd'],data_test_mep['w1mpro_ep'],color='red',label='W1')
#plot_secax(ax2,mi_interval=365,ma_interval=365*2)
#ax.scatter(data_test_neo['mjd'],data_test_neo['w2mpro'],color='blue',label='W2')
#ax.scatter(data_test_mep['mjd'],data_test_mep['w2mpro_ep'],color='blue',label='W2')
handles, labels = ax2.get_legend_handles_labels()
by_label = OrderedDict(zip(labels, handles))
ax2.legend(by_label.values(), by_label.keys(),bbox_to_anchor=(0.2, 0.99),ncol=1,
loc=2,fontsize=10)
ax1.invert_yaxis()
ax2.invert_yaxis()
ax2.set_xlabel('mjd')
ax1.set_ylabel(r'$mag$')
ax2.set_ylabel(r'$mag$')
ax1.set_title('%s'%(ztfsave_name))
save_lc_img_path=os.path.join(plotdir,'%s_mjd%d_%d.png'%(ztfsave_name,t0,t1))
#if not os.path.exists(save_lc_img_path):
#plt.savefig(save_lc_img_path,dpi=400, transparent=False, bbox_inches='tight')
#plt.close()
return fig,ax1,ax2
def plot_wise_ztf_lc_radec(name,ra,dec,ztfdir,wisedir,plotdir,marker=".", mec="0.7", ms=6, ecolor="0.7", ls="", zorder=4,**kwargs):
#filename_lamost='%d.fits'%obsid
prop = dict(marker=marker, mec=mec, ms=ms, ecolor=ecolor, ls=ls, zorder=zorder)
prop1 = dict(marker='*', mec=mec, ms=ms, ecolor='red', ls=ls, zorder=zorder)
ztfsave_name= 'ra%.5f_dec%.5f'%(ra,dec)
#data.to_csv('%s/%s.csv'%(pathztf,save_name),index=False)
#ztfsave_name= name.strip().replace(' ','')
ztf_file_path='%s/%s.csv'%(ztfdir,ztfsave_name)
data=pd.read_csv(ztf_file_path)
data_test_neo,data_test_mep=get_wise_RADEC_data(ra,dec,datawise_dir=wisedir)
data_rebin_neo,data_rebin_mep=get_rebinwise_data(data_test_neo,data_test_mep,sigma_outlier=100)
fig=plt.figure()
fig = plt.figure(figsize=(20,8))
fig.subplots_adjust(hspace=0.0, wspace = 0.2)
ax1 = fig.add_subplot(211)
lc_dataframe = data[data['filtercode']=='zg']
lc_dataframe=drop_index(lc_dataframe)
lc_dataframe_zr = data[data['filtercode']=='zr']
lc_dataframe_zr=drop_index(lc_dataframe_zr)
maxmjd= max(lc_dataframe['mjd'])
minmjd= min(lc_dataframe['mjd'])
delta_mjd=maxmjd-minmjd
if delta_mjd>750:
max_intervals= 365
min_intervals= 90
if delta_mjd<750:
max_intervals= delta_mjd//3
min_intervals= delta_mjd//6
ax1.errorbar(lc_dataframe['mjd'], lc_dataframe["mag"], yerr=lc_dataframe["magerr"],label='zg',
**{**prop,**kwargs})
ax1.errorbar(lc_dataframe_zr['mjd'], lc_dataframe_zr["mag"], yerr=lc_dataframe_zr["magerr"],label='zr',
**{**prop1,**kwargs})
ax1.set_xlabel("MJD")
plot_secax(ax1,mi_interval=min_intervals,ma_interval=max_intervals,rotation=30,)
#plot_secax(ax1,rotation=30,)
#ax1.axvline(mjd,color='black',ls=':',linewidth=1)
handles, labels = ax1.get_legend_handles_labels()
by_label = OrderedDict(zip(labels, handles))
ax1.legend(by_label.values(), by_label.keys(),bbox_to_anchor=(0.2, 0.99),ncol=1,
loc=2,fontsize=10)
#ax1.set_xlim(58100,59999)
ax2 = fig.add_subplot(212,sharex=ax1)
plot_rebinmag_w1(ax2,data_rebin_neo,data_rebin_mep)
plot_rebinmag_w2(ax2,data_rebin_neo,data_rebin_mep)
#ax.scatter(data_test_neo['mjd'],data_test_neo['w1mpro'],color='red',label='W1')
#ax.scatter(data_test_mep['mjd'],data_test_mep['w1mpro_ep'],color='red',label='W1')
#plot_secax(ax2,mi_interval=365,ma_interval=365*2)
#ax.scatter(data_test_neo['mjd'],data_test_neo['w2mpro'],color='blue',label='W2')
#ax.scatter(data_test_mep['mjd'],data_test_mep['w2mpro_ep'],color='blue',label='W2')
handles, labels = ax2.get_legend_handles_labels()
by_label = OrderedDict(zip(labels, handles))
ax2.legend(by_label.values(), by_label.keys(),bbox_to_anchor=(0.2, 0.99),ncol=1,
loc=2,fontsize=10)
ax1.invert_yaxis()
ax2.invert_yaxis()
ax2.set_xlabel('mjd')
ax1.set_ylabel(r'$mag$')
ax2.set_ylabel(r'$mag$')
ax1.set_title('%s'%(ztfsave_name))
save_lc_img_path=os.path.join(plotdir,'%s.png'%(ztfsave_name))
if not os.path.exists(save_lc_img_path):
plt.savefig(save_lc_img_path,dpi=400, transparent=False, bbox_inches='tight')
plt.close()
#return fig,ax,ax1,ax2
**T **
**T **
**T **
**T **
欢迎关注微信公众号：曜灵集
