Skip to content

Commit

Permalink
Add support for making a SkyWcs from cameraGeom
Browse files Browse the repository at this point in the history 
Add a new overload of `makeSkyWcs` that accepts a
PIXELS to FIELD_ANGLE transform, orientation, flipX and crval.
Include two tests: one with no distortion and one with distortion.
  • Loading branch information
@r-owen
r-owen committed Jul 10, 2018
1 parent a70beb8 commit fb6f4ce
Show file tree
Hide file tree
Showing 4 changed files with 95 additions and 0 deletions.
20 changes: 20 additions & 0 deletions include/lsst/afw/geom/SkyWcs.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -521,6 +521,26 @@ std::shared_ptr<SkyWcs> makeSkyWcs(daf::base::PropertySet &metadata, bool strip
std::shared_ptr<SkyWcs> makeSkyWcs(lsst::geom::Point2D const &crpix, lsst::geom::SpherePoint const &crval,
Eigen::Matrix2d const &cdMatrix, std::string const &projection = "TAN");

/**
* Construct a FITS SkyWcs from camera geometry
*
* @param[in] pixelsToFieldAngle Transformation from \ref afwCameraGeomPIXELS "pixels"
* to \ref afwCameraGeomFIELD_ANGLE "field angle" (in radians).
* @param[in] orientation Position angle of focal plane +Y, measured from N through E at crval.
* At 0 degrees, +Y is along N and +X is along W/E if flipX false/true.
* At 90 degrees, +Y is along E and +X is along N/S if flipX false/true.
* @param[in] flipX Fip x axis? See orientation for details.
* @param[in] crval Center of projection on the sky at field angle (0, 0).
* @param[in] projection The name of the FITS WCS projection, e.g. "TAN" or "STG".
*
* @note Unlike makeCdMatrix, `orientation` is with respect to the focal plane axes, not the CCD axes.
* This is because field angle is defined with respect to focal plane axes.
*/
std::shared_ptr<SkyWcs> makeSkyWcs(TransformPoint2ToPoint2 const &pixelsToFieldAngle,
lsst::geom::Angle const &orientation, bool flipX,
lsst::geom::SpherePoint const &crval,
std::string const &projection = "TAN");

/**
* Construct a TAN-SIP SkyWcs with forward SIP distortion terms and an iterative inverse.
*
Expand Down
4 changes: 4 additions & 0 deletions python/lsst/afw/geom/skyWcs/skyWcs.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -60,6 +60,10 @@ PYBIND11_PLUGIN(skyWcs) {
"crpix"_a, "crval"_a, "cdMatrix"_a, "projection"_a = "TAN");
mod.def("makeSkyWcs", (std::shared_ptr<SkyWcs>(*)(daf::base::PropertySet &, bool))makeSkyWcs,
"metadata"_a, "strip"_a = false);
mod.def("makeSkyWcs",
(std::shared_ptr<SkyWcs>(*)(TransformPoint2ToPoint2 const &, lsst::geom::Angle const &, bool,
lsst::geom::SpherePoint const &, std::string const &))makeSkyWcs,
"pixelsToFieldAngle"_a, "orientation"_a, "flipX"_a, "crval"_a, "projection"_a = "TAN");
mod.def("makeTanSipWcs",
(std::shared_ptr<SkyWcs>(*)(lsst::geom::Point2D const &, lsst::geom::SpherePoint const &,
Eigen::Matrix2d const &, Eigen::MatrixXd const &,
Expand Down
30 changes: 30 additions & 0 deletions src/geom/SkyWcs.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -106,6 +106,29 @@ std::string getSkyWcsPersistenceName() { return "SkyWcs"; }

SkyWcsFactory registration(getSkyWcsPersistenceName());

ast::FrameDict makeSkyWcsFrameDict(TransformPoint2ToPoint2 const& pixelsToFieldAngle,
lsst::geom::Angle const& orientation, bool flipX,
lsst::geom::SpherePoint const& crval,
std::string const& projection = "TAN") {
auto const orientationAndFlipXMatrix = makeCdMatrix(1 * lsst::geom::degrees, orientation, flipX);
auto const initialWcs = makeSkyWcs(lsst::geom::Point2D(0, 0), crval, orientationAndFlipXMatrix, projection);
auto const initialFrameDict = initialWcs->getFrameDict();
auto const iwcToSkyMap = initialFrameDict->getMapping("IWC", "SKY");
auto const pixelFrame = initialFrameDict->getFrame("PIXELS");
auto const iwcFrame = initialFrameDict->getFrame("IWC");
auto const skyFrame = initialFrameDict->getFrame("SKY");
// Field angle is in radians and is aligned to focal plane x and y;
// IWC is basically the same thing, but in degrees and with rotation and flipX applied
ndarray::Array<double, 2, 2> fieldAngleToIwcNdArray = ndarray::allocate(2, 2);
asEigenMatrix(fieldAngleToIwcNdArray) = orientationAndFlipXMatrix * 180.0 / lsst::geom::PI;
auto const pixelsToFieldAngleMap = pixelsToFieldAngle.getMapping();
auto const fieldAngleToIwcMap = ast::MatrixMap(fieldAngleToIwcNdArray);
auto const pixelsToIwcMap = pixelsToFieldAngleMap->then(fieldAngleToIwcMap);
auto finalFrameDict = ast::FrameDict(*pixelFrame, pixelsToIwcMap, *iwcFrame);
finalFrameDict.addFrame("IWC", *iwcToSkyMap, *skyFrame);
return finalFrameDict;
}

} // namespace

Eigen::Matrix2d makeCdMatrix(lsst::geom::Angle const& scale, lsst::geom::Angle const& orientation,
Expand Down Expand Up @@ -455,6 +478,13 @@ std::shared_ptr<SkyWcs> makeSkyWcs(lsst::geom::Point2D const& crpix, lsst::geom:
return std::make_shared<SkyWcs>(*metadata);
}

std::shared_ptr<SkyWcs> makeSkyWcs(TransformPoint2ToPoint2 const& pixelsToFieldAngle,
lsst::geom::Angle const& orientation, bool flipX,
lsst::geom::SpherePoint const& crval, std::string const& projection) {
auto frameDict = makeSkyWcsFrameDict(pixelsToFieldAngle, orientation, flipX, crval, projection);
return std::make_shared<SkyWcs>(frameDict);
}

std::shared_ptr<SkyWcs> makeTanSipWcs(lsst::geom::Point2D const& crpix, lsst::geom::SpherePoint const& crval,
Eigen::Matrix2d const& cdMatrix, Eigen::MatrixXd const& sipA,
Eigen::MatrixXd const& sipB) {
Expand Down
41 changes: 41 additions & 0 deletions tests/test_skyWcs.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -482,6 +482,47 @@ def testMakeModifiedWcsWithActualPixels(self):

self.checkNonFitsWcs(modifiedWcs)

def testMakeSkyWcsFromPixelsToFieldAngleNoDistortion(self):
"""Test makeSkyWcs from pixels to field angle with no distortion
"""
for crval, orientation, flipX, projection in itertools.product(self.crvalList, self.orientationList,
(False, True), ("TAN", "STG")):
with self.subTest(crval=crval, orientation=orientation):
cdMatrix = makeCdMatrix(scale=self.scale, orientation=orientation, flipX=flipX)
wcs1 = makeSkyWcs(crpix=self.crpix, crval=crval, cdMatrix=cdMatrix, projection=projection)
scaleRad = self.scale.asRadians()
shiftRad = lsst.geom.Extent2D(self.crpix)*-scaleRad
pixelsToFieldAngle = lsst.afw.geom.makeTransform(
lsst.geom.AffineTransform(lsst.geom.LinearTransform.makeScaling(scaleRad), shiftRad))
wcs2 = makeSkyWcs(pixelsToFieldAngle=pixelsToFieldAngle,
orientation=orientation,
flipX=flipX,
crval=crval,
projection=projection)
self.assertWcsAlmostEqualOverBBox(wcs1, wcs2, self.bbox)

def testMakeSkyWcsFromPixelsToFieldAngleWithDistortion(self):
"""Test makeSkyWcs from pixels to field angle with distortion
"""
# make an arbitrary but reasonable distortion transform
pixelTransform = makeRadialTransform([0.0, 1.0, 0.0, 0.0011])

for crval, orientation, flipX, projection in itertools.product(self.crvalList, self.orientationList,
(False, True), ("TAN", "STG")):
with self.subTest(crval=crval, orientation=orientation):
cdMatrix = makeCdMatrix(scale=self.scale, orientation=orientation, flipX=flipX)
wcs0 = makeSkyWcs(crpix=self.crpix, crval=crval, cdMatrix=cdMatrix, projection=projection)
# modifyActualPixels is ignored because there is no ACTUAL_PIXELS frame
wcs1 = makeModifiedWcs(pixelTransform=pixelTransform, wcs=wcs0, modifyActualPixels=False)
scaleRad = self.scale.asRadians()
shiftRad = lsst.geom.Extent2D(self.crpix)*-scaleRad
pixelsToFieldAngle = lsst.afw.geom.makeTransform(
lsst.geom.AffineTransform(lsst.geom.LinearTransform.makeScaling(scaleRad), shiftRad))
pixelsToFieldAngle = pixelTransform.then(pixelsToFieldAngle)
wcs2 = makeSkyWcs(pixelsToFieldAngle=pixelsToFieldAngle, orientation=orientation,
flipX=flipX, crval=crval, projection=projection)
self.assertWcsAlmostEqualOverBBox(wcs1, wcs2, self.bbox)

@unittest.skipIf(sys.version_info[0] < 3, "astropy.wcs rejects the header on py2")
def testAgainstAstropyWcs(self):
bbox = lsst.geom.Box2D(lsst.geom.Point2D(-1000, -1000), lsst.geom.Extent2D(2000, 2000))
Expand Down

0 comments on commit fb6f4ce

Please sign in to comment.