lsst::afw::geom::detail Namespace Reference

Functions
std::shared_ptr< ast::FrameSet >	readFitsWcs (daf::base::PropertySet &metadata, bool strip=true)
	Read a FITS convention WCS FrameSet from FITS metadata.
std::shared_ptr< ast::FrameDict >	readLsstSkyWcs (daf::base::PropertySet &metadata, bool strip=true)
	Read an LSST celestial WCS FrameDict from a FITS header.
void	stripWcsMetadata (daf::base::PropertySet &metadata)
	Strip all WCS metadata from a header.
std::shared_ptr< daf::base::PropertyList >	getPropertyListFromFitsChan (ast::FitsChan &fitsChan)
	Copy values from an AST FitsChan into a PropertyList.
ast::FitsChan	getFitsChanFromPropertyList (daf::base::PropertySet &metadata, std::set< std::string > const &excludeNames={}, std::string options="")
	Construct AST FitsChan from PropertyList.
template<class Transform >
std::shared_ptr< Transform >	readStream (std::istream &is)
	Deserialize a Transform from an input stream.
template<class Transform >
void	writeStream (Transform const &transform, std::ostream &os)
	Serialize a Transform to an output stream.
void	wrapFrameSetUtils (lsst::cpputils::python::WrapperCollection &)

Variables
constexpr int	serializationVersion = 1
	version of serialization used when writing (older versions may also be supported when reading)

Function Documentation

getFitsChanFromPropertyList()

ast::FitsChan lsst::afw::geom::detail::getFitsChanFromPropertyList	(	daf::base::PropertySet &	metadata,
		std::set< std::string > const &	excludeNames = {},
		std::string	options = "" )

Construct AST FitsChan from PropertyList.

Parameters

[in]	metadata	Metadata to transfer; if this is a PropertyList then the order of items is preserved.
[in]	excludeNames	Names of entries to exclude from the returned header string.
[in]	options	Options string to pass to ast::FitsChan constructor.

Returns

an ast::FitsChan representing this PropertyList

Definition at line 361 of file frameSetUtils.cc.

                                                             {
    // Create FitsChan to receive each of the parameters
    auto stream = ast::StringStream("");
    auto fc = ast::FitsChan(stream, options);
 
    // Get the parameter names (in order if necessary)
    daf::base::PropertyList const *pl = dynamic_cast<daf::base::PropertyList const *>(&metadata);
    std::vector<std::string> allParamNames;
    if (pl) {
        allParamNames = pl->getOrderedNames();
    } else {
        allParamNames = metadata.paramNames(false);
    }
 
    // Loop over the names and add them to the FitsChan if not excluded
    for (auto const &fullName : allParamNames) {
        if (excludeNames.count(fullName) == 0) {
            std::size_t lastPeriod = fullName.rfind(char('.'));
            auto name = (lastPeriod == std::string::npos) ? fullName : fullName.substr(lastPeriod + 1);
            std::type_info const &type = metadata.typeOf(name);
 
            if (name.size() > 8) {
                continue;  // The name is too long for a FITS keyword; skip this item
            }
 
            if (type == typeid(bool)) {
                fc.setFitsL(name, metadata.get<bool>(name));
            } else if (type == typeid(std::uint8_t)) {
                fc.setFitsI(name, static_cast<int>(metadata.get<std::uint8_t>(name)));
            } else if (type == typeid(int)) {
                fc.setFitsI(name, metadata.get<int>(name));
            } else if (type == typeid(double) || type == typeid(float)) {
                double value;
                if (type == typeid(double)) {
                    value = metadata.get<double>(name);
                } else {
                    value = static_cast<double>(metadata.get<float>(name));
                }
                // NaN is not allowed in a FitsChan (or in FITS)
                if (!std::isnan(value)) {
                    fc.setFitsF(name, value);
                } else {
                    // Treat it like an undefined value but warn about it
                    LOGLS_WARN("lsst.afw.geom.frameSetUtils",
                               boost::format("Found NaN in metadata item '%s'") % name);
                }
            } else if (type == typeid(std::string)) {
                std::string str = metadata.get<std::string>(name);
                // No support for long strings yet so skip those
                if (str.size() <= 68) {
                    fc.setFitsS(name, metadata.get<std::string>(name));
                }
            }
        }
    }
    // Rewind the channel
    fc.setCard(0);
    return fc;
}

