LSST Applications  21.0.0+04719a4bac,21.0.0-1-ga51b5d4+f5e6047307,21.0.0-11-g2b59f77+a9c1acf22d,21.0.0-11-ga42c5b2+86977b0b17,21.0.0-12-gf4ce030+76814010d2,21.0.0-13-g1721dae+760e7a6536,21.0.0-13-g3a573fe+768d78a30a,21.0.0-15-g5a7caf0+f21cbc5713,21.0.0-16-g0fb55c1+b60e2d390c,21.0.0-19-g4cded4ca+71a93a33c0,21.0.0-2-g103fe59+bb20972958,21.0.0-2-g45278ab+04719a4bac,21.0.0-2-g5242d73+3ad5d60fb1,21.0.0-2-g7f82c8f+8babb168e8,21.0.0-2-g8f08a60+06509c8b61,21.0.0-2-g8faa9b5+616205b9df,21.0.0-2-ga326454+8babb168e8,21.0.0-2-gde069b7+5e4aea9c2f,21.0.0-2-gecfae73+1d3a86e577,21.0.0-2-gfc62afb+3ad5d60fb1,21.0.0-25-g1d57be3cd+e73869a214,21.0.0-3-g357aad2+ed88757d29,21.0.0-3-g4a4ce7f+3ad5d60fb1,21.0.0-3-g4be5c26+3ad5d60fb1,21.0.0-3-g65f322c+e0b24896a3,21.0.0-3-g7d9da8d+616205b9df,21.0.0-3-ge02ed75+a9c1acf22d,21.0.0-4-g591bb35+a9c1acf22d,21.0.0-4-g65b4814+b60e2d390c,21.0.0-4-gccdca77+0de219a2bc,21.0.0-4-ge8a399c+6c55c39e83,21.0.0-5-gd00fb1e+05fce91b99,21.0.0-6-gc675373+3ad5d60fb1,21.0.0-64-g1122c245+4fb2b8f86e,21.0.0-7-g04766d7+cd19d05db2,21.0.0-7-gdf92d54+04719a4bac,21.0.0-8-g5674e7b+d1bd76f71f,master-gac4afde19b+a9c1acf22d,w.2021.13
LSST Data Management Base Package
frameSetUtils.cc
Go to the documentation of this file.
1 /*
2  * LSST Data Management System
3  * Copyright 2017 LSST Corporation.
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 <algorithm>
24 #include <exception>
25 #include <memory>
26 #include <ostream>
27 #include <set>
28 #include <sstream>
29 #include <vector>
30 #include <cmath>
31 
32 #include "boost/format.hpp"
33 
34 #include "astshim.h"
35 
37 #include "lsst/afw/fits.h"
38 #include "lsst/geom/Angle.h"
39 #include "lsst/geom/Point.h"
40 #include "lsst/afw/geom/wcsUtils.h"
42 #include "lsst/afw/image/ImageBase.h" // for wcsNameForXY0
45 #include "lsst/pex/exceptions.h"
46 #include "lsst/log/Log.h"
47 
48 namespace lsst {
49 namespace afw {
50 namespace geom {
51 namespace detail {
52 namespace {
53 
54 // destructively make a set of strings from a vector of strings
57 }
58 
59 /*
60  * Copy a FITS header card from a FitsChan to a PropertyList.
61  *
62  * Internal function for use by getPropertyListFromFitsChan.
63  *
64  * @param[in,out] metadata PropertyList to which to copy the value
65  * @param[in] name FITS header card name; used as the name for the new entry in `metadata`
66  * @param[in] foundValue Value and found flag returned by ast::FitsChan.getFits{X};
67  * foundValue.found must be true.
68  * @param[in] comment Card comment; if blank then no comment is written
69  *
70  * @throw lsst::pex::exceptions::LogicError if foundValue.found false.
71  */
72 template <typename T>
73 void setMetadataFromFoundValue(daf::base::PropertyList& metadata, std::string const& name,
74  ast::FoundValue<T> const& foundValue, std::string const& comment = "") {
75  if (!foundValue.found) {
76  throw LSST_EXCEPT(pex::exceptions::LogicError, "Bug! FitsChan card \"" + name + "\" not found");
77  }
78  if (comment.empty()) {
79  metadata.set(name, foundValue.value);
80  } else {
81  metadata.set(name, foundValue.value, comment);
82  }
83 }
84 
85 } // namespace
86 
88  // Exclude WCS A keywords because LSST uses them to store XY0
89  auto wcsANames = createTrivialWcsMetadata("A", lsst::geom::Point2I(0, 0))->names();
90  std::set<std::string> excludeNames(wcsANames.begin(), wcsANames.end());
91  // Ignore NAXIS1, NAXIS2 because if they are zero then AST will fail to read a WCS
92  // Ignore LTV1/2 because LSST adds it and this code should ignore it and not strip it
93  // Exclude comments and history to reduce clutter
94  std::set<std::string> moreNames{"NAXIS1", "NAXIS2", "LTV1", "LTV2", "COMMENT", "HISTORY"};
95  excludeNames.insert(moreNames.begin(), moreNames.end());
96 
97  // Replace RADECSYS with RADESYS if only the former is present
98  if (metadata.exists("RADECSYS") && !metadata.exists("RADESYS")) {
99  metadata.set("RADESYS", metadata.getAsString("RADECSYS"));
100  metadata.remove("RADECSYS");
101  }
102 
103  auto channel = getFitsChanFromPropertyList(metadata, excludeNames,
104  "Encoding=FITS-WCS, IWC=1, SipReplace=0, ReportLevel=3");
105  auto const initialNames = strip ? setFromVector(channel.getAllCardNames()) : std::set<std::string>();
107  try {
108  obj = channel.read();
109  } catch (std::runtime_error const&) {
110  throw LSST_EXCEPT(pex::exceptions::TypeError, "The metadata does not describe an AST object");
111  }
112 
113  auto frameSet = std::dynamic_pointer_cast<ast::FrameSet>(obj);
114  if (!frameSet) {
116  "metadata describes a " + obj->getClassName() + ", not a FrameSet");
117  }
118  if (strip) {
119  auto const finalNames = setFromVector(channel.getAllCardNames());
120 
121  // FITS keywords that FitsChan stripped
122  std::set<std::string> namesChannelStripped;
123  std::set_difference(initialNames.begin(), initialNames.end(), finalNames.begin(), finalNames.end(),
124  std::inserter(namesChannelStripped, namesChannelStripped.begin()));
125 
126  // FITS keywords that FitsChan may strip that we want to keep in `metadata`
127  std::set<std::string> const namesToKeep = {"DATE-OBS", "MJD-OBS"};
128 
129  std::set<std::string> namesToStrip; // names to strip from metadata
130  std::set_difference(namesChannelStripped.begin(), namesChannelStripped.end(), namesToKeep.begin(),
131  namesToKeep.end(), std::inserter(namesToStrip, namesToStrip.begin()));
132  for (auto const& name : namesToStrip) {
133  metadata.remove(name);
134  }
135  }
136  return frameSet;
137 }
138 
140  // Record CRPIX in GRID coordinates
141  // so we can compute CRVAL after standardizing the SkyFrame to ICRS
142  // (that standardization is why we don't simply save CRVAL now)
143  std::vector<double> crpixGrid(2);
144  try {
145  crpixGrid[0] = metadata.getAsDouble("CRPIX1");
146  crpixGrid[1] = metadata.getAsDouble("CRPIX2");
148  // std::string used because e.what() returns a C string and two C strings cannot be added
150  e.what() + std::string("; cannot read metadata as a SkyWcs"));
151  }
152 
153  auto rawFrameSet = readFitsWcs(metadata, strip);
154  auto const initialBaseIndex = rawFrameSet->getBase();
155 
156  // Find the GRID frame
157  auto gridIndex = ast::FrameSet::NOFRAME;
158  if (rawFrameSet->findFrame(ast::Frame(2, "Domain=GRID"))) {
159  gridIndex = rawFrameSet->getCurrent();
160  } else {
161  // No appropriate GRID frame found; if the original base frame is of type Frame
162  // with 2 axes and a blank domain then use that, else give up
163  auto const baseFrame = rawFrameSet->getFrame(initialBaseIndex, false);
164  auto const baseClassName = rawFrameSet->getClassName();
165  if (baseFrame->getClassName() != "Frame") {
167  "The base frame is of type " + baseFrame->getClassName() +
168  "instead of Frame; cannot read metadata as a SkyWcs");
169  }
170  if (baseFrame->getNAxes() != 2) {
172  "The base frame has " + std::to_string(baseFrame->getNAxes()) +
173  " axes instead of 2; cannot read metadata as a SkyWcs");
174  }
175  if (baseFrame->getDomain() != "") {
177  "The base frame has domain \"" + baseFrame->getDomain() +
178  "\" instead of blank or GRID; cannot read metadata as a SkyWcs");
179  }
180  // Original base frame has a blank Domain, is of type Frame, and has 2 axes, so
181  // Set its domain to GRID, and set some other potentially useful attributes.
182  gridIndex = initialBaseIndex;
183  }
184 
185  // Find the IWC frame
186  if (!rawFrameSet->findFrame(ast::Frame(2, "Domain=IWC"))) {
187  throw LSST_EXCEPT(pex::exceptions::TypeError, "No IWC frame found; cannot read metadata as a SkyWcs");
188  }
189  auto const iwcIndex = rawFrameSet->getCurrent();
190  auto const iwcFrame = rawFrameSet->getFrame(iwcIndex);
191 
192  // Create a standard sky frame: ICRS with axis order RA, Dec
193 
194  // Create the a template for the standard sky frame
195  auto const stdSkyFrameTemplate = ast::SkyFrame("System=ICRS");
196 
197  // Locate a Frame in the target FrameSet that looks like the template
198  // and hence can be used as the original sky frame.
199  // We ignore the frame set returned by findFrame because that goes from pixels to sky,
200  // and using it would add an unwanted extra branch to our WCS; instead, later on,
201  // we compute a mapping from the old sky frame to the new sky frame and add that.
202  if (!rawFrameSet->findFrame(stdSkyFrameTemplate)) {
204  "Could not find a SkyFrame; cannot read metadata as a SkyWcs");
205  }
206  auto initialSkyIndex = rawFrameSet->getCurrent();
207 
208  // Compute a frame set that maps from the original sky frame to our desired sky frame template;
209  // this contains the mapping and sky frame we will insert into the frame set.
210  // (Temporarily set the base frame to the sky frame, because findFrame
211  // produces a mapping from base to the found frame).
212  rawFrameSet->setBase(initialSkyIndex);
213  auto stdSkyFrameSet = rawFrameSet->findFrame(stdSkyFrameTemplate);
214  if (!stdSkyFrameSet) {
216  "Bug: found a SkyFrame the first time, but not the second time");
217  }
218 
219  // Add the new mapping into rawFrameSet, connecting it to the original SkyFrame.
220  // Note: we use stdSkyFrameSet as the new frame (meaning stdSkyFrameSet's current frame),
221  // because, unlike stdSkyFrameTemplate, stdSkyFrameSet's current frame has inherited some
222  // potentially useful attributes from the old sky frame, such as epoch.
223  rawFrameSet->addFrame(initialSkyIndex, *stdSkyFrameSet->getMapping()->simplified(),
224  *stdSkyFrameSet->getFrame(ast::FrameSet::CURRENT));
225  auto const stdSkyIndex = rawFrameSet->getCurrent();
226  auto const stdSkyFrame = rawFrameSet->getFrame(stdSkyIndex, false);
227 
228  // Compute a mapping from PIXELS (0-based in parent coordinates)
229  // to GRID (1-based in local coordinates)
231  std::vector<double> pixelToGridArray = {1.0 - xy0[0], 1.0 - xy0[1]}; // 1.0 for FITS vs LSST convention
232  auto pixelToGrid = ast::ShiftMap(pixelToGridArray);
233 
234  // Now construct the returned FrameDict
235  auto const gridToIwc = rawFrameSet->getMapping(gridIndex, iwcIndex)->simplified();
236  auto const pixelToIwc = pixelToGrid.then(*gridToIwc).simplified();
237  auto const iwcToStdSky = rawFrameSet->getMapping(iwcIndex, stdSkyIndex);
238 
239  auto frameDict = std::make_shared<ast::FrameDict>(ast::Frame(2, "Domain=PIXELS"), *pixelToIwc, *iwcFrame);
240  frameDict->addFrame("IWC", *iwcToStdSky, *stdSkyFrame);
241 
242  // Record CRVAL as SkyRef in the SkyFrame so it can easily be obtained later;
243  // set SkyRefIs = "Ignored" (the default) so SkyRef value is ignored instead of used as an offset
244  auto crpixPixels = pixelToGrid.applyInverse(crpixGrid);
245  auto crvalRad = frameDict->applyForward(crpixPixels);
246  auto skyFrame = std::dynamic_pointer_cast<ast::SkyFrame>(frameDict->getFrame("SKY", false));
247  if (!skyFrame) {
248  throw LSST_EXCEPT(pex::exceptions::LogicError, "SKY frame is not a SkyFrame");
249  }
250  skyFrame->setSkyRefIs("Ignored");
251  skyFrame->setSkyRef(crvalRad);
252 
253  return frameDict;
254 }
255 
257  int const numCards = fitsChan.getNCard();
258  auto metadata = std::make_shared<daf::base::PropertyList>();
259  for (int cardNum = 1; cardNum <= numCards; ++cardNum) {
260  fitsChan.setCard(cardNum);
261  auto const cardType = fitsChan.getCardType();
262  auto const cardName = fitsChan.getCardName();
263  auto const cardComment = fitsChan.getCardComm();
264  switch (cardType) {
265  case ast::CardType::FLOAT: {
266  auto foundValue = fitsChan.getFitsF();
267  setMetadataFromFoundValue(*metadata, cardName, foundValue, cardComment);
268  break;
269  }
270  case ast::CardType::INT: {
271  auto foundValue = fitsChan.getFitsI();
272  setMetadataFromFoundValue(*metadata, cardName, foundValue, cardComment);
273  break;
274  }
275  case ast::CardType::STRING: {
276  auto foundValue = fitsChan.getFitsS();
277  setMetadataFromFoundValue(*metadata, cardName, foundValue, cardComment);
278  break;
279  }
280  case ast::CardType::LOGICAL: {
281  auto foundValue = fitsChan.getFitsL();
282  setMetadataFromFoundValue(*metadata, cardName, foundValue, cardComment);
283  break;
284  }
285  case ast::CardType::UNDEF: {
286  if (cardComment.empty()) {
287  metadata->set(cardName, nullptr);
288  } else {
289  metadata->set(cardName, nullptr, cardComment);
290  }
291  break;
292  }
294  auto foundValue = fitsChan.getFitsCN();
295  setMetadataFromFoundValue(*metadata, cardName, foundValue, cardComment);
296  break;
297  }
299  // Drop HISTORY and COMMENT cards
300  break;
303  // PropertyList supports neither complex numbers nor cards with no value
305  os << "Card " << cardNum << " with name \"" << cardName << "\" has type "
306  << static_cast<int>(cardType) << ", which is not supported by PropertyList";
308  }
309  case ast::CardType::NOTYPE: {
310  // This should only occur if cardNum is invalid, and that should be impossible
312  os << "Bug! Card " << cardNum << " with name \"" << cardName
313  << "\" has type NOTYPE, which should not be possible";
315  }
316  }
317  }
318  return metadata;
319 }
320 
322  std::set<std::string> const& excludeNames,
323  std::string options) {
324  // Create FitsChan to receive each of the parameters
325  auto stream = ast::StringStream("");
326  auto fc = ast::FitsChan(stream, options);
327 
328  // Get the parameter names (in order if necessary)
329  daf::base::PropertyList const *pl = dynamic_cast<daf::base::PropertyList const *>(&metadata);
330  std::vector<std::string> allParamNames;
331  if (pl) {
332  allParamNames = pl->getOrderedNames();
333  } else {
334  allParamNames = metadata.paramNames(false);
335  }
336 
337  // Loop over the names and add them to the FitsChan if not excluded
338  for (auto const &fullName : allParamNames) {
339  if (excludeNames.count(fullName) == 0) {
340  std::size_t lastPeriod = fullName.rfind(char('.'));
341  auto name = (lastPeriod == std::string::npos) ? fullName : fullName.substr(lastPeriod + 1);
342  std::type_info const &type = metadata.typeOf(name);
343 
344  if (name.size() > 8) {
345  continue; // The name is too long for a FITS keyword; skip this item
346  }
347 
348  if (type == typeid(bool)) {
349  fc.setFitsL(name, metadata.get<bool>(name));
350  } else if (type == typeid(std::uint8_t)) {
351  fc.setFitsI(name, static_cast<int>(metadata.get<std::uint8_t>(name)));
352  } else if (type == typeid(int)) {
353  fc.setFitsI(name, metadata.get<int>(name));
354  } else if (type == typeid(double) || type == typeid(float)) {
355  double value;
356  if (type == typeid(double)) {
357  value = metadata.get<double>(name);
358  } else {
359  value = static_cast<double>(metadata.get<float>(name));
360  }
361  // NaN is not allowed in a FitsChan (or in FITS)
362  if (!std::isnan(value)) {
363  fc.setFitsF(name, value);
364  } else {
365  // Treat it like an undefined value but warn about it
366  LOGLS_WARN("afw.geom.frameSetUtils",
367  boost::format("Found NaN in metadata item '%s'") % name);
368  }
369  } else if (type == typeid(std::string)) {
370  std::string str = metadata.get<std::string>(name);
371  // No support for long strings yet so skip those
372  if (str.size() <= 68) {
373  fc.setFitsS(name, metadata.get<std::string>(name));
374  }
375  }
376  }
377  }
378  // Rewind the channel
379  fc.setCard(0);
380  return fc;
381 }
382 
383 } // namespace detail
384 } // namespace geom
385 } // namespace afw
386 } // namespace lsst
table::Key< std::string > name
Definition: Amplifier.cc:116
table::Key< int > type
Definition: Detector.cc:163
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition: Exception.h:48
LSST DM logging module built on log4cxx.
#define LOGLS_WARN(logger, message)
Log a warn-level message using an iostream-based interface.
Definition: Log.h:648
std::ostream * os
Definition: Schema.cc:746
table::Key< table::Array< std::uint8_t > > frameSet
T begin(T... args)
A specialized form of Channel which reads and writes FITS header cards.
Definition: FitsChan.h:202
std::string getCardComm() const
Get CardComm: the comment of the current card.
Definition: FitsChan.h:491
FoundValue< std::string > getFitsS(std::string const &name="", std::string defval="") const
Get the value of a string card.
Definition: FitsChan.cc:98
void setCard(int ind)
Set Card: the index of the current card, where 1 is the first card.
Definition: FitsChan.h:1036
FoundValue< bool > getFitsL(std::string const &name="", bool defval=false) const
Get the value of a bool card.
Definition: FitsChan.cc:91
std::string getCardName() const
Get CardName: the keyword name of the current card.
Definition: FitsChan.h:496
FoundValue< int > getFitsI(std::string const &name="", int defval=0) const
Get the value of a int card.
Definition: FitsChan.cc:84
CardType getCardType() const
Get CardType: data type of the current FITS card.
Definition: FitsChan.h:501
int getNCard() const
Get NCard: the number of cards.
Definition: FitsChan.h:555
FoundValue< double > getFitsF(std::string const &name="", double defval=0) const
Get the value of a double card.
Definition: FitsChan.cc:77
FoundValue< std::string > getFitsCN(std::string const &name="", std::string defval="") const
Get the value of a CONTINUE card.
Definition: FitsChan.cc:69
A value and associated validity flag.
Definition: FitsChan.h:69
T value
The found value; ignore if found is false.
Definition: FitsChan.h:82
bool found
Was the value found?
Definition: FitsChan.h:81
Frame is used to represent a coordinate system.
Definition: Frame.h:157
static constexpr int CURRENT
index of current frame
Definition: FrameSet.h:105
static constexpr int NOFRAME
an invalid frame index
Definition: FrameSet.h:106
ShiftMap is a linear Mapping which shifts each axis by a specified constant value.
Definition: ShiftMap.h:40
SkyFrame is a specialised form of Frame which describes celestial longitude/latitude coordinate syste...
Definition: SkyFrame.h:66
String-based source and sink for channels.
Definition: Stream.h:180
Class for storing ordered metadata with comments.
Definition: PropertyList.h:68
std::vector< std::string > getOrderedNames() const
Get the list of property names, in the order they were added.
Definition: PropertyList.cc:81
Class for storing generic metadata.
Definition: PropertySet.h:67
std::string getAsString(std::string const &name) const
Get the last value for a string property name (possibly hierarchical).
virtual void remove(std::string const &name)
Remove all values for a property name (possibly hierarchical).
bool exists(std::string const &name) const
Determine if a name (possibly hierarchical) exists.
std::vector< std::string > paramNames(bool topLevelOnly=true) const
A variant of names that excludes the names of subproperties.
void set(std::string const &name, T const &value)
Replace all values for a property name (possibly hierarchical) with a new scalar value.
double getAsDouble(std::string const &name) const
Get the last value for any arithmetic property name (possibly hierarchical).
std::type_info const & typeOf(std::string const &name) const
Get the type of values for a property name (possibly hierarchical).
T get(std::string const &name) const
Get the last value for a property name (possibly hierarchical).
virtual char const * what(void) const noexcept
Return a character string summarizing this exception.
Definition: Exception.cc:99
Reports errors in the logical structure of the program.
Definition: Runtime.h:46
Reports attempts to access elements using an invalid key.
Definition: Runtime.h:151
Reports errors from accepting an object of an unexpected or inappropriate type.
Definition: Runtime.h:167
T count(T... args)
T end(T... args)
bool strip
Definition: fits.cc:911
T insert(T... args)
T inserter(T... args)
T isnan(T... args)
T make_move_iterator(T... args)
@ NOTYPE
card does not exist (card number invalid)
@ CONTINUE
CONTINUE card.
@ LOGICAL
boolean
@ COMPLEXI
complex integer
@ INT
integer
@ COMPLEXF
complex floating point
@ STRING
string
@ UNDEF
card has no value
@ COMMENT
card is a comment-style card with no "=" (COMMENT, HISTORY, ...)
std::shared_ptr< ast::FrameSet > readFitsWcs(daf::base::PropertySet &metadata, bool strip=true)
Read a FITS convention WCS FrameSet from FITS metadata.
ast::FitsChan getFitsChanFromPropertyList(daf::base::PropertySet &metadata, std::set< std::string > const &excludeNames={}, std::string options="")
Construct AST FitsChan from PropertyList.
std::shared_ptr< ast::FrameDict > readLsstSkyWcs(daf::base::PropertySet &metadata, bool strip=true)
Read an LSST celestial WCS FrameDict from a FITS header.
std::shared_ptr< daf::base::PropertyList > getPropertyListFromFitsChan(ast::FitsChan &fitsChan)
Copy values from an AST FitsChan into a PropertyList.
std::shared_ptr< daf::base::PropertyList > createTrivialWcsMetadata(std::string const &wcsName, lsst::geom::Point2I const &xy0)
Definition: wcsUtils.cc:51
lsst::geom::Point2I getImageXY0FromMetadata(daf::base::PropertySet &metadata, std::string const &wcsName, bool strip=false)
Definition: wcsUtils.cc:98
std::string const wcsNameForXY0
Definition: ImageBase.h:70
def format(config, name=None, writeSourceLine=True, prefix="", verbose=False)
Definition: history.py:174
A base class for image defects.
T set_difference(T... args)
T size(T... args)
T to_string(T... args)