LSSTApplications  20.0.0
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 << ", 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  _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
lsst::afw::image
Backwards-compatibility support for depersisting the old Calib (FluxMag0/FluxMag0Err) objects.
Definition: imageAlgorithm.dox:1
LOG_LOGGER
#define LOG_LOGGER
Definition: Log.h:703
std::tan
T tan(T... args)
std::string
STL class.
std::shared_ptr< lsst::afw::geom::SkyWcs >
lsst.pipe.tasks.assembleCoadd.filter
filter
Definition: assembleCoadd.py:335
CcdImage.h
std::move
T move(T... args)
std::atan2
T atan2(T... args)
AstrometryTransform.h
std::pair< int, int >
lsst::jointcal::AstrometryTransformLinear
implements the linear transformations (6 real coefficients).
Definition: AstrometryTransform.h:426
wcs
table::Key< table::Array< std::uint8_t > > wcs
Definition: SkyWcs.cc:71
std::cos
T cos(T... args)
lsst::jointcal::TanRaDecToPixel
This one is the Tangent Plane (called gnomonic) projection (from celestial sphere to tangent plane)
Definition: AstrometryTransform.h:707
dec
double dec
Definition: Match.cc:41
std::stringstream
STL class.
LOG_GET
#define LOG_GET(logger)
Definition: Log.h:75
Image.h
crval
table::PointKey< double > crval
Definition: OldWcs.cc:130
lsst::jointcal::TanPixelToRaDec
The transformation that handles pixels to sideral transformations (Gnomonic, possibly with polynomial...
Definition: AstrometryTransform.h:627
CameraSys.h
std::sqrt
T sqrt(T... args)
std::list::clear
T clear(T... args)
lsst::jointcal::MeasuredStarList::setCcdImage
void setCcdImage(const CcdImage *_ccdImage)
Definition: MeasuredStar.cc:68
std::list::push_back
T push_back(T... args)
LOGLS_WARN
#define LOGLS_WARN(logger, message)
Definition: Log.h:648
std::isnan
T isnan(T... args)
Angle.h
image
afw::table::Key< afw::table::Array< ImagePixelT > > image
Definition: HeavyFootprint.cc:216
Point.h
std::ostream
STL class.
lsst::jointcal::operator<<
std::ostream & operator<<(std::ostream &stream, AstrometryMapping const &mapping)
Definition: AstrometryMapping.h:78
lsst::afw::table._source.SourceCatalog
Definition: _source.py:33
lsst::daf::base::DateTime::MJD
@ MJD
Definition: DateTime.h:66
lsst::jointcal
Definition: Associations.h:49
lsst
A base class for image defects.
Definition: imageAlgorithm.dox:1
std::sin
T sin(T... args)
photoCalib
Key< int > photoCalib
Definition: Exposure.cc:67
lsst::afw::cameraGeom::FOCAL_PLANE
CameraSys const FOCAL_PLANE
Focal plane coordinates: Position on a 2-d planar approximation to the focal plane (x,...
Definition: CameraSys.cc:30
PhotoCalib.h
Implementation of the Photometric Calibration class.
lsst::jointcal::CcdImage::countStars
std::pair< int, int > countStars() const
Count the number of valid measured and reference stars that fall within this ccdImage.
Definition: CcdImage.cc:159
detector
table::Key< int > detector
Definition: DetectorCollection.cc:172
key
Key< U > key
Definition: Schema.cc:281
lsst::geom::Point< double, 2 >
lsst::geom::SpherePoint::getLatitude
Angle getLatitude() const noexcept
The latitude of this point.
Definition: SpherePoint.h:190
lsst::geom::Box2I
An integer coordinate rectangle.
Definition: Box.h:55
lsst::jointcal::Point
A point in a plane.
Definition: Point.h:36
transform
table::Key< int > transform
Definition: TransformMap.cc:299
std::stringstream::str
T str(T... args)
std::make_pair
T make_pair(T... args)
Point.h
lsst::jointcal::CcdImageKey
For hashing a ccdImage: the pair of (visit, ccd) IDs should be unique to each ccdImage.
Definition: CcdImage.h:51
latitude
table::Key< lsst::geom::Angle > latitude
Definition: VisitInfo.cc:171
lsst::jointcal::compose
std::unique_ptr< AstrometryTransform > compose(AstrometryTransform const &left, AstrometryTransform const &right)
Returns a pointer to a composition of transforms, representing left(right()).
Definition: AstrometryTransform.cc:426
lsst.pipe.tasks.assembleCoadd.visit
visit
Definition: assembleCoadd.py:336
lsst::afw::cameraGeom::PIXELS
CameraSysPrefix const PIXELS
Pixel coordinates: Nominal position on the entry surface of a given detector (x, y unbinned pixels).
Definition: CameraSys.cc:34
lsst::jointcal::CcdImage::setCommonTangentPoint
void setCommonTangentPoint(Point const &commonTangentPoint)
Sets the common tangent point and computes necessary transforms.
Definition: CcdImage.cc:174
Log.h
LSST DM logging module built on log4cxx.
visitInfo
Key< int > visitInfo
Definition: Exposure.cc:70
lsst::jointcal::CcdImage::CcdImage
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
exceptions.h
lsst::jointcal::Frame
rectangle with sides parallel to axes.
Definition: Frame.h:38
bbox
AmpInfoBoxKey bbox
Definition: Amplifier.cc:117
std::pow
T pow(T... args)