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

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).

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?

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.

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.

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.

Jenkins run: https://ci.lsst.codes/blue/organizations/jenkins/stack-os-matrix/detail/stack-os-matrix/28985/pipeline

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

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.

Thanks for the review.

Merged and done.

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).

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.