48 std::shared_ptr<table::io::Persistable>
const&);
53int const SERIALIZATION_VERSION = 1;
61double const TIGHT_FITS_TOL = 0.0001;
63class SkyWcsPersistenceHelper {
68 static SkyWcsPersistenceHelper
const& get() {
69 static SkyWcsPersistenceHelper instance;
74 SkyWcsPersistenceHelper(
const SkyWcsPersistenceHelper&) =
delete;
75 SkyWcsPersistenceHelper& operator=(
const SkyWcsPersistenceHelper&) =
delete;
78 SkyWcsPersistenceHelper(SkyWcsPersistenceHelper&&) =
delete;
79 SkyWcsPersistenceHelper& operator=(SkyWcsPersistenceHelper&&) =
delete;
82 SkyWcsPersistenceHelper()
84 wcs(
schema.addField<
table::Array<std::uint8_t>>(
"wcs",
"wcs string representation",
"")) {}
89 explicit SkyWcsFactory(std::string
const& name) : table::io::PersistableFactory(
name) {}
91 std::shared_ptr<table::io::Persistable>
read(InputArchive
const& archive,
92 CatalogVector
const& catalogs)
const override {
93 SkyWcsPersistenceHelper
const&
keys = SkyWcsPersistenceHelper::get();
97 table::BaseRecord
const& record =
catalogs.front().front();
98 std::string stringRep = formatters::bytesToString(record.get(
keys.wcs));
99 return SkyWcs::readString(stringRep);
103std::string getSkyWcsPersistenceName() {
return "SkyWcs"; }
105SkyWcsFactory registration(getSkyWcsPersistenceName());
107ast::FrameDict makeSkyWcsFrameDict(TransformPoint2ToPoint2
const& pixelsToFieldAngle,
108 lsst::geom::Angle
const& orientation,
bool flipX,
109 lsst::geom::SpherePoint
const& crval,
110 std::string
const& projection =
"TAN") {
112 auto const initialWcs =
114 auto const initialFrameDict = initialWcs->getFrameDict();
115 auto const iwcToSkyMap = initialFrameDict->getMapping(
"IWC",
"SKY");
116 auto const pixelFrame = initialFrameDict->getFrame(
"PIXELS");
117 auto const iwcFrame = initialFrameDict->getFrame(
"IWC");
118 auto const skyFrame = initialFrameDict->getFrame(
"SKY");
121 ndarray::Array<double, 2, 2> fieldAngleToIwcNdArray = ndarray::allocate(2, 2);
122 asEigenMatrix(fieldAngleToIwcNdArray) = orientationAndFlipXMatrix * 180.0 /
lsst::geom::PI;
123 auto const pixelsToFieldAngleMap = pixelsToFieldAngle.getMapping();
124 auto const fieldAngleToIwcMap = ast::MatrixMap(fieldAngleToIwcNdArray);
125 auto const pixelsToIwcMap = pixelsToFieldAngleMap->then(fieldAngleToIwcMap);
126 auto finalFrameDict = ast::FrameDict(*pixelFrame, pixelsToIwcMap, *iwcFrame);
127 finalFrameDict.addFrame(
"IWC", *iwcToSkyMap, *skyFrame);
128 return finalFrameDict;
135 Eigen::Matrix2d cdMatrix;
136 double orientRad = orientation.asRadians();
137 double scaleDeg = scale.asDegrees();
138 double xmult = flipX ? 1.0 : -1.0;
139 cdMatrix(0, 0) =
std::cos(orientRad) * scaleDeg * xmult;
140 cdMatrix(0, 1) =
std::sin(orientRad) * scaleDeg;
141 cdMatrix(1, 0) = -
std::sin(orientRad) * scaleDeg * xmult;
142 cdMatrix(1, 1) =
std::cos(orientRad) * scaleDeg;
148 return src.getTransform()->then(*dstInverse);
161 double const side = 1.0;
171 auto skyLL = skyVec[0].getVector();
172 auto skyDx = skyVec[1].getVector() - skyLL;
173 auto skyDy = skyVec[2].getVector() - skyLL;
178 double skyAreaSq = skyDx.cross(skyDy).getSquaredNorm();
189 auto const crvalRad = skyFrame->getSkyRef();
196 return pixelToIwcTransform.getJacobian(pixel);
218 auto const iwcToSky = thisDict->getMapping(
"IWC",
"SKY");
219 auto const gridToSky = gridToPixel.
then(*pixelToIwc).
then(*iwcToSky);
224 os <<
"Encoding=FITS-WCS, CDMatrix=1, FitsAxisOrder=<copy>, FitsTol=" << TIGHT_FITS_TOL;
227 int const nObjectsWritten = fitsChan.
write(frameSet);
229 if (nObjectsWritten == 0) {
232 "Could not represent this SkyWcs using FITS-WCS metadata");
237 header = tanWcs->getFitsMetadata(
true);
244 header->remove(
"DATE-OBS");
245 header->remove(
"MJD-OBS");
248 bool const hasCd11 = header->exists(
"CD1_1");
249 bool const hasCd12 = header->exists(
"CD1_2");
250 bool const hasCd21 = header->exists(
"CD2_1");
251 bool const hasCd22 = header->exists(
"CD2_2");
252 if (hasCd11 || hasCd12 || hasCd21 || hasCd22) {
253 if (!hasCd11) header->set(
"CD1_1", 0.0,
"Transformation matrix element");
254 if (!hasCd12) header->set(
"CD1_2", 0.0,
"Transformation matrix element");
255 if (!hasCd21) header->set(
"CD2_1", 0.0,
"Transformation matrix element");
256 if (!hasCd22) header->set(
"CD2_2", 0.0,
"Transformation matrix element");
281 return _linearizePixelToSky(pix,
pixelToSky(pix), skyUnit);
291 return _linearizeSkyToPixel(pix,
pixelToSky(pix), skyUnit);
311 is >> shortClassName;
324 os <<
"The AST serialization was read as a " << astObjectPtr->getClassName()
325 <<
" instead of a FrameSet";
349 return std::make_unique<SkyWcs>(*
this);
355 os <<
"FITS standard SkyWcs:";
357 os <<
"Non-standard SkyWcs (Frames: ";
382 SkyWcsPersistenceHelper
const& keys = SkyWcsPersistenceHelper::get();
390 : _frameDict(frameDict), _transform(), _pixelOrigin(), _pixelScaleAtOrigin(0 *
lsst::
geom::radians) {
397 for (
auto const& domainName : domainNames) {
399 auto const frame = frameDict.
getFrame(domainName,
false);
400 if (frame->getNAxes() != 2) {
402 "Frame " + domainName +
" has " +
std::to_string(frame->getNAxes()) +
403 " axes instead of 2");
405 auto desiredClassName = domainName ==
"SKY" ?
"SkyFrame" :
"Frame";
406 if (frame->getClassName() != desiredClassName) {
408 "Frame " + domainName +
" is of type " + frame->getClassName() +
409 " instead of " + desiredClassName);
411 }
else if (domainName !=
"ACTUAL_PIXELS") {
413 throw LSST_EXCEPT(lsst::pex::exceptions::TypeError,
414 "No frame with domain " + domainName +
" found");
420 if (baseDomain !=
"ACTUAL_PIXELS" && baseDomain !=
"PIXELS") {
421 throw LSST_EXCEPT(lsst::pex::exceptions::TypeError,
422 "Base frame has domain " + baseDomain +
" instead of PIXELS or ACTUAL_PIXELS");
427 if (currentDomain !=
"SKY") {
428 throw LSST_EXCEPT(lsst::pex::exceptions::TypeError,
429 "Current frame has domain " + currentDomain +
" instead of SKY");
432 return frameDict.
copy();
436 lsst::geom::SpherePoint
const& coord,
437 lsst::geom::AngleUnit
const& skyUnit)
const {
441 const double side = 1.0;
447 m(0, 0) = dsky10.first.asAngularUnits(skyUnit) / side;
448 m(0, 1) = dsky01.first.asAngularUnits(skyUnit) / side;
449 m(1, 0) = dsky10.second.asAngularUnits(skyUnit) / side;
450 m(1, 1) = dsky01.second.asAngularUnits(skyUnit) / side;
452 Eigen::Vector2d sky00v;
453 sky00v << sky00.getX(), sky00.getY();
454 Eigen::Vector2d pix00v;
455 pix00v << pix00.getX(), pix00.getY();
457 return lsst::geom::AffineTransform(m, (sky00v - m * pix00v));
461 lsst::geom::SpherePoint
const& coord,
462 lsst::geom::AngleUnit
const& skyUnit)
const {
463 lsst::geom::AffineTransform inverse = _linearizePixelToSky(pix00, coord, skyUnit);
469 double const dx = 1000;
487 bool modifyActualPixels) {
488 auto const pixelMapping = pixelTransform.
getMapping();
489 auto oldFrameDict = wcs.getFrameDict();
490 bool const hasActualPixels = oldFrameDict->hasDomain(
"ACTUAL_PIXELS");
491 auto const pixelFrame = oldFrameDict->getFrame(
"PIXELS",
false);
492 auto const iwcFrame = oldFrameDict->getFrame(
"IWC",
false);
493 auto const skyFrame = oldFrameDict->getFrame(
"SKY",
false);
494 auto const oldPixelToIwc = oldFrameDict->getMapping(
"PIXELS",
"IWC");
495 auto const iwcToSky = oldFrameDict->getMapping(
"IWC",
"SKY");
499 if (hasActualPixels) {
500 auto const actualPixelFrame = oldFrameDict->getFrame(
"ACTUAL_PIXELS",
false);
501 auto const oldActualPixelToPixels = oldFrameDict->getMapping(
"ACTUAL_PIXELS",
"PIXELS");
503 if (modifyActualPixels) {
504 newActualPixelsToPixels = pixelMapping->then(*oldActualPixelToPixels).simplified();
505 newPixelToIwc = oldPixelToIwc;
507 newActualPixelsToPixels = oldActualPixelToPixels;
508 newPixelToIwc = pixelMapping->then(*oldPixelToIwc).simplified();
512 newFrameDict->addFrame(
"PIXELS", *newPixelToIwc, *iwcFrame);
514 newPixelToIwc = pixelMapping->then(*oldPixelToIwc).simplified();
517 newFrameDict->addFrame(
"IWC", *iwcToSky, *skyFrame);
526 Eigen::Matrix2d
const& cdMatrix,
std::string const& projection) {
534 auto frameDict = makeSkyWcsFrameDict(pixelsToFieldAngle, orientation, flipX, boresight, projection);
539 Eigen::Matrix2d
const& cdMatrix, Eigen::MatrixXd
const& sipA,
540 Eigen::MatrixXd
const& sipB) {
546 Eigen::Matrix2d
const& cdMatrix, Eigen::MatrixXd
const& sipA,
547 Eigen::MatrixXd
const& sipB, Eigen::MatrixXd
const& sipAp,
548 Eigen::MatrixXd
const& sipBp) {
555 auto iwcToSky = wcs.getFrameDict()->getMapping(
"IWC",
"SKY");
565 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.
int write(Object const &object)
Write an object to 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::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.
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.
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