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

Thomas Robitaille thomas.robitaille at gmail.com
Thu Dec 5 12:36:31 EST 2013

Hi Christoph and Mike,

Montage does exactly this and I started wrapping (with permission from
the Montage developers) the core routine in Python:


I didn't get very far due to lack of time, but you're very welcome to
check this out. You first need to read in an image with WCS, convert
the pixel corners to world coordinates, then use the routine I'm
attempting to wrap, and it should give the overlap area. Because I
didn't get very far yet, it may not return the area yet, but it
should. Feel free to work on the wrapper to improve it!


On 5 December 2013 18:11, Michael Droettboom <mdroe at stsci.edu> wrote:
> 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
> _______________________________________________
> AstroPy mailing list
> AstroPy at scipy.org
> http://mail.scipy.org/mailman/listinfo/astropy

More information about the AstroPy mailing list