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

[Christoph Deil]

Hi Mike and Tom,

thank you very much for your replies … very helpful!

Christoph

On 05 Dec 2013, at 18:36, Thomas Robitaille <thomas.robitaille at gmail.com> wrote:

> Hi Christoph and Mike,
> 
> Montage does exactly this and I started wrapping (with permission from
> the Montage developers) the core routine in Python:
> 
>  https://github.com/astrofrog/python-reprojection
> 
> 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!
> 
> Cheers,
> Tom
> 
> 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.

Is this library open source?

>> 
>> 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
> _______________________________________________
> AstroPy mailing list
> AstroPy at scipy.org
> http://mail.scipy.org/mailman/listinfo/astropy