Coverage for python / lsst / summit / utils / astrometry / plotting.py: 0%

1# This file is part of summit_utils. 

2# 

3# Developed for the LSST Data Management System. 

4# This product includes software developed by the LSST Project 

5# (https://www.lsst.org). 

6# See the COPYRIGHT file at the top-level directory of this distribution 

7# for details of code ownership. 

8# 

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

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

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

12# (at your option) any later version. 

13# 

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

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

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

17# GNU General Public License for more details. 

18# 

19# You should have received a copy of the GNU General Public License 

20# along with this program. If not, see <https://www.gnu.org/licenses/>. 

21 

22import copy 

23from typing import Any 

24 

25import astropy.units as u 

26import matplotlib 

27import matplotlib.pyplot as plt 

28import numpy as np 

29from astropy.coordinates import Angle 

30 

31import lsst.afw.image as afwImage 

32import lsst.afw.table as afwTable 

33from lsst.obs.lsst.translators.latiss import AUXTEL_LOCATION 

34from lsst.utils.plotting.figures import make_figure 

35 

36from .. import quickSmooth 

37 

38# TODO: Add some of Craig's nice overlay stuff here 

39 

40 

41def plot( 

42 exp: afwImage.Exposure, 

43 icSrc: afwTable.SourceCatalog = None, 

44 filteredSources: afwTable.SourceCatalog = None, 

45 saveAs: str | None = None, 

46 clipMin: float = 1, 

47 clipMax: float = 1000000, 

48 doSmooth: bool = True, 

49 fig: matplotlib.figure.Figure | None = None, 

50) -> None: 

51 """Plot an exposure, overlaying the selected sources and compass arrows. 

52 

53 Plots the exposure on a logNorm scale, with the brightest sources, as 

54 selected by the configuration, overlaid with an x. Plots compass arrows 

55 for both north/east and az/el. Optionally saves the output to a file if 

56 ``saveAs`` is supplied. 

57 

58 Parameters 

59 ---------- 

60 exp : `lsst.afw.image.Exposure` 

61 The exposure to get the astrometry for. 

62 icSrc : `lsst.afw.table.SourceCatalog` 

63 The source catalog for the exposure. 

64 filteredSources : `lsst.afw.table.SourceCatalog`, optional 

65 The filtered source catalog. If supplied, shows which sources were 

66 selected. 

67 saveAs : `str`, optional 

68 Saves the plot to this filename if specified. 

69 clipMin : `float` 

70 Clip values in the image below this value to ``clipMin``. 

71 clipMax : `float` 

72 Clip values in the image above this value to ``clipMax``. 

73 doSmooth : `bool`, optional 

74 Smooth the image slightly to improve the visability of low SNR sources? 

75 fig : `matplotlib.figure.Figure`, optional 

76 The figure to plot on. If not supplied, a new figure is created. 

77 """ 

78 if fig is None: 

79 fig = make_figure(figsize=(16, 16)) 

80 

81 fig.clear() 

82 ax = fig.gca() 

83 ax.clear() 

84 

85 data = copy.deepcopy(exp.image.array) 

86 data = np.clip(data, clipMin, clipMax) 

87 if doSmooth: 

88 data = quickSmooth(data) # smooth slightly to help visualize 

89 plt.imshow(np.arcsinh(data) / 10, interpolation="None", cmap="gray", origin="lower") 

90 

91 height, width = data.shape 

92 leftFraction = 0.15 # fraction into the image to start the N/E compass 

93 rightFraction = 0.225 # fraction into the image to start the az/el compass 

94 fontsize = 20 # for the compass labels 

95 compassSize = 500 

96 textDistance = 650 

97 compassCenter = (leftFraction * width, leftFraction * height) 

98 compassAzElCent = ((1 - rightFraction) * width, rightFraction * height) 

99 

100 vi = exp.getInfo().getVisitInfo() 