getPropertyListFromFitsChan()

std::shared_ptr< daf::base::PropertyList > lsst::afw::geom::detail::getPropertyListFromFitsChan

(

ast::FitsChan &

fitsChan

)

Copy values from an AST FitsChan into a PropertyList.

Warning

COMMENT and HISTORY cards are treated as string values

Exceptions

lsst::pex::exceptions::TypeError

if fitsChan contains cards whose type is not supported by PropertyList: complex numbers, or cards with no value

Definition at line 296 of file frameSetUtils.cc.

                                                                                      {
    int const numCards = fitsChan.getNCard();
    auto metadata = std::make_shared<daf::base::PropertyList>();
    for (int cardNum = 1; cardNum <= numCards; ++cardNum) {
        fitsChan.setCard(cardNum);
        auto const cardType = fitsChan.getCardType();
        auto const cardName = fitsChan.getCardName();
        auto const cardComment = fitsChan.getCardComm();
        switch (cardType) {
            case ast::CardType::FLOAT: {
                auto foundValue = fitsChan.getFitsF();
                setMetadataFromFoundValue(*metadata, cardName, foundValue, cardComment);
                break;
            }
            case ast::CardType::INT: {
                auto foundValue = fitsChan.getFitsI();
                setMetadataFromFoundValue(*metadata, cardName, foundValue, cardComment);
                break;
            }
            case ast::CardType::STRING: {
                auto foundValue = fitsChan.getFitsS();
                setMetadataFromFoundValue(*metadata, cardName, foundValue, cardComment);
                break;
            }
            case ast::CardType::LOGICAL: {
                auto foundValue = fitsChan.getFitsL();
                setMetadataFromFoundValue(*metadata, cardName, foundValue, cardComment);
                break;
            }
            case ast::CardType::UNDEF: {
                if (cardComment.empty()) {
                    metadata->set(cardName, nullptr);
                } else {
                    metadata->set(cardName, nullptr, cardComment);
                }
                break;
            }
            case ast::CardType::CONTINUE: {
                auto foundValue = fitsChan.getFitsCN();
                setMetadataFromFoundValue(*metadata, cardName, foundValue, cardComment);
                break;
            }
            case ast::CardType::COMMENT:
                // Drop HISTORY and COMMENT cards
                break;
            case ast::CardType::COMPLEXF:
            case ast::CardType::COMPLEXI: {
                // PropertyList supports neither complex numbers nor cards with no value
                std::ostringstream os;
                os << "Card " << cardNum << " with name \"" << cardName << "\" has type "
                   << static_cast<int>(cardType) << ", which is not supported by PropertyList";
                throw LSST_EXCEPT(lsst::pex::exceptions::TypeError, os.str());
            }
            case ast::CardType::NOTYPE: {
                // This should only occur if cardNum is invalid, and that should be impossible
                std::ostringstream os;
                os << "Bug! Card " << cardNum << " with name \"" << cardName
                   << "\" has type NOTYPE, which should not be possible";
                throw LSST_EXCEPT(lsst::pex::exceptions::TypeError, os.str());
            }
        }
    }
    return metadata;
}

readFitsWcs()

std::shared_ptr< ast::FrameSet > lsst::afw::geom::detail::readFitsWcs	(	daf::base::PropertySet &	metadata,
		bool	strip = true )

Read a FITS convention WCS FrameSet from FITS metadata.

The resulting FrameSet may be any kind of WCS supported by FITS; if it is a celestial WCS then 1,1 will be the lower left corner of the image (the FITS convention, not the LSST convention).

This routine replaces RADECSYS with RADESYS if the former is present and the latter is not, since that is a common misspelling in FITS headers.

The returned FrameSet will have an IWC (intermediate world coordinate system) frame.

Parameters

[in,out]	metadata	FITS header cards
[in]	strip	If true: strip items from metadata used to create the WCS, such as RADESYS, EQUINOX, CTYPE12, CRPIX12, CRVAL12, etc. Always keep keywords that might be wanted for other purpposes, including NAXIS12 and date-related keywords such as "DATE-OBS" and "TIMESYS" (but not "EQUINOX").

