astropy and pyrap.measures parallactic angle implementations out by 10 arc seconds

[sjperkins]

There's a test case for this here.

The coordinate systems don't quite match up for the different quantities:

• antenna positions. pyrap 'itrf', astropy 'itrs'
• field centre. pyrap 'J2000', astropy 'fk5'
• pole. pyrap 'unknown/default', astropy, 'fk5'

Need to figure out which coordinate systems are equivalent.

[o-smirnov]

AltAz is not a valid frame for antenna positions, must be something else...

[sjperkins]

I think I shouldn't have said AltAz. The antenna position code for the astropy case is here and uses EarthLocation.from_geocentric which specifies a position in ITRS (of which ITRF is a realisation)

Changed the initial comment.

[o-smirnov]

Can't you set the frame to 'j2000'? Does it recognize that?

[sjperkins]

Looks like the astropy equivalent of J2000 is ICRS (International Cannabinoid Research Society):

If you’re looking for “J2000” coordinates, and aren’t sure if you want to use this or FK5, you probably want to use ICRS. It’s more well-defined as a catalog coordinate and is an inertial system, and is very close (within tens of milliarcseconds) to J2000 equatorial.

I've tried setting the field centre and pole frames to 'icrs', but this doesn't change the numbers.

[sjperkins]

@ludwigschwardt @bmerry Any ideas on getting astropy and pyrap parallactic angles to agree? They're out by around 10 arcseconds

pyrap
astropy

[bmerry]

I also didn't get them to agree. Oleg put me in touch with some CASA person (forget who, but I can dig it out of my email if needed). From what I remember the short answer is that there are many things that affect it - so between two systems there is almost certainly something that one ignores, or for which they use different models. I got the impression that 10" is down in the noise for polarisation angles.

[ludwigschwardt]

Maybe use katpoint as the deciding vote 😀

[o-smirnov]

Looks like the astropy equivalent of J2000 is ICRS (International Cannabinoid Research Society):

But is J2000 offered as a separate option, subtly inferior as it may be? Radio astronomers are sticks-in-the-mud, they like their J2000 corrdinates.

[twillis449]

Anybody figure this one out? I think the problem is that astropy precessed positions seem to be incorrect. I compared the output of astropy, skyfield and CASA measures (see the attached script for 3C147). the outputs of skyfield and CASA measures are very close, but astropy output is way off, and consequently hour angle is also incorrect. Since parallactic angle depends on hour angle the astropy hour angle also disagrees with the other two. I'm attaching a script which shows the differences. If anyone can spot an error in the script, let me know. (As @o-smirnov knows, Wim Brouw is infaillable.) Its also interesting that despite the astropy precessed positions being wrong, the astropy azimuth and elevation agrees pretty well with the output from the other two packages so astropy would appear not to make use of the precessed positions to determine az and el.

positions_test_3C147.py.gz

[bmerry]

Something that looks odd to me from the output is that the errors you're getting in HA (over 35") are much bigger than the errors in az/el (which all agree to within a few arcsec). I would have expected any major issues with precession to show up in az/el as well.

My guess is it's the way you're computing apparent position in the astropy code: you're precessing to the current epoch, but that ignores all the other effects that affect apparent position (and which the conversion to AltAz will be handling), such as relativistic aberration. I think converting to CIRS will give you a better answer, but then hour angle has to be computed using Earth Rotation Angle (ERA) rather than sidereal time. I'm not sure how to get that from astropy, although there is a formula to get it from UT1.

Some other things to check:

• Are any of these packages doing refraction? I haven't used Skyfield, but I vaguely recall that PyEphem (the predecessor) does refraction more appropriate to optical astronomy unless it's disabled. Fiddling with the time a bit to increase elevation would reduce refraction effects.
• Is 3C147 near the Sun at the given time (within a few degrees)? There are general relativity effects that I don't know if all the packages do.

Reading over this it sounds like I understand this stuff better than I actually do - I have just enough knowledge to be dangerous, found by reading whatever Google turns up with no organised understanding. So please feel free to contradict me :-)

[twillis449]

You raise some of the issues I've thought about myself. Why does AZ/El agree pretty well with all three systems, while astrypy apparent RA/DEC is so out of wack? I believe that both skyfield and CASA measures include Earth wobbles / nutation etc when computing apparent positions as they can both update Earth orientation parameters (if you don't keep your CASA geodetic etc parameters updated it complains about accuracy problems.) astropy is also supposed to do so. I doubt that any of these packages are doing relativistic corrections. 3C147 is at declination +50 deg and the Sun never gets above declination 23 degrees so relativistic effects don't affect 3C147 (3C279 is the QSO that gets close to (eclipsed by?) the Sun once per year and is used to test general relativity, but the effect is only a few seconds of arc at most.) 3C147 does get close to zenith angle zero here so I'll try fiddling the time a bit to see if there's any refraction effects being included. I did a test run using CasA instead of 3C147 but similar differences appear.

[twillis449]

