Coverage for python/lsst/obs/hsc/makeTransmissionCurves.py: 95%

1#!/usr/bin/env python 

2# 

3# LSST Data Management System 

4# 

5# Copyright 2018 AURA/LSST. 

6# 

7# This product includes software developed by the 

8# LSST Project (http://www.lsst.org/). 

9# 

10# This program is free software: you can redistribute it and/or modify 

11# it under the terms of the GNU General Public License as published by 

12# the Free Software Foundation, either version 3 of the License, or 

13# (at your option) any later version. 

14# 

15# This program is distributed in the hope that it will be useful, 

16# but WITHOUT ANY WARRANTY; without even the implied warranty of 

17# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 

18# GNU General Public License for more details. 

19# 

20# You should have received a copy of the LSST License Statement and 

21# the GNU General Public License along with this program. If not, 

22# see <https://www.lsstcorp.org/LegalNotices/>. 

23# 

24import numpy as np 

25 

26from lsst.afw.image import TransmissionCurve 

27from lsst.resources import ResourcePath 

28 

29__all__ = ("getOpticsTransmission", "getSensorTransmission", "getAtmosphereTransmission", 

30 "getFilterTransmission",) 

31 

32DATA_DIR = ResourcePath("resource://lsst.obs.subaru/resources/hsc/transmission/", forceDirectory=True) 

33 

34HSC_BEGIN = "2012-12-18" # initial date for curves valid for entire lifetime of HSC 

35 

36# This maps the old names used in afw.image.Filter and various NAOJ-side things 

37# to the physical_filter names we've now standardized on. It includes only 

38# broad-band filters; the narrow-band filters are easier to transform without 

39# an explicit mapping (see getPhysicalFilterName). 

40OLD_FILTER_NAME_MAP = {b: f"HSC-{b.upper()}" for b in "grizy"} 

41OLD_FILTER_NAME_MAP.update(r2="HSC-R2", i2="HSC-I2") 

42 

43 

44def getPhysicalFilterName(short): 

45 """Return an HSC physical_filter name (e.g. 'HSC-R') that's usable as a 

46 data ID value from various names used in the transmission data files. 

47 """ 

48 if short.startswith("HSC"): 48 ↛ 49line 48 didn't jump to line 49 because the condition on line 48 was never true

49 return short 

50 # Sometimes narrow-band and intermediate-band filters start with NB or IB, 

51 # and we don't want to try to enumerate all of those if we can avoid it. 

52 if short.startswith("NB") or short.startswith("IB"): 

53 num = int(short[2:].lstrip("0")) 

54 return "%s%04d" % (short[:2], num) 

55 # Sometimes the narrow-band filters start with just N instead. 

56 if short.startswith("N"): 

57 num = int(short[1:].lstrip("0")) 

58 return "NB%04d" % (num,) 

59 return OLD_FILTER_NAME_MAP.get(short, short) 

60 

61 

62def readTransmissionCurveFromFile(filename, unit="angstrom", atMin=None, atMax=None): 

63 """Load a spatial TransmissionCurve from a text file with wavelengths and 

64 throughputs in columns. 

65 

66 Parameters 

67 ---------- 

68 filename : `str` 

69 Name of the file to read. 

70 unit : `str` 

71 Wavelength unit; one of "nm" or "angstrom". 

72 atMin : `float` 

73 Throughput to use at wavelengths below the tabulated minimum. If 

74 ``None``, the tabulated throughput at the mininum will be used. 

75 atMax : `float` 

76 Throughput to use at wavelengths above the tabulated maximum. If 

77 ``None``, the tabulated throughput at the maximum will be used. 

78 """ 

79 with DATA_DIR.join(filename).as_local() as local: 

80 wavelengths, throughput = np.loadtxt(local.ospath, usecols=[0, 1], unpack=True) 

81 i = np.argsort(wavelengths) 

82 wavelengths = wavelengths[i] 

83 throughput = throughput[i] 

84 if unit == "nm": 

85 wavelengths *= 10 

86 elif unit != "angstrom": 86 ↛ 87line 86 didn't jump to line 87 because the condition on line 86 was never true

87 raise ValueError("Invalid wavelength unit") 

88 if atMin is None: 

89 atMin = throughput[0] 

