SkyWcs.h

Go to the documentation of this file.

// -*- lsst-c++ -*-
/*
 * LSST Data Management System
 * Copyright 2017 AURA/LSST.
 *
 * This product includes software developed by the
 * LSST Project (http://www.lsst.org/).
 *
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the LSST License Statement and
 * the GNU General Public License along with this program.  If not,
 * see <http://www.lsstcorp.org/LegalNotices/>.
 */
 
#ifndef LSST_AFW_GEOM_SKYWCS_H
#define LSST_AFW_GEOM_SKYWCS_H
 
#include <memory>
#include <optional>
#include <utility>
#include <vector>
 
#include "Eigen/Core"
#include "astshim.h"
#include "ndarray.h"
 
#include "lsst/daf/base/PropertyList.h"
#include "lsst/geom/AffineTransform.h"
#include "lsst/geom/Angle.h"
#include "lsst/geom/Box.h"
#include "lsst/geom/Point.h"
#include "lsst/afw/geom/Endpoint.h"
#include "lsst/afw/geom/Transform.h"
#include "lsst/afw/table/io/Persistable.h"
#include "lsst/afw/typehandling/Storable.h"
 
namespace lsst {
namespace afw {
namespace geom {

Eigen::Matrix2d makeCdMatrix(lsst::geom::Angle const &scale,
                             lsst::geom::Angle const &orientation = 0 * lsst::geom::degrees,
                             bool flipX = false);

class SkyWcs final : public table::io::PersistableFacade<SkyWcs>, public typehandling::Storable {
public:
    SkyWcs(SkyWcs const &) = default;
    SkyWcs(SkyWcs &&) = default;
    SkyWcs &operator=(SkyWcs const &) = delete;
    SkyWcs &operator=(SkyWcs &&) = delete;
    ~SkyWcs() override = default;

    bool operator==(SkyWcs const &other) const;
    bool operator!=(SkyWcs const &other) const { return !(*this == other); }

    explicit SkyWcs(daf::base::PropertySet &metadata, bool strip = false);

    explicit SkyWcs(ast::FrameDict const &frameDict);

    std::shared_ptr<SkyWcs> copyAtShiftedPixelOrigin(lsst::geom::Extent2D const &shift) const;

    std::shared_ptr<daf::base::PropertyList> getFitsMetadata(
        bool precise = false,
        lsst::geom::Box2I const & bbox = lsst::geom::Box2I(
            lsst::geom::Point2I(0, 0), lsst::geom::Extent2I(100, 100)
        )
    ) const;

    lsst::geom::Angle getPixelScale(lsst::geom::Point2D const &pixel) const;

    lsst::geom::Angle getPixelScale() const { return _pixelScaleAtOrigin; };

    lsst::geom::Point2D getPixelOrigin() const { return _pixelOrigin; };

    lsst::geom::SpherePoint getSkyOrigin() const;

    Eigen::Matrix2d getCdMatrix(lsst::geom::Point2D const &pixel) const;

    Eigen::Matrix2d getCdMatrix() const;

    std::shared_ptr<SkyWcs> getTanWcs(lsst::geom::Point2D const &pixel) const;

    std::shared_ptr<const ast::FrameDict> getFrameDict() const;

    std::shared_ptr<const TransformPoint2ToSpherePoint> getTransform() const { return _transform; };

    bool isFlipped() const;
 
    bool isPersistable() const noexcept override { return true; }

    lsst::geom::AffineTransform linearizePixelToSky(lsst::geom::SpherePoint const &coord,
                                                    lsst::geom::AngleUnit const &skyUnit) const;

    lsst::geom::AffineTransform linearizePixelToSky(lsst::geom::Point2D const &pixel,
                                                    lsst::geom::AngleUnit const &skyUnit) const;

    lsst::geom::AffineTransform linearizeSkyToPixel(lsst::geom::SpherePoint const &coord,
                                                    lsst::geom::AngleUnit const &skyUnit) const;

    lsst::geom::AffineTransform linearizeSkyToPixel(lsst::geom::Point2D const &pixel,
                                                    lsst::geom::AngleUnit const &skyUnit) const;

    lsst::geom::SpherePoint pixelToSky(lsst::geom::Point2D const &pixel) const {
        return _transform->applyForward(pixel);
    }
    lsst::geom::SpherePoint pixelToSky(double x, double y) const {
        return pixelToSky(lsst::geom::Point2D(x, y));
    }
    std::vector<lsst::geom::SpherePoint> pixelToSky(std::vector<lsst::geom::Point2D> const &pixels) const {
        return _transform->applyForward(pixels);
    }


    lsst::geom::Point2D skyToPixel(lsst::geom::SpherePoint const &sky) const {
        return _transform->applyInverse(sky);
    }
    std::vector<lsst::geom::Point2D> skyToPixel(std::vector<lsst::geom::SpherePoint> const &sky) const {
        return _transform->applyInverse(sky);
    }

 
    static std::string getShortClassName();

    bool hasFitsApproximation() const { return bool(_fitsApproximation); }

    std::shared_ptr<SkyWcs> getFitsApproximation() const { return _fitsApproximation; }

