Repository navigation
celestial_pixel_scale() could return a "better" estimate of scale for non-square pixels #3113
Description
@astrofrog suggested to raise this issue with astropy. I believe that currently celestial_pixel_scale() returns the arithmetic mean of scales for each axis: scale=1/2*(scale_X+scale_Y).
However, current definition is not related in any way to other important quantities such as pixel area. If there are no strong reasons to keep this definition, e.g., if astropy does not try to mimic the behavior of some other package, I would like to propose to use a definition for the "average" scale such that scale**2=pixel area. This definition works for both square and non-square pixels. Of course this is just one possible definition and other choices may be possible, so it may be worth to discuss this in order to obtain more suggestions.
If proposed definition is adopted, then for non-square (but orthogonal) axes, one could use the geometric mean scale=sqrt(scale_X*scale_Y) - a formula used in the old pyregion package. For non-square and non-orthogonal axes, I would suggest using scale=sqrt(abs(det(CD))) as explained in astropy/pyregion#54 This latest formula is correct (I think) for square and non-square pixels (orthogonal and non-orthogonal) such that scale**2=pixel area.
Activity
nden commentedon Nov 17, 2014
In my experience a pixel scale is a somewhat controversial term. Would it be better if the code returns the pixel scale along each direction and let the application code decide what to do with them.
mcara commentedon Nov 17, 2014
@nden - that would definitely be one solution. It just depends on the original intention for introducing
celestial_pixel_scale(): do users want two scales or just one number? For the later I suggest to use a convention such thatscale**2=pixel area=sqrt(abs(det(CD)))**2. In my opinion this would be the most natural way of representing pixel scale through a single number (possibly there are other natural ways of representing scale through a single number). Maybe there should be two functions for the two usage models:celestial_pixel_scale()(oraverage_celestial_pixel_scale()) andcelestial_pixel_scales()(ends in s) that returns pixel scale along each axis?astrofrog commentedon Nov 17, 2014
The motivation was to have a 'simple' number that could be used in most cases, but of course the pixel scale technically varies over the image in addition to being potentially different along x and y.
mcara commentedon Nov 17, 2014
There is little that can be done about pixel variation across the image except, maybe, to introduce another function that will return pixel scale(s) at a given location. Or, we could average (integrate) over the entire image
abs(det(CD))which then will represent a more "averaged" pixel scale as it does now (only at CRPIX).However, if people are satisfied with current definition - I am fine with it too: in the end everyone can derive their own scale based on the WCS elements. What I want to say is that I do not want to insist on changing this function at all. But codes such as astropy/pyregion#54 will have to compute scale differently if they want to get results similar to DS9.
astrofrog commentedon Nov 18, 2014
I'll try and comment more on this later, but I do agree the geometric mean may be more appropriate here.
sargas commentedon Nov 18, 2014
I think some use case analysis would be useful here. Leaving everyone to derive their own methods could result in edge-case bugs, potentially ad-hoc motivations, misunderstandings of the properties of the technique they are applying, etc.
I think there is a good service in providing well-documented variants in
astropyof any reasonable way to do this. This doesn't include using the arithmetic mean if it has no justification over geometric.(As an aside, I'd like to propose an API of
celestial_pixel_scale[s](wcs, center_point='all'), wherecenter_point=='all'averages over the image somehow,center_point=='center'just uses the center (a fast option), andcenter_point==SkyCoord/BaseCoordinateFrameuses a specific point).mcara commentedon Jan 16, 2015
Based on @nden comment above (two scales for each axis) and on the discussion in #3230, I created PR #3304 that replaces
celestial_pixel_scalewithcelestial_pixel_scales(returns a tuple of pixel scales for each axis). Also #3304 introduces a newcelestial_pixel_areafunction 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.astrofrog commentedon Jan 24, 2015
Resolved by #3304