LSSTApplications  19.0.0-14-gb0260a2+72efe9b372,20.0.0+7927753e06,20.0.0+8829bf0056,20.0.0+995114c5d2,20.0.0+b6f4b2abd1,20.0.0+bddc4f4cbe,20.0.0-1-g253301a+8829bf0056,20.0.0-1-g2b7511a+0d71a2d77f,20.0.0-1-g5b95a8c+7461dd0434,20.0.0-12-g321c96ea+23efe4bbff,20.0.0-16-gfab17e72e+fdf35455f6,20.0.0-2-g0070d88+ba3ffc8f0b,20.0.0-2-g4dae9ad+ee58a624b3,20.0.0-2-g61b8584+5d3db074ba,20.0.0-2-gb780d76+d529cf1a41,20.0.0-2-ged6426c+226a441f5f,20.0.0-2-gf072044+8829bf0056,20.0.0-2-gf1f7952+ee58a624b3,20.0.0-20-geae50cf+e37fec0aee,20.0.0-25-g3dcad98+544a109665,20.0.0-25-g5eafb0f+ee58a624b3,20.0.0-27-g64178ef+f1f297b00a,20.0.0-3-g4cc78c6+e0676b0dc8,20.0.0-3-g8f21e14+4fd2c12c9a,20.0.0-3-gbd60e8c+187b78b4b8,20.0.0-3-gbecbe05+48431fa087,20.0.0-38-ge4adf513+a12e1f8e37,20.0.0-4-g97dc21a+544a109665,20.0.0-4-gb4befbc+087873070b,20.0.0-4-gf910f65+5d3db074ba,20.0.0-5-gdfe0fee+199202a608,20.0.0-5-gfbfe500+d529cf1a41,20.0.0-6-g64f541c+d529cf1a41,20.0.0-6-g9a5b7a1+a1cd37312e,20.0.0-68-ga3f3dda+5fca18c6a4,20.0.0-9-g4aef684+e18322736b,w.2020.45
LSSTDataManagementBasePackage
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 <sstream>
25 #include <cmath>
26 #include <cstring>
27 #include <limits>
28 
29 #include "lsst/daf/base.h"
30 #include "lsst/pex/exceptions.h"
32 #include "lsst/afw/geom/wcsUtils.h"
33 
34 namespace lsst {
35 namespace afw {
36 namespace geom {
37 namespace {
38 /*
39  * Get the name of a SIP matrix coefficient card
40  *
41  * This works for matrix coefficients of the form: <name>_i_j where i and j are zero-based.
42  * Thus it works for SIP matrix terms but not CD matrix terms (because CD terms use 1-based indexing
43  * and have no underscore between the name and the first index).
44  */
45 std::string getSipCoeffCardName(std::string const& name, int i, int j) {
46  return name + std::to_string(i) + "_" + std::to_string(j);
47 }
48 
49 } // namespace
50 
52  lsst::geom::Point2I const& xy0) {
54 
55  wcsMetaData->set("CTYPE1" + wcsName, "LINEAR", "Type of projection");
56  wcsMetaData->set("CTYPE2" + wcsName, "LINEAR", "Type of projection");
57  wcsMetaData->set("CRPIX1" + wcsName, static_cast<double>(1), "Column Pixel Coordinate of Reference");
58  wcsMetaData->set("CRPIX2" + wcsName, static_cast<double>(1), "Row Pixel Coordinate of Reference");
59  wcsMetaData->set("CRVAL1" + wcsName, static_cast<double>(xy0[0]), "Column pixel of Reference Pixel");
60  wcsMetaData->set("CRVAL2" + wcsName, static_cast<double>(xy0[1]), "Row pixel of Reference Pixel");
61  wcsMetaData->set("CUNIT1" + wcsName, "PIXEL", "Column unit");
62  wcsMetaData->set("CUNIT2" + wcsName, "PIXEL", "Row unit");
63 
64  return wcsMetaData;
65 }
66 
68  std::vector<std::string> const names = {"CRPIX1", "CRPIX2", "CRVAL1", "CRVAL2", "CTYPE1",
69  "CTYPE2", "CUNIT1", "CUNIT2", "CD1_1", "CD1_2",
70  "CD2_1", "CD2_2", "WCSAXES"};
71  for (auto const& name : names) {
72  if (metadata.exists(name + wcsName)) {
73  metadata.remove(name + wcsName);
74  }
75  }
76 }
77 
78 Eigen::Matrix2d getCdMatrixFromMetadata(daf::base::PropertySet& metadata) {
79  Eigen::Matrix2d matrix;
80  bool found{false};
81  for (int i = 0; i < 2; ++i) {
82  for (int j = 0; j < 2; ++j) {
83  std::string const cardName = "CD" + std::to_string(i + 1) + "_" + std::to_string(j + 1);
84  if (metadata.exists(cardName)) {
85  matrix(i, j) = metadata.getAsDouble(cardName);
86  found = true;
87  } else {
88  matrix(i, j) = 0.0;
89  }
90  }
91  }
92  if (!found) {
93  throw LSST_EXCEPT(pex::exceptions::TypeError, "No CD matrix coefficients found");
94  }
95  return matrix;
96 }
97 
99  bool strip) {
100  int x0 = 0; // Our value of X0
101  int y0 = 0; // Our value of Y0
102 
103  if (metadata.exists("CRPIX1" + wcsName) && metadata.exists("CRPIX2" + wcsName) &&
104  metadata.exists("CRVAL1" + wcsName) && metadata.exists("CRVAL2" + wcsName) &&
105  (metadata.getAsDouble("CRPIX1" + wcsName) == 1.0) &&
106  (metadata.getAsDouble("CRPIX2" + wcsName) == 1.0)) {
107  x0 = static_cast<int>(std::lround(metadata.getAsDouble("CRVAL1" + wcsName)));
108  y0 = static_cast<int>(std::lround(metadata.getAsDouble("CRVAL2" + wcsName)));
109  if (strip) {
110  deleteBasicWcsMetadata(metadata, wcsName);
111  }
112  }
113  return lsst::geom::Point2I(x0, y0);
114 }
115 
116 Eigen::MatrixXd getSipMatrixFromMetadata(daf::base::PropertySet const& metadata, std::string const& name) {
117  std::string cardName = name + "_ORDER";
118  if (!metadata.exists(cardName)) {
120  "metadata does not contain SIP matrix " + name + ": " + cardName + " not found");
121  };
122  int order = metadata.getAsInt(cardName);
123  if (order < 0) {
125  "matrix order = " + std::to_string(order) + " < 0");
126  }
127  Eigen::MatrixXd matrix(order + 1, order + 1);
128  auto const coeffName = name + "_";
129  for (int i = 0; i <= order; ++i) {
130  for (int j = 0; j <= order; ++j) {
131  std::string const cardName = getSipCoeffCardName(coeffName, i, j);
132  if (metadata.exists(cardName)) {
133  matrix(i, j) = metadata.getAsDouble(cardName);
134  } else {
135  matrix(i, j) = 0.0;
136  }
137  }
138  }
139  return matrix;
140 }
141 
142 bool hasSipMatrix(daf::base::PropertySet const& metadata, std::string const& name) {
143  std::string cardName = name + "_ORDER";
144  if (!metadata.exists(cardName)) {
145  return false;
146  }
147  return metadata.getAsInt(cardName) >= 0;
148 }
149 
151  std::string const& name) {
152  if (matrix.rows() != matrix.cols() || matrix.rows() < 1) {
154  "Matrix must be square and at least 1 x 1; dimensions = " +
155  std::to_string(matrix.rows()) + " x " + std::to_string(matrix.cols()));
156  }
157  int const order = matrix.rows() - 1;
158  auto metadata = std::make_shared<daf::base::PropertyList>();
159  std::string cardName = name + "_ORDER";
160  metadata->set(cardName, order);
161  auto const coeffName = name + "_";
162  for (int i = 0; i <= order; ++i) {
163  for (int j = 0; j <= order; ++j) {
164  if (matrix(i, j) != 0.0) {
165  metadata->set(getSipCoeffCardName(coeffName, i, j), matrix(i, j));
166  }
167  }
168  }
169  return metadata;
170 }
171 
174  Eigen::Matrix2d const& cdMatrix,
175  std::string const& projection) {
176  auto pl = std::make_shared<daf::base::PropertyList>();
177  pl->add("RADESYS", "ICRS");
178  pl->add("CTYPE1", "RA---" + projection);
179  pl->add("CTYPE2", "DEC--" + projection);
180  pl->add("CRPIX1", crpix[0] + 1);
181  pl->add("CRPIX2", crpix[1] + 1);
182  pl->add("CRVAL1", crval[0].asDegrees());
183  pl->add("CRVAL2", crval[1].asDegrees());
184  pl->add("CUNIT1", "deg");
185  pl->add("CUNIT2", "deg");
186  for (int i = 0; i < 2; ++i) {
187  for (int j = 0; j < 2; ++j) {
188  if (cdMatrix(i, j) != 0.0) {
189  std::string name = "CD" + std::to_string(i + 1) + "_" + std::to_string(j + 1);
190  pl->add(name, cdMatrix(i, j));
191  }
192  }
193  }
194  return pl;
195 }
196 
199  Eigen::Matrix2d const& cdMatrix,
200  Eigen::MatrixXd const& sipA,
201  Eigen::MatrixXd const& sipB) {
202  auto metadata = makeSimpleWcsMetadata(crpix, crval, cdMatrix, "TAN-SIP");
203  metadata->combine(makeSipMatrixMetadata(sipA, "A"));
204  metadata->combine(makeSipMatrixMetadata(sipB, "B"));
205  return metadata;
206 }
207 
210  Eigen::Matrix2d const& cdMatrix, Eigen::MatrixXd const& sipA, Eigen::MatrixXd const& sipB,
211  Eigen::MatrixXd const& sipAp, Eigen::MatrixXd const& sipBp) {
212  auto metadata = makeTanSipMetadata(crpix, crval, cdMatrix, sipA, sipB);
213  metadata->combine(makeSipMatrixMetadata(sipAp, "AP"));
214  metadata->combine(makeSipMatrixMetadata(sipBp, "BP"));
215  return metadata;
216 }
217 
218 } // namespace geom
219 } // namespace afw
220 } // namespace lsst
lsst::afw::geom::deleteBasicWcsMetadata
void deleteBasicWcsMetadata(daf::base::PropertySet &metadata, std::string const &wcsName)
Definition: wcsUtils.cc:67
lsst::daf::base::PropertySet::getAsInt
int getAsInt(std::string const &name) const
Get the last value for a bool/char/short/int property name (possibly hierarchical).
lsst::afw::geom::getSipMatrixFromMetadata
Eigen::MatrixXd getSipMatrixFromMetadata(daf::base::PropertySet const &metadata, std::string const &name)
Definition: wcsUtils.cc:116
std::string
STL class.
std::shared_ptr
STL class.
lsst::afw::geom::hasSipMatrix
bool hasSipMatrix(daf::base::PropertySet const &metadata, std::string const &name)
Definition: wcsUtils.cc:142
std::vector< std::string >
lsst::afw
Definition: imageAlgorithm.dox:1
lsst::daf::base::PropertyList
Class for storing ordered metadata with comments.
Definition: PropertyList.h:68
lsst::daf::base::PropertySet::exists
bool exists(std::string const &name) const
Determine if a name (possibly hierarchical) exists.
strip
bool strip
Definition: fits.cc:911
wcsUtils.h
lsst::afw::geom.transform.transformContinued.name
string name
Definition: transformContinued.py:32
crval
table::PointKey< double > crval
Definition: OldWcs.cc:130
base.h
lsst::afw::geom::createTrivialWcsMetadata
std::shared_ptr< daf::base::PropertyList > createTrivialWcsMetadata(std::string const &wcsName, lsst::geom::Point2I const &xy0)
Definition: wcsUtils.cc:51
crpix
table::PointKey< double > crpix
Definition: OldWcs.cc:131
std::to_string
T to_string(T... args)
lsst
A base class for image defects.
Definition: imageAlgorithm.dox:1
lsst::daf::base::PropertySet::remove
virtual void remove(std::string const &name)
Remove all values for a property name (possibly hierarchical).
LSST_EXCEPT
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition: Exception.h:48
lsst::daf::base::PropertySet::getAsDouble
double getAsDouble(std::string const &name) const
Get the last value for any arithmetic property name (possibly hierarchical).
lsst::geom
Definition: AffineTransform.h:36
std::lround
T lround(T... args)
lsst::geom::Point< int, 2 >
lsst::daf::base::PropertySet
Class for storing generic metadata.
Definition: PropertySet.h:67
lsst::afw::geom::makeTanSipMetadata
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:197
lsst::geom::Point2I
Point< int, 2 > Point2I
Definition: Point.h:321
lsst::afw::geom::getImageXY0FromMetadata
lsst::geom::Point2I getImageXY0FromMetadata(daf::base::PropertySet &metadata, std::string const &wcsName, bool strip=false)
Definition: wcsUtils.cc:98
lsst::afw::geom::makeSimpleWcsMetadata
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:172
lsst.pex::exceptions::TypeError
Reports errors from accepting an object of an unexpected or inappropriate type.
Definition: Runtime.h:167
lsst::geom::SpherePoint
Point in an unspecified spherical coordinate system.
Definition: SpherePoint.h:57
lsst::afw::geom::getCdMatrixFromMetadata
Eigen::Matrix2d getCdMatrixFromMetadata(daf::base::PropertySet &metadata)
Read a CD matrix from FITS WCS metadata.
Definition: wcsUtils.cc:78
frameSetUtils.h
exceptions.h
lsst::afw::geom::makeSipMatrixMetadata
std::shared_ptr< daf::base::PropertyList > makeSipMatrixMetadata(Eigen::MatrixXd const &matrix, std::string const &name)
Definition: wcsUtils.cc:150