101 az, _ = vi.boresightAzAlt 

102 _, dec = vi.boresightRaDec 

103 rotpa = vi.boresightRotAngle 

104 

105 az = Angle(az.asDegrees(), u.deg) 

106 dec = Angle(dec.asDegrees(), u.deg) 

107 rotpa = Angle(rotpa.asDegrees(), u.deg) 

108 

109 if icSrc: 

110 plt.scatter(icSrc["base_SdssCentroid_x"], icSrc["base_SdssCentroid_y"], color="red", marker="x") 

111 markerStyle: dict[str, Any] = dict() 

112 if filteredSources: 

113 markerStyle.update( 

114 dict( 

115 marker="o", 

116 linestyle="", 

117 markersize=20, 

118 linewidth=10, 

119 color="green", 

120 markeredgecolor="green", 

121 fillstyle="none", 

122 ) 

123 ) 

124 plt.plot( 

125 filteredSources["base_SdssCentroid_x"], filteredSources["base_SdssCentroid_y"], **markerStyle 

126 ) 

127 

128 if np.isfinite(rotpa): # need a rotation angle for the compass 

129 plt.arrow( 

130 compassCenter[0], 

131 compassCenter[1], 

132 -compassSize * np.sin(rotpa), 

133 compassSize * np.cos(rotpa), 

134 color="green", 

135 width=20, 

136 ) 

137 plt.text( 

138 compassCenter[0] - textDistance * np.sin(rotpa), 

139 compassCenter[1] + textDistance * np.cos(rotpa), 

140 "N", 

141 color="green", 

142 fontsize=fontsize, 

143 weight="bold", 

144 ) 

145 plt.arrow( 

146 compassCenter[0], 

147 compassCenter[1], 

148 compassSize * np.cos(rotpa), 

149 compassSize * np.sin(rotpa), 

150 color="green", 

151 width=20, 

152 ) 

153 plt.text( 

154 compassCenter[0] + textDistance * np.cos(rotpa), 

155 compassCenter[1] + textDistance * np.sin(rotpa), 

156 "E", 

157 color="green", 

158 fontsize=fontsize, 

159 weight="bold", 

160 ) 

161 

162 sinTheta = np.cos(AUXTEL_LOCATION.lat) / np.cos(dec) * np.sin(az) 

163 theta = Angle(np.arcsin(sinTheta)) 

164 rotAzEl = rotpa - theta - Angle(90.0 * u.deg) 

165 if np.isfinite(rotAzEl): # need a rotation angle for the compass 

166 plt.arrow( 

167 compassAzElCent[0], 

168 compassAzElCent[1], 

169 -compassSize * np.sin(rotAzEl), 

170 compassSize * np.cos(rotAzEl), 

171 color="cyan", 

172 width=20, 

173 ) 

174 plt.text( 

175 compassAzElCent[0] - textDistance * np.sin(rotAzEl), 

176 compassAzElCent[1] + textDistance * np.cos(rotAzEl), 

177 "EL", 

178 color="cyan", 

179 fontsize=fontsize, 

180 weight="bold", 

181 ) 

182 plt.arrow( 

183 compassAzElCent[0], 

184 compassAzElCent[1], 

185 compassSize * np.cos(rotAzEl), 

186 compassSize * np.sin(rotAzEl), 

187 color="cyan", 

188 width=20, 

189 ) 

190 plt.text( 

191 compassAzElCent[0] + textDistance * np.cos(rotAzEl), 

192 compassAzElCent[1] + textDistance * np.sin(rotAzEl), 

193 "AZ", 

194 color="cyan", 

195 fontsize=fontsize, 

196 weight="bold", 

197 ) 

198 

199 plt.ylim(0, height) 

200 plt.tight_layout() 

201 

202 if saveAs: 

203 fig.savefig(saveAs) 

204 plt.show() 

205 

206 plt.gcf().clear() 

207 del fig 

208 del data