Repository navigation
Cannot create GCRS frames with obsgeoloc shapes other than (3,) #4926
Description
After discussion with @eteq, the code that was merged for #4890 included the use of GCRS frames whose obsgeoloc and obsgeovel vectors could have shapes equal to (3, N), where N is the number of times at which the routine EarthLocation.get_gcrs_posvel calculates the position and velocity of an Earth-based observer.
However, in
, the shape of the obsgeoloc arrays are hardcoded to(3,).
This means that any attempt to call coordinates.get_body with a non-scalar time and a defined location will raise an error
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
<ipython-input-12-1102973f833a> in <module>()
----> 1 coord.get_body(tarr2, 'sun', location=wht)
/usr/local/lib/python3.5/site-packages/astropy/coordinates/solar_system.py in get_body(time, body, location)
252 gcrs = icrs.transform_to(GCRS(obstime=time,
253 obsgeoloc=obsgeoloc,
--> 254 obsgeovel=obsgeovel))
255 else:
256 gcrs = icrs.transform_to(GCRS(obstime=time))
/usr/local/lib/python3.5/site-packages/astropy/coordinates/baseframe.py in __init__(self, *args, **kwargs)
482
483 # Validate input by getting the attribute here.
--> 484 getattr(self, fnm)
485
486 args = list(args) # need to be able to pop them
/usr/local/lib/python3.5/site-packages/astropy/coordinates/baseframe.py in __get__(self, instance, frame_cls)
209 out = self.default
210
--> 211 out, converted = self.convert_input(out)
212 if instance is not None and converted:
213 setattr(instance, '_' + self.name, out)
/usr/local/lib/python3.5/site-packages/astropy/coordinates/baseframe.py in convert_input(self, value)
332 raise ValueError('The provided value has shape "{0}", but '
333 'should have shape "{1}"'.format(value.shape,
--> 334 self.shape))
335 if (oldvalue.unit == value.unit and hasattr(oldvalue, 'value') and
336 np.all(oldvalue.value == value.value)):
ValueError: The provided value has shape "(3, 2)", but should have shape "(3,)"
Activity
mhvk commentedon May 18, 2016
This particular problem might be more easily solved by making
geolocandgeovelbeCartesianRepresentationinstances with arbitrary shape.More generally, having seen the many
np.rollaxis(...xyz, 0, ...ndim)in the code, I wonder if rather than continuing to have quantities with a first dimension of 3, we shouldn't have a vector class. (And if we're going for that, I'd also have our ownMatrixclass -- numpy's is quite broken -- and start using__matmul__).eteq commentedon Jun 13, 2016
@mhvk:
I'm rather skeptical of this. In my experience with coordinates thus far, when I've experimented with doing this it always made things more confusing, not less. There's a lot of places where the broadcasting is non-trivial, involving one component broadcasting over another one in a way that needs to be aligned in some other way with the next operation, not necessarily in a 3-vector manner. Additionally, this would It adds one additional layer of abstraction over numpy arrays, which has to be documented/tested/etc, when this is exactly what numpy is meant for (and, more importantly, has a lot of potential devs who know it).
That said, I do see the potential value, and could believe that at some point adding something like this makes the transformation bits less complicated enough to justify it... But I guess I'm not convinced we're there yet.
mhvk commentedon Jun 13, 2016
Darn, it seems I am a bit like a broken record, having forgotten my comment here when creating #5088...
I think it is quite clear, though, that numpy is clear in what it expects: coordinates/matrices are the last (fastest varying) dimensions, not the first. The same holds for erfa, so there should be some benefit in making that the default on our side too.
Fortunately, I don't think there is a reason for a big change, especially if we stick to using things like the
CartesianRepresentationclass without exposing the internal API to such an extent that we're forever stuck on how the data is stored internally.StuartLittlefair commentedon Jun 13, 2016
I think I'm with @mhvk on this one. Whilst I'm much less familiar with the coordinates code than @eteq I believe that all the times I've had to confront this issue was to move the 3-vector to the last dimension for broadcasting and back again.
Added to this is the fact that many of the bugs we've found so far have been related to broadcasting not working correctly.
The simplest way I can imagine it working is a
Vectorclass, with a shape of(..., 3)and aMatrixclass with a shape of(..., 3, 3), and then convenience functions to make these objects from aCartesianRepresentation.eteq commentedon Jun 15, 2016
I'm willing (happy, even!) to believe that I'm wrong on this one, because it's definitely been a long source of frustration, as you say @StuartLittlefair ... But I think it needs to be rather carefully prototyped (with the possibility of it then being thrown out if it doesn't simplify things) to be sure.
StuartLittlefair commentedon Jun 15, 2016
👍 to that comment. There's little sense in introducing these classes unless they make dealing with non-scalar coordinates and times simple and robust.
in-progress band-aid fix to astropy#4926
eteq commentedon Jun 19, 2016
FYI, anyone who might be working on this (@mhvk or @StuartLittlefair), I have a regression test for this in 2ca07bb that any are welcome to grab
pllim commentedon Jul 26, 2016
Everyone, please review #5196.