Coverage for src/ufig/plugins/match_sextractor_catalog_multiband_read.py: 93%

1# Copyright (c) 2015 ETH Zurich, Institute of Astronomy, Claudio Bruderer 

2# <claudio.bruderer@phys.ethz.ch> 

3 

4""" 

5Created on Apr 12, 2015 

6@author: Claudio Bruderer 

7 

8""" 

9 

10from functools import reduce 

11 

12import numpy as np 

13from cosmic_toolbox import arraytools as at 

14from cosmic_toolbox import logger 

15from ivy.plugin.base_plugin import BasePlugin 

16from sklearn.neighbors import BallTree 

17 

18from ufig.plugins.match_sextractor_seg_catalog_multiband_read import NO_MATCH_VAL 

19 

20LOGGER = logger.get_logger(__file__) 

21 

22HDF5_COMPRESS = {"compression": "gzip", "compression_opts": 9, "shuffle": True} 

23 

24 

25def match(x_in, y_in, mag_in, x_out, y_out, mag_out, par): 

26 """ 

27 Match UFig input positions and magnitudes to SExtractor output positions and 

28 magnitudes. 

29 

30 :param x_in: Input x-coordinates (in SExtractor convention) 

31 :param y_in: Input y-coordinates (in SExtractor convention) 

32 :param mag_in: Input magnitudes 

33 :param x_out: Output x-coordinates 

34 :param y_out: Output y-coordinates 

35 :param mag_out: Output magnitudes 

36 :param par.max_radius: Maximum distance between an input and a detected objects to 

37 classify as a match. 

38 :param par.mag_diff: Maximum magnitude difference between an input and a detected 

39 object to classify as a match. 

40 

41 :return: Indices linking detected to input objects. 

42 """ 

43 

44 # use ball tree to find input objects close enough to detected objects 

45 tree = BallTree(np.stack((x_in, y_in), axis=-1), leaf_size=1000, metric="euclidean") 

46 inds_dist_match = tree.query_radius( 

47 np.stack((x_out, y_out), axis=-1), 

48 par.max_radius, 

49 return_distance=True, 

50 sort_results=True, 

51 )[0] 

52 

53 indices_link = np.full_like(x_out, NO_MATCH_VAL, dtype=int) 

54 

55 # compare input and output magnitudes, set closest object with sufficiently small 

56 # magnitude difference as match 

57 for c, ind_dist_match in enumerate(inds_dist_match): 

58 delta_mag = np.fabs(mag_out[c] - mag_in[ind_dist_match]) 

59 select_mag = delta_mag < par.mag_diff 

60 

61 if np.any(select_mag): 

62 indices_link[c] = ind_dist_match[select_mag][0] 

63 

64 return indices_link 

65 

66 

67class Plugin(BasePlugin): 

68 """ 

69 Matches a Source Extractor catalog to the input catalogs and adds all the input 

70 columns to the Source Extractor catalog. 

71 """ 

72 

73 def get_sexcat_paths(self, f): 

74 """ 

75 Get paths of SExtractor catalogs to be matched from list of plugins. 

76 

77 :param plugin_list: List of plugins 

78 :return: List of paths of SExtractor catalogs to be matched 

79 """ 

80 

81 paths = self.ctx.parameters.sextractor_forced_photo_catalog_name_dict[f] 

82 return paths 

83 

84 def get_ucat_paths(self, f, sg): 

85 if sg == "galaxy": 

86 paths = self.ctx.parameters.galaxy_catalog_name_dict[f] 

87 elif sg == "star": 87 ↛ 90line 87 didn't jump to line 90 because the condition on line 87 was always true

88 paths = self.ctx.parameters.star_catalog_name_dict[f] 

89 

90 return paths 

91 

92 def __call__(self): 

93 par = self.ctx.parameters 

94 

95 indices_link = {} 

96 sexcat_dict = {} 

97 gal_dict = {} 

98 star_dict = {} 

99 x_out_dict = {} 

100 y_out_dict = {} 

101 mag_out_dict = {} 

102 x_in = np.array([], float) 

103 y_in = np.array([], float) 

104 mag_in_dict = {} 

105 path_ucat = {} 

106 path_sexcat = {} 

107 indices_link_b = {} 

108 

109 for f in par.filters: 

110 if "galaxies" in self.ctx: 110 ↛ 117line 110 didn't jump to line 117 because the condition on line 110 was always true

111 path_ucat[f] = self.get_ucat_paths(f=f, sg="galaxy") 

112 gal_dict[f] = at.load_hdf(path_ucat[f]) 

113 # match 

114 x_in = gal_dict[f]["x"] 

115 y_in = gal_dict[f]["y"] 

116 mag_in_dict[f] = gal_dict[f]["mag"] 

117 if "stars" in self.ctx: 117 ↛ 124line 117 didn't jump to line 124 because the condition on line 117 was always true

118 path_ucat[f] = self.get_ucat_paths(f=f, sg="star") 