    std::shared_ptr<SkyWcs> copyWithFitsApproximation(std::shared_ptr<SkyWcs> fitsApproximation) const;

    bool isFits() const;

    static std::shared_ptr<SkyWcs> readStream(std::istream &is);

    static std::shared_ptr<SkyWcs> readString(std::string &str);

    void writeStream(std::ostream &os) const;

    std::string writeString() const;
 
    // Override methods required by afw::typehandling::Storable

    std::shared_ptr<typehandling::Storable> cloneStorable() const override;

    std::string toString() const override;

    bool equals(typehandling::Storable const &other) const noexcept override;
 
protected:
    // Override methods required by afw::table::io::Persistable
    std::string getPersistenceName() const override;
    std::string getPythonModule() const override;
    void write(OutputArchiveHandle &handle) const override;
 
private:
    /*
     * Construct a SkyWcs from a shared pointer to an ast::FrameDict
     *
     * The frameDict may be modified.
     */
    explicit SkyWcs(std::shared_ptr<ast::FrameDict> frameDict);
 
    /*
     * Check a FrameDict to see if it can safely be used for a SkyWcs
     * Return a copy so that it can be used as an argument to the SkyWcs(shared_ptr<FrameDict>) constructor
     */
    std::shared_ptr<ast::FrameDict> _checkFrameDict(ast::FrameDict const &frameDict) const;
 
    // the full FrameDict, for operations that need intermediate frames
    std::shared_ptr<const ast::FrameDict> _frameDict;
    // cached transform from _frameDict, for fast computation of pixels<->sky
    std::shared_ptr<const TransformPoint2ToSpherePoint> _transform;
    lsst::geom::Point2D _pixelOrigin;       // cached pixel origin
    lsst::geom::Angle _pixelScaleAtOrigin;  // cached pixel scale at pixel origin
 
    std::shared_ptr<SkyWcs> _fitsApproximation;
 
    /*
     * Implementation for the overloaded public linearizePixelToSky methods, requiring both a pixel coordinate
     * and the corresponding sky coordinate.
     */
    lsst::geom::AffineTransform _linearizePixelToSky(lsst::geom::Point2D const &pix,
                                                     lsst::geom::SpherePoint const &coord,
                                                     lsst::geom::AngleUnit const &skyUnit) const;
 
    /*
     * Implementation for the overloaded public linearizeSkyToPixel methods, requiring both a pixel coordinate
     * and the corresponding sky coordinate.
     */
    lsst::geom::AffineTransform _linearizeSkyToPixel(lsst::geom::Point2D const &pix,
                                                     lsst::geom::SpherePoint const &coord,
                                                     lsst::geom::AngleUnit const &skyUnit) const;
 
    /*
     * Compute the FITS representation of this WCS (without falling back to a FITS approximation) or return
     * a null pointer if it is not FITS-compatible.
    */
    std::shared_ptr<daf::base::PropertyList> _getDirectFitsMetadata(lsst::geom::Box2I const & bbox) const;

    void _computeCache() {
        _transform = std::make_shared<TransformPoint2ToSpherePoint>(*_frameDict->getMapping(), true);
        _pixelOrigin = skyToPixel(getSkyOrigin());
        _pixelScaleAtOrigin = getPixelScale(_pixelOrigin);
    }
};

std::shared_ptr<SkyWcs> makeFlippedWcs(SkyWcs const &wcs, bool flipLR, bool flipTB,
                                       lsst::geom::Point2D const &center);

std::shared_ptr<SkyWcs> makeModifiedWcs(TransformPoint2ToPoint2 const &pixelTransform, SkyWcs const &wcs,
                                        bool modifyActualPixels);

std::shared_ptr<SkyWcs> makeSkyWcs(daf::base::PropertySet &metadata, bool strip = false);

std::shared_ptr<SkyWcs> makeSkyWcs(lsst::geom::Point2D const &crpix, lsst::geom::SpherePoint const &crval,
                                   Eigen::Matrix2d const &cdMatrix, std::string const &projection = "TAN");

std::shared_ptr<SkyWcs> makeSkyWcs(TransformPoint2ToPoint2 const &pixelsToFieldAngle,
                                   lsst::geom::Angle const &orientation, bool flipX,
                                   lsst::geom::SpherePoint const &boresight,
                                   std::string const &projection = "TAN");

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);

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, Eigen::MatrixXd const &sipAp,
                                      Eigen::MatrixXd const &sipBp);
std::shared_ptr<TransformPoint2ToPoint2> makeWcsPairTransform(SkyWcs const &src, SkyWcs const &dst);

std::shared_ptr<TransformPoint2ToSpherePoint> getIntermediateWorldCoordsToSky(SkyWcs const &wcs,
                                                                              bool simplify = true);

std::shared_ptr<TransformPoint2ToPoint2> getPixelToIntermediateWorldCoords(SkyWcs const &wcs,
                                                                           bool simplify = true);

std::ostream &operator<<(std::ostream &os, SkyWcs const &wcs);
 
}  // namespace geom
}  // namespace afw
}  // namespace lsst
 
#endif