Exceptions

pex::exceptions::TypeError

if the metadata does not contain a FITS-WCS

Definition at line 82 of file frameSetUtils.cc.

                                                                                 {
    // Exclude WCS A keywords because LSST uses them to store XY0
    auto wcsANames = createTrivialWcsMetadata("A", lsst::geom::Point2I(0, 0))->names();
    std::set<std::string> excludeNames(wcsANames.begin(), wcsANames.end());
    // Ignore NAXIS1, NAXIS2 because if they are zero then AST will fail to read a WCS
    // Ignore LTV1/2 because LSST adds it and this code should ignore it and not strip it
    // Exclude comments and history to reduce clutter
    std::set<std::string> moreNames{"NAXIS1", "NAXIS2", "LTV1", "LTV2", "COMMENT", "HISTORY"};
    excludeNames.insert(moreNames.begin(), moreNames.end());
 
    // Replace RADECSYS with RADESYS if only the former is present
    if (metadata.exists("RADECSYS") && !metadata.exists("RADESYS")) {
        metadata.set("RADESYS", metadata.getAsString("RADECSYS"));
        metadata.remove("RADECSYS");
    }
 
    auto channel = getFitsChanFromPropertyList(metadata, excludeNames,
                                               "Encoding=FITS-WCS, IWC=1, SipReplace=0, ReportLevel=3");
    auto const initialNames = strip ? setFromVector(channel.getAllCardNames()) : std::set<std::string>();
    std::shared_ptr<ast::Object> obj;
    try {
        obj = channel.read();
    } catch (std::runtime_error const&) {
        throw LSST_EXCEPT(pex::exceptions::TypeError, "The metadata does not describe an AST object");
    }
 
    auto frameSet = std::dynamic_pointer_cast<ast::FrameSet>(obj);
    if (!frameSet) {
        throw LSST_EXCEPT(pex::exceptions::TypeError,
                          "metadata describes a " + obj->getClassName() + ", not a FrameSet");
    }
    if (strip) {
        auto const finalNames = setFromVector(channel.getAllCardNames());
 
        // FITS keywords that FitsChan stripped
        std::set<std::string> namesChannelStripped;
        std::set_difference(initialNames.begin(), initialNames.end(), finalNames.begin(), finalNames.end(),
                            std::inserter(namesChannelStripped, namesChannelStripped.begin()));
 
        // FITS keywords that FitsChan may strip that we want to keep in `metadata`
        std::set<std::string> const namesToKeep = {"DATE-OBS", "MJD-OBS"};
 
        std::set<std::string> namesToStrip;  // names to strip from metadata
        std::set_difference(namesChannelStripped.begin(), namesChannelStripped.end(), namesToKeep.begin(),
                            namesToKeep.end(), std::inserter(namesToStrip, namesToStrip.begin()));
        for (auto const& name : namesToStrip) {
            metadata.remove(name);
        }
    }
    return frameSet;
}

readLsstSkyWcs()

std::shared_ptr< ast::FrameDict > lsst::afw::geom::detail::readLsstSkyWcs	(	daf::base::PropertySet &	metadata,
		bool	strip = true )

Read an LSST celestial WCS FrameDict from a FITS header.

Calls getImageXY0FromMetadata to determine image XY0.

Saves CRVAL by setting the output SkyFrame's SkyRef to CRVAL and SkyRefIs="Ignore" (so SkyRef is not treated as an offset).

Warning

Does not compensate for the offset between a subimage and its parent image; callers must do that manually, e.g. by calling SkyWcs::copyAtShiftedPosition.

the general SkyWcs generated by LSST software cannot be exactly represented using standard WCS FITS cards, and so this function cannot read those. This function is intended for two purposes:

• Read the standard FITS WCS found in raw data and other images from non-LSST observatories and convert it to the LSST pixel convention.
• Read the approximate FITS WCS that LSST writes to FITS images (for use by non-LSST code).

The frames of the returned WCS will be as follows:

• "PIXELS" (base frame): pixel frame with 0,0 the lower left corner of the image (LSST convention)
• "IWC": FITS WCS intermediate world coordinate system
• "SKY" (current frame): an ast::SkyFrame with domain "SKY": ICRS RA, Dec

