LSST Applications g04e9c324dd+8c5ae1fdc5,g134cb467dc+1b3060144d,g18429d2f64+f642bf4753,g199a45376c+0ba108daf9,g1fd858c14a+2dcf163641,g262e1987ae+7b8c96d2ca,g29ae962dfc+3bd6ecb08a,g2cef7863aa+aef1011c0b,g35bb328faa+8c5ae1fdc5,g3fd5ace14f+53e1a9e7c5,g4595892280+fef73a337f,g47891489e3+2efcf17695,g4d44eb3520+642b70b07e,g53246c7159+8c5ae1fdc5,g67b6fd64d1+2efcf17695,g67fd3c3899+b70e05ef52,g74acd417e5+317eb4c7d4,g786e29fd12+668abc6043,g87389fa792+8856018cbb,g89139ef638+2efcf17695,g8d7436a09f+3be3c13596,g8ea07a8fe4+9f5ccc88ac,g90f42f885a+a4e7b16d9b,g97be763408+ad77d7208f,g9dd6db0277+b70e05ef52,ga681d05dcb+a3f46e7fff,gabf8522325+735880ea63,gac2eed3f23+2efcf17695,gb89ab40317+2efcf17695,gbf99507273+8c5ae1fdc5,gd8ff7fe66e+b70e05ef52,gdab6d2f7ff+317eb4c7d4,gdc713202bf+b70e05ef52,gdfd2d52018+b10e285e0f,ge365c994fd+310e8507c4,ge410e46f29+2efcf17695,geaed405ab2+562b3308c0,gffca2db377+8c5ae1fdc5,w.2025.35
LSST Data Management Base Package
Loading...
Searching...
No Matches
SkyWcs.cc
Go to the documentation of this file.
1/*
2 * LSST Data Management System
3 * Copyright 2017 AURA/LSST.
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#include <cmath>
24#include <cstdint>
25#include <exception>
26#include <memory>
27#include <vector>
28
29#include "astshim.h"
30
31#include "lsst/geom/Angle.h"
32#include "lsst/geom/Point.h"
41#include "lsst/pex/exceptions.h"
43
44namespace lsst {
45namespace afw {
46
47template std::shared_ptr<geom::SkyWcs> table::io::PersistableFacade<geom::SkyWcs>::dynamicCast(
48 std::shared_ptr<table::io::Persistable> const&);
49
50namespace geom {
51namespace {
52
53int const SERIALIZATION_VERSION = 1;
54
55// TIGHT_FITS_TOL is used by getFitsMetadata to determine if a WCS can accurately be represented as a FITS
56// WCS. It specifies the maximum departure from linearity (in pixels) allowed on either axis of the mapping
57// from pixel coordinates to Intermediate World Coordinates over a range of 100 x 100 pixels
58// (or the region specified by NAXIS[12], if provided, but we do not pass this to AST as of 2018-01-17).
59// For more information,
60// see FitsTol in the AST manual http://starlink.eao.hawaii.edu/devdocs/sun211.htx/sun211.html
61double const TIGHT_FITS_TOL = 0.0001;
62
63class SkyWcsPersistenceHelper {
64public:
68
69 // No copying
70 SkyWcsPersistenceHelper(const SkyWcsPersistenceHelper&) = delete;
71 SkyWcsPersistenceHelper& operator=(const SkyWcsPersistenceHelper&) = delete;
72
73 // No moving
74 SkyWcsPersistenceHelper(SkyWcsPersistenceHelper&&) = delete;
75 SkyWcsPersistenceHelper& operator=(SkyWcsPersistenceHelper&&) = delete;
76
77 explicit SkyWcsPersistenceHelper(bool hasFitsApproximation)
78 : schema(),
79 wcs(schema.addField<table::Array<std::uint8_t>>("wcs", "wcs string representation", "")) {
80 if (hasFitsApproximation) {
81 approx = schema.addField<table::Array<std::uint8_t>>(
82 "approx", "wcs string representation of FITS approximation", "");
83 }
84 }
85
86 explicit SkyWcsPersistenceHelper(table::Schema const& schema)
87 : schema(schema), wcs(schema["wcs"]), approx() {
88 try {
89 approx = schema["approx"];
90 } catch (pex::exceptions::NotFoundError&) {
91 }
92 }
93};
94
95class SkyWcsFactory : public table::io::PersistableFactory {
96public:
97 explicit SkyWcsFactory(std::string const& name) : table::io::PersistableFactory(name) {}
98
99 std::shared_ptr<table::io::Persistable> read(InputArchive const& archive,
100 CatalogVector const& catalogs) const override {
101 LSST_ARCHIVE_ASSERT(catalogs.size() == 1u);
102 LSST_ARCHIVE_ASSERT(catalogs.front().size() == 1u);
103 SkyWcsPersistenceHelper keys(catalogs.front().getSchema());
104 table::BaseRecord const& record = catalogs.front().front();
105 std::string stringRep = formatters::bytesToString(record.get(keys.wcs));
106 auto result = SkyWcs::readString(stringRep);
107 if (keys.approx.isValid()) {
108 auto bytes = record.get(keys.approx);
109 if (!bytes.isEmpty()) {
110 auto approxStringRep = formatters::bytesToString(bytes);
111 result = result->copyWithFitsApproximation(SkyWcs::readString(approxStringRep));
112 }
113 }
114 return result;
115 }
116};
117
118std::string getSkyWcsPersistenceName() { return "SkyWcs"; }
119
120SkyWcsFactory registration(getSkyWcsPersistenceName());
121
122ast::FrameDict makeSkyWcsFrameDict(TransformPoint2ToPoint2 const& pixelsToFieldAngle,
123 lsst::geom::Angle const& orientation, bool flipX,
124 lsst::geom::SpherePoint const& crval,
125 std::string const& projection = "TAN") {
126 auto const orientationAndFlipXMatrix = makeCdMatrix(1 * lsst::geom::degrees, orientation, flipX);
127 auto const initialWcs =
128 makeSkyWcs(lsst::geom::Point2D(0, 0), crval, orientationAndFlipXMatrix, projection);
129 auto const initialFrameDict = initialWcs->getFrameDict();
130 auto const iwcToSkyMap = initialFrameDict->getMapping("IWC", "SKY");
131 auto const pixelFrame = initialFrameDict->getFrame("PIXELS");
132 auto const iwcFrame = initialFrameDict->getFrame("IWC");
133 auto const skyFrame = initialFrameDict->getFrame("SKY");
134 // Field angle is in radians and is aligned to focal plane x and y;
135 // IWC is basically the same thing, but in degrees and with rotation and flipX applied
136 ndarray::Array<double, 2, 2> fieldAngleToIwcNdArray = ndarray::allocate(2, 2);
137 asEigenMatrix(fieldAngleToIwcNdArray) = orientationAndFlipXMatrix * 180.0 / lsst::geom::PI;
138 auto const pixelsToFieldAngleMap = pixelsToFieldAngle.getMapping();
139 auto const fieldAngleToIwcMap = ast::MatrixMap(fieldAngleToIwcNdArray);
140 auto const pixelsToIwcMap = pixelsToFieldAngleMap->then(fieldAngleToIwcMap);
141 auto finalFrameDict = ast::FrameDict(*pixelFrame, pixelsToIwcMap, *iwcFrame);
142 finalFrameDict.addFrame("IWC", *iwcToSkyMap, *skyFrame);
143 return finalFrameDict;
144}
145
146} // namespace
147
148Eigen::Matrix2d makeCdMatrix(lsst::geom::Angle const& scale, lsst::geom::Angle const& orientation,
149 bool flipX) {
150 Eigen::Matrix2d cdMatrix;
151 double orientRad = orientation.asRadians();
152 double scaleDeg = scale.asDegrees();
153 double xmult = flipX ? 1.0 : -1.0;
154 cdMatrix(0, 0) = std::cos(orientRad) * scaleDeg * xmult;
155 cdMatrix(0, 1) = std::sin(orientRad) * scaleDeg;
156 cdMatrix(1, 0) = -std::sin(orientRad) * scaleDeg * xmult;
157 cdMatrix(1, 1) = std::cos(orientRad) * scaleDeg;
158 return cdMatrix;
159}
160
162 auto const dstInverse = dst.getTransform()->inverted();
163 return src.getTransform()->then(*dstInverse);
164}
165
167 : SkyWcs(detail::readLsstSkyWcs(metadata, strip)) {}
168
169SkyWcs::SkyWcs(ast::FrameDict const& frameDict) : SkyWcs(_checkFrameDict(frameDict)) {}
170
171bool SkyWcs::operator==(SkyWcs const& other) const { return writeString() == other.writeString(); }
172
174 // Compute pixVec containing the pixel position and two nearby points
175 // (use a vector so all three points can be converted to sky in a single call)
176 double const side = 1.0;
178 pixel,
179 pixel + lsst::geom::Extent2D(side, 0),
180 pixel + lsst::geom::Extent2D(0, side),
181 };
182
183 auto skyVec = pixelToSky(pixVec);
184
185 // Work in 3-space to avoid RA wrapping and pole issues
186 auto skyLL = skyVec[0].getVector();
187 auto skyDx = skyVec[1].getVector() - skyLL;
188 auto skyDy = skyVec[2].getVector() - skyLL;
189
190 // Compute pixel scale in radians = sqrt(pixel area in radians^2)
191 // pixel area in radians^2 = area of parallelogram with sides skyDx, skyDy = |skyDx cross skyDy|
192 // Use squared norm to avoid two square roots
193 double skyAreaSq = skyDx.cross(skyDy).getSquaredNorm();
194 return (std::pow(skyAreaSq, 0.25) / side) * lsst::geom::radians;
195}
196
198 // CRVAL is stored as the SkyRef property of the sky frame (the current frame of the SkyWcs)
200 getFrameDict()->getFrame(ast::FrameDict::CURRENT, false)); // false: do not copy
201 if (!skyFrame) {
202 throw LSST_EXCEPT(pex::exceptions::LogicError, "Current frame is not a SkyFrame");
203 }
204 auto const crvalRad = skyFrame->getSkyRef();
205 return lsst::geom::SpherePoint(crvalRad[0] * lsst::geom::radians, crvalRad[1] * lsst::geom::radians);
206}
207
208Eigen::Matrix2d SkyWcs::getCdMatrix(lsst::geom::Point2D const& pixel) const {
209 auto const pixelToIwc = getFrameDict()->getMapping(ast::FrameSet::BASE, "IWC");
210 auto const pixelToIwcTransform = TransformPoint2ToPoint2(*pixelToIwc);
211 return pixelToIwcTransform.getJacobian(pixel);
212}
213
214Eigen::Matrix2d SkyWcs::getCdMatrix() const { return getCdMatrix(getPixelOrigin()); }
215
217 auto const crval = pixelToSky(pixel);
218 auto const cdMatrix = getCdMatrix(pixel);
219 auto metadata = makeSimpleWcsMetadata(pixel, crval, cdMatrix);
220 return std::make_shared<SkyWcs>(*metadata);
221}
222
224 auto newToOldPixel = TransformPoint2ToPoint2(ast::ShiftMap({-shift[0], -shift[1]}));
225 return makeModifiedWcs(newToOldPixel, *this, true);
226}
227
228std::shared_ptr<daf::base::PropertyList> SkyWcs::_getDirectFitsMetadata() const {
229 // Make a FrameSet that maps from GRID to SKY; GRID = the base frame (PIXELS or ACTUAL_PIXELS) + 1
230 auto const gridToPixel = ast::ShiftMap({-1.0, -1.0});
231 auto thisDict = getFrameDict();
232 auto const pixelToIwc = thisDict->getMapping(ast::FrameSet::BASE, "IWC");
233 auto const iwcToSky = thisDict->getMapping("IWC", "SKY");
234 auto const gridToSky = gridToPixel.then(*pixelToIwc).then(*iwcToSky);
235 ast::FrameSet frameSet(ast::Frame(2, "Domain=GRID"), gridToSky, *thisDict->getFrame("SKY", false));
236
237 // Write frameSet to a FitsChan and extract the metadata
239 os << "Encoding=FITS-WCS, CDMatrix=1, FitsAxisOrder=<copy>, FitsTol=" << TIGHT_FITS_TOL;
240 ast::StringStream strStream;
241 ast::FitsChan fitsChan(strStream, os.str());
242 int const nObjectsWritten = fitsChan.write(frameSet);
243 if (nObjectsWritten == 0) {
244 return nullptr;
245 }
247
248 // Remove DATE-OBS, MJD-OBS: AST writes these if the EQUINOX is set, but we set them via other mechanisms.
249 header->remove("DATE-OBS");
250 header->remove("MJD-OBS");
251
252 // If CD matrix is present, explicitly set any missing entries to zero, as a convenience to the user
253 bool const hasCd11 = header->exists("CD1_1");
254 bool const hasCd12 = header->exists("CD1_2");
255 bool const hasCd21 = header->exists("CD2_1");
256 bool const hasCd22 = header->exists("CD2_2");
257 if (hasCd11 || hasCd12 || hasCd21 || hasCd22) {
258 if (!hasCd11) header->set("CD1_1", 0.0, "Transformation matrix element");
259 if (!hasCd12) header->set("CD1_2", 0.0, "Transformation matrix element");
260 if (!hasCd21) header->set("CD2_1", 0.0, "Transformation matrix element");
261 if (!hasCd22) header->set("CD2_2", 0.0, "Transformation matrix element");
262 }
263
264 return header;
265}
266
268 if (!precise && hasFitsApproximation()) {
269 return getFitsApproximation()->_getDirectFitsMetadata();
270 }
271 auto result = _getDirectFitsMetadata();
272 if (!result) {
274 precise ? "WCS is not directly FITS-compatible."
275 : "WCS does not have an attached FITS approximation.");
276 }
277 return result;
278}
279
281
283 if (fitsApproximation && fitsApproximation->hasFitsApproximation()) {
285 "Cannot add a FITS approximation that itself already has a FITS approximation.");
286 }
287 auto result = std::make_shared<SkyWcs>(*this);
288 result->_fitsApproximation = fitsApproximation;
289 return result;
290}
291
292bool SkyWcs::isFits() const { return bool(_getDirectFitsMetadata()); }
293
295 lsst::geom::AngleUnit const& skyUnit) const {
296 return _linearizePixelToSky(skyToPixel(coord), coord, skyUnit);
297}
299 lsst::geom::AngleUnit const& skyUnit) const {
300 return _linearizePixelToSky(pix, pixelToSky(pix), skyUnit);
301}
302
304 lsst::geom::AngleUnit const& skyUnit) const {
305 return _linearizeSkyToPixel(skyToPixel(coord), coord, skyUnit);
306}
307
309 lsst::geom::AngleUnit const& skyUnit) const {
310 return _linearizeSkyToPixel(pix, pixelToSky(pix), skyUnit);
311}
312
314
315bool SkyWcs::isFlipped() const {
316 double det = getCdMatrix().determinant();
317 if (det == 0) {
318 throw(LSST_EXCEPT(pex::exceptions::RuntimeError, "CD matrix is singular"));
319 }
320 return (det > 0);
321}
322
324 int version;
325 is >> version;
326 if (version != 1) {
327 throw LSST_EXCEPT(pex::exceptions::TypeError, "Unsupported version " + std::to_string(version));
328 }
329 std::string shortClassName;
330 is >> shortClassName;
331 if (shortClassName != SkyWcs::getShortClassName()) {
333 os << "Class name in stream " << shortClassName << " != " << SkyWcs::getShortClassName();
335 }
338 auto astStream = ast::Stream(&is, nullptr);
339 auto astObjectPtr = ast::Channel(astStream).read();
340 auto frameSet = std::dynamic_pointer_cast<ast::FrameSet>(astObjectPtr);
341 if (!frameSet) {
343 os << "The AST serialization was read as a " << astObjectPtr->getClassName()
344 << " instead of a FrameSet";
346 }
347 ast::FrameDict frameDict(*frameSet);
348 return std::make_shared<SkyWcs>(frameDict);
349}
350
355
357 os << SERIALIZATION_VERSION << " " << SkyWcs::getShortClassName() << " " << hasFitsApproximation() << " ";
358 getFrameDict()->show(os, false); // false = do not write comments
359}
360
363 writeStream(os);
364 return os.str();
365}
366
368 return std::make_unique<SkyWcs>(*this);
369}
370
373 if (isFits()) {
374 os << "FITS standard SkyWcs:";
375 } else {
376 os << "Non-standard SkyWcs (Frames: ";
377 // Print the frames in index order (frames are numbered from 1).
378 std::string delimiter = "";
379 for (size_t i = 1; i <= getFrameDict()->getAllDomains().size(); ++i) {
380 os << delimiter << getFrameDict()->getFrame(i)->getDomain();
381 delimiter = ", ";
382 }
383 os << "): ";
384 }
385 std::string delimiter = "\n";
386 os << delimiter << "Sky Origin: " << getSkyOrigin();
387 os << delimiter << "Pixel Origin: " << getPixelOrigin();
388 os << delimiter << "Pixel Scale: " << getPixelScale().asArcseconds() << " arcsec/pixel";
389 return os.str();
390}
391
392bool SkyWcs::equals(typehandling::Storable const& other) const noexcept {
393 return singleClassEquals(*this, other);
394}
395
396std::string SkyWcs::getPersistenceName() const { return getSkyWcsPersistenceName(); }
397
398std::string SkyWcs::getPythonModule() const { return "lsst.afw.geom"; }
399
401 SkyWcsPersistenceHelper const keys(hasFitsApproximation());
402 table::BaseCatalog cat = handle.makeCatalog(keys.schema);
404 record->set(keys.wcs, formatters::stringToBytes(writeString()));
405 if (hasFitsApproximation()) {
406 record->set(keys.approx, formatters::stringToBytes(getFitsApproximation()->writeString()));
407 }
408 handle.saveCatalog(cat);
409}
410
412 : _frameDict(frameDict), _transform(), _pixelOrigin(), _pixelScaleAtOrigin(0 * lsst::geom::radians) {
413 _computeCache();
414};
415
416std::shared_ptr<ast::FrameDict> SkyWcs::_checkFrameDict(ast::FrameDict const& frameDict) const {
417 // Check that each frame is present and has the right type and number of axes
418 std::vector<std::string> const domainNames = {"ACTUAL_PIXELS", "PIXELS", "IWC", "SKY"};
419 for (auto const& domainName : domainNames) {
420 if (frameDict.hasDomain(domainName)) {
421 auto const frame = frameDict.getFrame(domainName, false);
422 if (frame->getNAxes() != 2) {
424 "Frame " + domainName + " has " + std::to_string(frame->getNAxes()) +
425 " axes instead of 2");
426 }
427 auto desiredClassName = domainName == "SKY" ? "SkyFrame" : "Frame";
428 if (frame->getClassName() != desiredClassName) {
430 "Frame " + domainName + " is of type " + frame->getClassName() +
431 " instead of " + desiredClassName);
432 }
433 } else if (domainName != "ACTUAL_PIXELS") {
434 // This is a required frame
435 throw LSST_EXCEPT(lsst::pex::exceptions::TypeError,
436 "No frame with domain " + domainName + " found");
437 }
438 }
439
440 // The base frame must have domain "PIXELS" or "ACTUAL_PIXELS"
441 auto baseDomain = frameDict.getFrame(ast::FrameSet::BASE, false)->getDomain();
442 if (baseDomain != "ACTUAL_PIXELS" && baseDomain != "PIXELS") {
443 throw LSST_EXCEPT(lsst::pex::exceptions::TypeError,
444 "Base frame has domain " + baseDomain + " instead of PIXELS or ACTUAL_PIXELS");
445 }
446
447 // The current frame must have domain "SKY"
448 auto currentDomain = frameDict.getFrame(ast::FrameSet::CURRENT, false)->getDomain();
449 if (currentDomain != "SKY") {
450 throw LSST_EXCEPT(lsst::pex::exceptions::TypeError,
451 "Current frame has domain " + currentDomain + " instead of SKY");
452 }
453
454 return frameDict.copy();
455}
456
457lsst::geom::AffineTransform SkyWcs::_linearizePixelToSky(lsst::geom::Point2D const& pix00,
458 lsst::geom::SpherePoint const& coord,
459 lsst::geom::AngleUnit const& skyUnit) const {
460 // Figure out the (0, 0), (0, 1), and (1, 0) ra/dec coordinates of the corners
461 // of a square drawn in pixel. It'd be better to center the square at sky00,
462 // but that would involve another conversion between sky and pixel coordinates
463 const double side = 1.0; // length of the square's sides in pixels
464 auto const sky00 = coord.getPosition(skyUnit);
465 auto const dsky10 = coord.getTangentPlaneOffset(pixelToSky(pix00 + lsst::geom::Extent2D(side, 0)));
466 auto const dsky01 = coord.getTangentPlaneOffset(pixelToSky(pix00 + lsst::geom::Extent2D(0, side)));
467
468 Eigen::Matrix2d m;
469 m(0, 0) = dsky10.first.asAngularUnits(skyUnit) / side;
470 m(0, 1) = dsky01.first.asAngularUnits(skyUnit) / side;
471 m(1, 0) = dsky10.second.asAngularUnits(skyUnit) / side;
472 m(1, 1) = dsky01.second.asAngularUnits(skyUnit) / side;
473
474 Eigen::Vector2d sky00v;
475 sky00v << sky00.getX(), sky00.getY();
476 Eigen::Vector2d pix00v;
477 pix00v << pix00.getX(), pix00.getY();
478 // return lsst::geom::AffineTransform(m, lsst::geom::Extent2D(sky00v - m * pix00v));
479 return lsst::geom::AffineTransform(m, (sky00v - m * pix00v));
480}
481
482lsst::geom::AffineTransform SkyWcs::_linearizeSkyToPixel(lsst::geom::Point2D const& pix00,
483 lsst::geom::SpherePoint const& coord,
484 lsst::geom::AngleUnit const& skyUnit) const {
485 lsst::geom::AffineTransform inverse = _linearizePixelToSky(pix00, coord, skyUnit);
486 return inverse.inverted();
487}
488
489std::shared_ptr<SkyWcs> makeFlippedWcs(SkyWcs const& wcs, bool flipLR, bool flipTB,
490 lsst::geom::Point2D const& center) {
491 double const dx = 1000; // any "reasonable" number of pixels will do
492 std::vector<double> inLL = {center[0] - dx, center[1] - dx};
493 std::vector<double> inUR = {center[0] + dx, center[1] + dx};
494 std::vector<double> outLL(inLL);
495 std::vector<double> outUR(inUR);
496 if (flipLR) {
497 outLL[0] = inUR[0];
498 outUR[0] = inLL[0];
499 }
500 if (flipTB) {
501 outLL[1] = inUR[1];
502 outUR[1] = inLL[1];
503 }
504 auto const flipPix = TransformPoint2ToPoint2(ast::WinMap(inLL, inUR, outLL, outUR));
505 return makeModifiedWcs(flipPix, wcs, true);
506}
507
509 bool modifyActualPixels) {
510 auto const pixelMapping = pixelTransform.getMapping();
511 auto oldFrameDict = wcs.getFrameDict();
512 bool const hasActualPixels = oldFrameDict->hasDomain("ACTUAL_PIXELS");
513 auto const pixelFrame = oldFrameDict->getFrame("PIXELS", false);
514 auto const iwcFrame = oldFrameDict->getFrame("IWC", false);
515 auto const skyFrame = oldFrameDict->getFrame("SKY", false);
516 auto const oldPixelToIwc = oldFrameDict->getMapping("PIXELS", "IWC");
517 auto const iwcToSky = oldFrameDict->getMapping("IWC", "SKY");
518
520 std::shared_ptr<ast::Mapping> newPixelToIwc;
521 if (hasActualPixels) {
522 auto const actualPixelFrame = oldFrameDict->getFrame("ACTUAL_PIXELS", false);
523 auto const oldActualPixelToPixels = oldFrameDict->getMapping("ACTUAL_PIXELS", "PIXELS");
524 std::shared_ptr<ast::Mapping> newActualPixelsToPixels;
525 if (modifyActualPixels) {
526 newActualPixelsToPixels = pixelMapping->then(*oldActualPixelToPixels).simplified();
527 newPixelToIwc = oldPixelToIwc;
528 } else {
529 newActualPixelsToPixels = oldActualPixelToPixels;
530 newPixelToIwc = pixelMapping->then(*oldPixelToIwc).simplified();
531 }
532 newFrameDict =
533 std::make_shared<ast::FrameDict>(*actualPixelFrame, *newActualPixelsToPixels, *pixelFrame);
534 newFrameDict->addFrame("PIXELS", *newPixelToIwc, *iwcFrame);
535 } else {
536 newPixelToIwc = pixelMapping->then(*oldPixelToIwc).simplified();
537 newFrameDict = std::make_shared<ast::FrameDict>(*pixelFrame, *newPixelToIwc, *iwcFrame);
538 }
539 newFrameDict->addFrame("IWC", *iwcToSky, *skyFrame);
540 return std::make_shared<SkyWcs>(*newFrameDict);
541}
542
544 return std::make_shared<SkyWcs>(metadata, strip);
545}
546
548 Eigen::Matrix2d const& cdMatrix, std::string const& projection) {
549 auto metadata = makeSimpleWcsMetadata(crpix, crval, cdMatrix, projection);
550 return std::make_shared<SkyWcs>(*metadata);
551}
552
554 lsst::geom::Angle const& orientation, bool flipX,
555 lsst::geom::SpherePoint const& boresight, std::string const& projection) {
556 auto frameDict = makeSkyWcsFrameDict(pixelsToFieldAngle, orientation, flipX, boresight, projection);
557 return std::make_shared<SkyWcs>(frameDict);
558}
559
561 Eigen::Matrix2d const& cdMatrix, Eigen::MatrixXd const& sipA,
562 Eigen::MatrixXd const& sipB) {
563 auto metadata = makeTanSipMetadata(crpix, crval, cdMatrix, sipA, sipB);
564 return std::make_shared<SkyWcs>(*metadata);
565}
566
568 Eigen::Matrix2d const& cdMatrix, Eigen::MatrixXd const& sipA,
569 Eigen::MatrixXd const& sipB, Eigen::MatrixXd const& sipAp,
570 Eigen::MatrixXd const& sipBp) {
571 auto metadata = makeTanSipMetadata(crpix, crval, cdMatrix, sipA, sipB, sipAp, sipBp);
572 return std::make_shared<SkyWcs>(*metadata);
573}
574
576 bool simplify) {
577 auto iwcToSky = wcs.getFrameDict()->getMapping("IWC", "SKY");
578 return std::make_shared<TransformPoint2ToSpherePoint>(*iwcToSky, simplify);
579}
580
582 auto pixelToIwc = wcs.getFrameDict()->getMapping(ast::FrameSet::BASE, "IWC");
583 return std::make_shared<TransformPoint2ToPoint2>(*pixelToIwc, simplify);
584}
585
587 os << wcs.toString();
588 return os;
589};
590
591} // namespace geom
592} // namespace afw
593} // namespace lsst
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition Exception.h:48
#define LSST_ARCHIVE_ASSERT(EXPR)
An assertion macro used to validate the structure of an InputArchive.
Definition Persistable.h:48
Channel provides input/output of AST objects.
Definition Channel.h:60
std::shared_ptr< Object > read()
Read an object from a channel.
Definition Channel.cc:52
A specialized form of Channel which reads and writes FITS header cards.
Definition FitsChan.h:202
A FrameSet whose frames can be referenced by domain name.
Definition FrameDict.h:67
bool hasDomain(std::string const &domain) const
Return True if a frame in this FrameDict has the specified domain.
Definition FrameDict.h:201
std::shared_ptr< Mapping > getMapping(int from, std::string const &to) const
Variant of FrameSet::getMapping with the second frame specified by domain.
Definition FrameDict.h:162
std::shared_ptr< FrameDict > copy() const
Return a deep copy of this object.
Definition FrameDict.h:123
std::set< std::string > getAllDomains() const
Get the domain names for all contained Frames (excluding frames with empty or defaulted domain names)...
Definition FrameDict.cc:45
std::shared_ptr< Frame > getFrame(std::string const &domain, bool copy=true) const
Variant of getFrame(int, bool) where the frame is specified by domain name.
Definition FrameDict.h:151
Frame is used to represent a coordinate system.
Definition Frame.h:157
A FrameSet consists of a set of one or more Frames (which describe coordinate systems),...
Definition FrameSet.h:99
static int constexpr CURRENT
index of current frame
Definition FrameSet.h:105
static int constexpr BASE
index of base frame
Definition FrameSet.h:104
SeriesMap then(Mapping const &next) const
Return a series compound mapping this(first(input)) containing shallow copies of the original.
Definition Mapping.cc:37
void show(std::ostream &os, bool showComments=true) const
Print a textual description the object to an ostream.
Definition Object.cc:158
ShiftMap is a linear Mapping which shifts each axis by a specified constant value.
Definition ShiftMap.h:40
A stream for ast::Channel.
Definition Stream.h:41
String-based source and sink for channels.
Definition Stream.h:180
A WinMap is a linear Mapping which transforms a rectangular window in one coordinate system into a si...
Definition WinMap.h:45
A 2-dimensional celestial WCS that transform pixels to ICRS RA/Dec, using the LSST standard for pixel...
Definition SkyWcs.h:118
std::shared_ptr< SkyWcs > copyAtShiftedPixelOrigin(lsst::geom::Extent2D const &shift) const
Return a copy of this SkyWcs with the pixel origin shifted by the specified amount.
Definition SkyWcs.cc:223
static std::shared_ptr< SkyWcs > readStream(std::istream &is)
Deserialize a SkyWcs from an input stream.
Definition SkyWcs.cc:323
std::shared_ptr< const ast::FrameDict > getFrameDict() const
Get the contained FrameDict.
Definition SkyWcs.cc:280
lsst::geom::SpherePoint pixelToSky(lsst::geom::Point2D const &pixel) const
Compute sky position(s) from pixel position(s)
Definition SkyWcs.h:334
lsst::geom::AffineTransform linearizeSkyToPixel(lsst::geom::SpherePoint const &coord, lsst::geom::AngleUnit const &skyUnit) const
Return the local linear approximation to skyToPixel at a point given in sky coordinates.
Definition SkyWcs.cc:303
std::string toString() const override
Create a string representation of this object.
Definition SkyWcs.cc:371
bool isFlipped() const
Does the WCS follow the convention of North=Up, East=Left?
Definition SkyWcs.cc:315
std::shared_ptr< typehandling::Storable > cloneStorable() const override
Create a new SkyWcs that is a copy of this one.
Definition SkyWcs.cc:367
lsst::geom::Angle getPixelScale() const
Get the pixel scale at the pixel origin.
Definition SkyWcs.h:208
std::shared_ptr< SkyWcs > getTanWcs(lsst::geom::Point2D const &pixel) const
Get a local TAN WCS approximation to this WCS at the specified pixel position.
Definition SkyWcs.cc:216
lsst::geom::Point2D getPixelOrigin() const
Get the pixel origin, in pixels, using the LSST convention.
Definition SkyWcs.h:215
void writeStream(std::ostream &os) const
Serialize this SkyWcs to an output stream.
Definition SkyWcs.cc:356
std::string writeString() const
Serialize this SkyWcs to a string, using the same format as writeStream.
Definition SkyWcs.cc:361
Eigen::Matrix2d getCdMatrix() const
Get the 2x2 CD matrix at the pixel origin.
Definition SkyWcs.cc:214
SkyWcs(SkyWcs const &)=default
lsst::geom::AffineTransform linearizePixelToSky(lsst::geom::SpherePoint const &coord, lsst::geom::AngleUnit const &skyUnit) const
Return the local linear approximation to pixelToSky at a point given in sky coordinates.
Definition SkyWcs.cc:294
std::shared_ptr< const TransformPoint2ToSpherePoint > getTransform() const
Get a TransformPoint2ToSpherePoint that transforms pixels to sky in the forward direction and sky to ...
Definition SkyWcs.h:257
lsst::geom::Point2D skyToPixel(lsst::geom::SpherePoint const &sky) const
Compute pixel position(s) from sky position(s)
Definition SkyWcs.h:349
static std::string getShortClassName()
Definition SkyWcs.cc:313
bool isFits() const
Return true getFitsMetadata(true) will succeed, false if not.
Definition SkyWcs.cc:292
bool hasFitsApproximation() const
Does this SkyWcs have an approximate SkyWcs that can be represented as standard FITS WCS?
Definition SkyWcs.h:362
std::shared_ptr< SkyWcs > getFitsApproximation() const
Return FITS SkyWcs that approximates this one.
Definition SkyWcs.h:369
std::string getPythonModule() const override
Return the fully-qualified Python module that should be imported to guarantee that its factory is reg...
Definition SkyWcs.cc:398
std::string getPersistenceName() const override
Return the unique name used to persist this object and look up its factory.
Definition SkyWcs.cc:396
lsst::geom::SpherePoint getSkyOrigin() const
Get the sky origin, the celestial fiducial point.
Definition SkyWcs.cc:197
std::shared_ptr< SkyWcs > copyWithFitsApproximation(std::shared_ptr< SkyWcs > fitsApproximation) const
Return a copy of this SkyWcs with the given FITS approximation.
Definition SkyWcs.cc:282
static std::shared_ptr< SkyWcs > readString(std::string &str)
Deserialize a SkyWcs from a string, using the same format as readStream.
Definition SkyWcs.cc:351
std::shared_ptr< daf::base::PropertyList > getFitsMetadata(bool precise=false) const
Return the WCS as FITS WCS metadata.
Definition SkyWcs.cc:267
void write(OutputArchiveHandle &handle) const override
Write the object to one or more catalogs.
Definition SkyWcs.cc:400
bool equals(typehandling::Storable const &other) const noexcept override
Compare this object to another Storable.
Definition SkyWcs.cc:392
bool operator==(SkyWcs const &other) const
Equality is based on the string representations being equal.
Definition SkyWcs.cc:171
std::shared_ptr< const ast::Mapping > getMapping() const
Get the contained mapping.
Definition Transform.h:135
Tag types used to declare specialized field types.
Definition misc.h:31
std::shared_ptr< RecordT > addNew()
Create a new record, add it to the end of the catalog, and return a pointer to it.
Definition Catalog.h:490
A class used as a handle to a particular field in a table.
Definition Key.h:53
Defines the fields and offsets for a table.
Definition Schema.h:51
void saveCatalog(BaseCatalog const &catalog)
Save a catalog in the archive.
BaseCatalog makeCatalog(Schema const &schema)
Return a new, empty catalog with the given schema.
static std::shared_ptr< T > dynamicCast(std::shared_ptr< Persistable > const &ptr)
Dynamically cast a shared_ptr.
A base class for factory classes used to reconstruct objects from records.
io::OutputArchiveHandle OutputArchiveHandle
Interface supporting iteration over heterogenous containers.
Definition Storable.h:58
static bool singleClassEquals(T const &lhs, Storable const &rhs)
Test if a Storable is of a particular class and equal to another object.
Definition Storable.h:151
Class for storing generic metadata.
Definition PropertySet.h:67
An affine coordinate transformation consisting of a linear transformation and an offset.
AffineTransform const inverted() const
Return the inverse transform.
A class representing an angle.
Definition Angle.h:128
constexpr double asArcseconds() const noexcept
Return an Angle's value in arcseconds.
Definition Angle.h:185
A class used to convert scalar POD types such as double to Angle.
Definition Angle.h:71
Point in an unspecified spherical coordinate system.
Definition SpherePoint.h:57
std::pair< Angle, Angle > getTangentPlaneOffset(SpherePoint const &other) const
Get the offset from a tangent plane centered at this point to another point.
Point2D getPosition(AngleUnit unit) const noexcept
Return longitude, latitude as a Point2D object.
Reports errors in the logical structure of the program.
Definition Runtime.h:46
Reports errors that are due to events beyond the control of the program.
Definition Runtime.h:104
Reports errors from accepting an object of an unexpected or inappropriate type.
Definition Runtime.h:167
T cos(T... args)
T make_shared(T... args)
ndarray::Array< std::uint8_t, 1, 1 > stringToBytes(std::string const &str)
Encode a std::string as a vector of uint8.
Definition Utils.cc:43
std::shared_ptr< daf::base::PropertyList > getPropertyListFromFitsChan(ast::FitsChan &fitsChan)
Copy values from an AST FitsChan into a PropertyList.
std::shared_ptr< SkyWcs > makeSkyWcs(daf::base::PropertySet &metadata, bool strip=false)
Construct a SkyWcs from FITS keywords.
Definition SkyWcs.cc:543
std::shared_ptr< daf::base::PropertyList > makeSimpleWcsMetadata(lsst::geom::Point2D const &crpix, lsst::geom::SpherePoint const &crval, Eigen::Matrix2d const &cdMatrix, std::string const &projection="TAN")
Make FITS metadata for a simple FITS WCS (one with no distortion).
Definition wcsUtils.cc:169
std::shared_ptr< TransformPoint2ToPoint2 > getPixelToIntermediateWorldCoords(SkyWcs const &wcs, bool simplify=true)
Return a transform from pixel coordinates to intermediate world coordinates.
Definition SkyWcs.cc:581
std::shared_ptr< SkyWcs > makeTanSipWcs(lsst::geom::Point2D const &crpix, lsst::geom::SpherePoint const &crval, Eigen::Matrix2d const &cdMatrix, Eigen::MatrixXd const &sipA, Eigen::MatrixXd const &sipB)
Construct a TAN-SIP SkyWcs with forward SIP distortion terms and an iterative inverse.
Definition SkyWcs.cc:560
std::ostream & operator<<(std::ostream &os, GenericEndpoint const &endpoint)
Print "GenericEndpoint(_n_)" to the ostream where _n_ is the number of axes, e.g. "GenericAxes(4)".
Definition Endpoint.cc:239
Eigen::Matrix2d makeCdMatrix(lsst::geom::Angle const &scale, lsst::geom::Angle const &orientation=0 *lsst::geom::degrees, bool flipX=false)
Make a WCS CD matrix.
Definition SkyWcs.cc:148
std::shared_ptr< TransformPoint2ToSpherePoint > getIntermediateWorldCoordsToSky(SkyWcs const &wcs, bool simplify=true)
Return a transform from intermediate world coordinates to sky.
Definition SkyWcs.cc:575
std::shared_ptr< daf::base::PropertyList > makeTanSipMetadata(lsst::geom::Point2D const &crpix, lsst::geom::SpherePoint const &crval, Eigen::Matrix2d const &cdMatrix, Eigen::MatrixXd const &sipA, Eigen::MatrixXd const &sipB)
Make metadata for a TAN-SIP WCS without inverse matrices.
Definition wcsUtils.cc:194
std::shared_ptr< SkyWcs > makeModifiedWcs(TransformPoint2ToPoint2 const &pixelTransform, SkyWcs const &wcs, bool modifyActualPixels)
Create a new SkyWcs whose pixels are transformed by pixelTransform, as described below.
Definition SkyWcs.cc:508
Transform< Point2Endpoint, Point2Endpoint > TransformPoint2ToPoint2
Definition Transform.h:300
std::shared_ptr< TransformPoint2ToPoint2 > makeWcsPairTransform(SkyWcs const &src, SkyWcs const &dst)
A Transform obtained by putting two SkyWcs objects "back to back".
Definition SkyWcs.cc:161
std::shared_ptr< SkyWcs > makeFlippedWcs(SkyWcs const &wcs, bool flipLR, bool flipTB, lsst::geom::Point2D const &center)
Return a copy of a FITS-WCS with pixel positions flipped around a specified center.
Definition SkyWcs.cc:489
CatalogT< BaseRecord > BaseCatalog
Definition fwd.h:72
AngleUnit constexpr degrees
constant with units of degrees
Definition Angle.h:110
Extent< double, 2 > Extent2D
Definition Extent.h:400
AngleUnit constexpr radians
constant with units of radians
Definition Angle.h:109
Point< double, 2 > Point2D
Definition Point.h:324
double constexpr PI
The ratio of a circle's circumference to diameter.
Definition Angle.h:40
T dynamic_pointer_cast(T... args)
T pow(T... args)
T sin(T... args)
T size(T... args)
T str(T... args)
T to_string(T... args)
std::shared_ptr< table::io::Persistable > read(table::io::InputArchive const &archive, table::io::CatalogVector const &catalogs) const override