It looks like J2000 and ICRS only differ by a few milli-arcseconds in positions at most whereas my astropy apparent positions are off by many arcsec so I don't think J2000 coordinates vs ICRS coordinates is the issue.

[bmerry]

Astropy and Skyfield both definitely account for aberration, and I would assume casacore does too. According to Wikipedia it can account for up to 20", so it could well be a contributor to the difference you're seeing.

[sjperkins]

@twillis449 (and @bmerry) Thanks very much for the script, the investigation and for making me aware of skyfield. It's a lighter dependency than astropy.

I haven't yet had a look at your discussion in depth, I aim to take a closer look next week.

[bmerry]

@twillis449 if I plug in these lines of code for the astropy analysis in place of the existing computation, then HA agrees with casacore to 0.01" and parallactic angle to 0.005". This code isn't even that precise because CIRS is a geocentric rather than topocentric system. From what I can see, the underlying ERFA calculations that astropy uses actually returns the hour angle, but astropy discards it.

ThreeC147_App = ThreeC147.transform_to(CIRS(obstime=t))
print 'Astropy time, app_ra app_dec:', t, ThreeC147_App.ra.deg, ThreeC147_App.dec.deg
print 'Astropy time, app_ra app_dec:', t, ThreeC147_App.ra.hms, ThreeC147_App.dec.dms

# Eqn (14.1) of Meeus' Astronomical Algorithms
era = (0.7790572732640 + (t.ut1.jd1 - 2451545.0 + t.ut1.jd2) * 1.00273781191135448) % 1 * 2 * math.pi
era = Longitude(era, u.rad).to(u.hourangle)
era_rel = Longitude(era + drao.lon, u.hourangle)
LST = t.sidereal_time('apparent')
print 'LST', LST, 'ERA', era, 'relative ERA', era_rel
Hd = (era_rel - ThreeC147_App.ra).deg
print 'astropy hour angle', Hd
print 'drao latitude ', drao.lat.deg

Hr = (era_rel - ThreeC147_App.ra).radian
astropy_pa = parallactic_angle(Hr, ThreeC147_App.dec.radian, drao.lat.radian)
print 'parallactic angle derived from astropy', astropy_pa

[twillis449]

@bmerry - thanks for this. :-) I'd never heard of CIRS before. You are correct re astropy EFRA calculations - that's why astropy returns a correct Az/El with ICRS coordinates. Of course the apparent ra is now the CIRS ra which differs from the 'standard' apparent ra, so I converted back to 'standard' ra as we now have LST and correct hour angle. I'm appending my complete updated script. Anyway I would say that CASA measures is still the best way to go as it seems to cover all bases pretty well. :-)

Reading the astropy docs a few things became clear - the ICRS system does not understand annual aberration so adjustment to FK5 as some particular equinox does not include annual aberration in the output and apparent positions will be wrong. I guess CIRS does include annual aberration. I also tried using GRCS, which also includes aberration. The apparent ra and dec output by GRCS still have about 3 arcsec errors with respect to CASA measures values, so CIRS is the clear winner here.

positions_test_3C147.py.gz

[sjperkins]

I've had a look at the script. It appears that skyfield is the one that's further out, rather than astropy?

Measures parallactic angle 99.1272843652
measures ha- astropy parallactic angle diff in arcsec 0.215725744749
measures ha - skyfield parallactic angle diff in arcsec -4.255417405
measures ha_app- astropy parallactic angle diff in arcsec 0.430822106949
measures ha_app - skyfield parallactic angle diff in arcsec -4.0403210428