All frames are instances of ast::Frame except the SKY frame. All have 2 axes.

Parameters

[in,out]	metadata	FITS header cards
[in]	strip	If true: strip items from metadata used to create the WCS, such as RADESYS, EQUINOX, CTYPE12, CRPIX12, CRVAL12, etc. Always keep keywords that might be wanted for other purpposes, including NAXIS12 and date-related keywords such as "DATE-OBS" and "TIMESYS" (but not "EQUINOX").

Exceptions

lsst::pex::exceptions::TypeError

if the metadata does not describe a celestial WCS.

Definition at line 179 of file frameSetUtils.cc.

                                                                                     {
    // Record CRPIX in GRID coordinates
    // so we can compute CRVAL after standardizing the SkyFrame to ICRS
    // (that standardization is why we don't simply save CRVAL now)
    std::vector<double> crpixGrid(2);
    try {
        crpixGrid[0] = metadata.getAsDouble("CRPIX1");
        crpixGrid[1] = metadata.getAsDouble("CRPIX2");
    } catch (lsst::pex::exceptions::NotFoundError& e) {
        // std::string used because e.what() returns a C string and two C strings cannot be added
        throw LSST_EXCEPT(lsst::pex::exceptions::TypeError,
                          e.what() + std::string("; cannot read metadata as a SkyWcs"));
    }
 
    auto rawFrameSet = readFitsWcs(metadata, strip);
    auto const initialBaseIndex = rawFrameSet->getBase();
 
    // Find the GRID frame
    auto gridIndex = ast::FrameSet::NOFRAME;
    if (rawFrameSet->findFrame(ast::Frame(2, "Domain=GRID"))) {
        gridIndex = rawFrameSet->getCurrent();
    } else {
        // No appropriate GRID frame found; if the original base frame is of type Frame
        // with 2 axes and a blank domain then use that, else give up
        auto const baseFrame = rawFrameSet->getFrame(initialBaseIndex, false);
        auto const baseClassName = rawFrameSet->getClassName();
        if (baseFrame->getClassName() != "Frame") {
            throw LSST_EXCEPT(pex::exceptions::TypeError,
                              "The base frame is of type " + baseFrame->getClassName() +
                                      "instead of Frame; cannot read metadata as a SkyWcs");
        }
        if (baseFrame->getNAxes() != 2) {
            throw LSST_EXCEPT(pex::exceptions::TypeError,
                              "The base frame has " + std::to_string(baseFrame->getNAxes()) +
                                      " axes instead of 2; cannot read metadata as a SkyWcs");
        }
        if (baseFrame->getDomain() != "") {
            throw LSST_EXCEPT(pex::exceptions::TypeError,
                              "The base frame has domain \"" + baseFrame->getDomain() +
                                      "\" instead of blank or GRID; cannot read metadata as a SkyWcs");
        }
        // Original base frame has a blank Domain, is of type Frame, and has 2 axes, so
        // Set its domain to GRID, and set some other potentially useful attributes.
        gridIndex = initialBaseIndex;
    }
 
    // Find the IWC frame
    if (!rawFrameSet->findFrame(ast::Frame(2, "Domain=IWC"))) {
        throw LSST_EXCEPT(pex::exceptions::TypeError, "No IWC frame found; cannot read metadata as a SkyWcs");
    }
    auto const iwcIndex = rawFrameSet->getCurrent();
    auto const iwcFrame = rawFrameSet->getFrame(iwcIndex);
 
    // Create a standard sky frame: ICRS with axis order RA, Dec
 
    // Create the a template for the standard sky frame
    auto const stdSkyFrameTemplate = ast::SkyFrame("System=ICRS");
 
    // Locate a Frame in the target FrameSet that looks like the template
    // and hence can be used as the original sky frame.
    // We ignore the frame set returned by findFrame because that goes from pixels to sky,
    // and using it would add an unwanted extra branch to our WCS; instead, later on,
    // we compute a mapping from the old sky frame to the new sky frame and add that.
    if (!rawFrameSet->findFrame(stdSkyFrameTemplate)) {
        throw LSST_EXCEPT(pex::exceptions::TypeError,
                          "Could not find a SkyFrame; cannot read metadata as a SkyWcs");
    }
    auto initialSkyIndex = rawFrameSet->getCurrent();
 
    // Compute a frame set that maps from the original sky frame to our desired sky frame template;
    // this contains the mapping and sky frame we will insert into the frame set.
    // (Temporarily set the base frame to the sky frame, because findFrame
    // produces a mapping from base to the found frame).
    rawFrameSet->setBase(initialSkyIndex);
    auto stdSkyFrameSet = rawFrameSet->findFrame(stdSkyFrameTemplate);
    if (!stdSkyFrameSet) {
        throw LSST_EXCEPT(pex::exceptions::LogicError,
                          "Bug: found a SkyFrame the first time, but not the second time");
    }
 
    // Add the new mapping into rawFrameSet, connecting it to the original SkyFrame.
    // Note: we use stdSkyFrameSet as the new frame (meaning stdSkyFrameSet's current frame),
    // because, unlike stdSkyFrameTemplate, stdSkyFrameSet's current frame has inherited some
    // potentially useful attributes from the old sky frame, such as epoch.
    rawFrameSet->addFrame(initialSkyIndex, *stdSkyFrameSet->getMapping()->simplified(),
                          *stdSkyFrameSet->getFrame(ast::FrameSet::CURRENT));
    auto const stdSkyIndex = rawFrameSet->getCurrent();
    auto const stdSkyFrame = rawFrameSet->getFrame(stdSkyIndex, false);
 
    // Compute a mapping from PIXELS (0-based in parent coordinates)
    // to GRID (1-based in local coordinates)
    auto xy0 = getImageXY0FromMetadata(metadata, image::detail::wcsNameForXY0, strip);
    std::vector<double> pixelToGridArray = {1.0 - xy0[0], 1.0 - xy0[1]};  // 1.0 for FITS vs LSST convention
    auto pixelToGrid = ast::ShiftMap(pixelToGridArray);
 
    // Now construct the returned FrameDict
    auto const gridToIwc = rawFrameSet->getMapping(gridIndex, iwcIndex)->simplified();
    auto const pixelToIwc = pixelToGrid.then(*gridToIwc).simplified();
    auto const iwcToStdSky = rawFrameSet->getMapping(iwcIndex, stdSkyIndex);
 
    auto frameDict = std::make_shared<ast::FrameDict>(ast::Frame(2, "Domain=PIXELS"), *pixelToIwc, *iwcFrame);
    frameDict->addFrame("IWC", *iwcToStdSky, *stdSkyFrame);
 
    // Record CRVAL as SkyRef in the SkyFrame so it can easily be obtained later;
    // set SkyRefIs = "Ignored" (the default) so SkyRef value is ignored instead of used as an offset
    auto crpixPixels = pixelToGrid.applyInverse(crpixGrid);
    auto crvalRad = frameDict->applyForward(crpixPixels);
    auto skyFrame = std::dynamic_pointer_cast<ast::SkyFrame>(frameDict->getFrame("SKY", false));
    if (!skyFrame) {
        throw LSST_EXCEPT(pex::exceptions::LogicError, "SKY frame is not a SkyFrame");
    }
    skyFrame->setSkyRefIs("Ignored");
    skyFrame->setSkyRef(crvalRad);
 
    return frameDict;
}

