# Add Jacobian BoundedField accessor to SkyWcs

XMLWordPrintable

## Details

• Type: Story
• Status: Done
• Resolution: Done
• Fix Version/s: None
• Component/s:
• Labels:
None
• Story Points:
6
• Sprint:
AP F18-4, AP F18-5, AP F18-6
• Team:

## Description

We've agreed that it should not be the responsibility of the PhotoCalib to include a Jacobian derived from the WCS, because we don't always want to apply that.  Instead, SkyWcs should have the ability (given a pixel bounding box) to return a BoundedField that evaluates to the determinant of the Jacobian of the pixelToSky transform (aka the pixel area in angular units) at every point in that bounding box efficiently.  This probably involves a new BoundedField class (though it's possible there's a way to do this with TransformBoundedField).

A reference implementation can be found in meas_mosaic's calculateJacobian functions, though a different approach to vectorization may be necessary to get good performance from AST.

## Attachments

1. wcsJacobian_muddling.py
3 kB

## Activity

Hide
John Parejko added a comment -

We may be able to use skyWcs.getTransform().getJacobian() to accomplish this. Whether it is performant is a different question.

Show
John Parejko added a comment - We may be able to use skyWcs.getTransform().getJacobian() to accomplish this. Whether it is performant is a different question.
Hide
John Parejko added a comment -

Summary of what I've learned so far regarding the determinant of the jacobian of the WCS (hereafter det(J)):

meas_mosaic.calculateJacobian(point) uses the square of the ratio of getPixelScale(point)/getPixelScale(origin). Conceptually, this is the value that we're interested in: how much a pixel is distorted at some pixel location vs. at the "center" of the WCS.

Alternate methods for computing this same value for a SkyWcs object (one would again take the ratio at the point vs. the pixel origin):

 1. np.linalg.det(wcs.getTransform().getJacobian(point))}}   2. linearMap = wcs.getFrameDict().getMapping("pixels", "iwc") t = lsst.afw.geom.TransformPoint2ToPoint2(linearMap) detJ = np.linalg.det(t.getJacobian(point))   3. linear = wcs.linearizePixelToSky(point, lsst.geom.radians) detJ = np.linalg.det(linear.getMatrix()) 

When compared with the pixel scale method, 1. gives results that can differ at the percent-ish level, 2. gives results that are closer, but slowly start to diverge as the point of interest is further from the origin, and 3. seems to stay within floating point precision even for points far from the origin. For example, using a skyWcs loaded from the testdata_jointcal, these are the ratios of the three methods above to the pixel scale calculation:

 In [37]: print(compare(skyWcs, point3, rawJacobianFull))   ...: print(compare(skyWcs, point3, rawJacobianIwc))   ...: print(compare(skyWcs, point3, rawJacobianLinearized))  1.0027808211807387 1.0000102561363005 0.9999999954468474 

I have been trying to come up with a comparable calculation using the gtransfo object in jointcal itself, but have not yet been able to produce any similar numbers. If I were to implement det(J) as an explicit non-fitted term in the photometric calibration, I would have done it similar to the "IWC" calculation above: take the partial derivatives of the AstrometryModel, which goes from pixel->tangent plane (where the tangent plane is the IWC in the resulting SkyWcs object). But this seems to produce a very different det(J).

Show
John Parejko added a comment - Summary of what I've learned so far regarding the determinant of the jacobian of the WCS (hereafter det(J)): meas_mosaic.calculateJacobian(point) uses the square of the ratio of getPixelScale(point)/getPixelScale(origin) . Conceptually, this is the value that we're interested in: how much a pixel is distorted at some pixel location vs. at the "center" of the WCS. Alternate methods for computing this same value for a SkyWcs object (one would again take the ratio at the point vs. the pixel origin): 1. np.linalg.det(wcs.getTransform().getJacobian(point))}}   2. linearMap = wcs.getFrameDict().getMapping("pixels", "iwc") t = lsst.afw.geom.TransformPoint2ToPoint2(linearMap) detJ = np.linalg.det(t.getJacobian(point))   3. linear = wcs.linearizePixelToSky(point, lsst.geom.radians) detJ = np.linalg.det(linear.getMatrix()) When compared with the pixel scale method, 1. gives results that can differ at the percent-ish level, 2. gives results that are closer, but slowly start to diverge as the point of interest is further from the origin, and 3. seems to stay within floating point precision even for points far from the origin. For example, using a skyWcs loaded from the testdata_jointcal, these are the ratios of the three methods above to the pixel scale calculation: In [37]: print(compare(skyWcs, point3, rawJacobianFull)) ...: print(compare(skyWcs, point3, rawJacobianIwc)) ...: print(compare(skyWcs, point3, rawJacobianLinearized)) 1.0027808211807387 1.0000102561363005 0.9999999954468474 I have been trying to come up with a comparable calculation using the gtransfo object in jointcal itself, but have not yet been able to produce any similar numbers. If I were to implement det(J) as an explicit non-fitted term in the photometric calibration, I would have done it similar to the "IWC" calculation above: take the partial derivatives of the AstrometryModel, which goes from pixel->tangent plane (where the tangent plane is the IWC in the resulting SkyWcs object). But this seems to produce a very different det(J).
Hide
Eli Rykoff added a comment -

