LSST Applications  21.0.0-147-g0e635eb1+1acddb5be5,22.0.0+052faf71bd,22.0.0+1ea9a8b2b2,22.0.0+6312710a6c,22.0.0+729191ecac,22.0.0+7589c3a021,22.0.0+9f079a9461,22.0.1-1-g7d6de66+b8044ec9de,22.0.1-1-g87000a6+536b1ee016,22.0.1-1-g8e32f31+6312710a6c,22.0.1-10-gd060f87+016f7cdc03,22.0.1-12-g9c3108e+df145f6f68,22.0.1-16-g314fa6d+c825727ab8,22.0.1-19-g93a5c75+d23f2fb6d8,22.0.1-19-gb93eaa13+aab3ef7709,22.0.1-2-g8ef0a89+b8044ec9de,22.0.1-2-g92698f7+9f079a9461,22.0.1-2-ga9b0f51+052faf71bd,22.0.1-2-gac51dbf+052faf71bd,22.0.1-2-gb66926d+6312710a6c,22.0.1-2-gcb770ba+09e3807989,22.0.1-20-g32debb5+b8044ec9de,22.0.1-23-gc2439a9a+fb0756638e,22.0.1-3-g496fd5d+09117f784f,22.0.1-3-g59f966b+1e6ba2c031,22.0.1-3-g849a1b8+f8b568069f,22.0.1-3-gaaec9c0+c5c846a8b1,22.0.1-32-g5ddfab5d3+60ce4897b0,22.0.1-4-g037fbe1+64e601228d,22.0.1-4-g8623105+b8044ec9de,22.0.1-5-g096abc9+d18c45d440,22.0.1-5-g15c806e+57f5c03693,22.0.1-7-gba73697+57f5c03693,master-g6e05de7fdc+c1283a92b8,master-g72cdda8301+729191ecac,w.2021.39
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 <vector>
29 #include <cmath>
30 
31 #include "boost/format.hpp"
32 
33 #include "astshim.h"
34 #include "lsst/geom/Point.h"
35 #include "lsst/afw/geom/wcsUtils.h"
37 #include "lsst/afw/image/ImageBase.h" // for wcsNameForXY0
40 #include "lsst/pex/exceptions.h"
41 #include "lsst/log/Log.h"
42 
43 namespace lsst {
44 namespace afw {
45 namespace geom {
46 namespace detail {
47 namespace {
48 
49 // destructively make a set of strings from a vector of strings
52 }
53 
54 /*
55  * Copy a FITS header card from a FitsChan to a PropertyList.
56  *
57  * Internal function for use by getPropertyListFromFitsChan.
58  *
59  * @param[in,out] metadata PropertyList to which to copy the value
60  * @param[in] name FITS header card name; used as the name for the new entry in `metadata`
61  * @param[in] foundValue Value and found flag returned by ast::FitsChan.getFits{X};
62  * foundValue.found must be true.
63  * @param[in] comment Card comment; if blank then no comment is written
64  *
65  * @throw lsst::pex::exceptions::LogicError if foundValue.found false.
66  */
67 template <typename T>
68 void setMetadataFromFoundValue(daf::base::PropertyList& metadata, std::string const& name,
69  ast::FoundValue<T> const& foundValue, std::string const& comment = "") {
70  if (!foundValue.found) {
71  throw LSST_EXCEPT(pex::exceptions::LogicError, "Bug! FitsChan card \"" + name + "\" not found");
72  }
73  if (comment.empty()) {
74  metadata.set(name, foundValue.value);
75  } else {
76  metadata.set(name, foundValue.value, comment);
77  }
78 }
79 
80 } // namespace
81 
83  // Exclude WCS A keywords because LSST uses them to store XY0
84  auto wcsANames = createTrivialWcsMetadata("A", lsst::geom::Point2I(0, 0))->names();
85  std::set<std::string> excludeNames(wcsANames.begin(), wcsANames.end());
86  // Ignore NAXIS1, NAXIS2 because if they are zero then AST will fail to read a WCS
87  // Ignore LTV1/2 because LSST adds it and this code should ignore it and not strip it
88  // Exclude comments and history to reduce clutter
89  std::set<std::string> moreNames{"NAXIS1", "NAXIS2", "LTV1", "LTV2", "COMMENT", "HISTORY"};
90  excludeNames.insert(moreNames.begin(), moreNames.end());
91 
92  // Replace RADECSYS with RADESYS if only the former is present
93  if (metadata.exists("RADECSYS") && !metadata.exists("RADESYS")) {
94  metadata.set("RADESYS", metadata.getAsString("RADECSYS"));
95  metadata.remove("RADECSYS");
96  }
97 
98  auto channel = getFitsChanFromPropertyList(metadata, excludeNames,
99  "Encoding=FITS-WCS, IWC=1, SipReplace=0, ReportLevel=3");
100  auto const initialNames = strip ? setFromVector(channel.getAllCardNames()) : std::set<std::string>();
102  try {
103  obj = channel.read();
104  } catch (std::runtime_error const&) {
105  throw LSST_EXCEPT(pex::exceptions::TypeError, "The metadata does not describe an AST object");
106  }
107 
108  auto frameSet = std::dynamic_pointer_cast<ast::FrameSet>(obj);
109  if (!frameSet) {
111  "metadata describes a " + obj->getClassName() + ", not a FrameSet");
112  }
113  if (strip) {
114  auto const finalNames = setFromVector(channel.getAllCardNames());
115 
116  // FITS keywords that FitsChan stripped
117  std::set<std::string> namesChannelStripped;
118  std::set_difference(initialNames.begin(), initialNames.end(), finalNames.begin(), finalNames.end(),
119  std::inserter(namesChannelStripped, namesChannelStripped.begin()));
120 
121  // FITS keywords that FitsChan may strip that we want to keep in `metadata`
122  std::set<std::string> const namesToKeep = {"DATE-OBS", "MJD-OBS"};
123 
124  std::set<std::string> namesToStrip; // names to strip from metadata
125  std::set_difference(namesChannelStripped.begin(), namesChannelStripped.end(), namesToKeep.begin(),
126  namesToKeep.end(), std::inserter(namesToStrip, namesToStrip.begin()));
127  for (auto const& name : namesToStrip) {
128  metadata.remove(name);
129  }
130  }
131  return frameSet;
132 }
133 
135  // Record CRPIX in GRID coordinates
136  // so we can compute CRVAL after standardizing the SkyFrame to ICRS
137  // (that standardization is why we don't simply save CRVAL now)
138  std::vector<double> crpixGrid(2);
139  try {
140  crpixGrid[0] = metadata.getAsDouble("CRPIX1");
141  crpixGrid[1] = metadata.getAsDouble("CRPIX2");
143  // std::string used because e.what() returns a C string and two C strings cannot be added
145  e.what() + std::string("; cannot read metadata as a SkyWcs"));
146  }
147 
148  auto rawFrameSet = readFitsWcs(metadata, strip);
149  auto const initialBaseIndex = rawFrameSet->getBase();
150 
151  // Find the GRID frame
152  auto gridIndex = ast::FrameSet::NOFRAME;
153  if (rawFrameSet->findFrame(ast::Frame(2, "Domain=GRID"))) {
154  gridIndex = rawFrameSet->getCurrent();
155  } else {
156  // No appropriate GRID frame found; if the original base frame is of type Frame
157  // with 2 axes and a blank domain then use that, else give up
158  auto const baseFrame = rawFrameSet->getFrame(initialBaseIndex, false);
159  auto const baseClassName = rawFrameSet->getClassName();
160  if (baseFrame->getClassName() != "Frame") {
162  "The base frame is of type " + baseFrame->getClassName() +
163  "instead of Frame; cannot read metadata as a SkyWcs");
164  }
165  if (baseFrame->getNAxes() != 2) {
167  "The base frame has " + std::to_string(baseFrame->getNAxes()) +
168  " axes instead of 2; cannot read metadata as a SkyWcs");
169  }
170  if (baseFrame->getDomain() != "") {
172  "The base frame has domain \"" + baseFrame->getDomain() +
173  "\" instead of blank or GRID; cannot read metadata as a SkyWcs");
174  }
175  // Original base frame has a blank Domain, is of type Frame, and has 2 axes, so
176  // Set its domain to GRID, and set some other potentially useful attributes.
177  gridIndex = initialBaseIndex;
178  }
179 
180  // Find the IWC frame
181  if (!rawFrameSet->findFrame(ast::Frame(2, "Domain=IWC"))) {
182  throw LSST_EXCEPT(pex::exceptions::TypeError, "No IWC frame found; cannot read metadata as a SkyWcs");
183  }
184  auto const iwcIndex = rawFrameSet->getCurrent();
185  auto const iwcFrame = rawFrameSet->getFrame(iwcIndex);
186 
187  // Create a standard sky frame: ICRS with axis order RA, Dec
188 
189  // Create the a template for the standard sky frame
190  auto const stdSkyFrameTemplate = ast::SkyFrame("System=ICRS");
191 
192  // Locate a Frame in the target FrameSet that looks like the template
193  // and hence can be used as the original sky frame.
194  // We ignore the frame set returned by findFrame because that goes from pixels to sky,
195  // and using it would add an unwanted extra branch to our WCS; instead, later on,
196  // we compute a mapping from the old sky frame to the new sky frame and add that.
197  if (!rawFrameSet->findFrame(stdSkyFrameTemplate)) {
199  "Could not find a SkyFrame; cannot read metadata as a SkyWcs");
200  }
201  auto initialSkyIndex = rawFrameSet->getCurrent();
202 
203  // Compute a frame set that maps from the original sky frame to our desired sky frame template;
204  // this contains the mapping and sky frame we will insert into the frame set.
205  // (Temporarily set the base frame to the sky frame, because findFrame
206  // produces a mapping from base to the found frame).
207  rawFrameSet->setBase(initialSkyIndex);
208  auto stdSkyFrameSet = rawFrameSet->findFrame(stdSkyFrameTemplate);
209  if (!stdSkyFrameSet) {
211  "Bug: found a SkyFrame the first time, but not the second time");
212  }
213 
214  // Add the new mapping into rawFrameSet, connecting it to the original SkyFrame.
215  // Note: we use stdSkyFrameSet as the new frame (meaning stdSkyFrameSet's current frame),
216  // because, unlike stdSkyFrameTemplate, stdSkyFrameSet's current frame has inherited some
217  // potentially useful attributes from the old sky frame, such as epoch.
218  rawFrameSet->addFrame(initialSkyIndex, *stdSkyFrameSet->getMapping()->simplified(),
219  *stdSkyFrameSet->getFrame(ast::FrameSet::CURRENT));
220  auto const stdSkyIndex = rawFrameSet->getCurrent();
221  auto const stdSkyFrame = rawFrameSet->getFrame(stdSkyIndex, false);
222 
223  // Compute a mapping from PIXELS (0-based in parent coordinates)
224  // to GRID (1-based in local coordinates)
226  std::vector<double> pixelToGridArray = {1.0 - xy0[0], 1.0 - xy0[1]}; // 1.0 for FITS vs LSST convention
227  auto pixelToGrid = ast::ShiftMap(pixelToGridArray);
228 
229  // Now construct the returned FrameDict
230  auto const gridToIwc = rawFrameSet->getMapping(gridIndex, iwcIndex)->simplified();
231  auto const pixelToIwc = pixelToGrid.then(*gridToIwc).simplified();
232  auto const iwcToStdSky = rawFrameSet->getMapping(iwcIndex, stdSkyIndex);
233 
234  auto frameDict = std::make_shared<ast::FrameDict>(ast::Frame(2, "Domain=PIXELS"), *pixelToIwc, *iwcFrame);
235  frameDict->addFrame("IWC", *iwcToStdSky, *stdSkyFrame);
236 
237  // Record CRVAL as SkyRef in the SkyFrame so it can easily be obtained later;
238  // set SkyRefIs = "Ignored" (the default) so SkyRef value is ignored instead of used as an offset
239  auto crpixPixels = pixelToGrid.applyInverse(crpixGrid);
240  auto crvalRad = frameDict->applyForward(crpixPixels);
241  auto skyFrame = std::dynamic_pointer_cast<ast::SkyFrame>(frameDict->getFrame("SKY", false));
242  if (!skyFrame) {
243  throw LSST_EXCEPT(pex::exceptions::LogicError, "SKY frame is not a SkyFrame");
244  }
245  skyFrame->setSkyRefIs("Ignored");
246  skyFrame->setSkyRef(crvalRad);
247 
248  return frameDict;
249 }
250 
252  int const numCards = fitsChan.getNCard();
253  auto metadata = std::make_shared<daf::base::PropertyList>();
254  for (int cardNum = 1; cardNum <= numCards; ++cardNum) {
255  fitsChan.setCard(cardNum);
256  auto const cardType = fitsChan.getCardType();
257  auto const cardName = fitsChan.getCardName();
258  auto const cardComment = fitsChan.getCardComm();
259  switch (cardType) {
260  case ast::CardType::FLOAT: {
261  auto foundValue = fitsChan.getFitsF();
262  setMetadataFromFoundValue(*metadata, cardName, foundValue, cardComment);
263  break;
264  }
265  case ast::CardType::INT: {
266  auto foundValue = fitsChan.getFitsI();
267  setMetadataFromFoundValue(*metadata, cardName, foundValue, cardComment);
268  break;
269  }
270  case ast::CardType::STRING: {
271  auto foundValue = fitsChan.getFitsS();
272  setMetadataFromFoundValue(*metadata, cardName, foundValue, cardComment);
273  break;
274  }
275  case ast::CardType::LOGICAL: {
276  auto foundValue = fitsChan.getFitsL();
277  setMetadataFromFoundValue(*metadata, cardName, foundValue, cardComment);
278  break;
279  }
280  case ast::CardType::UNDEF: {
281  if (cardComment.empty()) {
282  metadata->set(cardName, nullptr);
283  } else {
284  metadata->set(cardName, nullptr, cardComment);
285  }
286  break;
287  }
289  auto foundValue = fitsChan.getFitsCN();
290  setMetadataFromFoundValue(*metadata, cardName, foundValue, cardComment);
291  break;
292  }
294  // Drop HISTORY and COMMENT cards
295  break;
298  // PropertyList supports neither complex numbers nor cards with no value
300  os << "Card " << cardNum << " with name \"" << cardName << "\" has type "
301  << static_cast<int>(cardType) << ", which is not supported by PropertyList";
303  }
304  case ast::CardType::NOTYPE: {
305  // This should only occur if cardNum is invalid, and that should be impossible
307  os << "Bug! Card " << cardNum << " with name \"" << cardName
308  << "\" has type NOTYPE, which should not be possible";
310  }
311  }
312  }
313  return metadata;
314 }
315 
317  std::set<std::string> const& excludeNames,
318  std::string options) {
319  // Create FitsChan to receive each of the parameters
320  auto stream = ast::StringStream("");
321  auto fc = ast::FitsChan(stream, options);
322 
323  // Get the parameter names (in order if necessary)
324  daf::base::PropertyList const *pl = dynamic_cast<daf::base::PropertyList const *>(&metadata);
325  std::vector<std::string> allParamNames;
326  if (pl) {
327  allParamNames = pl->getOrderedNames();
328  } else {
329  allParamNames = metadata.paramNames(false);
330  }
331 
332  // Loop over the names and add them to the FitsChan if not excluded
333  for (auto const &fullName : allParamNames) {
334  if (excludeNames.count(fullName) == 0) {
335  std::size_t lastPeriod = fullName.rfind(char('.'));
336  auto name = (lastPeriod == std::string::npos) ? fullName : fullName.substr(lastPeriod + 1);
337  std::type_info const &type = metadata.typeOf(name);
338 
339  if (name.size() > 8) {
340  continue; // The name is too long for a FITS keyword; skip this item
341  }
342 
343  if (type == typeid(bool)) {
344  fc.setFitsL(name, metadata.get<bool>(name));
345  } else if (type == typeid(std::uint8_t)) {
346  fc.setFitsI(name, static_cast<int>(metadata.get<std::uint8_t>(name)));
347  } else if (type == typeid(int)) {
348  fc.setFitsI(name, metadata.get<int>(name));
349  } else if (type == typeid(double) || type == typeid(float)) {
350  double value;
351  if (type == typeid(double)) {
352  value = metadata.get<double>(name);
353  } else {
354  value = static_cast<double>(metadata.get<float>(name));
355  }
356  // NaN is not allowed in a FitsChan (or in FITS)
357  if (!std::isnan(value)) {
358  fc.setFitsF(name, value);
359  } else {
360  // Treat it like an undefined value but warn about it
361  LOGLS_WARN("afw.geom.frameSetUtils",
362  boost::format("Found NaN in metadata item '%s'") % name);
363  }
364  } else if (type == typeid(std::string)) {
365  std::string str = metadata.get<std::string>(name);
366  // No support for long strings yet so skip those
367  if (str.size() <= 68) {
368  fc.setFitsS(name, metadata.get<std::string>(name));
369  }
370  }
371  }
372  }
373  // Rewind the channel
374  fc.setCard(0);
375  return fc;
376 }
377 
378 } // namespace detail
379 } // namespace geom
380 } // namespace afw
381 } // 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:659
std::ostream * os
Definition: Schema.cc:557
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:82
Class for storing generic metadata.
Definition: PropertySet.h:66
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:48
lsst::geom::Point2I getImageXY0FromMetadata(daf::base::PropertySet &metadata, std::string const &wcsName, bool strip=false)
Definition: wcsUtils.cc:95
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)