Renewed attempt to add the InterferometricVisibility frame - #10737

Hello @Joshuaalbert 👋! Thanks for opening this pull request, we are very grateful for your contribution! I'm a friendly 🤖 that checks for style issues in this pull request, since this project follows the PEP8 style guidelines. I've listed some small issues I found below, but please don't hesitate to ask if any of them are unclear!

• In the file astropy/coordinates/interferometric_visibility.py:

Line 22:5: W292 no newline at end of file

Hello. I converted your PR to a draft PR due to #10706 .

Thanks very much - the description of the tasks looks good. Agreed also to start with just getting as-exact-as-possible solutions, even if they are expensive, and worry about interpolation later (though I can resist noting that we just merged a PR that allows it for the current transforms - #10646; even though that focuses on erfa routines, I think the infrastructure would be useful here too).

I'd also encourage consideration of an approximate solution that is appropriate for connected interferometers. Computational intensity is one concern. It would be good to know if the SkyOffset approach will allow a fast, vectorized calculation in time (use case: calculate uvw for 8 hours of data with 1 second samples).

I'm going to look at the SkyOffset code to get acquainted with it.

Ok after looking at SkyOffset we may not use it for the exact solution. From SkyOffset docs:

    Coordinates in this frame are both centered on the position specified by the
    ``origin`` coordinate, *and* they are oriented in the same manner as the
    ``origin`` frame.  E.g., if ``origin`` is `~astropy.coordinates.ICRS`, this
    object's ``lat`` will be pointed in the direction of Dec, while ``lon``
    will point in the direction of RA.

This won't work for two antennas in different reference frames.

@Joshuaalbert - have you checked the rotation argument to the initializer? Also, if it is close, I think it would be fine to adapt or subclass SkyOffsetFrame. The question may be the extent to which actual transformations are similar to what we need:

p.s. Just to be sure: I think a conclusion from the previous discussion was that the SkyOffsetFrame idea was quite similar to what's needed for InterferometricVisibility, in things being based around a particular pointing direction. If it is not obvious and easy, it would seem totally fine to start with a set of classes that work and verify they reproduce calc, before worrying too much about whether SkyOffsetFrame is directly usable itself.

Hey @mhvk I mean due to relativistic effects the notion of a pointing offset is incorrect in the exact formulation suitable for VLBI. In the exact solution the UVW coordinates are a function of both the antenna position and the direction of the radio source. The w coordinate is the speed of like times the group delay,

and u and v are the partial derivatives with respect to the first two components of s in a Cartesian system.

p.s. that's not exact still because it's lacking some terms, but is still very accurate. The full group delay is in section 11.4.1 of https://www.iers.org/SharedDocs/Publikationen/EN/IERS/Publications/tn/TechnNote36/tn36.pdf

Ah, thanks! (I fear you had in fact told me this before...) So, it sounds like, if anything, SkyOffsetFrame would be a special case. The idea of SkyOffsetFrame to set oneself up with a meta-class might still be a good one, though.

To get an idea of the importance of the relativistic kinematic (non gravitational) corrections here I show the difference in delay between including and not including them:

Baseline=1.00 km, max diff=0.33ns, mean diff=0.17ns
Baseline=2.78 km, max diff=0.93ns, mean diff=0.47ns
Baseline=7.74 km, max diff=2.58ns, mean diff=1.30ns
Baseline=21.54 km, max diff=7.18ns, mean diff=3.62ns
Baseline=59.95 km, max diff=19.99ns, mean diff=10.07ns
Baseline=166.81 km, max diff=55.63ns, mean diff=28.02ns
Baseline=464.16 km, max diff=154.80ns, mean diff=77.96ns
Baseline=1291.55 km, max diff=430.73ns, mean diff=216.93ns
Baseline=3593.81 km, max diff=1198.52ns, mean diff=603.62ns
Baseline=10000.00 km, max diff=3334.96ns, mean diff=1679.61ns

What we can do is provide an option that computes all components with delays above a certain threshold. E.g. Most telescopes will get away with 0.1ns precision.

Hello, given the draft status of this PR and the impending feature freeze tomorrow, I am moving the milestone. Thank you for your understanding.

I have created a python package that wraps CALC11 (https://github.com/casangi/radio_telescope_delay_model/blob/main/README.md). Currently, it only wraps the almacalc function which calls the main CALC function aDRIVR. The Fortran files of interest are calc11/src/almacalc.f and calc11/src/adrvr.f

Here is a link to the example notebook that can be run on Google Colab (https://colab.research.google.com/github/casangi/radio_telescope_delay_model/blob/main/example_notebooks/Radio_Telescope_Delay_Model_Example.ipynb). It compares CALC, CASA, Astropy, and data from a measurement set. The Astropy and CASA calculations are done using code from @demorest. There is no exact agreement and the ms uvw coordinates show a strange difference pattern with the other methods. If I understand correctly the uvw coordinates in an Alma ms (converted from asdm) are done in https://open-bitbucket.nrao.edu/projects/CASA/repos/casa6/browse/casa5/code/alma/apps/asdm2MS/UvwCoords.cc . I have not yet done a careful comparison between UvwCoords.cc and CASA code in the notebook. In the "Run and Compare Functions" section of the notebook, you can change the pointing direction by changing use_user_pointing and user_pointing_direction.

Let me know if I should expose any of the other CALC functions.

@Jan-Willem - thanks so much! Will need some time to really go through it all...

@Jan-Willem great job! I'm taking a look at the notebook now.

Quite strange the difference between MS and the others. Perhaps something related to precision? Thanks for making this notebook. I think it'll be good for comparing against the astropy version which should be as precise or better than CALC.

Hi all, after reviewing some of the older discussion, I'm still not sure a new frame is needed. It still seems to me that the existing SkyOffsetFrame does the job of defining the "l" and "m" directions in the usual way relative to a specified sky position. It would be great to have VLBI-level UVW calculation implemented using astropy, possibly using the derivative method as mentioned above, but I'm not seeing how this needs a new coordinate frame or what that frame would be.

Thanks @Jan-Willem for posting the new comparisons! When I looked at this before, it seemed that the difference between my astropy and CASA UVWs was a small rotation that got larger for sky positions near the poles. But I never really tracked down the cause. Also I looked a bit at CASA code but I think you have much more experience with it than I do.. if UvwCoords.cc is doing something different that would be good to know.

It seems like it would be easy to implement a high-level function that could call either the "astropy-like" function or "calc-like" function, based on the application and precision required. If we took that approach, we could create a PR now that uses the astropy-like function and document the fact that it has known variances from CASA and calc.
I agree that we haven't actually fully characterized those variances, though.

Indeed, there is precedence for that idea of allowing users to chose in how one sets, e.g., the solar system ephemeris. But as you note, we first have to be sure we understand the variances.

Hi @demorest, the point is that for the precise formulation, the UVW coordinates are a function of position, direction (of the geodesic), and velocity with respect to the solar system barycenter. I agree SkyOffsetFrame could be useful if delay is linearised in direction, as it just does an arbitrary rotation of a right-handed-system:

origin = skyoffset_frame.origin.spherical
mat1 = rotation_matrix(-skyoffset_frame.rotation, 'x')
mat2 = rotation_matrix(-origin.lat, 'y')
mat3 = rotation_matrix(origin.lon, 'z')

However, even in that case the origin and rotation arguments would need to be calculated per antenna, direction, and time. I also agree that we can provide varying levels of precision. I personally, don't see why it would be difficult to provide the VLBI-precise transform. Astropy has all the ingredients already.

Slightly randomly but I started another way to access CALC11, one that uses numpy.f2py so one can access all subroutines, over at https://github.com/mhvk/pycalc11. For now, mostly just tests.

This happened to be picked up by @aelanman who mentioned his own efforts to see how consistent astropy and CALC11 were (and more interface routines at https://github.com/aelanman/difxcalc_astropy). He also mentioned the https://github.com/RadioAstronomySoftwareGroup/pyuvdata package that does non-relativistic UVW (and apparently uses astropy throughout).

No chance yet to look into any of this in detail...

Hi humans 👋 - this pull request hasn't had any new commits for approximately 11 months. I plan to close this in 30 days if the pull request doesn't have any new commits by then.

In lieu of a stalled pull request, please consider closing this and open an issue instead if a reminder is needed to revisit in the future. Maintainers may also choose to add keep-open label to keep this PR open but it is discouraged unless absolutely necessary.

If this PR still needs to be reviewed, as an author, you can rebase it to reset the clock.

If you believe I commented on this pull request incorrectly, please report this here.

Hi humans 👋 - this pull request hasn't had any new commits for approximately 11 months. I plan to close this in 30 days if the pull request doesn't have any new commits by then.

In lieu of a stalled pull request, please consider closing this and open an issue instead if a reminder is needed to revisit in the future. Maintainers may also choose to add keep-open label to keep this PR open but it is discouraged unless absolutely necessary.

If this PR still needs to be reviewed, as an author, you can rebase it to reset the clock.

If you believe I commented on this pull request incorrectly, please report this here.

Ops, sorry for the double bot post. Hopefully I fixed the bug upstream.

I'm going to close this pull request as per my previous message. If you think what is being added/fixed here is still important, please remember to open an issue to keep track of it. Thanks!

If this is the first time I am commenting on this issue, or if you believe I closed this issue incorrectly, please report this here.

@caseyjlaw I have written a JAX+astropy based VLBI-accurate UVW engine here https://github.com/Joshuaalbert/DSA2000-Cal/blob/main/dsa2000_cal/dsa2000_cal/common/uvw_utils.py for the DSA2000 project and wonder if there is any desire still to have this functionality in astropy. #9869

Definitely still would like to have UVW in astropy (though not a dependency on JAX; is it optional?)

The implementation relies on the auto-diff functionality of JAX to compute the I've coordinates. Without it derivatives would have to be written out by hand.

This is the core part that relies on auto diff https://github.com/Joshuaalbert/DSA2000-Cal/blob/4dfd9d2d85c215850de2ef01fc6efef940f3ce7f/dsa2000_cal/dsa2000_cal/common/uvw_utils.py#L312

We do have transformations that do automatic differentation too -- FunctionTransformWithFiniteDifference -- so hopefully this is not a showstopper!

Nice to see this work coming back into astropy.
@Joshuaalbert Would you suggest we open a PR to start discussing the question of dependencies, etc. with astropy developers?

Sure we can open another PR for that. @mhvk that might be an option, do you have some examples of how that might work?

I think I wrote a bit too quickly, FunctionTransformWithFiniteDifference does the finite difference only to transform velocities, so that is not directly relevant here. But I think we really should discuss in #7766, since this is a closed PR, and hence will not be easily found/tracked!!

I started an entry at #7766 (comment) to transfer us over to #7766.