Simplified celestial_pixel_scale (now re-named to proj_plane_pixel_scale) and new proj_plane_pixel_area - #3304

@astrofrog The tests should be passing, I believe... Travis failed probably because of some network error. I do not know what it wrong with AppVeyor.

EDIT: Just found the following message in AppVeyor: "Build execution time has reached the maximum allowed time for your plan (30 minutes)." So, I think the code is just fine but I do not know how to restart the tests.

@mcara - thanks! I've restarted the tests.

@astrofrog AppVeyor still fails with message: "Build execution time has reached the maximum allowed time for your plan (30 minutes)." Is it possible to increase the allocated time for the test, or we just disregard this particular test?

Also, what you you think of the proposed changes? CC: @nden

I am ready to submit a new PR with a function that will test in a correct way if pixels are non-square (that is, if they contain distortions other than projection effects - projection effects are normally understood to be present with any kind of non-linear mapping, except Cartesian).

It would have been nice (=easier for me) if this PR was merged and then I would just modify the new master instead of an older version of the file. (I also understood from #3230 that there is some urgency in getting this done ASAP).

I opened an issue with appveyor: http://help.appveyor.com/discussions/problems/1371-build-being-interrupted-after-30-minutes

In the mean time, maybe @keflavich (who wrote the original function) and @mdboom and/or @nden could review this?

@mcara - hopefully this will get merged in the next couple of days, so that you can open the other PR early next week.

I think it is better if name of the functions indicates that the pixel scales (pixel area) are computed in the projection plane from CD (PC, CDELT) and not on the "celestial sphere" (as the old name might have suggested - "celestial" in the old name was referring to "celestial axes" as opposite to, e.g., spectral). This is to clearly indicate that the computed scales ignore distortions due to non-linear nature of most projections.

Also, I have modified proj_plane_pixel_scales to compute the scales of the original input WCS and not of the astropy.wcs.WCS.celestial WCS object as it was done in the old celestial_pixel_scale. It seems to me that it is perfectly fine to compute the scale for non-celestial axes.

The PR seems fine to me, I just nitpicked one small thing. Of course, this PR removes the convenience function of deriving an approximate pixel scale, but the pixel area ought to be good enough.

The previous version also returned a Quantity. I think it's fine and sensible to return just a float and let users inspect the WCS to figure out what the units are.

My only use case involves using separate scales so I'm happy with the way this is written.

@nden your suggestion to return a numpy array is actually what I was thinking of doing (see my answer to @keflavich below). The reason I have used tuple(map(float,...)) was that when I was testing possible uses for the returned values, I got an error because the elements of the tuple were np.float64 (instead of Python floats) - I do not remember exactly how I got this (a few days ago) and I cannot reproduce it anymore.

@keflavich [Issue Nr. 1]: "Of course, this PR removes the convenience function of deriving an approximate pixel scale, but the pixel area ought to be good enough." It is true that this function removes the convenience factor, but it seems that (after discussing this with @nden and also by looking at various ways these scales are being used in the STScI software and Pyregion, etc.) there is no single good "average" pixel scale that is suitable for a large number of possible use cases. For this reason it seems better to let users decide what they want when they need a single number. It should be pretty easy to say proj_plane_pixel_scales(wcs).mean() if they need to, or sx,sy=proj_plane_pixel_scales(wcs) to get two scales as, e.g., @nden wants, etc.

I also thought that adding the which parameter (as well as allow_nonsquare) this early may hinder the ability to expand this function in the future. Specifically, I have the following "vision" (please forgive apparent but not intended grandioseness of the word "vision") for this function. At this moment we compute the scales and area only at CRPIX location ignoring SIP and table lookup distortions. However, in the future, I was thinking (assuming there will be interest in doing this) to add (x,y) parameters to this function to compute scales (area) at other locations in the image. Essentially, these two functions could be made to behave similarly to wcs.all_pix2world or wcs.pix2foc. That is, we could pass either arrays of [X], [Y] coordinates (arrays of length N), or arrays of form [[x1,y1],[x2,y2],...] of size Nx2. Adding such ability after the optional parameters which and allow_nonsquare, I think, will be difficult.

@keflavich [Issue Nr. 2]: "The previous version also returned a Quantity. I think it's fine and sensible to return just a float and let users inspect the WCS to figure out what the units are." Again it is true that I remove the "Quantity"/units from being returned. It seems fine to return a quantity for two scales (or a single pixel area). However, if in the future we will return an array of pixel areas and if this function will have input parameters (x,y) or just "xy", then it seems to me that it will be more difficult to deal with quantities for the following reasons:

1. This function will behave differently from similar functions such as wcs.all_pix2world, wcs.pix2foc, etc. in the sense that it will return Quantities instead of values;
2. Really, when (x,y) is provided then the units of x(y) will be important too. It is clear what to do when input units are 'pixels' but it is not clear to me what to do if units are 'mm', etc. Of course we can throw an exception if units are != 'pixels' - it is a solution. Except I do not think I have time to figure this out right now: I think it would be better to have some open discussion about this. However, if we do not want to replicate the behavior of wcs.all_pix2world or wcs.pix2foc then of course it is possible to return Quantities.

Please let me know if there is a strong desire to return Quantities and I can put them back (I do not have any strong opinion on this except that I think that by returning just values it allows more options for future expansion).

So, the vision (for the future) is that a function like proj_plane_pixel_area(wcs, x, y) could be used to produce "pixel area" maps.

@mcara - I'm happy with all the changes you made, I was just noting them. You have a reasonable idea of what this function is supposed to be; I had just implemented something that worked for one use case at the time. So, this is a resounding 👍 for the PR.

@keflavich Thank you very much for your comments! I really appreciate them and I just want to make clear that I do not have a "strong opinion" on how these functions should be (it is just my own "vision" based on my own limited experience) and I am really open to all suggestions as to how to make these functions more useful, clear, etc. to the users. Personally in my work I would use "pixel area" (and not "scale") so any suggestions/comments on "scale" are welcome... Except my understanding is that this should be done before the v1.0 release.

@astrofrog I have no idea what is wrong with the tests (travis): I have not even modified anything in the astropy/io/votable/tests/util_test.py or whatever test follows it, or why the following error appears in another failing test: ImportError: No module named 'pkg_resources'. AppVeyor test are being "queued" for half a day now...

I believe all tests should pass as they passed yesterday and the changes that I made to the code are too small to affect the tests.

@mcara - I think the failures were unrelated indeed and I've restarted Travis.

@astrofrog @mdboom If anyone could merge this, then I would be able to create another PR for checking if pixels are non-square.

This is great, thanks!

If I understand correctly I'm pretty sure this is affects-dev. But none of this work is mentioned in the changelog and I'm wondering if it should be, and to what extent. Are these new interfaces public, or just for internal use?

It seems the addition of celestial_pixel_scale was mentioned in the changelog, but now that it has been renamed that changelog entry should also be updated, and should include a link to this PR as well as the original one.

Yes, this needs an updated changelog entry (i.e. we need to change the existing one). I'll do that.

Done in a654376

The function that can perform a more complete check for pixels being non-square is proposed in #3329. Users can call that function to see whether pixels are non-square and then decide what to do about the two scales returned by proj_plane_pixel_scale.