Support for InterferometricVisibility frame

[caseyjlaw]

Description

First step reviving this interferometric coordinate calculation, as described in Issue #7766. This work was started as AAS Hack Together Day! 🎉 This is not complete, but hopefully helps start a conversation.

I modified the original work by @Joshuaalbert to use a time and earth location in the frame definition. Then a given RA, Dec can be projected into that frame to get "uvw" coordinates. An example of how to calculate a single baseline with two antennas looking at the same location and time:

uvw_frame0 = coordinates.InterferometricVisibility(location=coordinates.EarthLocation(lon=-u.deg*(107+37/60+6.9/3600), lat=u.deg*(33+54/60+10.3/3600)), obstime=now)
uvw_frame1 = coordinates.InterferometricVisibility(location=coordinates.EarthLocation(lon=-u.deg*(107+36/60+58.4/3600), lat=u.deg*(33+53/60+58.8/3600)), obstime=now)
uvw0 = coo.transform_to(uvw_frame0)
uvw1= coo.transform_to(uvw_frame1)
print(uvw1.u-uvw0.u, uvw1.v-uvw0.v, uvw1.w-uvw0.w)                                
--> -320.8834414416924 -294.0936500374228 356.1691571404226

This result was confirmed (once) by comparing to the code used by the EHT (https://github.com/achael/eht-imaging/blob/master/ehtim/observing/obs_helpers.py#L55). More checks would be helpful. For example, I tried to compare to a memo by CASA developers (https://casa.nrao.edu/Memos/CoordConvention.pdf), but could not get identical uvw values.

closes #7766

[pep8speaks]

Hello @caseyjlaw 👋! It looks like you've made some changes in your pull request, so I've checked the code again for style.

• In the file astropy/coordinates/builtin_frames/intvis.py:

Line 29:71: W291 trailing whitespace
Line 90:65: W291 trailing whitespace
Line 107:34: W291 trailing whitespace

• In the file astropy/coordinates/tests/test_frames.py:

Line 1241:1: E302 expected 2 blank lines, found 1

Comment last updated at 2020-01-24 18:37:39 UTC

[Joshuaalbert]

@caseyjlaw Thanks for taking this up! I have two things. To clarify, the tests defined in here, which should reflect (https://casa.nrao.edu/Memos/CoordConvention.pdf), are failing in all cases, or just some edge cases? It might be that the compVectors function I defined is insufficient, perhaps simply replace with np.all(np.isclose(a,b)).

Secondly, for those tests I make use of the East-North-Up coordinate frame that I also co-committed during that PR, but I don't think it was accepted because it is equivalent to the AltAz frame with the x and y axes swapped. Does this effect the unit tests at all?

[caseyjlaw]

@Joshuaalbert Indeed, I am removing reference to ENU in the code (just pushed revised tests). Your tests make sense to me, but it seems that I have an error somewhere. Could you check the logic of the way I use uvw_frame in the latest definition in test_frames.py?

[mhvk]

Thanks for the revival! I do worry this is not quite accurate; e.g., the definition of "North" (for RA, Dec) needs to include polar motion - it is only approximately in the direction of the Earth's polar axis (which of course precesses...). I think that one will need to follow what is done in AltAz though omitting the pressure and temperature parts (which are relevant for optical only). Indeed, AltAz should give one the right unit vectors to help project against.

p.s. If this is only confusing, I can look into it in more detail later, but it would in any case help greatly to include the test cases that are discussed above - we should really reproduce what other packages do (or understand why they are wrong ;l-)

[caseyjlaw]

@mhvk Glad to help!
I think there is some more fundamental error in the current version. I started from Thompson, Moran and Swenson to get the xyz to uvw coord transform. But it seems rotated relative to the tests suggested in the NRAO doc. For example, the first test is to confirm that xyz aligns with uvw when LST=0 and Dec=+90. However, this is what I get:

> print(coo.transform_to(uvw_frame).u.value, coo.transform_to(uvw_frame).v.value, coo.transform_to(uvw_frame).w.value)                                       
1090.8357691960427 -6186.437066030218 1100.2485477353619
> print(uvw_frame.location.x.value, uvw_frame.location.y.value, uvw_frame.location.z.value)                                                                  
6186.437066030218 1090.8357691960427 1100.2485477353614

[mhvk]

@caseyjlaw, @Joshuaalbert - a more basic question occurred to me: are we doing this the right way around? Normally, our coordinate frames are representations of the same sky coordinates of actual objects, but here UVW are really just representations of the positions of the observatory, right? Given that, might it make more sense to have it be a method on either EarthLocation or SkyCoord (latter in analogy with radial_velocity_correction, which is similar in that it really just calculates a velocity of the observatory projected along the line of sight to an object)? cc @eteq who is one of the few who really understood the coordinate transformation graph...

p.s. this frame of reference question might also explain the sign change...

[Joshuaalbert]

@caseyjlaw You have to be careful using the ITRS frame to compare with xyz as defined in the NRAO doc. They are not the same. That is why I defined the East-North-Up system. Take a close look at the ITRS documentation and the coordinate frame talked about in the NRAO doc. Use your right hand, and perhaps draw on an empty coconut shell, and you'll see.

@mhvk I don't see why the fact that astrophysical objects wouldn't commonly be written in the UVW frame should mean that it's not a coordinate frame. What about the ITRS for example? Fundamentally, astronomy is an optical science, and there are many ways to represent optical objects. We use UVW to represent the locations, in the plane of the sky, of point sources. Thus, I think we are doing this the right way. It's just that our frame is representing the conjugate coordinates of physical locations at infinity.

[mhvk]

@Joshuaalbert - thanks, that does clarify it!

I do still think one will need to go through something like CIRS to get the coordinates relative to the real North pole, not the ICRS one - what is there now can, I think, only be approximately right. Or, looking at that documentation linked to above, the observatory coordinates have to be transformed into an XYZ system where Z points at the ITRS north celestial pole (dec(J2000)=90), not the Earths (lat=90).

[Joshuaalbert]

@mhvk Good point. Okay, so one way to resolve this question is to look at the motivation behind the construction of east-west baseline radio antennas that exist like WSRT. These telescopes all fall on a line such that as the earth rotates a source at DEC=0 has zero V coordinate always, and a source at DEC=90 traces a circle in the UV plane. That was the intuition behind them at least. Of course because the celestial pole is not the same as Earth's pole, it's only an approximation that's tied to to closest epoch coordinate standard e.g. J2000. So I think we want to indeed fix the north pole to the CIRS north (nutation-precession taken into account) pole of Earth, and not the dec=90 pole on the J2000 or similar celestial sphere.

[mhvk]

Makes sense. We just have to be sure we follow whatever is the actual convention, and do that accurately. Probably the tightest matches we need to make would be from VLBI, where typical precisions in telescope positions (and thus UVW) are at the sub-cm level.

[caseyjlaw]

Based on comments above, I've pushed a change to transform to CIRS to improve accuracy.

@Joshuaalbert I can see that the xyz from ICRS is different than the xyz in Thompson, Moran, and Swenson, but the transformation is not clear to me. I've looked over your east-north-up transformation, but I don't think that's what I need to get to the TMS xyz system. If you have a rotation matrix in your head that I can implement, I'm happy to do it.

[mhvk]

@caseyjlaw - for the transformations, note that you can do quite a bit within the representation classes themselves - see http://docs.astropy.org/en/latest/coordinates/representations.html#vector-arithmetic. In particular, I think you don't really need to calculate most angles and do the explicit projections.

But that is an implementation detail. Most important really is to add tests against results from trusted software!

[caseyjlaw]

@mhvk I didn't know about those transformations. I'm happy to use them, although I'm still unclear on how to go from the "xyz" defined by ICRS and "xyz" in TMS.
I believe the revised tests in tests_frames.py cover the comparison of the code to the NRAO description.

[Joshuaalbert]

@mhvk agree. Get the tests well defined first.

@caseyjlaw Currently, I don't see test_frames.py failing. I do see

astropy/coordinates/builtin_frames/intvis.py F                           [  9%]

[caseyjlaw]

The docstring test in intvis.py wasn't written properly, but I just fixed it. It should pass now.
The test in test_frames.py is not being run by travis because it seems to skip those that require remote_data. Testing locally, I can see it fails because the xyz definition is wrong, as we've discussed here.

[StuartLittlefair]

This isn't quite right. Hour angle is E + L - RA, where E is the earth rotation angle, L is the longitude (east positive) and RA is the CIRS RA (see 2.2 of this paper). E can be found using erfa.era00

[StuartLittlefair]

I think we also need to be careful about the transform to CIRS. The icrs_location may not have an obstime attached to it, and the frame you want to transform to is the CIRS frame at the obstime of the uvw_frame - e.g

cirs_frame = CIRS(obstime=uvw_frame.obstime)
cirs_loc = icrs_loc.transform_to(cirs_frame)

[StuartLittlefair]

If you put this together, I think the correct code should be (you'll need to import erfa and get_jd12, which is is coordinates.builtin_frames.utils):

cirs_frame = CIRS(obstime=uvw_frame.obstime)
cirs_loc = icrs_loc.transform_to(cirs_frame)
era = erfa.era00(*get_jd12(uvw_frame.obstime, 'ut1'))* u.rad
ha = era + lon - cirs_loc.ra  # unit-ful calculation
dec = icrs_loc.transform_to(CIRS).dec  # no need for transformation of units

[caseyjlaw]

@StuartLittlefair Thanks. I'm working that into my local branch.
After that, the uvw is calculated by the dot product with xyz in a local coordinate system. We start with xyz as EarthLocation, but we need to transform to xyz as defined in TMS (East-North-Up on the equator).
Is that transformation clear to you?

[caseyjlaw]

Ugh, now I see that the TMS definition for uvw is different from that of the NRAO! I think that is part of my confusion.

[StuartLittlefair]

TBH I’m still trying to wrap my head around the orientation of the XYZ coordinate frame, let alone the transformation to UVW.

[caseyjlaw]

@mhvk I think we are stuck on some fundamental confusion on the coordinate systems. Part of the problem is that I wanted to base everything on the TMS textbook transform. But they use a coordinate system that I don't understand in the astropy context. Can you give us any tips on that interpretation?

[mhvk]

Agreed. What I think we need is someone who can point us to more definitive definitions. So, I wrote Adam Deller. From a quick glance myself, it would seem that we really should be looking inside the veritable calc/solve software...

[mhvk]

Here is Adam Deller's reply:

Sometimes, it seems like the more fundamental the issue, the less likely it
is to be well documented! I suspect that for a lot of radio
interferometers they don't care at the level you're going to want to dig
down to, which is an explanation if not an excuse.

For DiFX-correlated data, we do indeed use CALC (and SPICE for orbital
antennas, but that's an even bigger can of worms.) At different sites and
different times, CALC8, CALC9, and CALC11 have all been supported. And of
course anyone is free to update their local site data to have whatever
position for whatever antenna they like. For normal VLBI operations using
DiFX, the site positions are taken from the vex file, so it is ultimately
sched and its catalog that is setting the positions of the antennas.
Sched has the following to say about the x,y,z positions:

X: Station X coordinate in meters. This is in the direction of the
Greenwich meridian. Can be put in the locations catalog.
Y: Station Y coordinate in meters. This makes a right handed coordinate
system with X and Z. Can be put in the locations catalog.
Z: Station Z coordinate in meters. This is in the direction of the north
pole. Can be put in the locations catalog.

And the locations.dat file itself (at least, my one) has the following
entry in the catalog info for all of the well-known radio antennas:

FRAME='GSFC Feb. 26, 2015: 2015a'

What does that exactly mean? I suppose Craig Walker must know, or Walter
Brisken would probably know too. I assume it is a given realisation of an
ITRF produced by GSFC.

The problem with station positions is that the goedesists don't care (they
work with total delay, so they don't care what the initial model was) and
for pretty much anyone else the differences between different (modern, at
least) realisations of the ITRF so small relative to other unknown
propagation effects that you really have to go digging to find a case where
it matters. (Now I'm making excuses for myself!!)

Anyway, that is how the delays are derived. On to the uvws. w is of
course delay multiplied by the speed of light. The u and v are derived by
taking the derivative of delay with respect to the l and m coordinates -
see Section 9.2 of
www.atnf.csiro.au/vlbi/dokuwiki/lib/exe/fetch.php/difx/difxuserguide.pdf
for a short description. I should note also that for the purpose of
calculation uvw this way, only the pure free space delays are used - for
the delay model used to align the data we also include all of the wet/dry
troposphere delays, etc, that CALC can estimate, but using these would bias
the uvw coordinates, so they are all stripped out from the delays used to
calculate the uvws. Said differently, the affects of atmospheric
refraction are ignored when calculating the uvws, as they should be.

[mhvk]

Thinking about Adam's reply, I again wonder what is actually the best description of what UVW are. It is really relative to a "delay centre", which means not like a normal coordinate frame but rather like a SkyOffsetFrame. I think this is helpful because it separates the concerns of which direction that phase centre is relative to the telescopes (which should be a calculation very similar to what is already done for AltAz) with how the directions of UVW are defined (which is already done by SkyOffsetFrame, and depends on the system the phase centre is defined in, i.e., can be SkyOffsetICRS, SkyOffsetFK5, etc.).

It also seems we may have to get a CALC server running to write tests against...

[mhvk]

Further on my last comment, I asked Adam about how to get to the source code of CALC:

CALC lives in several places. The easiest one to work with (at least, the
one I'm most familiar with) is probably the CALC 11 version bundled with
DiFX: svn.atnf.csiro.au/trac/difx/browser/applications/difxcalc11/trunk if
you have access to the DiFX trac pages. If not, you can still just check
out the code from svn as described on
www.atnf.csiro.au/vlbi/dokuwiki/doku.php/difx/installation (use the latest
stable tag, which is DiFX-2.6.1 as described on that page), and then the
code itself is under applications/difxcalc11/trunk/. You should be able to
compile all of DiFX and run vexdifx + difxcalc on a sample vex file to
produce a .im file, which has the delay stored in polynomial form for each
antenna, should then be easy enough to check against. Somewhat oddly,
no-one has ever written a convenient python interface to the im file,
although it would be pretty trivial adn could be added to the parseDiFX
library under libraries/parseDiFX/.

[ludwigschwardt]

Wow, talk about a timely issue... 😄

TL;DR: I want to help!

I'm busy converting our telescope coordinate library (katpoint) to use Astropy instead of PyEphem. When I created it 12 years ago, PyEphem seemed like the easiest backend, and Astropy did not even exist. It brought us much joy and pretty images, but now it is time to upgrade. I want to replace as much homegrown functionality with Astropy equivalents, and this PR is certainly welcome.

I implemented UVW coordinates the same way @Joshuaalbert did, via AltAz and an ENU local coordinate system. I quickly realised that this is a murky part of radio astronomy, with sparse documentation and inconsistent implementations. I checked it against DiFX's CALC back in 2013 and it was surprisingly accurate, even with single-precision floats and without proper leap second support, as long as you stick to a connected interferometer.

We revisited this as part of SKA prototyping in the past 3 months, and with some help I now have a Pythonised version of CALC 11 running for new benchmarks. You can imagine my delight when I stumbled upon this thread... 😂

[mhvk]

@ludwigschwardt - help is definitely appreciated! And in some sense good to know that you found little documentation - then at least it is not just us. To me, it continues to make sense to then try to match the de-facto standard that is calc.

Have you thought about how this might work in detail in the astropy context? Specifically, I think that since UVW is always relative to some pointing centre, which is different for each observation, the SkyOffsetFrame idea might be a good model to work from.

p.s. I'd be quite interested in just having access to calc11 via python - do you have it somewhere accessible?

[Joshuaalbert]

@ludwigschwardt Welcome! and glad to know that our approach is pretty close to CALC back in 2013. @mhvk I am not sure what the SkyOffsetFrame is, can you explain better why that would be useful? In my implementation I just include a ICRS coordinate to act as the phase tracking center, similar to how you can set an observing time and reference location.

If we were to take to same approach as DiFX, then we should look firstly at:

This above is similar to using an AltAz and ENU frame. I think that is the direction of attack that we should continue to look into. That is, compute uhat, vhat, and what in any convienient coordinate frame. what is easily found by transforming the pointing in ICRS to the choosen coordinate frame. vhat is tricky. They say delta_J2000=90 (so I guess they mean the mean pole), but I think it verified if it should be fixed to the time of observation. From this uhat follows, and you can transform any antenna coordinate to UVW.

Then, consider that a radio user will want to calculate uvw's for millions of times. In this case, we should look for when the obs_time is a vector and if it is then apply DiFX's trick. They use numeric differentiation and fit a 5-th order polynomial to d tau/dl and d tau/dm every 120 seconds. Each model point requires 3 calls to CALC which computes the d tau/dl and d tau/dm, where tau is w/c. This is explained here:

Crucially, what ever you choose for defining the celestial point vhat, should agree with the definition of l and m.

As I'm Canadian I'm obliged to say sorry for the long rant!

[mhvk]

@Joshuaalbert - see link in #9869 (comment); it seems custom made for UVW, since those also are coordinates relative to some reference centre. Let me cc @eteq, who was more involved in SkyOffsetFrame to see what he thinks about the analogy. To me, given that UVW involves both an observatory and a reference position, it certainly is much more like SkyOffsetFrame than, say, Galactic. Though of course an AltAz aspect does come in.

But most important remains to have some good test cases...

[ludwigschwardt]

p.s. I'd be quite interested in just having access to calc11 via python - do you have it somewhere accessible?

@mhvk, I'm depending on ALMA's software but I'm still trying to uncover its license. Let me get back to you. I'll also have a look at SkyOffsetFrame in the meantime (I'm an Astropy novice, which is embarrassing to admit 😂).

That is, compute uhat, vhat, and what in any convienient coordinate frame. what is easily found by transforming the pointing in ICRS to the choosen coordinate frame. vhat is tricky. They say delta_J2000=90 (so I guess they mean the mean pole), but I think it verified if it should be fixed to the time of observation. From this uhat follows, and you can transform any antenna coordinate to UVW.

@Joshuaalbert, this is exactly how we started out, and as you say, vhat is the tricky part. Then we realised (quite recently!) that vhat=NCP_J2000 is ever so slightly wrong sub-optimal, since the astrometric-to-topocentric conversion doesn't simply rotate the celestial sphere, but also distorts it. It's therefore better to use a local North and not the global NCP. See more details in ska-sa/katpoint@a99cdf1.

Each model point requires 3 calls to CALC which computes the d tau/dl and d tau/dm, where tau is w/c.

To elaborate: CALC basically produces tau for (ra, dec). To get UVW the "DiFX way", you need 3 delays: tau(ra, dec), tau(ra + dx, dec) and tau(ra, dec + dy). The small steps in the l and m directions are of the order of 0.0001 rad (20"), and the differences in tau are used to compute UVW.

[mhvk]

@Joshuaalbert - thanks! Definitely don't worry about not knowing yet about SkyOffsetFrame - nor about being a novice - it is after all what every one of us was before.

[mhvk]

Still super-important but also sadly still not entirely clear how to get it right, so changing milestone to 4.2...

[demorest]

Hi all, I've been looking into this recently. I agree with @mhvk that SkyOffsetFrame sounds like it should be essentially the same thing as the usual UVW coordinate system. I think at least in principle there shouldn't be much/any extra code needed in astropy to support this, except maybe some convenience functions if people want those.

I've gone through a test UVW calculation using this method, you can see the code here:
https://gist.github.com/demorest/8c8bca4ac5860796593ca07006cc3df6

This uses VLA B-configuration antenna positions (~10-km maximum baseline). It also includes the same calculation using CASA for comparison. The astropy results agree with CASA to ~10s of cm, which is probably good enough for some purposes, but might not scale well out to VLBI baselines. I don't know whether CASA is expected to work in the VLBI regime either, so a comparison to CALC as suggested earlier in this thread is still a good idea.

I suspect the discrepancy with CASA is related to some GCRS vs ICRS subtlety but I haven't fully got my head around this yet. Also CASA is basically a black box here, I don't really know what it is doing at this level of detail. One clue may be that in astropy, the SkyOffsetFrame can be set up using either ICRS or GCRS as a base. I naively expected ICRS would be the right choice, but using GCRS seems to agree better with CASA. When GCRS is used, the W values agree very well, to about 1 cm, and the differences in U,V seem to be a small rotation about the W-axis. With ICRS, the differences don't have an obvious easy interpretation but don't seem random (ie I don't think it's a numerical precision thing for example).

I have not yet tried the "delay derivatives" method for computing U, V for comparison but that shouldn't be too hard so I will give that a look soon. Any other comments or questions are welcome!

[Joshuaalbert]

Thanks @demorest for doing this comparison. Ideally astropy be as precise as possible to be usable for VLBI also. Much work has gone into making astropy's coordinate system incredibly precise so it would be shame not to carry that through all the way to UVW coordinates. I'd be curious to see how the delay coordinate derivative method compares. Also, I think you are right to use GCRS.

[ludwigschwardt]

Hi all, sorry for the absence... I've spent the past two months benchmarking my package (katpoint) against CALC and Astropy for calculating interferometric delays. I tested it on our telescope layout, with baselines planned out to 20 km. The results also suggest that Astropy should do well, although I used a simplistic ICRS -> ITRF path length difference calculation.

I could share the report with anyone that's interested. @demorest, I'll take your script as a cue and add CASA to the mix :-) I'll also try out GCRS to see if it closes the small gap I still have.

I'm currently stuck converting katpoint to Astropy but would really like to help this effort, since katpoint also needs (U, V, W) :-)

[ludwigschwardt]

A perk of my study is an incrementally better understanding of CALC... 😂

The code describes the relevant coordinate systems thus:

!     The basic coordinate system is referenced to the Epoch of 2000.0 and is a
!     right-handed Cartesian system oriented to the mean celestial pole and mean
!     equator of that epoch. The nominal origin is the solar system barycenter.

This sounds like ICRF?

!     There is also an Earth fixed coordinate system which is a right-handed
!     Cartesian system oriented to the mean geographic pole of 1900-1906 and the
!     Greenwich Meridian. The nominal origin is the Earth's center of mass.

Almost ITRF, with the older but still close CIO pole? (That's Conventional International Origin, not SOFA's Celestial Intermediate Origin 😄)

There are several options for UVW:

• "uncorrected"
• "aberrated"
• delay derivative version, potentially including aberration and/or atmospheric refraction

All of them operate on J2000.0 source vectors and "J2000.0 geocentric" baseline vectors. The latter sounds like GCRS.

[demorest]

After thinking a bit more, computing UVW baselines in GCRS coords probably does make sense. These will eventually be used to calculate phases. If the radio freq (or wavelength) values used for phase calculation are topocentric (or maybe geocentric for VLBI?), the baselines should be in a consistent reference frame.

I mentioned a rotation between GCRS SkyOffsetFrame and CASA UV. The magnitude of this increases as a function of declination. Near the equator it's arcsec-level, at dec=80deg it's about an arcmin, and seems to blow up at the pole (these are rotations of the image field about its center, so shifts in source positions will be much smaller depending on FoV). Maybe this comes back to the "which north pole" question.

Was anyone able to track down (or implement) a CALC python interface? This seems like something that would be useful generally..

[ludwigschwardt]

If the radio freq (or wavelength) values used for phase calculation are topocentric (or maybe geocentric for VLBI?), the baselines should be in a consistent reference frame.

Yes. The UVW coordinates essentially have a Fourier relationship with the sky coordinates, so it makes sense that they are locked to the sky. The Earth is just incidental in the process :-)

Maybe this comes back to the "which north pole" question.

Or maybe no fixed (global) North Pole... We got rid of this issue (blow-up at the poles) by doing something like the delay derivative approach in ska-sa/katpoint#46. We find a local North direction by adding a small offset in declination towards the equator, which seems to give better results.

Was anyone able to track down (or implement) a CALC python interface?

Indeed. I've been using a combination of ALMA's CALC11 and a Python wrapper developed as part of SKA work. I was unsure if there would be a licensing issue but that seems OK, so I'm hoping to release this on PyPI in the near future.

[astropy-bot]

Hi humans 👋 - this pull request hasn't had any new commits for approximately 5 months. I plan to close this in a month 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.

[astropy-bot]

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.