readStream()

template<class Transform >

std::shared_ptr< Transform > lsst::afw::geom::detail::readStream

(

std::istream &

is

)

Deserialize a Transform from an input stream.

Template Parameters

Transform

the Transform class; can be Transform<FromEndpoint, ToEndpoint>, SkyWcs, or any other compatible class, i.e. it must support the following (see Transform.h for details): a constructor that takes an ast::FrameSet static method getShortClassName method getMapping

Parameters

[in]

is

input stream from which to deserialize this Transform

Definition at line 79 of file transformUtils.h.

                                                    {
    int version;
    is >> version;
    if (version != 1) {
        throw LSST_EXCEPT(pex::exceptions::InvalidParameterError,
                          "Unsupported version " + std::to_string(version));
    }
    std::string shortClassName;
    is >> shortClassName;
    if (shortClassName != Transform::getShortClassName()) {
        std::ostringstream os;
        os << "Class name in stream " << shortClassName << " != " << Transform::getShortClassName();
        throw LSST_EXCEPT(pex::exceptions::InvalidParameterError, os.str());
    }
    auto astStream = ast::Stream(&is, nullptr);
    auto astObjectPtr = ast::Channel(astStream).read();
    auto mapping = std::dynamic_pointer_cast<ast::Mapping>(astObjectPtr);
    if (!mapping) {
        std::ostringstream os;
        os << "The AST serialization was read as a " << astObjectPtr->getClassName()
           << " instead of a Mapping";
        throw LSST_EXCEPT(pex::exceptions::InvalidParameterError, os.str());
    }
 
    return std::make_shared<Transform>(*mapping);
}

