SkyOffsetFrame with longitude / inclination instead of longitude / latitude

[maxnoe]

Description

To express quantities that are expected to be radially in the field of view, a version of the SkyOffsetFrame that
works with longitude / inclination, instead of longitude / latitude as the current implementation does, would be very helpful.

This frame would rotate so that the origin is then at inclination = 0, phi=0 (arbritrary really, since it is the pole).
In small angle approximation, this would be equivalent to having polar coordinates around the origin.

I would be interested in contributing this frame.

Any pointers on how to achieve this, would be great.
I assume this needs a new representation class that uses inclination instead of latitude?

EDIT by @mhvk: link to how the frames are hoped to work: https://gamma-astro-data-formats.readthedocs.io/en/latest/general/coordinates.html#field-of-view

[adrn]

Hi @maxnoe --

Could you say a little more about the idea and some use cases? I'm a little unclear on the motivation, but not opposed to the idea! FYI we do have a representation PhysicsSphericalRepresentation that expects longitude/colatitude instead of longitude/latitude, so I think that would work for this case too.

[maxnoe]

Could you say a little more about the idea and some use cases?

Sure! The SkyOffsetFrames are really helpful for us to express detector coordinate systems.
Our telescopes have large fields of view (up to ten degrees diameter), so treating this as a spherical system is required.

We are using the current implementation of SkyOffsetFrame to have the telescope pointing position in lon = 0, lat=0, which creates a something close to a cartesian system.

For some applications, having a polar coordinate system around the pointing position would be more appropriate. This would be a SkyOffsetFrame that rotates so, that the pole is the new origin, so lon = 0, colat = 0 is equal to origin.

FYI we do have a representation PhysicsSphericalRepresentation

Ah, nice, so that's already there.

[nstarman]

Given PhysicsSphericalRepresentation, is this issue more about initializing a SkyOffsetFrame with longitude / colatitude?

[mhvk]

Discussed this in coordinates telecon: right now, the location of the origin in the offset frame is hardcoded to be along the X axis (ra=dec=0). It should be possible to add another attribute that tells where it should be instead (we have the rotation attribute already). Would need a change in the transformation, done right now as

astropy/astropy/coordinates/builtin_frames/skyoffset.py

Lines 68 to 78 in a371a59

	@frame_transform_graph.transform(DynamicMatrixTransform, framecls, _SkyOffsetFramecls)
	def reference_to_skyoffset(reference_frame, skyoffset_frame):
	"""Convert a reference coordinate to an sky offset frame."""
	# Define rotation matrices along the position angle vector, and
	# relative to the origin.
	origin = skyoffset_frame.origin.spherical
	mat1 = rotation_matrix(-skyoffset_frame.rotation, 'x')
	mat2 = rotation_matrix(-origin.lat, 'y')
	mat3 = rotation_matrix(origin.lon, 'z')
	return matrix_product(mat1, mat2, mat3)

[mhvk]

Though since the representation would have different attributes, not d_lon, d_lat, but rather d_theta and d_phi.

Maybe a new class is needed, though one could also let SkyOffsetFrame have an __init__ method that overrides _default_representation based on what is input/wanted.

[maxnoe]

Today, I was looking into this again, as we'd like to standardize our field of view coordinate systems soon and need also the version with "polar coordinates" in the FoV.

I had the idea that it might be possible to already use the existing SkyOffsetFrame with a 90° rotated origin and the PhysicsSphericalRepresentation (which btw. does not seem to have a Unit... equivalent):

from astropy.coordinates import SkyCoord, SkyOffsetFrame, PhysicsSphericalRepresentation, Angle
import numpy as np
import astropy.units as u

crab = SkyCoord.from_name("Crab")
zeta_tauri = SkyCoord.from_name("Zeta Tauri")


offset_frame = SkyOffsetFrame(origin=crab)
crab_offset = crab.transform_to(offset_frame)
zeta_tauri_offset = zeta_tauri.transform_to(offset_frame)

# find the origin that would put crab into the north pole
polar_frame = SkyOffsetFrame(origin=crab.directional_offset_by(0 * u.deg, 90 * u.deg))
crab_polar = crab.transform_to(polar_frame)
zeta_tauri_polar = zeta_tauri.transform_to(polar_frame)


theta_ticks = np.arange(0, 2, 0.25)

fig = plt.figure(figsize=(8, 4))

ax1 = fig.add_subplot(1, 2, 1)
ax1.set_aspect(1)


for s, label in zip([crab_offset, zeta_tauri_offset], ['Crab', 'Zeta Tauri']):
    ax1.plot(s.lon.deg, s.lat.deg, 'o', label=label)

    
ax1.legend()
ax1.set_xlim(-2, 2)
ax1.set_ylim(-2, 2)
ax1.set_xlabel('offset lon / deg')
ax1.set_ylabel('offset lat / deg')

for theta in theta_ticks:
    ax1.add_artist(plt.Circle(xy=(0, 0), radius=theta, fill=False, linewidth=0.5, color='gray'))
    
for y1, y2 in [(-1, 1), (1, -1)]: 
    ax1.axline([-1, y1], [1, y2], color='gray', lw=0.5)

    
ax1.grid()
ax1.set_title('SkyOffsetFrame(origin)')



ax2 = fig.add_subplot(1, 2, 2, projection='polar')
ax2.set_title('SkyOffsetFrame(origin + 90°) with using theta/phi')
# ax2.set_theta_zero_location('N')


for s, label in zip([crab_polar, zeta_tauri_polar], ['Crab', 'Zeta Tauri']):
    
    r = s.represent_as(PhysicsSphericalRepresentation)
    ax2.plot(r.phi.rad, r.theta.deg, 'o', label=label)

ax2.set_yticks(theta_ticks)  
ax2.set_ylim(0, 2)
ax2.legend()

It results in a rather exactl 90° degrees rotation with respect to the "naive" polar angle calculated from the lat/lon in the normal sky offset frame:

>>> zeta_tauri_polar.represent_as(PhysicsSphericalRepresentation).phi.to(u.deg)
<Longitude 219.83057441 deg>
>>>  np.arctan2(zeta_tauri_offset.lat, zeta_tauri_offset.lon).wrap_at(360 * u.deg).to(u.deg)
<Angle 309.83349484 deg>

which I don't really get...

Plot:

[maxnoe]

Adding also 114 Tauri shows it's clearly just a rotation:

[maxnoe]

I experimented with adding a polar=True attribute to the SkyOffsetFrame. Will open a PR, but at the moment I still have a different 90° rotation 🙈