Frame attributes in sky coordinates do not obey Numpy broadcasting rules #11571
Description
Description
When creating a SkyCoord, the dimensions of the representation components and the frame attributes do not respect Numpy broadcasting rules.
In the example below, I am trying to apply a grid of sky offset frames to calculate the footprints of a telescope on a grid in right ascension, declination, and roll. I had to resort to manually broadcasting the frame attributes using np.broadcast_to for this to work.
Expected behavior
Expected output:
fov_lon.shape=(4,) fov_lon.shape=(4,)
center.shape=(360, 181)
frame.shape=(360, 181, 90)
new_shape=(360, 181, 90, 4)
with manual broadcasting, result.shape=(360, 181, 90, 4)
with automatic broadcasting, result.shape=(360, 181, 90, 4)
Actual behavior
Actual output:
fov_lon.shape=(4,) fov_lon.shape=(4,)
center.shape=(360, 181)
frame.shape=(360, 181, 90)
new_shape=(360, 181, 90, 4)
with manual broadcasting, result.shape=(360, 181, 90, 4)
Traceback (most recent call last):
File "/opt/local/Library/Frameworks/Python.framework/Versions/3.8/lib/python3.8/site-packages/astropy/coordinates/attributes.py", line 114, in __get__
out = np.broadcast_to(out, instance_shape, subok=True)
File "<__array_function__ internals>", line 5, in broadcast_to
File "/opt/local/Library/Frameworks/Python.framework/Versions/3.8/lib/python3.8/site-packages/astropy/units/quantity.py", line 1512, in __array_function__
return super().__array_function__(function, types, args, kwargs)
File "/opt/local/Library/Frameworks/Python.framework/Versions/3.8/lib/python3.8/site-packages/numpy/lib/stride_tricks.py", line 411, in broadcast_to
return _broadcast_to(array, shape, subok=subok, readonly=True)
File "/opt/local/Library/Frameworks/Python.framework/Versions/3.8/lib/python3.8/site-packages/numpy/lib/stride_tricks.py", line 348, in _broadcast_to
it = np.nditer(
ValueError: input operand has more dimensions than allowed by the axis remapping
During handling of the above exception, another exception occurred:
Traceback (most recent call last):
File "./test.py", line 37, in <module>
result = SkyCoord(fov_lon, fov_lat, frame=frame[..., np.newaxis]) # KABOOM
File "/opt/local/Library/Frameworks/Python.framework/Versions/3.8/lib/python3.8/site-packages/astropy/coordinates/sky_coordinate.py", line 330, in __init__
self._sky_coord_frame = frame_cls(copy=copy, **frame_kwargs)
File "/opt/local/Library/Frameworks/Python.framework/Versions/3.8/lib/python3.8/site-packages/astropy/coordinates/builtin_frames/skyoffset.py", line 178, in __init__
super().__init__(*args, **kwargs)
File "/opt/local/Library/Frameworks/Python.framework/Versions/3.8/lib/python3.8/site-packages/astropy/coordinates/baseframe.py", line 601, in __init__
values[fnm] = getattr(self, fnm)
File "/opt/local/Library/Frameworks/Python.framework/Versions/3.8/lib/python3.8/site-packages/astropy/coordinates/attributes.py", line 117, in __get__
raise ValueError(
ValueError: attribute rotation should be scalar or have shape (4,), but is has shape (360, 181, 90, 1) and could not be broadcast.
Steps to Reproduce
Sample code:
#!/usr/bin/env python
from astropy.coordinates import SkyCoord
from astropy import units as u
import numpy as np
# Square field of view of telescope
fov_width = 50 * u.deg
fov_lon = np.asarray([0.5, 0.5, -0.5, -0.5]) * fov_width
fov_lat = np.asarray([0.5, -0.5, -0.5, 0.5]) * fov_width
print(f'{fov_lon.shape=} {fov_lon.shape=}')
# RA, Dec, and roll grid
ra = np.arange(0, 360) * u.deg
dec = np.arange(-90, 91) * u.deg
roll = np.arange(90) * u.deg
# Telescope boresight grid
center = SkyCoord(ra[:, np.newaxis], dec[np.newaxis, :])
print(f'{center.shape=}')
# Frame grid with boresights and roll angles
frame = center[..., np.newaxis].skyoffset_frame(roll)
print(f'{frame.shape=}')
# Expected output shape
new_shape = np.broadcast_shapes(frame[..., np.newaxis].shape, fov_lon.shape)
print(f'{new_shape=}')
# Broadcast manually. This works fine!
frame_bcast = np.broadcast_to(frame[..., np.newaxis], new_shape, subok=True)
fov_lon_bcast = np.broadcast_to(fov_lon, new_shape, subok=True)
fov_lat_bcast = np.broadcast_to(fov_lat, new_shape, subok=True)
result = SkyCoord(fov_lon_bcast, fov_lat_bcast, frame=frame_bcast)
print(f'with manual broadcasting, {result.shape=}')
# Let Numpy do the broadcasting for us.
result = SkyCoord(fov_lon, fov_lat, frame=frame[..., np.newaxis]) # KABOOM
print(f'with automatic broadcasting, {result.shape=}')System Details
>>> import platform; print(platform.platform())
macOS-10.15.7-x86_64-i386-64bit
>>> import sys; print("Python", sys.version)
Python 3.8.9 (default, Apr 3 2021, 06:24:49)
[Clang 11.0.3 (clang-1103.0.32.62)]
>>> import numpy; print("Numpy", numpy.__version__)
Numpy 1.20.2
>>> import astropy; print("astropy", astropy.__version__)
astropy 4.2.1
>>> import scipy; print("Scipy", scipy.__version__)
Scipy 1.6.2
>>> import matplotlib; print("Matplotlib", matplotlib.__version__)
Matplotlib 3.4.0
Activity
pllim commentedon Apr 15, 2021
Just for completeness, have you tried with the latest dev version of
astropy?lpsinger commentedon Apr 15, 2021
It's still broken on main. I get the same error message with these package versions:
adrn commentedon Apr 15, 2021
Ah yes. I think this is also likely related to #8812
pllim commentedon Apr 15, 2021
@adrn , is this issue a duplicate of #8812 then?
mhvk commentedon Apr 15, 2021
Agreed it is the same issue as #8812, but this is a particularly good example for why we need to solve it!
lpsinger commentedon Apr 15, 2021
Yes, this is a duplicate, but the title of this issue makes it more clear that the problem is general to any frame attribute, not just times.
pllim commentedon Apr 15, 2021
I updated the title of the older issue and closing this one. Thanks, everyone!