Add support for specifying the axis when obtaining a wcs pixel scale - #3230

I'm fine with merging this, but I'd like to have feedback from @mcara and/or @nden first.

I fixed the docstrings - although I couldn't get the double backtick to be parsed as literal code in the Raises section.

I like that both scales are now available since this can be used now in legacy code we have. I left a minor comment about units. @mcara may have other comments but 👍 from me.

I thought about this issue, i.e., what do we call "pixel scale", and I think it is not simple and that the answer depends on how this pixel scale is intended to be used. Now, I have a really bad habit of being verbose, so I want to apologize for this ahead of time. Also, I will split my comments into thee parts: 1) I would like to discuss how is "pixel scale" being used; 2) I will discuss issues related to the current function as it is (assuming you insist on keeping this function); and 3) Propose two separate functions instead of the current function.

Edit: went from 2 parts -> 3 parts (due to verbosity) 😃

With regard to the use of the "pixel scale"... I can imagine people using the pixel scale in the following ways:

1. Pixel scale may be used to compute pixel area so that we can convert flux to brightness. This is how I use it. Of course, if distortions are present, then pixel area will change across the image, but this is a different issue which could be addressed in the future. Often knowing a "representative" pixel area may be good enough.
2. Recently I have realized that for problems such as computing the linear extent of regions (as in pyregion), it may be more useful to have a representation for the "average" linear size. For example, if you want to estimate the extent (in angle) of a side of a rectangular region then a linear scale may be more useful.

For the first use case above, I believe, as explained in #3113, that we should be using sqrt(|det(CD)|) to compute average pixel scale (thus pscale**2=pixel area). On the other hand, for the average linear scale, it would be better to use pixel perimeter/ number of sides. Essentially, it would be better to take the average of the lengths of each side of the pixel (as it was done in the past).

Now, if users were to get two scales for each axis, how hard would it be to compute their mean? Well, I was thinking that CDELT could give us exactly the lengths of the two sides. However, @nden has mentioned to me that CDELT=(1,1) for WCS that have CD matrix. So essentially we have to re-compute CDELT to get these dimensions. For this reason I will propose below to split these functions. according to the intended use case.

Possibly, there may be other uses for an "average" pixel scale besides the ones mentioned above, so if anyone is aware of those, please contribute to this list.

Assuming you want to keep this function, then I would like to make the following comments about it:

1. I think the "average" pixel scale should be computed as sqrt(|det(CD)|) instead of geometric mean (geometric mean was proposed in celestial_pixel_scale() could return a "better" estimate of scale for non-square pixels #3113 as a sort of "compromise"). The determinant will give the correct area for any king of parallelogram: it does not have to be rectangular or equilateral. The geometric mean of two sides will always overestimate pixel scale (except for rectangles for which sqrt(scaleX*scaleY)=sqrt(|det(CD)|)).

2. So, I am proposing to use sqrt(|det(CD)|) when a single representative pixel scale is required. Those who need to estimate the average linear pixel scale as explained in a previous comment, they could set which='both' and then compute the arithmetic mean of the two values of pixels scales (which is easy to do). All this should be documented in the docstring of the function.

3. I think if you do want to introduce the new parameter which then it should be added after the allow_nonsquare parameter so that this change is backward compatible with the case when users might have been calling celestial_pixel_scale as (for example): celestial_pixel_scale(inwcs,False) in the old code.

4. I would suggest dropping allow_nonsquare parameter and removing any warnings about pixels being "non-square". I cannot imagine a typical use in which knowing that a pixel is non-square is of any importance. For undistorted images (that is, the images containing distortions) we already know that pixels are not square. For a distortion-corrected image (e.g., obtained after running it through AstroDrizzle if it is an HST image) it would be surprising not to have square pixels, etc. In the rare instances in which a user needs to know if the pixel is a square pixel, they could ask to get both pixel scales for each axis and compare them as it is done (though incorrectly - see next point) in this function right now.

5. Let's assume that we have a pixel with sides of length a and b. Testing for np.allclose(a,b) to conclude that a pixel is square is not enough: one can have a rhomb that has all sides equal (a=b) and still it is not a square. So, what I want to say here, is that even if you decide to keep allow_nonsquare parameter, this test for "non-square" pixels should be modified as right now it tests only for non-equilateralness. (Alternatively, you can rename allow_nonsquare to allow_nonequilateral.)

6. I think the line:

   if np.allclose(scale[0], scale[1]) or allow_nonsquare:

   should be modified as right now this function will warn users that pixels are "non-square" even if np.allclose(scale[0], scale[1]) is True...

Because of the different use of "pixel scale" by different users, I would suggest having two different functions:

1. celestial_pixel_scale1 (or linear_celestial_pixel_scale, linear_celestial_pscale, or something like that - so please propose other names) that will always return two pixel scales (one for each axis). Then users can take their mean if they need a single number to be used for problems related to linear sizes of regions (for best results both pixel scales should be used instead of the average). No need to warn users if a!=b - they can check it themselves. NOTE: Actually these two scales should be identical to CDELT (modulo unit conversion), however (see my previous comment above) CDELT does not have the "right" values when CD matrix is defined.
2. celestial_pixel_scale2 (or areal_celestial_pixel_scale, areal_celestial_pscale, or something like that) that will always return sqrt(|det(CD|). This is the pixel scale that when being squared it will give the area of the reference pixel. Again no need to warn users about pixels being non-square (or non-equilateral) - they can investigate this themselves.

These functions (personally I do not like that their names are soooo looong hence "pixel_scale" -> "pscale") should take a single argument - the WCS for which pixel scale needs to be computed.

In the future, if need (or interest) be, we can add (x,y) arguments to both functions so that pixel scale(s) can be computed in any part of the image based on distortion model in the WCS. For now, I think a number (or two) at CRPIX location is what most people want.

EDIT: maybe call these functions celest_lpscale(wcs) (returns a tuple of two pixel scales) and celest_apscale(wcs) (returns sqrt(|det(CD|))?

@mcara - this all sounds reasonable to me, but I wouldn't consider myself a pixel scale power-user (I usually just need a 'representative'. In the interest of time, since it would be good to resolve this before the 1.0 release candidate (in a few days) maybe you could open a new pull request with your suggested changes since it will make it easier to review? If possible, you could include @sargas's commits in the branch so as to give him credit for his work here?

@astrofrog which one of my two suggestions do you want me to implement: 1) splitting pixel scale into two functions; or 2) keeping the old function and applying the suggested (6) fixes?

Also, can anyone instruct on how to do this: "If possible, you could include @sargas's commits in the branch so as to give him credit for his work here?"

Based on @nden comment in #3113 (two scales for each axis) and on the discussion in this PR, I created PR #3304 that replaces celestial_pixel_scale with celestial_pixel_scales (returns a tuple of pixel scales for each axis). Also #3304 introduces a new celestial_pixel_area function that can be used for conversion of fluxes to brightness or for computing a representative pixel scale based as the scale of a square pixel of the same area as the area of the input WCS pixel.

@mcara - can this now be closed?

@astrofrog I believe that #3304 essentially solved the issue of what do we call a "pixel scale(s)"... So, from my point of view, I think this issue can be closed.

@mcara - thanks!

@sargas - thanks for working on this PR, and sorry for not including it in the end!