What is the relative compute time of these four methods?

And maybe I'm missing something in why you can't just use method 3 inside jointcal.  Doesn't this only have to be done once to compute these corrections for each star observation, once you've computed the astrometric solution?

Show
Eli Rykoff added a comment - What is the relative compute time of these four methods?  And maybe I'm missing something in why you can't just use method 3 inside jointcal.  Doesn't this only have to be done once to compute these corrections for each star observation, once you've computed the astrometric solution?
Hide
Russell Owen added a comment -

Note that SkyWcs already offers getPixelScale(pixel) which returns the square root of the pixel area (as an Angle so you can pick your favorite units for the area). Whether or not it is presently efficient enough is a valid question, but if it's not, it can be the basis of an interpolation of some kind.

Show
Russell Owen added a comment - Note that SkyWcs already offers getPixelScale(pixel) which returns the square root of the pixel area (as an Angle so you can pick your favorite units for the area). Whether or not it is presently efficient enough is a valid question, but if it's not, it can be the basis of an interpolation of some kind.
Hide
John Parejko added a comment -

Jim Bosch: do you think I need to implement the persistence parts of this BoundedField? I feel like it will only be used in-memory.

Show
John Parejko added a comment - Jim Bosch : do you think I need to implement the persistence parts of this BoundedField? I feel like it will only be used in-memory.
Hide
Jim Bosch added a comment -

I don't think it needs to be done immediately, and whether we need it in the end depends on:

• whether we decide to make the Jointcal photometric fit include it explicitly as a fixed model component;
• how we want to persist that composite model.

So I'm happy at least deferring until we know the answers to those questions.

Show
Jim Bosch added a comment - I don't think it needs to be done immediately, and whether we need it in the end depends on: whether we decide to make the Jointcal photometric fit include it explicitly as a fixed model component; how we want to persist that composite model. So I'm happy at least deferring until we know the answers to those questions.
Hide
John Parejko added a comment - - edited
Show
John Parejko added a comment - - edited Jenkins run: https://ci.lsst.codes/blue/organizations/jenkins/stack-os-matrix/detail/stack-os-matrix/28985/pipeline
Hide
John Parejko added a comment -

Paul Price: do you mind taking on this medium sized review (~350 new lines)?

The handful of new lines in jointcal were from my initial testing, but I figured they were useful enough to be added to the python layer.

If Jira doesn't pick up the afw PR, it's here: https://github.com/lsst/afw/pull/411

Show
John Parejko added a comment - Paul Price : do you mind taking on this medium sized review (~350 new lines)? The handful of new lines in jointcal were from my initial testing, but I figured they were useful enough to be added to the python layer. If Jira doesn't pick up the afw PR, it's here: https://github.com/lsst/afw/pull/411
Hide
Paul Price added a comment -

Reviews completed on the GitHub PRs. Minor comments in afw that should be addressed before merging. I don't understand how the changes in jointcal relate, since they don't use the new functionality in afw, but you say it's useful so I won't stand in your way; nevertheless, the commit messages should be fixed.

Show
Paul Price added a comment - Reviews completed on the GitHub PRs. Minor comments in afw that should be addressed before merging. I don't understand how the changes in jointcal relate, since they don't use the new functionality in afw, but you say it's useful so I won't stand in your way; nevertheless, the commit messages should be fixed.
Hide
John Parejko added a comment -

Thanks for the review.

Merged and done.

Show
John Parejko added a comment - Thanks for the review. Merged and done.
Hide
John Parejko added a comment -

While reviewing fgcm, I discovered SingleFrameJacobianPlugin in meas_base plugins.py. It uses linearizePixelToSky().determinant (my version 2 above), so it should be very close to the pixel area calculation implemented here. However, it does result in the question of whether we need both of these approaches (the plugin itself doesn't appear to be used anywhere).

Show
John Parejko added a comment - While reviewing fgcm, I discovered SingleFrameJacobianPlugin in meas_base plugins.py . It uses linearizePixelToSky().determinant (my version 2 above), so it should be very close to the pixel area calculation implemented here. However, it does result in the question of whether we need both of these approaches (the plugin itself doesn't appear to be used anywhere).
Hide
Jim Bosch added a comment -

I suspect the plugin as it exists now would be better off using something based on getPixelScale() (and perhaps renaming it to make that clear).   For completely different reasons (computing shapes and centroids in sky coordinates from pixel-coordinate measurements) we're going to also want to start saving all of the matrix elements returned from linearizePixelToSky independently, and I imagine we'll do that in the same plugin.

Show
Jim Bosch added a comment - I suspect the plugin as it exists now would be better off using something based on getPixelScale() (and perhaps renaming it to make that clear).   For completely different reasons (computing shapes and centroids in sky coordinates from pixel-coordinate measurements) we're going to also want to start saving all of the matrix elements returned from linearizePixelToSky independently, and I imagine we'll do that in the same plugin.

## People

• Assignee:
John Parejko
Reporter:
Jim Bosch
Reviewers:
Paul Price
Watchers:
Eli Rykoff, Jim Bosch, John Parejko, John Swinbank, Paul Price, Russell Owen