-
Notifications
You must be signed in to change notification settings - Fork 1
/
get_data_hsla.py
158 lines (140 loc) · 6.94 KB
/
get_data_hsla.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
149
150
151
152
153
154
155
156
157
158
"""
.. module:: get_data_hsla
:synopsis: Returns coadded or exposure-level spectra from the
Hubble Spectroscopic Legacy Archive.
.. moduleauthor:: Scott W. Fleming <fleming@stsci.edu>
"""
import collections
import os
import sys
from astropy.io import fits
from data_series import DataSeries
from parse_obsid_hsla import parse_obsid_hsla
#--------------------
def get_data_hsla(obsid, targ):
"""
Given an HSLA observation ID, returns the spectral data. If a
coadd-level spectrum, must supply the target name via the 'targ'
parameter.
:param obsid: The HSLA grism observation ID to retrieve the data from.
:type obsid: str
:param targ: The name of the target, if a coadd-level spectrum.
:type targ: str
:returns: JSON -- The spectral data for this observation ID.
Error codes:
From parse_obsid_hsla_grism:
0 = No error parsing observation ID.
1 = Directory not found.
2 = Extracted spectra FITS file not found.
From this module:
3 = Could not open one or more FITS file for reading.
"""
# This defines the maximum (approximate) size in bytes allowed to be
# returned as a DataSeries.
max_size = 50000000.
# This error code will be used unless there's a problem reading any
# of the FITS files in the list.
errcode = 0
# This defines a data point for a DataSeries object as a namedtuple.
data_point = collections.namedtuple('DataPoint', ['x', 'y'])
# For HSLA grisms, this defines the x-axis and y-axis units as a string.
hsla_xunit = "Angstroms"
hsla_yunit = "ergs/cm^2/s/Angstrom"
# Parse the obsID string to determine the paths+files to read.
parsed_files_result = parse_obsid_hsla(obsid, targ)
errcode = parsed_files_result.errcode
# We create a list of return DataSeries for each segment.
all_data_series = []
# For each file, read in the contents and create a return JSON object.
if errcode == 0:
# This keeps track of the (approximate) total size of the return
# object so that it doesn't get too large (in bytes).
total_size = 0.0
for sfile in parsed_files_result.specfiles:
try:
with fits.open(sfile) as hdulist:
if obsid.lower().strip() != "hsla_coadd":
# Get the segments, which can be ['FUVA' or 'FUVB']
segments = hdulist[1].data['segment']
# Wavelengths and fluxes are stored as a list for
# each seg.
segment_wls = hdulist[1].data['wavelength']
segment_fls = hdulist[1].data["flux"]
segment_flerrs = hdulist[1].data["error"]
else:
# Then this is a coadd, so there aren't multiple segs.
wls = hdulist[1].data['wave']
fls = hdulist[1].data["fluxwgt"]
flerrs = hdulist[1].data["fluxwgt_err"]
except IOError:
errcode = 4
all_data_series.append(DataSeries(
'hsla', obsid, [], [''], [''], [''], errcode,
is_ancillary=[1]))
else:
# Create DataSeries if this is an exposure-level spectrum.
if obsid.lower().strip() != "hsla_coadd":
for this_seg, this_wl, this_fl, this_flerr in zip(
segments, segment_wls, segment_fls, segment_flerrs):
wls = [float(x) for x in this_wl]
fls = [float(x) for x in this_fl]
flerrs = [float(x) for x in this_flerr]
wlfls = [x for x in zip(wls, fls)]
wlfls_err = [x for x in zip(wls, flerrs)]
# Append the wl-fl DataSeries for this segment.
all_data_series.append(DataSeries(
'hsla', obsid,
[[data_point(x=float("{0:.8e}".format(x)),
y=float("{0:.8e}".format(y)))
for x, y in wlfls]],
[obsid+'_'+this_seg], [hsla_xunit], [hsla_yunit],
errcode, is_ancillary=[0]))
this_dataseries = DataSeries(
'hsla', obsid,
[[data_point(x=float("{0:.8e}".format(x)),
y=float("{0:.8e}".format(y)))
for x, y in wlfls_err]],
[obsid+'_'+this_seg+'_ERR'], [hsla_xunit],
[hsla_yunit], errcode, is_ancillary=[1])
this_dataseries_size = sys.getsizeof(
this_dataseries.plot_series[0])
# Append the wl-flerr DataSeries for this segment.
if total_size + this_dataseries_size <= max_size:
all_data_series.append(this_dataseries)
total_size = total_size + this_dataseries_size
else:
# Create DataSeries if this is a coadd-level spectrum.
wls = [float(x) for x in wls]
fls = [float(x) for x in fls]
flerrs = [float(x) for x in flerrs]
wlfls = [x for x in zip(wls, fls)]
wlfls_err = [x for x in zip(wls, flerrs)]
# Only plot the coadd across lifetime positions by default.
if '_all.fits.gz' in os.path.basename(sfile):
is_anc = [0]
else:
is_anc = [1]
# Append the wl-fl DataSeries for this segment.
this_dataseries = DataSeries(
'hsla', obsid,
[[data_point(x=float("{0:.8e}".format(x)),
y=float("{0:.8e}".format(y)))
for x, y in wlfls]],
[os.path.basename(sfile).strip(".fits.gz")],
[hsla_xunit], [hsla_yunit], errcode,
is_ancillary=is_anc)
this_dataseries_size = sys.getsizeof(
this_dataseries.plot_series[0])
# Append the wl-flerr DataSeries for this segment.
if total_size + this_dataseries_size <= max_size:
all_data_series.append(this_dataseries)
total_size = total_size + this_dataseries_size
else:
# This is where an error DataSeries object would be returned.
all_data_series.append(DataSeries(
'hsla', obsid, [], [], [], [], errcode, is_ancillary=[1]))
# Return the DataSeries object back to the calling module.
if len(all_data_series) == 1:
return all_data_series[0]
return all_data_series
#--------------------