LSST Applications  21.0.0-147-g0e635eb1+1acddb5be5,22.0.0+052faf71bd,22.0.0+1ea9a8b2b2,22.0.0+6312710a6c,22.0.0+729191ecac,22.0.0+7589c3a021,22.0.0+9f079a9461,22.0.1-1-g7d6de66+b8044ec9de,22.0.1-1-g87000a6+536b1ee016,22.0.1-1-g8e32f31+6312710a6c,22.0.1-10-gd060f87+016f7cdc03,22.0.1-12-g9c3108e+df145f6f68,22.0.1-16-g314fa6d+c825727ab8,22.0.1-19-g93a5c75+d23f2fb6d8,22.0.1-19-gb93eaa13+aab3ef7709,22.0.1-2-g8ef0a89+b8044ec9de,22.0.1-2-g92698f7+9f079a9461,22.0.1-2-ga9b0f51+052faf71bd,22.0.1-2-gac51dbf+052faf71bd,22.0.1-2-gb66926d+6312710a6c,22.0.1-2-gcb770ba+09e3807989,22.0.1-20-g32debb5+b8044ec9de,22.0.1-23-gc2439a9a+fb0756638e,22.0.1-3-g496fd5d+09117f784f,22.0.1-3-g59f966b+1e6ba2c031,22.0.1-3-g849a1b8+f8b568069f,22.0.1-3-gaaec9c0+c5c846a8b1,22.0.1-32-g5ddfab5d3+60ce4897b0,22.0.1-4-g037fbe1+64e601228d,22.0.1-4-g8623105+b8044ec9de,22.0.1-5-g096abc9+d18c45d440,22.0.1-5-g15c806e+57f5c03693,22.0.1-7-gba73697+57f5c03693,master-g6e05de7fdc+c1283a92b8,master-g72cdda8301+729191ecac,w.2021.39
LSST Data Management Base Package
CcdImage.cc
Go to the documentation of this file.
1 // -*- LSST-C++ -*-
2 /*
3  * This file is part of jointcal.
4  *
5  * Developed for the LSST Data Management System.
6  * This product includes software developed by the LSST Project
7  * (https://www.lsst.org).
8  * See the COPYRIGHT file at the top-level directory of this distribution
9  * for details of code ownership.
10  *
11  * This program is free software: you can redistribute it and/or modify
12  * it under the terms of the GNU General Public License as published by
13  * the Free Software Foundation, either version 3 of the License, or
14  * (at your option) any later version.
15  *
16  * This program is distributed in the hope that it will be useful,
17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
19  * GNU General Public License for more details.
20  *
21  * You should have received a copy of the GNU General Public License
22  * along with this program. If not, see <https://www.gnu.org/licenses/>.
23  */
24 
25 #include <assert.h>
26 #include <string>
27 #include <sstream>
28 #include <cmath>
29 
31 #include "lsst/pex/exceptions.h"
32 #include "lsst/afw/image/Image.h"
34 #include "lsst/geom/Angle.h"
35 #include "lsst/geom/Point.h"
36 
37 #include "lsst/log/Log.h"
38 #include "lsst/jointcal/CcdImage.h"
40 #include "lsst/jointcal/Point.h"
41 
42 namespace jointcal = lsst::jointcal;
43 namespace afwImg = lsst::afw::image;
44 
45 namespace {
46 LOG_LOGGER _log = LOG_GET("jointcal.CcdImage");
47 }
48 
49 namespace lsst {
50 namespace jointcal {
51 
53  out << "(visit: " << key.visit << ", detector: " << key.ccd << ")";
54  return out;
55 }
56 
57 void CcdImage::loadCatalog(afw::table::SourceCatalog const &catalog, std::string const &fluxField) {
58  auto xKey = catalog.getSchema().find<double>("slot_Centroid_x").key;
59  auto yKey = catalog.getSchema().find<double>("slot_Centroid_y").key;
60  auto xsKey = catalog.getSchema().find<float>("slot_Centroid_xErr").key;
61  auto ysKey = catalog.getSchema().find<float>("slot_Centroid_yErr").key;
62  auto mxxKey = catalog.getSchema().find<double>("slot_Shape_xx").key;
63  auto myyKey = catalog.getSchema().find<double>("slot_Shape_yy").key;
64  auto mxyKey = catalog.getSchema().find<double>("slot_Shape_xy").key;
65  auto instFluxKey = catalog.getSchema().find<double>(fluxField + "_instFlux").key;
66  auto instFluxErrKey = catalog.getSchema().find<double>(fluxField + "_instFluxErr").key;
67 
68  auto transform = _detector->getTransform(afw::cameraGeom::PIXELS, afw::cameraGeom::FOCAL_PLANE);
69 
70  _wholeCatalog.clear();
71  for (auto const &record : catalog) {
72  auto ms = std::make_shared<MeasuredStar>();
73  ms->setId(record.getId());
74  ms->x = record.get(xKey);
75  ms->y = record.get(yKey);
76  ms->vx = std::pow(record.get(xsKey), 2);
77  ms->vy = std::pow(record.get(ysKey), 2);
78  auto pointFocal = transform->applyForward(record.getCentroid());
79  ms->setXFocal(pointFocal.getX());
80  ms->setYFocal(pointFocal.getY());
81  /* the xy covariance is not provided in the input catalog: we
82  cook it up from the x and y position variance and the shape
83  measurements: */
84  double mxx = record.get(mxxKey);
85  double myy = record.get(myyKey);
86  double mxy = record.get(mxyKey);
87  ms->vxy = mxy * (ms->vx + ms->vy) / (mxx + myy);
88  if (std::isnan(ms->vxy) || ms->vx < 0 || ms->vy < 0 || (ms->vxy * ms->vxy) > (ms->vx * ms->vy)) {
89  LOGLS_WARN(_log, "Bad source detected during loadCatalog id: "
90  << ms->getId() << " with vx,vy: " << ms->vx << "," << ms->vy
91  << " vxy^2: " << ms->vxy * ms->vxy << " vx*vy: " << ms->vx * ms->vy);
92  continue;
93  }
94  ms->setInstFluxAndErr(record.get(instFluxKey), record.get(instFluxErrKey));
95  // TODO: the below lines will be less clumsy once DM-4044 is cleaned up and we can say:
96  // TODO: instFluxToNanojansky(ms->getInstFlux(), ms) (because ms will be derived from
97  // geom::Point).
98  geom::Point<double, 2> point(ms->x, ms->y);
99  auto flux = _photoCalib->instFluxToNanojansky(ms->getInstFlux(), ms->getInstFluxErr(), point);
100  ms->setFlux(flux.value);
101  ms->setFluxErr(flux.error);
102  auto mag = _photoCalib->instFluxToMagnitude(ms->getInstFlux(), ms->getInstFluxErr(), point);
103  ms->getMag() = mag.value;
104  ms->setMagErr(mag.error);
105  ms->setCcdImage(this);
106  _wholeCatalog.push_back(std::move(ms));
107  }
108  _wholeCatalog.setCcdImage(this);
109 }
110 
115  std::string const &fluxField)
116  : _ccdId(ccdId), _visit(visit), _photoCalib(photoCalib), _detector(detector), _filter(filter) {
117  loadCatalog(catalog, fluxField);
118 
119  Point lowerLeft(bbox.getMinX(), bbox.getMinY());
120  Point upperRight(bbox.getMaxX(), bbox.getMaxY());
121  _imageFrame = Frame(lowerLeft, upperRight);
122 
123  _readWcs = std::make_shared<AstrometryTransformSkyWcs>(wcs);
124 
125  std::stringstream out;
126  out << visit << "_" << ccdId;
127  _name = out.str();
128 
129  _boresightRaDec = visitInfo->getBoresightRaDec();
130  _airMass = visitInfo->getBoresightAirmass();
131  _epoch = visitInfo->getDate().get(lsst::daf::base::DateTime::EPOCH);
132  _lstObs = visitInfo->getEra();
133  _hourAngle = visitInfo->getBoresightHourAngle();
134 
135  // Some cameras don't manage ERA (and thus Hour Angle) properly, so it's going to be NaN.
136  // Because we need the refraction vector later, go with 0 HA to prevent crashes on that NaN.
137  if (std::isnan(_hourAngle) == true) {
138  _hourAngle = 0;
139  }
140 
141  if (_airMass == 1)
142  _sinEta = _cosEta = _tanZ = 0;
143  else {
144  double cosz = 1. / _airMass;
145  double sinz = std::sqrt(1 - cosz * cosz); // astronomers usually observe above the horizon
146  _tanZ = sinz / cosz;
147  lsst::geom::Angle eta = visitInfo->getBoresightParAngle();
148  _sinEta = std::sin(eta.asRadians());
149  _cosEta = std::cos(eta.asRadians());
150  }
151 }
152 
154  int measuredStars = 0;
155  int refStars = 0;
156  for (auto const &measuredStar : _catalogForFit) {
157  if (measuredStar->isValid()) {
158  measuredStars++;
159  }
160  if ((measuredStar->getFittedStar() != nullptr) &&
161  (measuredStar->getFittedStar()->getRefStar() != nullptr)) {
162  refStars++;
163  }
164  }
165  return std::make_pair(measuredStars, refStars);
166 }
167 
168 void CcdImage::setCommonTangentPoint(Point const &commonTangentPoint) {
169  _commonTangentPoint = commonTangentPoint;
170 
171  auto const crval = _readWcs->getSkyWcs()->getSkyOrigin();
172  jointcal::Point tangentPoint(crval[0].asDegrees(), crval[1].asDegrees());
173 
174  /* we don't assume here that we know the internals of TanPixelToRaDec:
175  to construct pix->TP, we do pix->sky->TP, although pix->sky
176  actually goes through TP */
177  AstrometryTransformLinear identity;
178  TanRaDecToPixel raDecToTangentPlane(identity, tangentPoint);
179  _pixelToTangentPlane = compose(raDecToTangentPlane, *_readWcs);
180  TanPixelToRaDec CommonTangentPlane2RaDec(identity, commonTangentPoint);
181  _commonTangentPlaneToTangentPlane = compose(raDecToTangentPlane, CommonTangentPlane2RaDec);
182 
183  // jump from one TP to an other:
184  TanRaDecToPixel raDecToCommonTangentPlane(identity, commonTangentPoint);
185  TanPixelToRaDec TangentPlaneToRaDec(identity, tangentPoint);
186  _tangentPlaneToCommonTangentPlane = compose(raDecToCommonTangentPlane, TangentPlaneToRaDec);
187  _skyToTangentPlane.reset(new TanRaDecToPixel(identity, tangentPoint));
188 
189  // this one is needed for matches :
190  _pixelToCommonTangentPlane = compose(raDecToCommonTangentPlane, *_readWcs);
191 }
192 } // namespace jointcal
193 } // namespace lsst
AmpInfoBoxKey bbox
Definition: Amplifier.cc:117
table::Key< int > detector
afw::table::Key< afw::table::Array< ImagePixelT > > image
LSST DM logging module built on log4cxx.
#define LOGLS_WARN(logger, message)
Log a warn-level message using an iostream-based interface.
Definition: Log.h:659
#define LOG_GET(logger)
Returns a Log object associated with logger.
Definition: Log.h:75
#define LOG_LOGGER
Definition: Log.h:714
table::PointKey< double > crval
Definition: OldWcs.cc:128
Implementation of the Photometric Calibration class.
table::Key< table::Array< std::uint8_t > > wcs
Definition: SkyWcs.cc:66
table::Key< int > transform
A class representing an angle.
Definition: Angle.h:127
constexpr double asRadians() const noexcept
Return an Angle's value in radians.
Definition: Angle.h:166
An integer coordinate rectangle.
Definition: Box.h:55
implements the linear transformations (6 real coefficients).
void setCommonTangentPoint(Point const &commonTangentPoint)
Sets the common tangent point and computes necessary transforms.
Definition: CcdImage.cc:168
CcdImage(afw::table::SourceCatalog &record, std::shared_ptr< lsst::afw::geom::SkyWcs > wcs, std::shared_ptr< lsst::afw::image::VisitInfo > visitInfo, lsst::geom::Box2I const &bbox, std::string const &filter, std::shared_ptr< afw::image::PhotoCalib > photoCalib, std::shared_ptr< afw::cameraGeom::Detector > detector, int visit, int ccd, std::string const &fluxField)
Definition: CcdImage.cc:111
std::pair< int, int > countStars() const
Count the number of valid measured and reference stars that fall within this ccdImage.
Definition: CcdImage.cc:153
rectangle with sides parallel to axes.
Definition: Frame.h:38
void setCcdImage(const CcdImage *_ccdImage)
Definition: MeasuredStar.cc:68
A point in a plane.
Definition: Point.h:37
The transformation that handles pixels to sideral transformations (Gnomonic, possibly with polynomial...
This one is the Tangent Plane (called gnomonic) projection (from celestial sphere to tangent plane)
T clear(T... args)
T cos(T... args)
T isnan(T... args)
T make_pair(T... args)
T move(T... args)
CameraSys const FOCAL_PLANE
Focal plane coordinates: Position on a 2-d planar approximation to the focal plane (x,...
Definition: CameraSys.cc:30
CameraSysPrefix const PIXELS
Pixel coordinates: Nominal position on the entry surface of a given detector (x, y unbinned pixels).
Definition: CameraSys.cc:34
Backwards-compatibility support for depersisting the old Calib (FluxMag0/FluxMag0Err) objects.
std::unique_ptr< AstrometryTransform > compose(AstrometryTransform const &left, AstrometryTransform const &right)
Returns a pointer to a composition of transforms, representing left(right()).
std::ostream & operator<<(std::ostream &stream, AstrometryMapping const &mapping)
A base class for image defects.
T pow(T... args)
T push_back(T... args)
T sin(T... args)
T sqrt(T... args)
T str(T... args)
For hashing a ccdImage: the pair of (visit, ccd) IDs should be unique to each ccdImage.
Definition: CcdImage.h:51
Key< int > visitInfo
Definition: Exposure.cc:70
Key< int > photoCalib
Definition: Exposure.cc:67