@twillis449 Would you consider opening an issue on skyfield (or astropy if I've got the wrong end of the stick completely) with your script?

In any case, I've reverted to CASA as the VOICE OF AUTHORITY in #49, but kept the test case illustrating the difference between casa and astropy.

I'm keeping this issue open.

[bmerry]

I've recently realised there is a small problem in the way I compute parallactic angles in katsdpimager, which I think was the basis for @sjperkins' implementation, and also affects katpoint: the transformation from astrometric (ra, dec) to topocentric (az, el) isn't just a rotation of the celestial sphere, it's also a distortion (I think mostly due to aberration, but refraction would also play a role if one modelled it). The result is that the direction "north" of a point on the sky (i.e. derivative of position with respect to dec) is not the same as the direction of a great circle passing through the target and the north pole.

I don't recall exactly how much the distortion is, but I became aware of this issue because there was a very similar problem affecting katpoint's UVW calculations (that caused the UV plane to be slightly rotated because V wasn't aligned with north) that was significant enough to be noticeable in images.

[sjperkins]

@bmerry Thanks, I assume you're referring to the changes you made in ska-sa/katpoint#46?

[bmerry]

@bmerry Thanks, I assume you're referring to the changes you made in ska-sa/katpoint#46?

Exactly.

[bmerry]

I've taken another look at this to ensure that katsdpimager does the right thing. From what I can see, there are a number of differences between implementations (referring to the explicit parallactic angle functionality in katsdpimager, casacore and katpoint plus the assorted test scripts that @twillis449 and @sjperkins have pointed me at):

• How the latitude is interpreted. The casacore derived columns work in AZEL (which treats it as geocentric i.e. as if the Earth was spherical), whereas AZELGEO surely makes more sense. I think everything not based on casacore is effectively AZELGEO.
• Which definition of "north" is used. Some code (including katpoint, and casacore derived columns) use the apparent pole i.e. the position that will appear fixed as the Earth rotates. That might be useful for building polar-mount telescopes but doesn't seem that useful for polarisation work (converting between antenna H/V and Stokes parameters) for which the ICRS/J2000 pole seems more useful. This can make almost arbitrarily large differences for targets near the pole.
• How the angle is measured. Pretty much all the code I've seen just takes a position angle of some pole relative to the target in some coordinate system (mostly converting the north pole into AzEl, but I also saw some code that converts Zenith back into J2000). That's fine when there is no change in velocity (e.g. when casacore derived columns convert between HADEC and AZEL), but when there is, the distortion of great circles skews the result. In a quick test (that doesn't explore the parameter space much) I've seen deviations of up to 2 arcmin between using the north and south poles in the calculation, for example. A better solution (which I plan to implement in katsdpimager) is to use a point infinitesimally close to the target but offset towards north (in practice, some finite amount to avoid problems with numeric stability).
• The coordinate system in which the angle is measured. My gut feel is that the effects here are too small to worry about, but I suspect the conversion between astrometric and topocentric frames is not conformal (i.e., won't preserve angles between vectors).

Of course, a lot of the effects might not matter since they'll be swallowed up in calibration (as long as one is consistent).

[twillis449]

@bmerry As I recall you are correct in that casacore (still) uses AZEL as the default specification for the derivation of Azel from a latitude and longitude. I have no idea why that is the case as most latitudes and longitudes are specified on the oblate earth geoid so if AZELGEO is not specified by the user during the conversion process then the derived elevation is definitely wrong. The discrepancy becomes greatest around absolute latitude 45 degrees and decreases to zero at the equator and poles. Since my telescope is located at a geodetic latitude of 49 degrees its important that I get this conversion right. I filed a ticket with the NRAO casa group about this issue (see https://help.nrao.edu/index.php?/Tickets/Ticket/View/10739) but it was really never resolved in a satisfactory way. @tammojan might be able to comment on this issue.

[bennahugo]

This is an interesting discussion. From our analysis between CASA and cubical the parallactic angle derotation seems to be very similar which points to CASA itself using spherical AZEL frames when computing angles (as we currently do) within the various gain and polcal steps. However, the parallactic angle derotation do flatten the system response in Q and U as a function of time, so I'm not sure whether the induced angular errors in the TOPO frame here are a significantly large contributing factor? I have a feeling the quoted 2' offset is negligible compared to the possible system calibration performance. My feeling is that the imager itself need not do this derotation (no commonly-used imager that I know of does) - the derotation of the corrected data for imaging purposes and rotation of the model data into sky relative frame for self calibration is the job of the calibration procedure in my opinion and is better implemented within katsdpcal? It is far more telling to study the polarimetric response of the calibrated data in visibility space before imaging a target. From long term studies I can recover the ~33 degree angle of 3C286 and the crosshand phase solutions remain stable to a few degrees (well within calibration SNR error bars) over long spans of the observation - something that won't happen if the first order effect of parallactic angle is not correctly accounted for. The much larger systematic imaging polarimetric errors I think lie in the direction dependent calibration off boresight (> 5') which remains a unsolved problem. The proof of the pudding so to speak would be to check the angles against the moon after accounting for ionospheric RM. Cheers,

…

On Thu, Dec 12, 2019 at 12:07 AM Tony Willis ***@***.***> wrote: @bmerry <https://github.com/bmerry> As I recall you are correct in that casacore (still) uses AZEL as the default specification for the derivation of Azel from a latitude and longitude. I have no idea why that is the case as most latitudes and longitudes are specified on the oblate earth geoid so if AZELGEO is not specified by the user during the conversion process then the derived elevation is definitely wrong. The discrepancy becomes greatest around absolute latitude 45 degrees and decreases to zero at the equator and poles. Since my telescope is located at a geodetic latitude of 49 degrees its important that I get this conversion right. I filed a ticket with the NRAO casa group about this issue (see https://help.nrao.edu/index.php?/Tickets/Ticket/View/10739) but it was really never resolved in a satisfactory way. — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub <#21>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AB4RE6WUFTDQRNANTWQODTDQYFQD3ANCNFSM4E4ACCFA> .

-- -- Benjamin Hugo PhD. student, Centre for Radio Astronomy Techniques and Technologies Department of Physics and Electronics Rhodes University Junior software developer Radio Astronomy Research Group South African Radio Astronomy Observatory Black River Business Park Observatory Cape Town

[twillis449]

This issue has also cropped up in ratt-ru/meqtrees#877 where there is a plot of the difference between geodetic and geocentric elevation as a function of observer latitude.