LSSTApplications  17.0+124,17.0+14,17.0+73,18.0.0+37,18.0.0+80,18.0.0-4-g68ffd23+4,18.1.0-1-g0001055+12,18.1.0-1-g03d53ef+5,18.1.0-1-g1349e88+55,18.1.0-1-g2505f39+44,18.1.0-1-g5315e5e+4,18.1.0-1-g5e4b7ea+14,18.1.0-1-g7e8fceb+4,18.1.0-1-g85f8cd4+48,18.1.0-1-g8ff0b9f+4,18.1.0-1-ga2c679d+1,18.1.0-1-gd55f500+35,18.1.0-10-gb58edde+2,18.1.0-11-g0997b02+4,18.1.0-13-gfe4edf0b+12,18.1.0-14-g259bd21+21,18.1.0-19-gdb69f3f+2,18.1.0-2-g5f9922c+24,18.1.0-2-gd3b74e5+11,18.1.0-2-gfbf3545+32,18.1.0-26-g728bddb4+5,18.1.0-27-g6ff7ca9+2,18.1.0-3-g52aa583+25,18.1.0-3-g8ea57af+9,18.1.0-3-gb69f684+42,18.1.0-3-gfcaddf3+6,18.1.0-32-gd8786685a,18.1.0-4-gf3f9b77+6,18.1.0-5-g1dd662b+2,18.1.0-5-g6dbcb01+41,18.1.0-6-gae77429+3,18.1.0-7-g9d75d83+9,18.1.0-7-gae09a6d+30,18.1.0-9-gc381ef5+4,w.2019.45
LSSTDataManagementBasePackage
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 << ", ccd: " << 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  _mjd = visitInfo->getDate().get(lsst::daf::base::DateTime::MJD);
132  double latitude = visitInfo->getObservatory().getLatitude();
133  _lstObs = visitInfo->getEra();
134  _hourAngle = visitInfo->getBoresightHourAngle();
135 
136  // Some cameras don't manage ERA (and thus Hour Angle) properly, so it's going to be NaN.
137  // Because we need the refraction vector later, go with 0 HA to prevent crashes on that NaN.
138  if (std::isnan(_hourAngle) == true) {
139  _hourAngle = 0;
140  }
141 
142  if (_airMass == 1)
143  _sinEta = _cosEta = _tanZ = 0;
144  else {
145  double cosz = 1. / _airMass;
146  double sinz = std::sqrt(1 - cosz * cosz); // astronomers usually observe above the horizon
147  _tanZ = sinz / cosz;
148  // TODO: as part of DM-12473, we can remove all of this and just call _visitInfo.getParallacticAngle()
149  double dec = _boresightRaDec.getLatitude();
150  // x/y components of refraction angle, eta.]
151  double yEta = std::sin(_hourAngle);
152  double xEta = std::cos(dec) * std::tan(latitude) - std::sin(dec) * std::cos(_hourAngle);
153  double eta = std::atan2(yEta, xEta);
154  _sinEta = std::sin(eta);
155  _cosEta = std::cos(eta);
156  }
157 }
158 
160  int measuredStars = 0;
161  int refStars = 0;
162  for (auto const &measuredStar : _catalogForFit) {
163  if (measuredStar->isValid()) {
164  measuredStars++;
165  }
166  if ((measuredStar->getFittedStar() != nullptr) &&
167  (measuredStar->getFittedStar()->getRefStar() != nullptr)) {
168  refStars++;
169  }
170  }
171  return std::make_pair(measuredStars, refStars);
172 }
173 
174 void CcdImage::setCommonTangentPoint(Point const &commonTangentPoint) {
175  _commonTangentPoint = commonTangentPoint;
176 
177  auto const crval = _readWcs->getSkyWcs()->getSkyOrigin();
178  jointcal::Point tangentPoint(crval[0].asDegrees(), crval[1].asDegrees());
179 
180  /* we don't assume here that we know the internals of TanPixelToRaDec:
181  to construct pix->TP, we do pix->sky->TP, although pix->sky
182  actually goes through TP */
183  AstrometryTransformLinear identity;
184  TanRaDecToPixel raDecToTangentPlane(identity, tangentPoint);
185  _pixelToTangentPlane = compose(raDecToTangentPlane, *_readWcs);
186  TanPixelToRaDec CommonTangentPlane2RaDec(identity, commonTangentPoint);
187  _commonTangentPlaneToTangentPlane = compose(raDecToTangentPlane, CommonTangentPlane2RaDec);
188 
189  // jump from one TP to an other:
190  TanRaDecToPixel raDecToCommonTangentPlane(identity, commonTangentPoint);
191  TanPixelToRaDec TangentPlaneToRaDec(identity, tangentPoint);
192  _tangentPlaneToCommonTangentPlane = compose(raDecToCommonTangentPlane, TangentPlaneToRaDec);
193  _skyToTangentPlane.reset(new TanRaDecToPixel(identity, tangentPoint));
194 
195  // this one is needed for matches :
196  _pixelToCommonTangentPlane = compose(raDecToCommonTangentPlane, *_readWcs);
197 }
198 } // namespace jointcal
199 } // namespace lsst
AmpInfoBoxKey bbox
Definition: Amplifier.cc:117
#define LOGLS_WARN(logger, message)
Log a warn-level message using an iostream-based interface.
Definition: Log.h:648
The transformation that handles pixels to sideral transformations (Gnomonic, possibly with polynomial...
double dec
Definition: Match.cc:41
table::Key< lsst::geom::Angle > latitude
Definition: VisitInfo.cc:171
#define LOG_LOGGER
Definition: Log.h:703
afw::table::Key< afw::table::Array< ImagePixelT > > image
implements the linear transformations (6 real coefficients).
A point in a plane.
Definition: Point.h:36
CameraSysPrefix const PIXELS
Pixel coordinates: Nominal position on the entry surface of a given detector (x, y unbinned pixels)...
Definition: CameraSys.cc:34
For hashing a ccdImage: the pair of (visit, ccd) IDs should be unique to each ccdImage.
Definition: CcdImage.h:51
table::PointKey< double > crval
Definition: OldWcs.cc:130
table::Key< table::Array< std::uint8_t > > wcs
Definition: SkyWcs.cc:71
std::ostream & operator<<(std::ostream &os, CameraSysPrefix const &detSysPrefix)
Definition: CameraSys.cc:47
T atan2(T... args)
STL class.
LSST DM logging module built on log4cxx.
T sin(T... args)
rectangle with sides parallel to axes.
Definition: Frame.h:38
A base class for image defects.
int getMaxY() const noexcept
Definition: Box.h:162
std::unique_ptr< AstrometryTransform > compose(AstrometryTransform const &left, AstrometryTransform const &right)
Returns a pointer to a composition of transforms, representing left(right()).
Key< U > key
Definition: Schema.cc:281
Angle getLatitude() const noexcept
The latitude of this point.
Definition: SpherePoint.h:190
T make_pair(T... args)
This one is the Tangent Plane (called gnomonic) projection (from celestial sphere to tangent plane) ...
T cos(T... args)
table::Key< int > detector
T move(T... args)
int getMaxX() const noexcept
Definition: Box.h:161
T get(T... args)
int getMinX() const noexcept
Definition: Box.h:157
T pow(T... args)
T tan(T... args)
T isnan(T... args)
T sqrt(T... args)
Backwards-compatibility support for depersisting the old Calib (FluxMag0/FluxMag0Err) objects...
#define LOG_GET(logger)
Returns a Log object associated with logger.
Definition: Log.h:75
An integer coordinate rectangle.
Definition: Box.h:55
CameraSys const FOCAL_PLANE
Focal plane coordinates: Position on a 2-d planar approximation to the focal plane (x...
Definition: CameraSys.cc:30
STL class.
Implementation of the Photometric Calibration class.
int getMinY() const noexcept
Definition: Box.h:158
void setCommonTangentPoint(Point const &commonTangentPoint)
Sets the common tangent point and computes necessary transforms.
Definition: CcdImage.cc:174
std::pair< int, int > countStars() const
Count the number of valid measured and reference stars that fall within this ccdImage.
Definition: CcdImage.cc:159