[AstroPy] General method to compute pixel solid angle for any WCS projection?

[Michael Droettboom]

On 12/05/2013 09:07 AM, Christoph Deil wrote:
> Hi,
>
> does anyone know a good general method to compute pixel solid angles for arbitrary WCS FITS maps?
>
> I know there is HEALPIX and for some WCS projections there are analytical formulae for pixel area, but I was wondering if there is a standard robust method to implement it for any WCS.
>
> I did find methods to compute the area of polygons on the sphere:
> http://stackoverflow.com/questions/4681737/how-to-calculate-the-area-of-a-polygon-on-the-earths-surface-using-python
> Would representing each pixel as a 4-corner polygon (or subdividing the pixel sides into parts if the pixels are large and precision is needed) and using that method be the best option?

Yes, I was going to suggest something similar.  From a quick glance at 
the stackoverflow page, it seems some of those methods may have 
distortions or discontinuities at the poles, and this is the shortcoming 
of many geography-based methods.

I would do something like:

1) transform the four corners of the pixel using wcs.all_pix2sky

2) project the Ra and Dec into Cartesian point on the unit sphere:

         x = cos(ra) * cos(dec),
         y = sin(ra) * cos(dec),
         z = sin(dec)

3) The area of a well-behaved polygon on the unit sphere is 
"sum(interior_angles) - (n - 2)pi", where n is the number of angles.  So 
first we calculate the interior angles:

Given points A, B, C, the angle at B is:

     acos((B x (A x B)) . (B x (C x B)))

(x for cross product, . for dot product).

I have a spherical geometry library that *somewhat* does this, but given 
that it tries to handle arbitrary polygons (including concave and 
self-crossing ones), it's a lot more complex and unfortunately still a 
little prone to blowing up.  The simple case of a quadrilateral is much 
simpler -- and I should probably add that to my library.

This is still an approximation, of course, as it assumes the edges of 
the pixel are a straight shot along on the unit sphere rather than 
taking into account any distortion at the subpixel level. Usually this 
doesn't matter, of course.

Mike


>
> Christoph
> _______________________________________________
> AstroPy mailing list
> AstroPy at scipy.org
> http://mail.scipy.org/mailman/listinfo/astropy


-- 
                    _
|\/|o _|_  _. _ | | \.__  __|__|_|_  _  _ ._ _
|  ||(_| |(_|(/_| |_/|(_)(/_|_ |_|_)(_)(_)| | |

http://www.droettboom.com