LSST Applications  21.0.0-172-gfb10e10a+18fedfabac,22.0.0+297cba6710,22.0.0+80564b0ff1,22.0.0+8d77f4f51a,22.0.0+a28f4c53b1,22.0.0+dcf3732eb2,22.0.1-1-g7d6de66+2a20fdde0d,22.0.1-1-g8e32f31+297cba6710,22.0.1-1-geca5380+7fa3b7d9b6,22.0.1-12-g44dc1dc+2a20fdde0d,22.0.1-15-g6a90155+515f58c32b,22.0.1-16-g9282f48+790f5f2caa,22.0.1-2-g92698f7+dcf3732eb2,22.0.1-2-ga9b0f51+7fa3b7d9b6,22.0.1-2-gd1925c9+bf4f0e694f,22.0.1-24-g1ad7a390+a9625a72a8,22.0.1-25-g5bf6245+3ad8ecd50b,22.0.1-25-gb120d7b+8b5510f75f,22.0.1-27-g97737f7+2a20fdde0d,22.0.1-32-gf62ce7b1+aa4237961e,22.0.1-4-g0b3f228+2a20fdde0d,22.0.1-4-g243d05b+871c1b8305,22.0.1-4-g3a563be+32dcf1063f,22.0.1-4-g44f2e3d+9e4ab0f4fa,22.0.1-42-gca6935d93+ba5e5ca3eb,22.0.1-5-g15c806e+85460ae5f3,22.0.1-5-g58711c4+611d128589,22.0.1-5-g75bb458+99c117b92f,22.0.1-6-g1c63a23+7fa3b7d9b6,22.0.1-6-g50866e6+84ff5a128b,22.0.1-6-g8d3140d+720564cf76,22.0.1-6-gd805d02+cc5644f571,22.0.1-8-ge5750ce+85460ae5f3,master-g6e05de7fdc+babf819c66,master-g99da0e417a+8d77f4f51a,w.2021.48
LSST Data Management Base Package
wcsUtils.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 <iostream>
24 #include <cmath>
25 
26 #include "lsst/daf/base.h"
27 #include "lsst/pex/exceptions.h"
29 #include "lsst/afw/geom/wcsUtils.h"
30 
31 namespace lsst {
32 namespace afw {
33 namespace geom {
34 namespace {
35 /*
36  * Get the name of a SIP matrix coefficient card
37  *
38  * This works for matrix coefficients of the form: <name>_i_j where i and j are zero-based.
39  * Thus it works for SIP matrix terms but not CD matrix terms (because CD terms use 1-based indexing
40  * and have no underscore between the name and the first index).
41  */
42 std::string getSipCoeffCardName(std::string const& name, int i, int j) {
43  return name + std::to_string(i) + "_" + std::to_string(j);
44 }
45 
46 } // namespace
47 
49  lsst::geom::Point2I const& xy0) {
51 
52  wcsMetaData->set("CTYPE1" + wcsName, "LINEAR", "Type of projection");
53  wcsMetaData->set("CTYPE2" + wcsName, "LINEAR", "Type of projection");
54  wcsMetaData->set("CRPIX1" + wcsName, static_cast<double>(1), "Column Pixel Coordinate of Reference");
55  wcsMetaData->set("CRPIX2" + wcsName, static_cast<double>(1), "Row Pixel Coordinate of Reference");
56  wcsMetaData->set("CRVAL1" + wcsName, static_cast<double>(xy0[0]), "Column pixel of Reference Pixel");
57  wcsMetaData->set("CRVAL2" + wcsName, static_cast<double>(xy0[1]), "Row pixel of Reference Pixel");
58  wcsMetaData->set("CUNIT1" + wcsName, "PIXEL", "Column unit");
59  wcsMetaData->set("CUNIT2" + wcsName, "PIXEL", "Row unit");
60 
61  return wcsMetaData;
62 }
63 
65  std::vector<std::string> const names = {"CRPIX1", "CRPIX2", "CRVAL1", "CRVAL2", "CTYPE1",
66  "CTYPE2", "CUNIT1", "CUNIT2", "CD1_1", "CD1_2",
67  "CD2_1", "CD2_2", "WCSAXES"};
68  for (auto const& name : names) {
69  if (metadata.exists(name + wcsName)) {
70  metadata.remove(name + wcsName);
71  }
72  }
73 }
74 
75 Eigen::Matrix2d getCdMatrixFromMetadata(daf::base::PropertySet& metadata) {
76  Eigen::Matrix2d matrix;
77  bool found{false};
78  for (int i = 0; i < 2; ++i) {
79  for (int j = 0; j < 2; ++j) {
80  std::string const cardName = "CD" + std::to_string(i + 1) + "_" + std::to_string(j + 1);
81  if (metadata.exists(cardName)) {
82  matrix(i, j) = metadata.getAsDouble(cardName);
83  found = true;
84  } else {
85  matrix(i, j) = 0.0;
86  }
87  }
88  }
89  if (!found) {
90  throw LSST_EXCEPT(pex::exceptions::TypeError, "No CD matrix coefficients found");
91  }
92  return matrix;
93 }
94 
96  bool strip) {
97  int x0 = 0; // Our value of X0
98  int y0 = 0; // Our value of Y0
99 
100  if (metadata.exists("CRPIX1" + wcsName) && metadata.exists("CRPIX2" + wcsName) &&
101  metadata.exists("CRVAL1" + wcsName) && metadata.exists("CRVAL2" + wcsName) &&
102  (metadata.getAsDouble("CRPIX1" + wcsName) == 1.0) &&
103  (metadata.getAsDouble("CRPIX2" + wcsName) == 1.0)) {
104  x0 = static_cast<int>(std::lround(metadata.getAsDouble("CRVAL1" + wcsName)));
105  y0 = static_cast<int>(std::lround(metadata.getAsDouble("CRVAL2" + wcsName)));
106  if (strip) {
107  deleteBasicWcsMetadata(metadata, wcsName);
108  }
109  }
110  return lsst::geom::Point2I(x0, y0);
111 }
112 
113 Eigen::MatrixXd getSipMatrixFromMetadata(daf::base::PropertySet const& metadata, std::string const& name) {
114  std::string cardName = name + "_ORDER";
115  if (!metadata.exists(cardName)) {
117  "metadata does not contain SIP matrix " + name + ": " + cardName + " not found");
118  };
119  int order = metadata.getAsInt(cardName);
120  if (order < 0) {
122  "matrix order = " + std::to_string(order) + " < 0");
123  }
124  Eigen::MatrixXd matrix(order + 1, order + 1);
125  auto const coeffName = name + "_";
126  for (int i = 0; i <= order; ++i) {
127  for (int j = 0; j <= order; ++j) {
128  std::string const cardName = getSipCoeffCardName(coeffName, i, j);
129  if (metadata.exists(cardName)) {
130  matrix(i, j) = metadata.getAsDouble(cardName);
131  } else {
132  matrix(i, j) = 0.0;
133  }
134  }
135  }
136  return matrix;
137 }
138 
139 bool hasSipMatrix(daf::base::PropertySet const& metadata, std::string const& name) {
140  std::string cardName = name + "_ORDER";
141  if (!metadata.exists(cardName)) {
142  return false;
143  }
144  return metadata.getAsInt(cardName) >= 0;
145 }
146 
148  std::string const& name) {
149  if (matrix.rows() != matrix.cols() || matrix.rows() < 1) {
151  "Matrix must be square and at least 1 x 1; dimensions = " +
152  std::to_string(matrix.rows()) + " x " + std::to_string(matrix.cols()));
153  }
154  int const order = matrix.rows() - 1;
155  auto metadata = std::make_shared<daf::base::PropertyList>();
156  std::string cardName = name + "_ORDER";
157  metadata->set(cardName, order);
158  auto const coeffName = name + "_";
159  for (int i = 0; i <= order; ++i) {
160  for (int j = 0; j <= order; ++j) {
161  if (matrix(i, j) != 0.0) {
162  metadata->set(getSipCoeffCardName(coeffName, i, j), matrix(i, j));
163  }
164  }
165  }
166  return metadata;
167 }
168 
171  Eigen::Matrix2d const& cdMatrix,
172  std::string const& projection) {
173  auto pl = std::make_shared<daf::base::PropertyList>();
174  pl->add("RADESYS", "ICRS");
175  pl->add("CTYPE1", "RA---" + projection);
176  pl->add("CTYPE2", "DEC--" + projection);
177  pl->add("CRPIX1", crpix[0] + 1);
178  pl->add("CRPIX2", crpix[1] + 1);
179  pl->add("CRVAL1", crval[0].asDegrees());
180  pl->add("CRVAL2", crval[1].asDegrees());
181  pl->add("CUNIT1", "deg");
182  pl->add("CUNIT2", "deg");
183  for (int i = 0; i < 2; ++i) {
184  for (int j = 0; j < 2; ++j) {
185  if (cdMatrix(i, j) != 0.0) {
186  std::string name = "CD" + std::to_string(i + 1) + "_" + std::to_string(j + 1);
187  pl->add(name, cdMatrix(i, j));
188  }
189  }
190  }
191  return pl;
192 }
193 
196  Eigen::Matrix2d const& cdMatrix,
197  Eigen::MatrixXd const& sipA,
198  Eigen::MatrixXd const& sipB) {
199  auto metadata = makeSimpleWcsMetadata(crpix, crval, cdMatrix, "TAN-SIP");
200  metadata->combine(makeSipMatrixMetadata(sipA, "A"));
201  metadata->combine(makeSipMatrixMetadata(sipB, "B"));
202  return metadata;
203 }
204 
207  Eigen::Matrix2d const& cdMatrix, Eigen::MatrixXd const& sipA, Eigen::MatrixXd const& sipB,
208  Eigen::MatrixXd const& sipAp, Eigen::MatrixXd const& sipBp) {
209  auto metadata = makeTanSipMetadata(crpix, crval, cdMatrix, sipA, sipB);
210  metadata->combine(makeSipMatrixMetadata(sipAp, "AP"));
211  metadata->combine(makeSipMatrixMetadata(sipBp, "BP"));
212  return metadata;
213 }
214 
215 } // namespace geom
216 } // namespace afw
217 } // 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::PointKey< double > crpix
Definition: OldWcs.cc:129
table::PointKey< double > crval
Definition: OldWcs.cc:128
Class for storing ordered metadata with comments.
Definition: PropertyList.h:68
Class for storing generic metadata.
Definition: PropertySet.h:66
int getAsInt(std::string const &name) const
Get the last value for a bool/char/short/int 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.
double getAsDouble(std::string const &name) const
Get the last value for any arithmetic property name (possibly hierarchical).
Point in an unspecified spherical coordinate system.
Definition: SpherePoint.h:57
Reports errors from accepting an object of an unexpected or inappropriate type.
Definition: Runtime.h:167
bool strip
Definition: fits.cc:911
std::shared_ptr< daf::base::PropertyList > makeSimpleWcsMetadata(lsst::geom::Point2D const &crpix, lsst::geom::SpherePoint const &crval, Eigen::Matrix2d const &cdMatrix, std::string const &projection="TAN")
Make FITS metadata for a simple FITS WCS (one with no distortion).
Definition: wcsUtils.cc:169
Eigen::MatrixXd getSipMatrixFromMetadata(daf::base::PropertySet const &metadata, std::string const &name)
Definition: wcsUtils.cc:113
bool hasSipMatrix(daf::base::PropertySet const &metadata, std::string const &name)
Definition: wcsUtils.cc:139
std::shared_ptr< daf::base::PropertyList > makeTanSipMetadata(lsst::geom::Point2D const &crpix, lsst::geom::SpherePoint const &crval, Eigen::Matrix2d const &cdMatrix, Eigen::MatrixXd const &sipA, Eigen::MatrixXd const &sipB)
Make metadata for a TAN-SIP WCS without inverse matrices.
Definition: wcsUtils.cc:194
Eigen::Matrix2d getCdMatrixFromMetadata(daf::base::PropertySet &metadata)
Read a CD matrix from FITS WCS metadata.
Definition: wcsUtils.cc:75
std::shared_ptr< daf::base::PropertyList > createTrivialWcsMetadata(std::string const &wcsName, lsst::geom::Point2I const &xy0)
Definition: wcsUtils.cc:48
std::shared_ptr< daf::base::PropertyList > makeSipMatrixMetadata(Eigen::MatrixXd const &matrix, std::string const &name)
Definition: wcsUtils.cc:147
lsst::geom::Point2I getImageXY0FromMetadata(daf::base::PropertySet &metadata, std::string const &wcsName, bool strip=false)
Definition: wcsUtils.cc:95
void deleteBasicWcsMetadata(daf::base::PropertySet &metadata, std::string const &wcsName)
Definition: wcsUtils.cc:64
Point< int, 2 > Point2I
Definition: Point.h:321
A base class for image defects.
T lround(T... args)
T to_string(T... args)
table::Key< int > order