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
OldWcs.cc
Go to the documentation of this file.
1 /*
2  * LSST Data Management System
3  * Copyright 2017 AURA/LSST.
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 
27 #include "astshim.h"
28 
35 #include "lsst/afw/geom/wcsUtils.h"
36 #include "lsst/afw/geom/SkyWcs.h"
37 
38 namespace lsst {
39 namespace afw {
40 namespace geom {
41 namespace {
42 
43 /*
44  * Does one string end with another?
45  *
46  * From https://stackoverflow.com/a/2072890
47  */
48 inline bool endsWith(std::string const& value, std::string const& suffix) {
49  if (suffix.size() > value.size()) {
50  return false;
51  }
52  return std::equal(suffix.rbegin(), suffix.rend(), value.rbegin());
53 }
54 
55 /*
56  * @internal Get FITS WCS metadata from a record for an lsst::afw::image::Wcs
57  * or the non-SIP terms of an lsst::afw::image::TanWcs
58  *
59  * @param[in] record Record holding data for an old lsst::afw::image::Wcs
60  * @return FITS WCS metadata
61  */
62 std::shared_ptr<daf::base::PropertyList> getOldWcsMetadata(table::BaseRecord const& record);
63 
64 /*
65  * @internal Get FITS WCS metadata for the SIP terms of a TAN-SIP WCS
66  * saved by lsst::afw::image::TanWcs
67  *
68  * @param[in] record Record holding SIP terms for an old lsst::afw::image::TanWcs
69  * @return FITS WCS metadata
70  */
71 std::shared_ptr<daf::base::PropertyList> getOldSipMetadata(table::BaseRecord const& record);
72 
73 // Unpersist a SkyWcs from table data for an afw::image::Wcs
74 class OldWcsFactory : public table::io::PersistableFactory {
75 public:
76  explicit OldWcsFactory(std::string const& name) : table::io::PersistableFactory(name) {}
77 
78  std::shared_ptr<table::io::Persistable> read(InputArchive const& archive,
79  CatalogVector const& catalogs) const override {
80  LSST_ARCHIVE_ASSERT(catalogs.size() >= 1u);
81  LSST_ARCHIVE_ASSERT(catalogs.front().size() == 1u);
82  auto const& record = catalogs.front().front();
83  auto const metadata = getOldWcsMetadata(record);
84  return std::make_shared<SkyWcs>(*metadata);
85  }
86 };
87 
88 OldWcsFactory registerWcs("Wcs");
89 
90 // Unpersist a SkyWcs from table data for an afw::image::TanWcs, which may have a record for SIP terms
91 class OldTanWcsFactory : public table::io::PersistableFactory {
92 public:
93  explicit OldTanWcsFactory(std::string const& name) : table::io::PersistableFactory(name) {}
94 
95  std::shared_ptr<table::io::Persistable> read(InputArchive const& archive,
96  CatalogVector const& catalogs) const override {
97  LSST_ARCHIVE_ASSERT(catalogs.size() >= 1u);
98  auto const& record = catalogs.front().front();
99  auto const metadata = getOldWcsMetadata(record);
100 
101  if (catalogs.size() > 1u) {
102  LSST_ARCHIVE_ASSERT(catalogs.size() == 2u);
103  LSST_ARCHIVE_ASSERT(catalogs.front().size() == 1u);
104  LSST_ARCHIVE_ASSERT(catalogs.back().size() == 1u);
105  std::shared_ptr<table::BaseRecord const> sipRecord = catalogs.back().begin();
106 
107  // if CTYPE1 or CTYPE2 does not end with -SIP (and likely it never does), append it
108  std::string const ctypeSuffix = "-SIP";
109  for (auto i = 1; i <= 2; ++i) {
110  std::string const ctypeName = "CTYPE" + std::to_string(i);
111  auto const ctypeValue = metadata->getAsString(ctypeName);
112  if (!endsWith(ctypeValue, ctypeSuffix)) {
113  metadata->set(ctypeName, ctypeValue + ctypeSuffix);
114  }
115  }
116 
117  auto sipMetadata = getOldSipMetadata(*sipRecord);
118  metadata->combine(sipMetadata);
119  }
120  return std::make_shared<SkyWcs>(*metadata);
121  }
122 };
123 
124 OldTanWcsFactory registerTanWcs("TanWcs");
125 
126 // Read-only singleton struct containing the schema and keys that lsst::afw::image::Wcs
127 // and the non-SIP portion of lsst::afw::image::TanWcs were mapped to.
128 struct OldWcsPersistenceHelper {
129  table::Schema schema;
130  table::PointKey<double> crval;
131  table::PointKey<double> crpix;
132  table::Key<table::Array<double>> cd;
133  table::Key<std::string> ctype1;
134  table::Key<std::string> ctype2;
135  table::Key<double> equinox;
136  table::Key<std::string> radesys;
137  table::Key<std::string> cunit1;
138  table::Key<std::string> cunit2;
139 
140  static OldWcsPersistenceHelper const& get() {
141  static OldWcsPersistenceHelper instance;
142  return instance;
143  };
144 
145  // No copying
146  OldWcsPersistenceHelper(const OldWcsPersistenceHelper&) = delete;
147  OldWcsPersistenceHelper& operator=(const OldWcsPersistenceHelper&) = delete;
148 
149  // No moving
150  OldWcsPersistenceHelper(OldWcsPersistenceHelper&&) = delete;
151  OldWcsPersistenceHelper& operator=(OldWcsPersistenceHelper&&) = delete;
152 
153 private:
154  OldWcsPersistenceHelper()
155  : schema(),
156  crval(table::PointKey<double>::addFields(schema, "crval", "celestial reference point", "deg")),
157  crpix(table::PointKey<double>::addFields(schema, "crpix", "pixel reference point", "pixel")),
158  cd(schema.addField<table::Array<double>>(
159  "cd", "linear transform matrix, ordered (1_1, 2_1, 1_2, 2_2)", 4)),
160  ctype1(schema.addField<std::string>("ctype1", "coordinate type", 72)),
161  ctype2(schema.addField<std::string>("ctype2", "coordinate type", 72)),
162  equinox(schema.addField<double>("equinox", "equinox of coordinates")),
163  radesys(schema.addField<std::string>("radesys", "coordinate system for equinox", 72)),
164  cunit1(schema.addField<std::string>("cunit1", "coordinate units", 72)),
165  cunit2(schema.addField<std::string>("cunit2", "coordinate units", 72)) {}
166 };
167 
168 std::shared_ptr<daf::base::PropertyList> getOldWcsMetadata(table::BaseRecord const& record) {
169  auto const& keys = OldWcsPersistenceHelper::get();
170  LSST_ARCHIVE_ASSERT(record.getSchema() == keys.schema);
171 
172  auto metadata = std::make_shared<daf::base::PropertyList>();
173 
174  auto crvalDeg = record.get(keys.crval);
175  auto crpix = record.get(keys.crpix);
176  auto cd = record.get(keys.cd);
177  metadata->set("CRVAL1", crvalDeg[0]);
178  metadata->set("CRVAL2", crvalDeg[1]);
179  // Add 1 to CRPIX because the field was saved using the LSST standard: 0,0 is lower left pixel
180  metadata->set("CRPIX1", crpix[0] + 1);
181  metadata->set("CRPIX2", crpix[1] + 1);
182  metadata->set("CD1_1", cd[0]);
183  metadata->set("CD2_1", cd[1]);
184  metadata->set("CD1_2", cd[2]);
185  metadata->set("CD2_2", cd[3]);
186  metadata->set("CTYPE1", record.get(keys.ctype1));
187  metadata->set("CTYPE2", record.get(keys.ctype2));
188  metadata->set("EQUINOX", record.get(keys.equinox));
189  metadata->set("RADESYS", record.get(keys.radesys));
190  metadata->set("CUNIT1", record.get(keys.cunit1));
191  metadata->set("CUNIT2", record.get(keys.cunit2));
192  return metadata;
193 }
194 
195 std::shared_ptr<daf::base::PropertyList> getOldSipMetadata(table::BaseRecord const& record) {
196  // Cannot use a PersistenceHelper for the SIP terms because the length of each SIP array
197  // must be read from the schema
198  afw::table::Key<table::Array<double>> kA;
199  afw::table::Key<table::Array<double>> kB;
200  afw::table::Key<table::Array<double>> kAp;
201  afw::table::Key<table::Array<double>> kBp;
202  try {
203  kA = record.getSchema()["A"];
204  kB = record.getSchema()["B"];
205  kAp = record.getSchema()["Ap"];
206  kBp = record.getSchema()["Bp"];
207  } catch (...) {
208  throw LSST_EXCEPT(afw::table::io::MalformedArchiveError,
209  "Incorrect schema for TanWcs distortion terms");
210  }
211 
212  // Adding 0.5 and truncating the result guarantees we'll get the right answer
213  // for small ints even when round-off error is involved.
214  int nA = static_cast<int>(std::sqrt(kA.getSize() + 0.5));
215  int nB = static_cast<int>(std::sqrt(kB.getSize() + 0.5));
216  int nAp = static_cast<int>(std::sqrt(kAp.getSize() + 0.5));
217  int nBp = static_cast<int>(std::sqrt(kBp.getSize() + 0.5));
218  if (nA * nA != kA.getSize()) {
219  throw LSST_EXCEPT(table::io::MalformedArchiveError, "Forward X SIP matrix is not square.");
220  }
221  if (nB * nB != kB.getSize()) {
222  throw LSST_EXCEPT(table::io::MalformedArchiveError, "Forward Y SIP matrix is not square.");
223  }
224  if (nAp * nAp != kAp.getSize()) {
225  throw LSST_EXCEPT(table::io::MalformedArchiveError, "Reverse X SIP matrix is not square.");
226  }
227  if (nBp * nBp != kBp.getSize()) {
228  throw LSST_EXCEPT(table::io::MalformedArchiveError, "Reverse Y SIP matrix is not square.");
229  }
230  Eigen::Map<Eigen::MatrixXd const> mapA((record)[kA].getData(), nA, nA);
231  Eigen::Map<Eigen::MatrixXd const> mapB((record)[kB].getData(), nB, nB);
232  Eigen::Map<Eigen::MatrixXd const> mapAp((record)[kAp].getData(), nAp, nAp);
233  Eigen::Map<Eigen::MatrixXd const> mapBp((record)[kBp].getData(), nBp, nBp);
234 
235  auto metadata = makeSipMatrixMetadata(mapA, "A");
236  metadata->combine(makeSipMatrixMetadata(mapB, "B"));
237  metadata->combine(makeSipMatrixMetadata(mapAp, "AP"));
238  metadata->combine(makeSipMatrixMetadata(mapBp, "BP"));
239  return metadata;
240 }
241 
242 } // namespace
243 } // namespace geom
244 } // namespace afw
245 } // namespace lsst
table::Key< std::string > name
Definition: Amplifier.cc:116
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition: Exception.h:48
table::Key< std::string > ctype1
Definition: OldWcs.cc:133
table::PointKey< double > crpix
Definition: OldWcs.cc:131
table::Key< std::string > radesys
Definition: OldWcs.cc:136
table::Schema schema
Definition: OldWcs.cc:129
table::Key< table::Array< double > > cd
Definition: OldWcs.cc:132
table::Key< std::string > cunit2
Definition: OldWcs.cc:138
table::Key< std::string > cunit1
Definition: OldWcs.cc:137
table::Key< std::string > ctype2
Definition: OldWcs.cc:134
table::Key< double > equinox
Definition: OldWcs.cc:135
table::PointKey< double > crval
Definition: OldWcs.cc:130
#define LSST_ARCHIVE_ASSERT(EXPR)
An assertion macro used to validate the structure of an InputArchive.
Definition: Persistable.h:48
T equal(T... args)
std::shared_ptr< daf::base::PropertyList > makeSipMatrixMetadata(Eigen::MatrixXd const &matrix, std::string const &name)
Definition: wcsUtils.cc:150
FilterProperty & operator=(FilterProperty const &)=default
A base class for image defects.
STL namespace.
T rbegin(T... args)
T rend(T... args)
T size(T... args)
T sqrt(T... args)
T to_string(T... args)