stripWcsMetadata()

void lsst::afw::geom::detail::stripWcsMetadata

(

daf::base::PropertySet &

metadata

)

Strip all WCS metadata from a header.

Parameters

[in,out]

metadata

FITS header cards

Definition at line 134 of file frameSetUtils.cc.

                                                    {
    // Exclude WCS A keywords because LSST uses them to store XY0
    auto wcsANames = createTrivialWcsMetadata("A", lsst::geom::Point2I(0, 0))->names();
    std::set<std::string> excludeNames(wcsANames.begin(), wcsANames.end());
    // Ignore NAXIS1, NAXIS2 because if they are zero then AST will fail to read a WCS
    // Ignore LTV1/2 because LSST adds it and this code should ignore it and not strip it
    // Exclude comments and history to reduce clutter
    std::set<std::string> moreNames{"NAXIS1", "NAXIS2", "LTV1", "LTV2", "COMMENT", "HISTORY"};
    excludeNames.insert(moreNames.begin(), moreNames.end());
 
    // Replace RADECSYS with RADESYS if only the former is present
    if (metadata.exists("RADECSYS") && !metadata.exists("RADESYS")) {
        metadata.set("RADESYS", metadata.getAsString("RADECSYS"));
        metadata.remove("RADECSYS");
    }
 
    auto channel = getFitsChanFromPropertyList(metadata, excludeNames,
                                               "Encoding=FITS-WCS, IWC=1, SipReplace=0, ReportLevel=3");
    auto const initialNames = setFromVector(channel.getAllCardNames());
 
    std::shared_ptr<ast::Object> obj;
    try {
        obj = channel.read();
    } catch (std::runtime_error const&) {
        return;
    }
 
    auto const finalNames = setFromVector(channel.getAllCardNames());
 
    // FITS keywords that FitsChan stripped
    std::set<std::string> namesChannelStripped;
    std::set_difference(initialNames.begin(), initialNames.end(), finalNames.begin(), finalNames.end(),
                        std::inserter(namesChannelStripped, namesChannelStripped.begin()));
 
    // FITS keywords that FitsChan may strip that we want to keep in `metadata`
    std::set<std::string> const namesToKeep = {"DATE-OBS", "MJD-OBS"};
 
    std::set<std::string> namesToStrip;  // names to strip from metadata
    std::set_difference(namesChannelStripped.begin(), namesChannelStripped.end(), namesToKeep.begin(),
                        namesToKeep.end(), std::inserter(namesToStrip, namesToStrip.begin()));
    for (auto const& name : namesToStrip) {
        metadata.remove(name);
    }
}

wrapFrameSetUtils()

void lsst::afw::geom::detail::wrapFrameSetUtils

(

lsst::cpputils::python::WrapperCollection &

)

writeStream()

template<class Transform >

void lsst::afw::geom::detail::writeStream	(	Transform const &	transform,
		std::ostream &	os )

Serialize a Transform to an output stream.

Version 1 format is as follows:

• The version number (an integer)
• A space
• The short class name, as obtained from getShortClassName
• A space
• The contained ast::FrameSet written using FrameSet.show(os, false)

Parameters

[out]	os	output stream to which to serialize this Transform
[in]	transform	Transform to serialize

Definition at line 107 of file transformUtils.h.

                                                             {
    os << serializationVersion << " " << Transform::getShortClassName();
    transform.getMapping()->show(os, false);  // false = do not write comments
}

Variable Documentation

serializationVersion

constexpr int lsst::afw::geom::detail::serializationVersion = 1

constexpr

version of serialization used when writing (older versions may also be supported when reading)

Definition at line 42 of file transformUtils.h.