On-sky rotations of coordinates with NaN distance are NaN #11423

Activity

1. pllim

   added

   coordinates

   on Mar 25, 2021
2. adrn commented on Mar 25, 2021

   

   adrn

   on Mar 25, 2021

   Member

   More actions

   Hi @nstarman! Thanks for bringing this up. As you may already know from looking at the internals, this is because frame data gets turned into cartesian 3-vectors, so that rotation matrices can be applied simply to the coordinate data to do the transformation. NaN's are greedy, so anything that touches the NaN distance value (i.e. all 3 cartesian components) get turned into NaN values as well, which then get re-represented to spherical and...hence the NaN's appearing in the angle components of the output.

   It's hard to think of a generic solution to this that doesn't require a lot of special casing: It's tempting to want to do something like "if transformation is a rotation, if distance is not finite, put the source on the unit sphere rather than using distance information", but maybe the example transform here is just one in a chain and the next transform requires a distance (e.g., Galactocentric). @eteq @mhvk @StuartLittlefair @astrojuanlu ?
3. mhvk commented on Mar 25, 2021

   

   mhvk

   on Mar 25, 2021

   Contributor

   More actions

   Working around this may be possible, but is there an actual use case? I could see having a catalogue of objects in which some have distances and others have not. In that case, one could have something more like the marker @nstarman suggested, something like "arbitrary distance".
4. nstarman commented on Mar 25, 2021

   

   nstarman

   on Mar 25, 2021

   Last edited by nstarman

   MemberAuthor

   More actions

   @mhvk. That's my use case exactly. I had a table of coordinates from 2 different catalogs, outer vstacked, and tried to make a SkyOffsetFrame with a rotation. I like the concept of "arbitrary distance". I recall but can't find a numpy proposal to have two separate types of NaN, one that indicated a lack of information, the other problematic values. https://numpy.org/neps/nep-0012-missing-data.html

   @adrn Thanks for the explanation. I've been looking at the transforms code. If we want to go the special-case route then I don't think we need to drop the distance information, just restore it after the transformation (see example kludge). This doesn't apply to AffineTransform, but for Static and DynamicMatrix, it should work. This might also solve the edge case of rotations at the origin, which currently lose their angular information.
5. mhvk commented on Mar 25, 2021

   

   mhvk

   on Mar 25, 2021

   Contributor

   More actions

   @nstarman - I can see a case for trying to ensure that if it is a pure rotation matrix, it doesn't touch the distances. This could either be in the general transformation machinery, but another option might be to make .transform() a method that doesn't exist only on CartesianRepresentation. Instead, by default it would go via cartesian, but is overridden for spherical/unitspherical and do not use distance (if the matrix is a rotation matrix - might need some keyword argument or so to indicate that).
6. nstarman commented on Mar 25, 2021

   

   nstarman

   on Mar 25, 2021

   Last edited by nstarman

   MemberAuthor

   More actions

   This could either be in the general transformation machinery

   Looking in the code, two spots in the transformation machinery stand out to me as places to a) check for a pure rotation matrix and avoid routing through cartesian coordinates and b) expand the special casing from UnitSpherical to just Spherical. I haven't started experimenting with fixes. I need to first sit down and better understand the implementation of the transformation machinery. I think part of such a solution would mandate that spherical representations default to being returned in a spherical representation, whereas they now default to cartesian.

   astropy/astropy/coordinates/transformations.py

   Lines 1105 to 1110 in 3b751a0

   	# Convert the representation and differentials to cartesian without
   	# having them attached to a frame
   	rep = data.to_cartesian()
   	diffs = dict([(k, diff.represent_as(CartesianDifferential, data))
   	for k, diff in data.differentials.items()])
   	rep = rep.with_differentials(diffs)

   astropy/astropy/coordinates/transformations.py

   Lines 1136 to 1139 in 3b751a0

   	if isinstance(fromcoord.data, UnitSphericalRepresentation):
   	# Special-case this because otherwise the return object will think
   	# it has a valid distance with the default return (a
   	# CartesianRepresentation instance)

   
   

   another option might be to make .transform() a method that doesn't exist only on CartesianRepresentation

   I like this idea as well. My only concern is that custom representations, such as a log-polar representation, would not immediately benefit.
7. mhvk commented on Mar 25, 2021

   

   mhvk

   on Mar 25, 2021

   Contributor

   More actions

   For custom representations, I don't see how one can very easily benefit - the transformation machinery would also need to know whether the representation is a form of spherical. And if it is really similar, one could inherit from SphericalRepresentation rather than RepresentationBase and just have the .tranform() method. In any case, probably best first to think how most easily to solve the common case of this already quite uncommon example!
8. nstarman commented on Mar 27, 2021

   

   nstarman

   on Mar 27, 2021

   Last edited by nstarman

   MemberAuthor

   More actions

   I think I figured out a recipe for making the spherical .transform method.

   1.) if the matrix is given in a spherical basis this is very easy. Rotation matrices don't have any off-diagonal elements on the 3rd column. To deal with nan, we just need to do the (2x2) angular matrix and the distance row separately.

           _p = np.stack([self.lon.value, self.lat.value, self.distance.value],
                 axis=-1)
           newp = np.empty_like(_p)
           if matrix.shape == (2, 2):
               newp[:, :-1] = np.matmul(matrix, _p[:, :-1].T).T
               newp[:, -1] = _p[:, -1]
           # is it a rotation matrix?
           elif matrix.shape == (3,3) and np.allclose(matrix[:-1, -1], 0):
               # deal with NaN by multiplying the 2x2 first and then the last row
               newp[:, :-1] = np.matmul(matrix[:-1,:-1], _p[:, :-1].T).T
               newp[:, -1] = np.matmul(matrix[-1], _p.T).T
           else:
               newp = erfa_ufunc.rxp(matrix, _p)
           return self.__class__(
               **{k: v * u for v, (k, u) in zip(newp.T, self._units.items())},
               copy=False, differentials=None)
   

   2. For a matrix in a Cartesian basis it's a bit more tricky because the change-of-basis matrix from spherical to cartesian and back will have NaN. One possible solution is to recast the problem into complex numbers with Nan -> 0+1j. Now any use of NaN can be tracked non-destructively by the imaginary component. IDK. I'm having some difficulties implementing this.
9. mhvk commented on Mar 27, 2021

   

   mhvk

   on Mar 27, 2021

   Contributor

   More actions

   The general idea makes sense, but unfortunately the transformation routines all use matrices in a cartesian basis. That said, in those places, we know whether the matrix is a rotation matrix or not, so one can simply set an appropriate keyword argument.
10. mhvk

    closed this as completedin #11444on Apr 20, 2021