90 if atMax is None: 

91 atMax = throughput[-1] 

92 return TransmissionCurve.makeSpatiallyConstant(throughput=throughput, wavelengths=wavelengths, 

93 throughputAtMin=atMin, throughputAtMax=atMax) 

94 

95 

96def getOpticsTransmission(): 

97 """Return a dictionary of TransmissionCurves describing the combined 

98 throughput of HSC and the Subaru primary mirror. 

99 

100 Dictionary keys are string dates (YYYY-MM-DD) indicating the beginning of 

101 the validity period for the curve stored as the associated dictionary 

102 value. If the curve is spatially varying, it will be defined in focal 

103 plane coordinates. 

104 

105 Dictionary values may be None to indicate that no TransmissionCurve is 

106 valid after the date provided in the key. 

107 """ 

108 mirror2010 = readTransmissionCurveFromFile("M1-2010s.txt", unit="nm") 

109 camera = readTransmissionCurveFromFile("throughput_popt2.txt") 

110 camera *= readTransmissionCurveFromFile("throughput_win.txt") 

111 return { 

112 HSC_BEGIN: mirror2010*camera, 

113 "2017-10-01": None # mirror recoating begins, approximately 

114 } 

115 

116 

117def getSensorTransmission(): 

118 """Return a nested dictionary of TransmissionCurves describing the 

119 throughput of each sensor. 

120 

121 Outer directionary keys are string dates (YYYY-MM-DD), with values 

122 a dictionary mapping CCD ID to TransmissionCurve. If the curve 

123 is spatially varying, it will be defined in pixel coordinates. 

124 

125 Outer dictionary values may be None to indicate that no TransmissionCurve 

126 is valid after the date provided in the key. 

127 """ 

128 qe = readTransmissionCurveFromFile("qe_ccd_HSC.txt", atMin=0.0, atMax=0.0) 

129 return {HSC_BEGIN: {n: qe for n in range(112)}} 

130 

131 

132def getAtmosphereTransmission(): 

133 """Return a dictionary of TransmissionCurves describing the atmospheric 

134 throughput at Mauna Kea. 

135 

136 Dictionary keys are string dates (YYYY-MM-DD) indicating the beginning of 

137 the validity period for the curve stored as the associated dictionary 

138 value. The curve is guaranteed not to be spatially-varying. 

139 

140 Dictionary values may be None to indicate that no TransmissionCurve is 

141 valid after the date provided in the key. 

142 """ 

143 average = readTransmissionCurveFromFile("modtran_maunakea_am12_pwv15_binned10ang.dat") 

144 return {HSC_BEGIN: average} 

145 

146 

147def getFilterTransmission(): 

148 """Return a nested dictionary of TransmissionCurves describing the 

149 throughput of each HSC filter. 

150 

151 Outer directionary keys are string dates (YYYY-MM-DD), with values 

152 a dictionary mapping filter name to TransmissionCurve. If the curve 

153 is spatially varying, it will be defined in pixel coordinates. 

154 

155 Filter curve names are in the long form used as data ID values (e.g. 

156 'HSC-I'). 

157 

158 Outer dictionary values may be None to indicate that no TransmissionCurve 

159 is valid after the date provided in the key. 

160 """ 

161 module = {} 

162 filterTraces = DATA_DIR.join("filterTraces.py") 

163 exec(compile(filterTraces.read().decode(), str(filterTraces), mode='exec'), module) 

164 result = {} 

165 for band, data in module["FILTER_DATA"].items(): 

166 result[getPhysicalFilterName(band)] = TransmissionCurve.makeRadial( 

167 throughput=data["T"], wavelengths=data["lam"]*10, 

168 radii=data['radius']/module["PIXEL_SIZE"], 

169 throughputAtMin=0.0, throughputAtMax=0.0 

170 ) 

171 for uri in ResourcePath.findFileResources([DATA_DIR], file_filter=r"wHSC-.*\.txt$"): 

172 fname = uri.basename() 

173 band = getPhysicalFilterName(fname[len("wHSC-"): -len(".txt")]) 

174 if band not in result: 

175 result[band] = readTransmissionCurveFromFile(fname, atMin=0.0, atMax=0.0) 

176 return {HSC_BEGIN: result}