LSSTApplications  10.0-2-g4f67435,11.0.rc2+1,11.0.rc2+12,11.0.rc2+3,11.0.rc2+4,11.0.rc2+5,11.0.rc2+6,11.0.rc2+7,11.0.rc2+8
LSSTDataManagementBasePackage
makeWcs.cc
Go to the documentation of this file.
1 // -*- lsst-c++ -*-
2 
3 /*
4  * LSST Data Management System
5  * Copyright 2008, 2009, 2010 LSST Corporation.
6  *
7  * This product includes software developed by the
8  * LSST Project (http://www.lsst.org/).
9  *
10  * This program is free software: you can redistribute it and/or modify
11  * it under the terms of the GNU General Public License as published by
12  * the Free Software Foundation, either version 3 of the License, or
13  * (at your option) any later version.
14  *
15  * This program is distributed in the hope that it will be useful,
16  * but WITHOUT ANY WARRANTY; without even the implied warranty of
17  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18  * GNU General Public License for more details.
19  *
20  * You should have received a copy of the LSST License Statement and
21  * the GNU General Public License along with this program. If not,
22  * see <http://www.lsstcorp.org/LegalNotices/>.
23  */
24 
25 #include "boost/make_shared.hpp"
26 #include "Eigen/Core"
27 #include "lsst/pex/logging.h"
28 #include "lsst/afw/image/Wcs.h"
29 #include "lsst/afw/image/TanWcs.h"
30 
31 namespace except = lsst::pex::exceptions;
32 namespace afwImg = lsst::afw::image;
33 
42  bool stripMetadata
43  )
44 {
45  //
46  // _metadata is not const (it is probably meant to be), but we don't want to modify it.
47  //
48  auto metadata = _metadata; // we'll make a copy and modify metadata if needs be
49  auto modifyable = false; // ... and set this variable to say that we did
50 
51  std::string ctype1, ctype2;
52  if (metadata->exists("CTYPE1") && metadata->exists("CTYPE2")) {
53  ctype1 = metadata->getAsString("CTYPE1");
54  ctype2 = metadata->getAsString("CTYPE2");
55  } else {
56  return PTR(Wcs)();
57  }
58  //
59  // SCAMP used to use PVi_j keys with a CTYPE of TAN to specify a "TPV" projection
60  // (cf. https://github.com/astropy/astropy/issues/299
61  // and the discussion from Dave Berry in https://jira.lsstcorp.org/browse/DM-2883)
62  //
63  // Follow Dave's AST and switch TAN to TPV
64  //
65  using pex::logging::Log;
66  auto log = Log(Log::getDefaultLog(), "makeWcs");
67 
68  if (ctype1.substr(5, 3) == "TAN" &&
69  (metadata->exists("PV1_5") || metadata->exists("PV2_1"))) {
70  log.log(Log::INFO, str(boost::format("Interpreting %s/%s + PVi_j as TPV") % ctype1 % ctype2));
71 
72  if (!modifyable) {
73  metadata = _metadata->deepCopy();
74  modifyable = true;
75  }
76 
77  ctype1.replace(5, 3, "TPV");
78  metadata->set<std::string>("CTYPE1", ctype1);
79 
80  ctype2.replace(5, 3, "TPV");
81  metadata->set<std::string>("CTYPE2", ctype2);
82  }
83 
84  afwImg::Wcs::Ptr wcs; // we can't use make_shared as ctor is private
85  if (ctype1.substr(5, 3) == "TAN") {
86  wcs = afwImg::Wcs::Ptr(new afwImg::TanWcs(metadata));
87  } else if (ctype1.substr(5, 3) == "TPV") { // unfortunately we don't support TPV
88  if (!modifyable) {
89  metadata = _metadata->deepCopy();
90  modifyable = true;
91  }
92 
93  log.log(Log::WARN, str(boost::format("Stripping PVi_j keys from projection %s/%s") % ctype1 % ctype2));
94 
95  metadata->set<std::string>("CTYPE1", "RA---TAN");
96  metadata->set<std::string>("CTYPE2", "DEC--TAN");
97  metadata->set<bool>("TPV_WCS", true);
98  // PV1_[1-4] are in principle legal, although Swarp reuses them as part of the TPV parameterisation.
99  // It turns out that leaving PV1_[1-4] in the header breaks wcslib 4.14's ability to read
100  // DECam headers (DM-3196) so we'll delete them all for now.
101  //
102  // John Swinbank points out the maximum value of j in TVi_j is 39;
103  // http://fits.gsfc.nasa.gov/registry/tpvwcs/tpv.html
104  for (int i = 0; i != 2; ++i) {
105  for (int j = 1; j <= 39; ++j) { // 39's the max in the TPV standard
106  char pvName[8];
107  sprintf(pvName, "PV%d_%d", i, j);
108  if (metadata->exists(pvName)) {
109  metadata->remove(pvName);
110  }
111  }
112  }
113 
114  wcs = afwImg::Wcs::Ptr(new afwImg::TanWcs(metadata));
115  } else {
116  wcs = afwImg::Wcs::Ptr(new afwImg::Wcs(metadata));
117  }
118 
119  //If keywords LTV[1,2] are present, the image on disk is already a subimage, so
120  //we should shift the wcs to allow for this.
121  std::string key = "LTV1";
122  if (metadata->exists(key)) {
123  wcs->shiftReferencePixel(-metadata->getAsDouble(key), 0);
124  }
125 
126  key = "LTV2";
127  if (metadata->exists(key) ) {
128  wcs->shiftReferencePixel(0, -metadata->getAsDouble(key));
129  }
130 
131  if (stripMetadata) {
132  afwImg::detail::stripWcsKeywords(metadata, wcs);
133  }
134 
135  return wcs;
136 }
137 
144  double CD11,
145  double CD12,
146  double CD21,
147  double CD22
148  ) {
149  Eigen::Matrix2d CD;
150  CD << CD11, CD12, CD21, CD22;
151  lsst::afw::geom::Point2D crvalTmp;
152  crvalTmp[0] = crval.toIcrs().getLongitude().asDegrees();
153  crvalTmp[1] = crval.toIcrs().getLatitude().asDegrees();
154  return afwImg::Wcs::Ptr(new lsst::afw::image::TanWcs(crvalTmp, crpix, CD));
155 }
#define PTR(...)
Definition: base.h:41
boost::shared_ptr< Wcs > Ptr
Definition: Wcs.h:113
table::PointKey< double > crval
Definition: Wcs.cc:1035
tbl::Key< int > wcs
Implementation of the WCS standard for a any projection.
Definition: Wcs.h:107
table::Key< std::string > ctype2
Definition: Wcs.cc:1039
a place to record messages and descriptions of the state of processing.
Definition: Log.h:154
double asDegrees() const
Definition: Angle.h:124
table::Key< table::Array< Kernel::Pixel > > image
Definition: FixedKernel.cc:117
def log
Definition: log.py:85
Wcs::Ptr makeWcs(boost::shared_ptr< lsst::daf::base::PropertySet > const &fitsMetadata, bool stripMetadata=false)
Definition: makeWcs.cc:40
Implementation of the WCS standard for the special case of the Gnomonic (tangent plane) projection...
Definition: TanWcs.h:68
int INFO
Definition: log.py:37
lsst::afw::geom::Angle getLongitude() const
The main access method for the longitudinal coordinate.
Definition: Coord.h:113
Class for storing generic metadata.
Definition: PropertySet.h:82
virtual IcrsCoord toIcrs() const
Convert ourself to ICRS: RA, Dec (basically J2000)
Definition: Coord.cc:792
virtual Ptr deepCopy(void) const
Definition: PropertySet.cc:70
int WARN
Definition: log.py:38
int stripWcsKeywords(boost::shared_ptr< lsst::daf::base::PropertySet > const &metadata, boost::shared_ptr< Wcs const > const &wcs)
Definition: Wcs.cc:1253
daf::base::PropertySet & _metadata
Definition: fits_io_mpl.h:79
lsst::afw::geom::Angle getLatitude() const
The main access method for the latitudinal coordinate.
Definition: Coord.h:120
table::Key< std::string > ctype1
Definition: Wcs.cc:1038
table::PointKey< double > crpix
Definition: Wcs.cc:1036