-
Notifications
You must be signed in to change notification settings - Fork 0
/
mage_multipanel-spectrum-ids.py
148 lines (131 loc) · 6.53 KB
/
mage_multipanel-spectrum-ids.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
'''This python script plots a MagE/Magellan echelle spectrum that spans
many panels and pages. Written jrigby 4/2012, edited 10/2013.
**********************************************************************
UPDATE 9/2016: THIS SCRIPT IS DEPRECATED!!!!
I WROTE IT A LONG TIME AGO, BEFORE I KNEW ANY PYTHON.
I REWROTE THE ALGORITHM INTO A FUNCTION, NOW IN jrr.plot.echelle_spectrum().
HERE ARE TWO EXAMPLES OF CALLING THAT FUNCTION:
JRR_Code/multipanel-stacks.py and JRR_Code/multipanel-redo_as_function.py
**********************************************************************
#
# This script needs a concatenated linelist, with redshifts customized for the spectrum.
# BEFORE YOU RUN THIS SCRIPT, you must use concat-linelist.pl to generate such a linelist.
#
# Required arguments of this script are:
# 1) The spectrum to plot. That's it.
#
#example: python multipanel-spectrum-ids.py "RCS0327/KnotE/rcs0327-knotE-allres-comb.txt"
'''
from __future__ import print_function
from future import standard_library
standard_library.install_aliases()
from builtins import str
from builtins import range
import jrr
import matplotlib.pyplot as plt
from numpy import ones_like, zeros_like, median
from astropy.io import ascii
from subprocess import getoutput
import sys
from re import split, sub, search
#import lineid_plot
from mpl_toolkits.axes_grid1 import host_subplot
from matplotlib.ticker import MultipleLocator, FormatStrFormatter, MaxNLocator
from matplotlib.backends.backend_pdf import PdfPages
mage_mode = "reduction"
#mage_mode = "released"
#print "DEBUG, arguments received were: ", sys.argv, len(sys.argv)
if(len(sys.argv) == 2):
this_script = sys.argv[0] # echos the name of this python program
inspec = sys.argv[1] # spectrum file to plot.
else:
print("ERROR! One argument required: spectrum to plot")
sys.exit(1)
# Set up plotting
plotsize = (11,16)
Npanels = 24 # N panels total
Nfigs = 4 # divided into this many pages of plots
lincolor = 'k'
errcolor = '0.55'
contcolor = '0.75'
def onclick(event): # Setting up interactive clicking. Right now, just prints location. Need to add fitting.
print('button=%d, x=%d, y=%d, xdata=%f, ydata=%f'%(
event.button, event.x, event.y, event.xdata, event.ydata))
# Handle filenames with regexp
rootname = split(".txt", (split("/", inspec)[-1]))[0] # Grab the last thing after the last slash. Strip out .txt
the_pdf = "PDF_Out/" + rootname + ".pdf"
pp = PdfPages(the_pdf) # output
linelist = "Lines/" + rootname + ".linelist" # Should have been generated by concat-linelist.pl
(LL, zz) = jrr.mage.get_linelist(linelist) # Load the linelist into Pandas data frame, and get systemic redshift
print("Rest-frame wavelength axis will assume redshift of", zz)
if inspec == "stacked" : # This is the stacked spectrum!
(wave, fnu, X_clipavg, X_median, sig, X_jack_std, Ngal) = jrr.mage.open_stacked_spectrum(mage_mode)
cont = ones_like(sig)
cont_u = zeros_like(sig)
factor = 1.0
skipcont = False
else :
factor = 1E29 # for plotting purposes, multiply by scale factor to bring to order unity
(sp, resoln, dresoln) = jrr.mage.open_spectrum(inspec, zz, mage_mode) # sp is the spectrum in a pandas data frame
wave = sp.wave
fnu = sp.fnu
sig = sp.fnu_u
if 'cont' in sp.columns :
cont = sp.fnu_cont
cont_u = sp.fnu_cont_u
else :
skipcont = True
fnu *= factor
sig *= factor
if not skipcont :
cont *= factor
cont_u *= factor
length = len(wave)
for j in range(0,Nfigs):
print("Drawing page ", j, " of ", Nfigs-1)
fig = plt.figure(num=j+1, figsize=plotsize) #from matplotlib
for k in range(0, Npanels/Nfigs):
format = (Npanels/Nfigs)*100 + 11 # makes label: 311 for 3 rows, 1 col, start at 1
subit = host_subplot(format + k)
start = j*(length/Nfigs) + k*(length/Npanels)
end = start + length/Npanels
top = median(fnu[start:end])*1.1 + jrr.util.mad(fnu[start:end])*5
plt.ylim(0, top) # trying to fix weird autoscaling from bad pixels
# plt.autoscale(enable=False, axis='y', tight=True)
plt.plot(wave[start:end], fnu[start:end], lincolor, linestyle='steps') # plot spectrum
plt.plot(wave[start:end], sig[start:end], errcolor, linestyle='steps') # plot 1 sigma uncert spectrum
plt.autoscale(axis='x', tight=True)
if not skipcont :
plt.plot(wave[start:end], cont[start:end], contcolor, linestyle='steps', zorder=1, linewidth=2) # plot the continuum
jrr.mage.plot_linelist(LL, zz) # Plot the line IDs.
upper = subit.twiny() # make upper x-axis in rest wave (systemic, or of abs lines)
upper.set_xlim(wave[start]/(1.0+zz), wave[end]/(1.0+zz))
majorLocator = MultipleLocator(50) # put subticks on lower x-axis
minorLocator = MultipleLocator(10)
subit.xaxis.set_major_locator(majorLocator)
subit.xaxis.set_minor_locator(minorLocator)
subit.xaxis.tick_bottom() # don't let lower ticks be mirrored on upper axis
subit.yaxis.set_major_locator(MaxNLocator(nbins=4, integer=False, Symmetric=False))
if(k == Npanels/Nfigs-1):
subit.set_xlabel(ur"observed-frame vacuum wavelength (\u00c5)")
plt.ylabel('fnu') # fnu in cgs units: erg/s/cm^2/Hz
if(k == 0):
upper.set_xlabel(ur"rest-frame vacuum wavelength (\u00c5)")
plt.annotate("f_nu * "+str(factor) + " erg/s/cm^2/Hz", xy=(0.05,0.055), color="black", xycoords="figure fraction")
plt.annotate("Line color coding:", xy=(0.05,0.043), color="black", xycoords="figure fraction")
plt.annotate("Photosphere", xy=(0.05,0.03), color="blue", xycoords="figure fraction")
plt.annotate("Emission", xy=(0.17,0.03), color="red", xycoords="figure fraction")
plt.annotate("ISM", xy=(0.25,0.03), color="green", xycoords="figure fraction")
plt.annotate("Wind", xy=(0.32,0.03), color="cyan", xycoords="figure fraction")
plt.annotate("Fine structure", xy=(0.4,0.03), color="purple", xycoords="figure fraction")
plt.annotate("Intervening", xy=(0.55,0.03), color="orange", xycoords="figure fraction")
plt.suptitle(rootname + " z=" + str(zz)) # global title
pp.savefig()
fig.canvas.draw()
fig.canvas.mpl_connect('button_press_event', onclick) # This is called an event connection.
pp.close()
print("Ran ", this_script)
print(" Plotted spectrum: ", inspec)
print(" Used linelist: ", linelist)
print(" Generated PDF: ", the_pdf)
print(" DONE.")