48 std::shared_ptr<table::io::Persistable>
const&);
53int const SERIALIZATION_VERSION = 1;
61double const TIGHT_FITS_TOL = 0.0001;
63class SkyWcsPersistenceHelper {
70 SkyWcsPersistenceHelper(
const SkyWcsPersistenceHelper&) =
delete;
71 SkyWcsPersistenceHelper& operator=(
const SkyWcsPersistenceHelper&) =
delete;
74 SkyWcsPersistenceHelper(SkyWcsPersistenceHelper&&) =
delete;
75 SkyWcsPersistenceHelper& operator=(SkyWcsPersistenceHelper&&) =
delete;
77 explicit SkyWcsPersistenceHelper(
bool hasFitsApproximation)
79 wcs(
schema.addField<
table::Array<std::uint8_t>>(
"wcs",
"wcs string representation",
"")) {
80 if (hasFitsApproximation) {
82 "approx",
"wcs string representation of FITS approximation",
"");
90 }
catch (pex::exceptions::NotFoundError&) {
97 explicit SkyWcsFactory(std::string
const& name) : table::io::PersistableFactory(
name) {}
99 std::shared_ptr<table::io::Persistable>
read(InputArchive
const& archive,
100 CatalogVector
const& catalogs)
const override {
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));
118std::string getSkyWcsPersistenceName() {
return "SkyWcs"; }
120SkyWcsFactory registration(getSkyWcsPersistenceName());
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") {
127 auto const initialWcs =
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");
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;
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;
163 return src.getTransform()->then(*dstInverse);
176 double const side = 1.0;
186 auto skyLL = skyVec[0].getVector();
187 auto skyDx = skyVec[1].getVector() - skyLL;
188 auto skyDy = skyVec[2].getVector() - skyLL;
193 double skyAreaSq = skyDx.cross(skyDy).getSquaredNorm();
204 auto const crvalRad = skyFrame->getSkyRef();
211 return pixelToIwcTransform.getJacobian(pixel);
233 auto const iwcToSky = thisDict->getMapping(
"IWC",
"SKY");
234 auto const gridToSky = gridToPixel.
then(*pixelToIwc).
then(*iwcToSky);
239 os <<
"Encoding=FITS-WCS, CDMatrix=1, FitsAxisOrder=<copy>, FitsTol=" << TIGHT_FITS_TOL;
242 int const nObjectsWritten = fitsChan.write(frameSet);
243 if (nObjectsWritten == 0) {
249 header->remove(
"DATE-OBS");
250 header->remove(
"MJD-OBS");
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");
271 auto result = _getDirectFitsMetadata();
274 precise ?
"WCS is not directly FITS-compatible."
275 :
"WCS does not have an attached FITS approximation.");
283 if (fitsApproximation && fitsApproximation->hasFitsApproximation()) {
285 "Cannot add a FITS approximation that itself already has a FITS approximation.");
288 result->_fitsApproximation = fitsApproximation;
300 return _linearizePixelToSky(pix,
pixelToSky(pix), skyUnit);
310 return _linearizeSkyToPixel(pix,
pixelToSky(pix), skyUnit);
330 is >> shortClassName;
343 os <<
"The AST serialization was read as a " << astObjectPtr->getClassName()
344 <<
" instead of a FrameSet";
368 return std::make_unique<SkyWcs>(*
this);
374 os <<
"FITS standard SkyWcs:";
376 os <<
"Non-standard SkyWcs (Frames: ";
412 : _frameDict(frameDict), _transform(), _pixelOrigin(), _pixelScaleAtOrigin(0 *
lsst::
geom::radians) {
419 for (
auto const& domainName : domainNames) {
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");
427 auto desiredClassName = domainName ==
"SKY" ?
"SkyFrame" :
"Frame";
428 if (frame->getClassName() != desiredClassName) {
430 "Frame " + domainName +
" is of type " + frame->getClassName() +
431 " instead of " + desiredClassName);
433 }
else if (domainName !=
"ACTUAL_PIXELS") {
435 throw LSST_EXCEPT(lsst::pex::exceptions::TypeError,
436 "No frame with domain " + domainName +
" found");
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");
449 if (currentDomain !=
"SKY") {
450 throw LSST_EXCEPT(lsst::pex::exceptions::TypeError,
451 "Current frame has domain " + currentDomain +
" instead of SKY");
454 return frameDict.
copy();
458 lsst::geom::SpherePoint
const& coord,
459 lsst::geom::AngleUnit
const& skyUnit)
const {
463 const double side = 1.0;
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;
474 Eigen::Vector2d sky00v;
475 sky00v << sky00.getX(), sky00.getY();
476 Eigen::Vector2d pix00v;
477 pix00v << pix00.getX(), pix00.getY();
479 return lsst::geom::AffineTransform(m, (sky00v - m * pix00v));
483 lsst::geom::SpherePoint
const& coord,
484 lsst::geom::AngleUnit
const& skyUnit)
const {
485 lsst::geom::AffineTransform inverse = _linearizePixelToSky(pix00, coord, skyUnit);
491 double const dx = 1000;
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");
521 if (hasActualPixels) {
522 auto const actualPixelFrame = oldFrameDict->getFrame(
"ACTUAL_PIXELS",
false);
523 auto const oldActualPixelToPixels = oldFrameDict->getMapping(
"ACTUAL_PIXELS",
"PIXELS");
525 if (modifyActualPixels) {
526 newActualPixelsToPixels = pixelMapping->then(*oldActualPixelToPixels).simplified();
527 newPixelToIwc = oldPixelToIwc;
529 newActualPixelsToPixels = oldActualPixelToPixels;
530 newPixelToIwc = pixelMapping->then(*oldPixelToIwc).simplified();
534 newFrameDict->addFrame(
"PIXELS", *newPixelToIwc, *iwcFrame);
536 newPixelToIwc = pixelMapping->then(*oldPixelToIwc).simplified();
539 newFrameDict->addFrame(
"IWC", *iwcToSky, *skyFrame);
548 Eigen::Matrix2d
const& cdMatrix,
std::string const& projection) {
556 auto frameDict = makeSkyWcsFrameDict(pixelsToFieldAngle, orientation, flipX, boresight, projection);
561 Eigen::Matrix2d
const& cdMatrix, Eigen::MatrixXd
const& sipA,
562 Eigen::MatrixXd
const& sipB) {
568 Eigen::Matrix2d
const& cdMatrix, Eigen::MatrixXd
const& sipA,
569 Eigen::MatrixXd
const& sipB, Eigen::MatrixXd
const& sipAp,
570 Eigen::MatrixXd
const& sipBp) {
577 auto iwcToSky = wcs.getFrameDict()->getMapping(
"IWC",
"SKY");
587 os << wcs.toString();
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
#define LSST_ARCHIVE_ASSERT(EXPR)
An assertion macro used to validate the structure of an InputArchive.
Channel provides input/output of AST objects.
std::shared_ptr< Object > read()
Read an object from a channel.
A specialized form of Channel which reads and writes FITS header cards.
A FrameSet whose frames can be referenced by domain name.
bool hasDomain(std::string const &domain) const
Return True if a frame in this FrameDict has the specified domain.
std::shared_ptr< Mapping > getMapping(int from, std::string const &to) const
Variant of FrameSet::getMapping with the second frame specified by domain.
std::shared_ptr< FrameDict > copy() const
Return a deep copy of this object.
std::set< std::string > getAllDomains() const
Get the domain names for all contained Frames (excluding frames with empty or defaulted domain names)...
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.
Frame is used to represent a coordinate system.
A FrameSet consists of a set of one or more Frames (which describe coordinate systems),...
static int constexpr CURRENT
index of current frame
static int constexpr BASE
index of base frame
SeriesMap then(Mapping const &next) const
Return a series compound mapping this(first(input)) containing shallow copies of the original.
void show(std::ostream &os, bool showComments=true) const
Print a textual description the object to an ostream.
ShiftMap is a linear Mapping which shifts each axis by a specified constant value.
A stream for ast::Channel.
String-based source and sink for channels.
A WinMap is a linear Mapping which transforms a rectangular window in one coordinate system into a si...
A 2-dimensional celestial WCS that transform pixels to ICRS RA/Dec, using the LSST standard for pixel...
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.
static std::shared_ptr< SkyWcs > readStream(std::istream &is)
Deserialize a SkyWcs from an input stream.
std::shared_ptr< const ast::FrameDict > getFrameDict() const
Get the contained FrameDict.
lsst::geom::SpherePoint pixelToSky(lsst::geom::Point2D const &pixel) const
Compute sky position(s) from pixel position(s)
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.
std::string toString() const override
Create a string representation of this object.
bool isFlipped() const
Does the WCS follow the convention of North=Up, East=Left?
std::shared_ptr< typehandling::Storable > cloneStorable() const override
Create a new SkyWcs that is a copy of this one.
lsst::geom::Angle getPixelScale() const
Get the pixel scale at the pixel origin.
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.
lsst::geom::Point2D getPixelOrigin() const
Get the pixel origin, in pixels, using the LSST convention.
void writeStream(std::ostream &os) const
Serialize this SkyWcs to an output stream.
std::string writeString() const
Serialize this SkyWcs to a string, using the same format as writeStream.
Eigen::Matrix2d getCdMatrix() const
Get the 2x2 CD matrix at the pixel origin.
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.
std::shared_ptr< const TransformPoint2ToSpherePoint > getTransform() const
Get a TransformPoint2ToSpherePoint that transforms pixels to sky in the forward direction and sky to ...
lsst::geom::Point2D skyToPixel(lsst::geom::SpherePoint const &sky) const
Compute pixel position(s) from sky position(s)
static std::string getShortClassName()
bool isFits() const
Return true getFitsMetadata(true) will succeed, false if not.
bool hasFitsApproximation() const
Does this SkyWcs have an approximate SkyWcs that can be represented as standard FITS WCS?
std::shared_ptr< SkyWcs > getFitsApproximation() const
Return FITS SkyWcs that approximates this one.
std::string getPythonModule() const override
Return the fully-qualified Python module that should be imported to guarantee that its factory is reg...
std::string getPersistenceName() const override
Return the unique name used to persist this object and look up its factory.
lsst::geom::SpherePoint getSkyOrigin() const
Get the sky origin, the celestial fiducial point.
std::shared_ptr< SkyWcs > copyWithFitsApproximation(std::shared_ptr< SkyWcs > fitsApproximation) const
Return a copy of this SkyWcs with the given FITS approximation.
static std::shared_ptr< SkyWcs > readString(std::string &str)
Deserialize a SkyWcs from a string, using the same format as readStream.
std::shared_ptr< daf::base::PropertyList > getFitsMetadata(bool precise=false) const
Return the WCS as FITS WCS metadata.
void write(OutputArchiveHandle &handle) const override
Write the object to one or more catalogs.
bool equals(typehandling::Storable const &other) const noexcept override
Compare this object to another Storable.
bool operator==(SkyWcs const &other) const
Equality is based on the string representations being equal.
Tag types used to declare specialized field types.
std::shared_ptr< RecordT > addNew()
Create a new record, add it to the end of the catalog, and return a pointer to it.
A class used as a handle to a particular field in a table.
Defines the fields and offsets for a table.
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.
static bool singleClassEquals(T const &lhs, Storable const &rhs)
Test if a Storable is of a particular class and equal to another object.
Class for storing generic metadata.
A class representing an angle.
constexpr double asArcseconds() const noexcept
Return an Angle's value in arcseconds.
A class used to convert scalar POD types such as double to Angle.
Point in an unspecified spherical coordinate system.
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.
Reports errors that are due to events beyond the control of the program.
Reports errors from accepting an object of an unexpected or inappropriate type.
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.
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).
std::shared_ptr< TransformPoint2ToPoint2 > getPixelToIntermediateWorldCoords(SkyWcs const &wcs, bool simplify=true)
Return a transform from pixel coordinates to intermediate world coordinates.
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.
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)".
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.
std::shared_ptr< TransformPoint2ToSpherePoint > getIntermediateWorldCoordsToSky(SkyWcs const &wcs, bool simplify=true)
Return a transform from intermediate world coordinates to sky.
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.
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.
Transform< Point2Endpoint, Point2Endpoint > TransformPoint2ToPoint2
std::shared_ptr< TransformPoint2ToPoint2 > makeWcsPairTransform(SkyWcs const &src, SkyWcs const &dst)
A Transform obtained by putting two SkyWcs objects "back to back".
std::shared_ptr< SkyWcs > makeFlippedWcs(SkyWcs const &wcs, bool flipLR, bool flipTB, lsst::geom::Point2D const ¢er)
Return a copy of a FITS-WCS with pixel positions flipped around a specified center.
CatalogT< BaseRecord > BaseCatalog
AngleUnit constexpr degrees
constant with units of degrees
Extent< double, 2 > Extent2D
AngleUnit constexpr radians
constant with units of radians
Point< double, 2 > Point2D
double constexpr PI
The ratio of a circle's circumference to diameter.
T dynamic_pointer_cast(T... args)
std::shared_ptr< table::io::Persistable > read(table::io::InputArchive const &archive, table::io::CatalogVector const &catalogs) const override