119 star_dict[f] = at.load_hdf(path_ucat[f]) 

120 x_in = np.append(x_in, star_dict[f]["x"]) 

121 y_in = np.append(y_in, star_dict[f]["y"]) 

122 mag_in_dict[f] = np.append(mag_in_dict[f], star_dict[f]["mag"]) 

123 # 0.5-shift due to different pixel conventions in UFig and Sextractor 

124 x_in += 0.5 

125 y_in += 0.5 

126 path_sexcat[f] = self.get_sexcat_paths(f=f) 

127 # load catalog 

128 sexcat_dict[f] = at.load_hdf_cols(path_sexcat[f]) 

129 x_out_dict[f] = sexcat_dict[f][par.matching_x] 

130 y_out_dict[f] = sexcat_dict[f][par.matching_y] 

131 mag_out_dict[f] = sexcat_dict[f][par.matching_mag] 

132 

133 indices_link_b[f] = match( 

134 x_in, 

135 y_in, 

136 mag_in_dict[f], 

137 x_out_dict[f], 

138 y_out_dict[f], 

139 mag_out_dict[f], 

140 par, 

141 ) 

142 

143 indices_link = reduce( 

144 lambda x, y: np.where(x == y, x, NO_MATCH_VAL), indices_link_b.values() 

145 ) 

146 

147 # check that there is no multiple counting 

148 inds, counts = np.unique( 

149 indices_link[indices_link != NO_MATCH_VAL], return_counts=True 

150 ) 

151 multi = sum(np.where(counts != 1, True, False)) 

152 LOGGER.info(f"There are {multi} matches to the same index") 

153 LOGGER.info( 

154 f"This means {sum(counts[counts != 1])} galaxies instead of {multi}" 

155 ) 

156 

157 # split matches into stars and galaxies 

158 galaxy_mask = (indices_link < self.ctx.numgalaxies) & ( 

159 indices_link != NO_MATCH_VAL 

160 ) 

161 star_mask = (indices_link >= self.ctx.numgalaxies) & ( 

162 indices_link != NO_MATCH_VAL 

163 ) 

164 galaxy_indices = indices_link[galaxy_mask] 

165 star_indices = indices_link[star_mask] - self.ctx.numgalaxies 

166 

167 def append_(obj, col): 

168 new_columns.append(f"{col}:{obj[col].dtype.str}") 

169 new_column_shapes.append(obj[col].shape[1:]) 

170 new_columns_names.append(col) 

171 

172 for f in par.filters: 

173 # add columns 

174 new_columns = [] 

175 new_column_shapes = [] 

176 new_columns_names = [] 

177 

178 for col in gal_dict[f].dtype.names: 

179 if col not in sexcat_dict[f].dtype.names: 179 ↛ 178line 179 didn't jump to line 178 because the condition on line 179 was always true

180 append_(gal_dict[f], col) 

181 

182 for col in star_dict[f].dtype.names: 

183 if ( 183 ↛ 187line 183 didn't jump to line 187

184 col not in sexcat_dict[f].dtype.names 

185 and col not in gal_dict[f].dtype.names 

186 ): 

187 append_(star_dict[f], col) 

188 

189 if "star_gal" not in sexcat_dict[f].dtype.names: 189 ↛ 194line 189 didn't jump to line 194 because the condition on line 189 was always true

190 new_columns.append(f"star_gal:{np.array([NO_MATCH_VAL]).dtype.str}") 

191 new_column_shapes.append(()) 

192 new_columns_names.append("star_gal") 

193 

194 sexcat_dict[f] = at.add_cols( 

195 sexcat_dict[f], new_columns, shapes=new_column_shapes, data=NO_MATCH_VAL 

196 ) 

197 

198 # add matched values 

199 if "star_gal" in new_columns_names: 199 ↛ 204line 199 didn't jump to line 204 because the condition on line 199 was always true

200 sexcat_dict[f]["star_gal"][galaxy_mask] = 0 

201 sexcat_dict[f]["star_gal"][star_mask] = 1 

202 new_columns_names.remove("star_gal") 

203 

204 for new_col in new_columns_names: 

205 if new_col in gal_dict[f].dtype.names: 205 ↛ 210line 205 didn't jump to line 210 because the condition on line 205 was always true

206 sexcat_dict[f][new_col][galaxy_mask] = gal_dict[f][new_col][ 

207 galaxy_indices 

208 ] 

209 

210 if new_col in star_dict[f].dtype.names: 210 ↛ 204line 210 didn't jump to line 204 because the condition on line 210 was always true

211 sexcat_dict[f][new_col][star_mask] = star_dict[f][new_col][ 

212 star_indices 

213 ] 

214 

215 # write extended catalog to disk 

216 at.save_hdf_cols(path_sexcat[f], sexcat_dict[f], compression=HDF5_COMPRESS) 

217 

218 def __str__(self): 

219 return "match catalog to input"