LSST Applications 27.0.0,g0265f82a02+469cd937ee,g02d81e74bb+21ad69e7e1,g1470d8bcf6+cbe83ee85a,g2079a07aa2+e67c6346a6,g212a7c68fe+04a9158687,g2305ad1205+94392ce272,g295015adf3+81dd352a9d,g2bbee38e9b+469cd937ee,g337abbeb29+469cd937ee,g3939d97d7f+72a9f7b576,g487adcacf7+71499e7cba,g50ff169b8f+5929b3527e,g52b1c1532d+a6fc98d2e7,g591dd9f2cf+df404f777f,g5a732f18d5+be83d3ecdb,g64a986408d+21ad69e7e1,g858d7b2824+21ad69e7e1,g8a8a8dda67+a6fc98d2e7,g99cad8db69+f62e5b0af5,g9ddcbc5298+d4bad12328,ga1e77700b3+9c366c4306,ga8c6da7877+71e4819109,gb0e22166c9+25ba2f69a1,gb6a65358fc+469cd937ee,gbb8dafda3b+69d3c0e320,gc07e1c2157+a98bf949bb,gc120e1dc64+615ec43309,gc28159a63d+469cd937ee,gcf0d15dbbd+72a9f7b576,gdaeeff99f8+a38ce5ea23,ge6526c86ff+3a7c1ac5f1,ge79ae78c31+469cd937ee,gee10cc3b42+a6fc98d2e7,gf1cff7945b+21ad69e7e1,gfbcc870c63+9a11dc8c8f
LSST Data Management Base Package
Loading...
Searching...
No Matches
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 <cassert>
26#include <cmath>
27#include <sstream>
28#include <string>
29#include <utility>
30
32#include "lsst/pex/exceptions.h"
35#include "lsst/geom/Angle.h"
36#include "lsst/geom/Point.h"
37
38#include "lsst/log/Log.h"
41#include "lsst/jointcal/Point.h"
42
43namespace jointcal = lsst::jointcal;
44namespace afwImg = lsst::afw::image;
45
46namespace {
47LOG_LOGGER _log = LOG_GET("lsst.jointcal.CcdImage");
48}
49
50namespace lsst {
51namespace jointcal {
52
54 out << "(visit: " << key.visit << ", detector: " << key.ccd << ")";
55 return out;
56}
57
58void CcdImage::loadCatalog(afw::table::SourceCatalog const &catalog, std::string const &fluxField) {
59 auto xKey = catalog.getSchema().find<double>("slot_Centroid_x").key;
60 auto yKey = catalog.getSchema().find<double>("slot_Centroid_y").key;
61 auto xsKey = catalog.getSchema().find<float>("slot_Centroid_xErr").key;
62 auto ysKey = catalog.getSchema().find<float>("slot_Centroid_yErr").key;
63 auto mxxKey = catalog.getSchema().find<double>("slot_Shape_xx").key;
64 auto myyKey = catalog.getSchema().find<double>("slot_Shape_yy").key;
65 auto mxyKey = catalog.getSchema().find<double>("slot_Shape_xy").key;
66 auto instFluxKey = catalog.getSchema().find<double>(fluxField + "_instFlux").key;
67 auto instFluxErrKey = catalog.getSchema().find<double>(fluxField + "_instFluxErr").key;
68
69 auto transform = _detector->getTransform(afw::cameraGeom::PIXELS, afw::cameraGeom::FOCAL_PLANE);
70
71 _wholeCatalog.clear();
72 for (auto const &record : catalog) {
73 auto ms = std::make_shared<MeasuredStar>();
74 ms->setId(record.getId());
75 ms->x = record.get(xKey);
76 ms->y = record.get(yKey);
77 ms->vx = std::pow(record.get(xsKey), 2);
78 ms->vy = std::pow(record.get(ysKey), 2);
79 auto pointFocal = transform->applyForward(record.getCentroid());
80 ms->setXFocal(pointFocal.getX());
81 ms->setYFocal(pointFocal.getY());
82 /* the xy covariance is not provided in the input catalog: we
83 cook it up from the x and y position variance and the shape
84 measurements: */
85 double mxx = record.get(mxxKey);
86 double myy = record.get(myyKey);
87 double mxy = record.get(mxyKey);
88 ms->vxy = mxy * (ms->vx + ms->vy) / (mxx + myy);
89 if (std::isnan(ms->vxy) || ms->vx < 0 || ms->vy < 0 || (ms->vxy * ms->vxy) > (ms->vx * ms->vy)) {
90 LOGLS_WARN(_log, "Bad source detected during loadCatalog, "
91 << "detector=" << _ccdId << ", visit=" << _visit
92 << " id: " << ms->getId() << " with vx,vy: " << ms->vx << "," << ms->vy
93 << " vxy^2: " << ms->vxy * ms->vxy << " vx*vy: " << ms->vx * ms->vy);
94 continue;
95 }
96 ms->setInstFluxAndErr(record.get(instFluxKey), record.get(instFluxErrKey));
97 // TODO: the below lines will be less clumsy once DM-4044 is cleaned up and we can say:
98 // TODO: instFluxToNanojansky(ms->getInstFlux(), ms) (because ms will be derived from
99 // geom::Point).
100 geom::Point<double, 2> point(ms->x, ms->y);
101 auto flux = _photoCalib->instFluxToNanojansky(ms->getInstFlux(), ms->getInstFluxErr(), point);
102 ms->setFlux(flux.value);
103 ms->setFluxErr(flux.error);
104 auto mag = _photoCalib->instFluxToMagnitude(ms->getInstFlux(), ms->getInstFluxErr(), point);
105 ms->getMag() = mag.value;
106 ms->setMagErr(mag.error);
107 ms->setCcdImage(this);
108 _wholeCatalog.push_back(std::move(ms));
109 }
110 _wholeCatalog.setCcdImage(this);
111}
112
116 std::shared_ptr<afw::cameraGeom::Detector> detector, int visit, int ccdId,
117 std::string const &fluxField)
118 : _ccdId(ccdId), _visit(visit), _photoCalib(std::move(photoCalib)), _detector(std::move(detector)), _filter(std::move(filter)) {
119 loadCatalog(catalog, fluxField);
120
121 Point lowerLeft(bbox.getMinX(), bbox.getMinY());
122 Point upperRight(bbox.getMaxX(), bbox.getMaxY());
123 _imageFrame = Frame(lowerLeft, upperRight);
124
125 _readWcs = std::make_shared<AstrometryTransformSkyWcs>(wcs);
126
128 out << visit << "_" << ccdId;
129 _name = out.str();
130
131 _boresightRaDec = visitInfo->getBoresightRaDec();
132 _airMass = visitInfo->getBoresightAirmass();
133 _epoch = visitInfo->getDate().get(lsst::daf::base::DateTime::EPOCH);
134 _lstObs = visitInfo->getEra();
135 _hourAngle = visitInfo->getBoresightHourAngle();
136
137 // Some cameras don't manage ERA (and thus Hour Angle) properly, so it's going to be NaN.
138 // Because we need the refraction vector later, go with 0 HA to prevent crashes on that NaN.
139 if (std::isnan(_hourAngle) == true) {
140 _hourAngle = 0;
141 }
142
143 if (_airMass == 1)
144 _sinEta = _cosEta = _tanZ = 0;
145 else {
146 double cosz = 1. / _airMass;
147 double sinz = std::sqrt(1 - cosz * cosz); // astronomers usually observe above the horizon
148 _tanZ = sinz / cosz;
149 lsst::geom::Angle eta = visitInfo->getBoresightParAngle();
150 _sinEta = std::sin(eta.asRadians());
151 _cosEta = std::cos(eta.asRadians());
152 }
153}
154
156 int measuredStars = 0;
157 int refStars = 0;
158 for (auto const &measuredStar : _catalogForFit) {
159 if (measuredStar->isValid()) {
160 measuredStars++;
161 }
162 if ((measuredStar->getFittedStar() != nullptr) &&
163 (measuredStar->getFittedStar()->getRefStar() != nullptr)) {
164 refStars++;
165 }
166 }
167 return std::make_pair(measuredStars, refStars);
168}
169
170void CcdImage::setCommonTangentPoint(Point const &commonTangentPoint) {
171 _commonTangentPoint = commonTangentPoint;
172
173 auto const crval = _readWcs->getSkyWcs()->getSkyOrigin();
174 jointcal::Point tangentPoint(crval[0].asDegrees(), crval[1].asDegrees());
175
176 /* we don't assume here that we know the internals of TanPixelToRaDec:
177 to construct pix->TP, we do pix->sky->TP, although pix->sky
178 actually goes through TP */
180 TanRaDecToPixel raDecToTangentPlane(identity, tangentPoint);
181 _pixelToTangentPlane = compose(raDecToTangentPlane, *_readWcs);
182 TanPixelToRaDec CommonTangentPlane2RaDec(identity, commonTangentPoint);
183 _commonTangentPlaneToTangentPlane = compose(raDecToTangentPlane, CommonTangentPlane2RaDec);
184
185 // jump from one TP to an other:
186 TanRaDecToPixel raDecToCommonTangentPlane(identity, commonTangentPoint);
187 TanPixelToRaDec TangentPlaneToRaDec(identity, tangentPoint);
188 _tangentPlaneToCommonTangentPlane = compose(raDecToCommonTangentPlane, TangentPlaneToRaDec);
189 _skyToTangentPlane.reset(new TanRaDecToPixel(identity, tangentPoint));
190
191 // this one is needed for matches :
192 _pixelToCommonTangentPlane = compose(raDecToCommonTangentPlane, *_readWcs);
193}
194} // namespace jointcal
195} // namespace lsst
AmpInfoBoxKey bbox
Definition Amplifier.cc:117
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.
A class representing an angle.
Definition Angle.h:128
constexpr double asRadians() const noexcept
Return an Angle's value in radians.
Definition Angle.h:173
An integer coordinate rectangle.
Definition Box.h:55
A coordinate class intended to represent absolute positions.
Definition Point.h:169
implements the linear transformations (6 real coefficients).
void setCommonTangentPoint(Point const &commonTangentPoint)
Sets the common tangent point and computes necessary transforms.
Definition CcdImage.cc:170
std::pair< int, int > countStars() const
Count the number of valid measured and reference stars that fall within this ccdImage.
Definition CcdImage.cc:155
CcdImage(afw::table::SourceCatalog &record, const std::shared_ptr< lsst::afw::geom::SkyWcs > &wcs, const std::shared_ptr< lsst::afw::image::VisitInfo > &visitInfo, lsst::geom::Box2I const &bbox, std::string 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:113
rectangle with sides parallel to axes.
Definition Frame.h:38
void setCcdImage(const CcdImage *_ccdImage)
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)
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)
STL namespace.
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