search

casacore / casacore

Suite of C++ libraries for radio astronomy data processing
http://casacore.github.io/casacore
Other
118 stars 87 forks source link

Difference in ITRF to ICRS conversion values between casacore and astropy. #1107

Open dasingh92 opened 3 years ago

dasingh92 commented 3 years ago

I'm trying to convert antenna positions in ITRF to ICRS frame using casacore. I use the following code for the conversion.

// Antenna Positions as a vector.
    casacore::Vector<casacore::Double> position{-4752405.23838012, 2790372.78676016, -3200482.26399622};

    // Observation time as given by ms file
    casacore::MEpoch tbm{casacore::Quantity{4810734844.928455, "s"}, casacore::MEpoch::Ref(casacore::MEpoch::UTC)};

    // Creating position with reference system
    casacore::MPosition pos{casacore::MVPosition{position}, casacore::MPosition::ITRF};

    // Time and position as a joint reference frame
    casacore::MeasFrame frame{tbm, pos};

    std::cout << " Current frame: " << frame << std::endl;
    // ITRF direction frame
    casacore::MDirection dir{casacore::MVDirection{pos.getAngle()}, casacore::MDirection::ITRF};

    // Conversion object
    casacore::MDirection::Convert conv{casacore::MDirection::ITRF, casacore::MDirection::Ref{casacore::MDirection::ICRS, frame}};

    printf("Direction value before conversion: %.8f, %.8f, %.8f\n",dir.getValue().getVector()[0], dir.getValue().getVector()[1], dir.getValue().getVector()[2]);

    casacore::MDirection conv_dir = conv(dir);

    // double norm = 6372960.26670721;
    casacore::Vector<casacore::Double> nDir = conv_dir.getValue().getVector();
    printf("%.8f\t", nDir[0]);
    printf("%.8f\t", nDir[1]);
    printf("%.8f\n", nDir[2]);

The problem I'm facing is that the values returned by casacore and astropy are different. From casacore I get the output : 0.35615424 -0.78773723 -0.50261736 and astropy returns the output as : 0.35615498, -0.78773578, -0.50261912. Both casacore and astropy give the same values before conversion (the ITRF direction values (variable dir) that I print in the code above are the same in astropy as well).

I would appreciate if someone can let me know if I'm doing something wrong or why this might be happening.

dpetry commented 3 years ago

Convert terrestrial coordinates to celestial ones? Please explain the purpose ...

dasingh92 commented 3 years ago

