Coverage for python / lsst / meas / astrom / verifyWcs.py: 0%

1# 

2# LSST Data Management System 

3# Copyright 2008, 2009, 2010 LSST Corporation. 

4# 

5# This product includes software developed by the 

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

7# 

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

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

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

11# (at your option) any later version. 

12# 

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

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

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

16# GNU General Public License for more details. 

17# 

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

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

20# see <http://www.lsstcorp.org/LegalNotices/>. 

21# 

22 

23__all__ = ["checkMatches"] 

24 

25 

26import numpy as np 

27import logging 

28 

29import lsst.geom 

30import lsst.afw.detection as afwDetection 

31import lsst.afw.math as afwMath 

32import lsst.meas.algorithms as measAlg 

33 

34_LOG = logging.getLogger(__name__) 

35 

36 

37def checkMatches(srcMatchSet, exposure, log=None): 

38 """Check astrometric matches and assess Wcs quality by computing statics 

39 over spacial cells in the image. 

40 

41 Parameters 

42 ---------- 

43 srcMatchSet : `list` of `lsst.afw.table.ReferenceMatch` 

44 List of matched sources to a reference catalog. 

45 exposure : `lsst.afw.image.Exposure` 

46 Image the sources in srcMatchSet were detected/measured in. 

47 log : `lsst.log.Log` or `logging.Logger` 

48 Logger object. 

49 

50 Returns 

51 ------- 

52 values : `dict` 

53 Result dictionary with fields: 

54 

55 - ``minObjectsPerCell`` : (`int`) 

56 - ``maxObjectsPerCell`` : (`int`) 

57 - ``meanObjectsPerCell`` : (`float`) 

58 - ``stdObjectsPerCell`` : (`float`) 

59 """ 

60 if not exposure: 

61 return {} 

62 

63 if log is None: 

64 log = _LOG 

65 

66 im = exposure.getMaskedImage().getImage() 

67 width, height = im.getWidth(), im.getHeight() 

68 nx, ny = 3, 3 

69 w, h = width//nx, height//ny 

70 

71 if w == 0: 

72 w = 1 

73 while nx*w < width: 

74 w += 1 

75 

76 if h == 0: 

77 h = 1 

78 while ny*h < height: 

79 h += 1 

80 

81 cellSet = afwMath.SpatialCellSet( 

82 lsst.geom.Box2I(lsst.geom.Point2I(0, 0), lsst.geom.Extent2I(width, height)), w, h) 

83 # 

84 # Populate cellSet 

85 # 

86 i = -1 

87 for m in srcMatchSet: 

88 i += 1 

89 

90 src = m.second 

91 csrc = afwDetection.Source() 

92 csrc.setId(i) 

93 csrc.setXAstrom(src.getXAstrom()) 

94 csrc.setYAstrom(src.getYAstrom()) 

95 

96 try: 

97 cellSet.insertCandidate(measAlg.PsfCandidateF(csrc, exposure.getMaskedImage())) 

98 except Exception as e: 

99 log.warning("%s", e) 

100 

101 ncell = len(cellSet.getCellList()) 

102 nobj = np.ndarray(ncell, dtype='i') 

103 

104 for i in range(ncell): 

105 cell = cellSet.getCellList()[i] 

106 

107 nobj[i] = cell.size() 

108 

109 dx = np.ndarray(cell.size()) 

110 dy = np.ndarray(cell.size()) 

111 

112 j = 0 

113 for cand in cell: 

114 # 

115 # Swig doesn't know that we're a SpatialCellImageCandidate; all it knows is that we have 

116 # a SpatialCellCandidate so we need an explicit (dynamic) cast 

117 # 

118 mid = cand.getSource().getId() 

119 dx[j] = srcMatchSet[mid].first.getXAstrom() - srcMatchSet[mid].second.getXAstrom() 

120 dy[j] = srcMatchSet[mid].first.getYAstrom() - srcMatchSet[mid].second.getYAstrom() 

121 

122 j += 1 

123 

124 log.debug("%s %-30s %8s dx,dy = %5.2f,%5.2f rms_x,y = %5.2f,%5.2f", 

125 cell.getLabel(), cell.getBBox(), ("nobj=%d" % cell.size()), 

126 dx.mean(), dy.mean(), dx.std(), dy.std()) 

127 

128 nobj.sort() 

129 

130 values = {} 

131 values["minObjectsPerCell"] = int(nobj[0]) # otherwise it's a numpy integral type 

132 values["maxObjectsPerCell"] = int(nobj[-1]) 

133 values["meanObjectsPerCell"] = nobj.mean() 

134 values["stdObjectsPerCell"] = nobj.std() 

135 

136 return values