I was following along this SunPy example to create a helioprojective WCS from a RA/DEC one. I noticed in the P angle calculation code it uses ITRS to define the geocentric frame. Is this the only valid choice? Would GCRS also work?
I ask because for PUNCH, we are including both a celestial and helio WCS in our products. In our tests to make sure our WCSes actually match, we found offsets of a few pixels. Chris Lowder asked some questions in the SunPy chat a while back related to this. (For context, the test in particular was going from a pixel coordinate in an image to the celestial coordinate to the helio coordinate and then converting back to the pixel coordinate. If the WCSes matched, we were expecting that our input pixel coordinates and output pixel coordinates would agree. That’s where we saw the few pixel mismatch.) However, if we change that P angle code to use GCRS for the geocentric instead of ITRS, our errors greatly reduced, from ~3 pixels to hundredths of a pixel. We wondered if our change made sense because in the example GCRS is used for everything else. I can share a code snippet of our test. I just wanted to check if what we’re doing is ludicrous before going down the rabbit hole. We may be misunderstanding and need more guidance on coordinate systems because there are a lot of details one can get wrong.
If using GCRS for that calculation is okay, I can open a PR that allows the user to pick between frames.
Indeed, there are multiple valid choices for the definition of “north” for an equatorial (i.e., RA/dec) coordinate system. Glossing over some details:
North of ICRF axes, referenced to distant stars: use Astropy’s GCRS
Then adding frame bias: use Astropy’s PrecessedGeocentric with equinox=J2000
Then further adding precession: use Astropy’s PrecessedGeocentric with equinox set to the observation time
Then even further adding nutation: (I don’t think there is a convenient class for this)
Then at last adding polar motion, which gets you to actual geographic north: use Astropy’s ITRS
That is, using GCRS for the calculation of P angle means that you are not including the effects of frame bias, precession, nutation, and polar motion. Whether not including those effects is the correct thing to do depends on the situation.
The reason that the function P() uses ITRS is because it is checked against The Astronomical Almanac, which is a resource primarily for Earth-based observers, and thus the “right” north is the geographic north. There could certainly be user options for other norths, but they should be presented in the API as including/excluding effects.
On the other hand, if you are cross-checking a WCS defined in GCRS, then you need the P angle to be calculated against GCRS as well. As you noted, the example is internally inconsistent. There’s a sentence in the example text that acknowledges that the result will be slightly inaccurate because the angle returned by P() is not quite correct for how it’s going to be used.
I think we want to use GCRS. For now, we’ve implemented a work around in our pipeline where we calculate the P angle using the underlying _sun_north_angle_to_z directly with GCRS instead of the implemented P() function.