I need to calculate the position of earth bound telescope in ICRS coordinates as these coordinates are used later for calculating beamforming matrix and gram matrices. (I'm trying to implementing an imaging algorithm in C++).

dpetry commented 3 years ago

Which code are you using in astropy?

dasingh92 commented 3 years ago

The astropy code that was used is : 

import numpy as np
import astropy.coordinates as coord
import astropy.time as atime

# This is MJD equivalent of the time in C++ example
time = atime.Time(55679.801445931196, scale="utc", format = "mjd")

layout = np.array(-4752405.23838012, 2790372.78676016, -3200482.26399622)
itrs_layout = coord.CartesianRepresentation(layout)
itrs_position = coord.SkyCoord(itrs_layout, obstime=time, frame="itrs")
icrs_position = itrs_position.transform_to("icrs").cartesian
bennahugo commented 3 years ago

It may be slightly unrelated but I think there may be precision issues in the astropy codebase. Even with something as simple as elevation computation casacore is closer to pyephem than it is to astropy when I turn off precession corrections from UTC to UT.

I'm not sure if the same issue is biting you here (looking at the variations on the least significant digits.

What remains to be done is a cross comparison to CALC?

On Tue, May 18, 2021 at 11:16 AM Dewan Arun Singh @.***> wrote:

The astropy code that was used is :

import numpy as np import astropy.coordinates as coord import astropy.time as atime

This is MJD equivalent of the time in C++ example

time = atime.Time(55679.801445931196, scale="utc", format = "mjd")

layout = np.array(-4752405.23838012, 2790372.78676016, -3200482.26399622) itrs_layout = coord.CartesianRepresentation(layout) itrs_position = coord.SkyCoord(itrs_layout, obstime=time, frame="itrs") icrs_position = itrs_position.transform_to("icrs").cartesian

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/casacore/casacore/issues/1107#issuecomment-843002922, or unsubscribe https://github.com/notifications/unsubscribe-auth/AB4RE6SHH63R7EHFWY3W75TTOIV5PANCNFSM44JDKHLQ .

--

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

bennahugo commented 3 years ago

The precision issues is a fraction of an arcsecond I forgot to mention.

On Tue, May 18, 2021 at 11:27 AM Benna Hugo @.***> wrote:

It may be slightly unrelated but I think there may be precision issues in the astropy codebase. Even with something as simple as elevation computation casacore is closer to pyephem than it is to astropy when I turn off precession corrections from UTC to UT.

I'm not sure if the same issue is biting you here (looking at the variations on the least significant digits.

What remains to be done is a cross comparison to CALC?

On Tue, May 18, 2021 at 11:16 AM Dewan Arun Singh < @.***> wrote:

The astropy code that was used is :

import numpy as np import astropy.coordinates as coord import astropy.time as atime

This is MJD equivalent of the time in C++ example

time = atime.Time(55679.801445931196, scale="utc", format = "mjd")

layout = np.array(-4752405.23838012, 2790372.78676016, -3200482.26399622) itrs_layout = coord.CartesianRepresentation(layout) itrs_position = coord.SkyCoord(itrs_layout, obstime=time, frame="itrs") icrs_position = itrs_position.transform_to("icrs").cartesian

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/casacore/casacore/issues/1107#issuecomment-843002922, or unsubscribe https://github.com/notifications/unsubscribe-auth/AB4RE6SHH63R7EHFWY3W75TTOIV5PANCNFSM44JDKHLQ .

--

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

--

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

dpetry commented 3 years ago

In order to test the precision, you might try converting a given position from ITRF to ICRS and back again and see if you get the same value as the original. You could also repeat the conversion to see if it is reproducible.

dasingh92 commented 3 years ago

That works in astropy, but I haven't tried in casacore (I just assumed it would !). I'll try doing that. Thanks!!

tammojan commented 3 years ago

Astropy uses erfa under the hood. This is a copy of sofa, the IAU implementation. Casacore internally tests against sofa, and there are indeed a few (very small) differences. This test is in tIAU2000, which is only built if SOFA is installed. Example output:

Test IAU2000 coordinate conversions ...
---------------------------------------------------------------------------
A default date of 2003/03/01 is used (JD 2452699.5, MJD 52699)
---------------------------------------------------------------------------
Precession
---------------------------------------------------------------------------
SOFA precession 2000A and bias:
Axis Lengths: [3, 3]  (NB: Matrix in Row/Column order)
[0.999999703039, 0.000706855546336, 0.000307044232973
 -0.000706855534985, 0.999999750178, -1.45484959268e-07
 -0.000307044259103, -7.15509994777e-08, 0.999999952862]

Casacore precession and bias:
Axis Lengths: [3, 3]  (NB: Matrix in Row/Column order)
[0.999999703039, 0.000706855546329, 0.00030704423298
 -0.000706855534978, 0.999999750178, -1.45484959345e-07
 -0.00030704425911, -7.15509994235e-08, 0.999999952862]

SOFA precession     (no bias):
Axis Lengths: [3, 3]  (NB: Matrix in Row/Column order)
[0.999999703065, 0.00070678476355, 0.000307124795151
 -0.000706784762363, 0.999999750228, -1.12402812835e-07
 -0.000307124797884, -1.04668345878e-07, 0.999999952837]

Casacore precession (no bias):
Axis Lengths: [3, 3]  (NB: Matrix in Row/Column order)
[0.999999703065, 0.000706784763543, 0.000307124795157
 -0.000706784762356, 0.999999750228, -1.12402812912e-07
 -0.000307124797891, -1.04668345823e-07, 0.999999952837]

SOFA precession corrections:     -4.59190465063e-08, -3.86783492013e-09
Casacore precession corrections: -4.59190465063e-08, -3.86783492013e-09

SOFA 2000A bias matrix:
Axis Lengths: [3, 3]  (NB: Matrix in Row/Column order)
[1, 7.07827947786e-08, -8.05621738099e-08
 -7.0782797442e-08, 1, -3.30604088398e-08
 8.05621714698e-08, 3.30604145422e-08, 1]

Casacore 2000A bias matrix:
Axis Lengths: [3, 3]  (NB: Matrix in Row/Column order)
[1, 7.07827947786e-08, -8.05621738099e-08
 -7.0782797442e-08, 1, -3.30604088398e-08
 8.05621714698e-08, 3.30604145422e-08, 1]

Difference (arcsec) SOFA 2000A - Casacore nobias:  0.000000002
Difference (arcsec) SOFA 2000A - Casacore   bias:  0.000000002
Difference (arcsec) SOFA 2000A - Casacore nobias:  0.000000002
(The above must all 3 be .002 uas)
---------------------------------------------------------------------------
Nutation
---------------------------------------------------------------------------
Difference (arcsec) SOFA 2000A - SOFA 2000B:         0.000511578
Difference (arcsec) Casacore 2000A - Casacore 2000B: 0.000511973
(The above two differences are 0.512mas (note next note))
Difference (arcsec) SOFA 2000A - Casacore 2000A:   0.000000000
Difference (arcsec) SOFA 2000B - Casacore 2000B:   0.000000422
(This should be 0.422 uas, due to different definition of
    the fundamental arguments)
Equation of equinoxes (IAU2000A)(SOFA, Casacore, diff):
    -0.000062314, -0.000062314, -0.000000001
---------------------------------------------------------------------------
SOFA method comparisons ...
---------------------------------------------------------------------------
Difference (arcsec) SOFA CEO and Equinox based in 1935:   0.000000119
Difference (arcsec) SOFA CEO based and old J2000 in 1935: 0.065034218
(The above should be 0.119 uas and 65.034 mas respectively)
---------------------------------------------------------------------------
IAU2000A/B comparisons ...
---------------------------------------------------------------------------
Registrations old: 3, 4
New J2000 0, J2000A 0
Direction: [0.557172304, 0.320952009, 0.765864761]
New J2000 1, J2000A 0
Direction: [0.557172269, 0.320952031, 0.765864777]
Difference: [-0.006987056, 0.004427103, 0.003227863]
New J2000 1, J2000A 1
Direction: [0.557172268, 0.320952033, 0.765864777]
Difference: [-0.007155719, 0.004781831, 0.003201910]
---------------------------------------------------------------------------
Test of some details ...
---------------------------------------------------------------------------
UT: -23741.250746343, -23741.250746343
TT: -23741.250000000, -23741.250000000
s'   (Sofa, Casacore, diff): 0.000000000, 0.000000000, 0.000000000
ERA  (Sofa, Casacore, diff): 3.325261764, 3.325261764, 0.000000000
GMST (Sofa, Casacore, diff): 3.310730455, 3.310730455, -0.000000000
GMST82 (GMST82, GMST00, diff): 34451.526918968, 34451.526919117, -0.000000149
EqEqCT00 (Sofa, Casacore, diff): -0.000000011, -0.000000011, 0.000000000
---------------------------------------------------------------------------
Test Aipsrc value cross talk ...
---------------------------------------------------------------------------
Registrations now: 3, 4
New J2000 0, J2000A 0
---------------------------------------------------------------------------
OK

There are some very small differences, which due to differences in the equation of equinoxes. I must confess I don't understand this yet, but Wim Brouw (who implemented this) does.

Are the differences you're seeing of this order of magnitude?

dasingh92 commented 3 years ago

Yes they are. The differences usually fluctuate between 10e-6 & 10e-8.