Skip to content

Astropy 4.3+ Produces Inaccurate Range Predictions via AltAz Calc #13319

New issue
New issue
Open
Open
Astropy 4.3+ Produces Inaccurate Range Predictions via AltAz Calc#13319
Labels
Bugcoordinatestime

Description

@mstevensonll
mstevensonll
opened on Jun 9, 2022

Description

We regularly compare our radar Range measurements to ILRS ephems. Upon upgrading to Astropy 5.1, we started seeing very large residuals. We've tried different versions and found astropy 4.2 gives expected agreement, while astropy 4.3 and later give the unexpected behavior.

We believe that the astropy 4.3+ predictions are in error because we have internal residual analysis and comparison to other data sources (not using astropy) that do not show this issue.

Expected behavior

When we run it with astropy 4.2, we get nice small residuals like this:

(astropy5) mas@ip-192-168-4-110 scripts % python astropy_snippet_20220608.py
	Calculated measurement residual
		range: 0.1m
	Calculated measurement residual
		range: 12.4m
	Calculated measurement residual
		range: 17.8m

Actual behavior

When we run it with astropy 4.3, we get great big residuals like this:

(astropy5) mas@ip-192-168-4-110 scripts % python astropy_snippet_20220608.py
	Calculated measurement residual
		range: -112.2m
	Calculated measurement residual
		range: -99.9m
	Calculated measurement residual
		range: -94.4m

Steps to Reproduce

  1. Get some radar Range measurements
  2. Get ILRS Ephem predictions
  3. Create an EarthLocation for the radar location.
  4. Interpolate the ILRS ephem for every measurement timestamp.
  5. Create a SkyCoord object for the interpolated position.
  6. Do an AltAz transformation to get the Range prediction.
import numpy as np
import scipy.interpolate
import datetime
import dateutil.parser

import astropy.time, astropy.coordinates
import astropy.units as u
from astropy.utils import iers
iers.conf.iers_auto_url='https://maia.usno.navy.mil/ser7/finals2000A.all'

# Instrument:
inst = {
    'longitude': -85.479833,
    'latitude': 10.413333,
    'altitude': 23.0,
    }

measurements = [
    {
        'timestamp_measurement': dateutil.parser.parse('2022-06-08 08:18:45.103596+00'),
        'range_value': 1813723.24341158,
    },
    {
        'timestamp_measurement': dateutil.parser.parse('2022-06-08 08:18:44.983595+00'),
        'range_value':  1813538.3470108,
    },
    {
        'timestamp_measurement': dateutil.parser.parse('2022-06-08 08:18:44.926049+00'),
        'range_value':  1813449.3366306,
    },
]

ilrs_mjd = 59738
ilrs_secs = np.array([
    27900.00000, 28080.00000, 28260.00000, 28440.00000, 28620.00000,
    28800.00000, 28980.00000, 29160.00000, 29340.00000, 29520.00000,
    29700.00000, 29880.00000, 30060.00000, 30240.00000, 30420.00000,
    30600.00000, 30780.00000, 30960.00000, 31140.00000, 31320.00000,
    31500.00000, 31680.00000, 31860.00000, 32040.00000, 32220.00000,
    ])

ilrs_data = np.array([
    [-2741944.579,   -235689.755,  -7325387.829],
    [-2711657.122,  -1474364.598,  -7189654.207],
    [-2641154.833,  -2675083.522,  -6860757.255],
    [-2530855.876,  -3806852.586,  -6347386.337],
    [-2381441.653,  -4840435.015,  -5663181.834],
    [-2193903.565,  -5749124.122,  -4826398.133],
    [-1969601.044,  -6509460.381,  -3859440.304],
    [-1710326.586,  -7101871.185,  -2788282.056],
    [-1418372.400,  -7511212.894,  -1641776.504],
    [-1096593.106,  -7727194.734,   -450875.156],
    [ -748458.835,  -7744665.880,    752222.831],
    [ -378094.716,  -7563753.660,   1934963.517],
    [    9698.097,  -7189847.764,   3065344.786],
    [  409443.848,  -6633429.388,   4112807.145],
    [  815017.796,  -5909751.857,   5049077.091],
    [ 1219697.263,  -5038388.983,   5848936.440],
    [ 1616238.610,  -4042669.245,   6490897.504],
    [ 1996980.100,  -2949016.560,   6957766.909],
    [ 2353970.392,  -1786221.732,   7237086.985],
    [ 2679123.360,   -584665.510,   7321450.016],
    [ 2964396.995,    624488.873,   7208680.076],
    [ 3201992.320,   1810113.002,   6901878.965],
    [ 3384567.960,   2941919.905,   6409337.507],
    [ 3505465.893,   3991242.320,   5744316.504],
    [ 3558942.459,   4931767.385,   4924700.825],
    ])

ilrs_time = astropy.time.Time(ilrs_mjd + ilrs_secs / 86400, format='mjd', scale='utc')
x = ilrs_data[:,0]
y = ilrs_data[:,1]
z = ilrs_data[:,2]

# instrument location
topo_loc = astropy.coordinates.EarthLocation(
    lat=inst['latitude']*u.deg, lon=inst['longitude']*u.deg, height=inst['altitude']*u.m
)

for meas in measurements:

    # measurements
    time_meas = meas['timestamp_measurement']
    ttime = astropy.time.Time(time_meas, format='datetime', scale='utc')
    trange = meas['range_value']

    # find closest point in time
    imin = abs(ilrs_time - ttime).argmin()

    # interpolate, 9 point
    tinterp = range(imin - 4, imin + 5)
    t0 = ilrs_time[tinterp][0]
    tl = (ilrs_time[tinterp] - t0).sec  # ilrs times
    teval = (ttime - t0).sec + np.array([-0.05, 0, 0.05])  # times to evaluate at
    px = scipy.interpolate.lagrange(tl, x[tinterp])
    py = scipy.interpolate.lagrange(tl, y[tinterp])
    pz = scipy.interpolate.lagrange(tl, z[tinterp])
    posi = np.array([px(teval), py(teval), pz(teval)])

    # position in ITRS
    dpos = astropy.coordinates.SkyCoord(
        x=posi[0, :], y=posi[1, :], z=posi[2, :],
        representation_type='cartesian',
        frame="itrs", unit="meter"
    )

    # position relative to site
    dpos_altaz = dpos.transform_to(
        astropy.coordinates.AltAz(location=topo_loc)
    )
    rs = dpos_altaz.distance.m

    # residuals
    r_resid = trange - rs[1]

    print('\tCalculated measurement residual')
    print('\t\trange: %2.1fm' % r_resid)

System Details

Please note that we have reproduced this in several environments.

We were able to reproduce it with the following environment:

Python 3.8.13 (default, Mar 16 2022, 20:38:02) 
[Clang 13.0.0 (clang-1300.0.29.30)] on darwin
Type "help", "copyright", "credits" or "license" for more information.
>>> import platform; print(platform.platform())
macOS-11.6-x86_64-i386-64bit
>>> import sys; print("Python", sys.version)
Python 3.8.13 (default, Mar 16 2022, 20:38:02) 
[Clang 13.0.0 (clang-1300.0.29.30)]
>>> import numpy; print("Numpy", numpy.__version__)
Numpy 1.22.4
>>> import erfa; print("pyerfa", erfa.__version__)
pyerfa 2.0.0.1
>>> import astropy; print("astropy", astropy.__version__)
astropy 4.3
>>> import scipy; print("Scipy", scipy.__version__)
Scipy 1.8.0
>>> import matplotlib; print("Matplotlib", matplotlib.__version__)
Matplotlib 3.5.2
>>>

Activity

  1. mstevensonll
    added
    Bug
    on Jun 9, 2022

  2. github-actions commented on Jun 9, 2022

    @github-actions
    github-actionsBot
    on Jun 9, 2022
    Contributor

    Welcome to Astropy 👋 and thank you for your first issue!

    A project member will respond to you as soon as possible; in the meantime, please double-check the guidelines for submitting issues and make sure you've provided the requested details.

    GitHub issues in the Astropy repository are used to track bug reports and feature requests; If your issue poses a question about how to use Astropy, please instead raise your question in the Astropy Discourse user forum and close this issue.

    If you feel that this issue has not been responded to in a timely manner, please leave a comment mentioning our software support engineer @embray, or send a message directly to the development mailing list. If the issue is urgent or sensitive in nature (e.g., a security vulnerability) please send an e-mail directly to the private e-mail [email protected].

  3. pllim
    added
    coordinates
    time
    on Jun 9, 2022

  4. mkbrewer commented on Jun 16, 2022

    @mkbrewer
    mkbrewer
    on Jun 16, 2022
    Contributor

    This behavior is due to a bug that was fixed in 4.3. Previously, the AltAz transform ignored diurnal aberration. The way that Astropy works when transforming an ITRS position to AltAz is rather round about. It first transforms the ITRS position to ICRS removing geocentric annual aberration in the process. Then it transforms back to AltAz adding in both geocentric annual and topocentric diurnal aberration. It is this addition of diurnal aberration that is causing your problem. There is a way to get around this though. What you need to do is transform both the target and instrument positions to geocentric CIRS. Then you can difference these directly to form a topocentric CIRS position that has no aberration included. Transforming this to AltAz will yield your desired result. I have modified your script to do that with all credit due to Stuart Littlefair for the added code from his test_straight_overhead() function from astropy.coordinates.tests.test_intermediate_transformations.

    import numpy as np
import scipy.interpolate
import datetime
import dateutil.parser

import astropy.time, astropy.coordinates
from astropy.coordinates import CIRS
import astropy.units as u
from astropy.utils import iers
iers.conf.iers_auto_url='https://maia.usno.navy.mil/ser7/finals2000A.all'


# Instrument:
inst = {
    'longitude': -85.479833,
    'latitude': 10.413333,
    'altitude': 23.0,
    }

measurements = [
    {
        'timestamp_measurement': dateutil.parser.parse('2022-06-08 08:18:45.103596+00'),
        'range_value': 1813723.24341158,
    },
    {
        'timestamp_measurement': dateutil.parser.parse('2022-06-08 08:18:44.983595+00'),
        'range_value':  1813538.3470108,
    },
    {
        'timestamp_measurement': dateutil.parser.parse('2022-06-08 08:18:44.926049+00'),
        'range_value':  1813449.3366306,
    },
]

ilrs_mjd = 59738
ilrs_secs = np.array([
    27900.00000, 28080.00000, 28260.00000, 28440.00000, 28620.00000,
    28800.00000, 28980.00000, 29160.00000, 29340.00000, 29520.00000,
    29700.00000, 29880.00000, 30060.00000, 30240.00000, 30420.00000,
    30600.00000, 30780.00000, 30960.00000, 31140.00000, 31320.00000,
    31500.00000, 31680.00000, 31860.00000, 32040.00000, 32220.00000,
    ])

ilrs_data = np.array([
    [-2741944.579,   -235689.755,  -7325387.829],
    [-2711657.122,  -1474364.598,  -7189654.207],
    [-2641154.833,  -2675083.522,  -6860757.255],
    [-2530855.876,  -3806852.586,  -6347386.337],
    [-2381441.653,  -4840435.015,  -5663181.834],
    [-2193903.565,  -5749124.122,  -4826398.133],
    [-1969601.044,  -6509460.381,  -3859440.304],
    [-1710326.586,  -7101871.185,  -2788282.056],
    [-1418372.400,  -7511212.894,  -1641776.504],
    [-1096593.106,  -7727194.734,   -450875.156],
    [ -748458.835,  -7744665.880,    752222.831],
    [ -378094.716,  -7563753.660,   1934963.517],
    [    9698.097,  -7189847.764,   3065344.786],
    [  409443.848,  -6633429.388,   4112807.145],
    [  815017.796,  -5909751.857,   5049077.091],
    [ 1219697.263,  -5038388.983,   5848936.440],
    [ 1616238.610,  -4042669.245,   6490897.504],
    [ 1996980.100,  -2949016.560,   6957766.909],
    [ 2353970.392,  -1786221.732,   7237086.985],
    [ 2679123.360,   -584665.510,   7321450.016],
    [ 2964396.995,    624488.873,   7208680.076],
    [ 3201992.320,   1810113.002,   6901878.965],
    [ 3384567.960,   2941919.905,   6409337.507],
    [ 3505465.893,   3991242.320,   5744316.504],
    [ 3558942.459,   4931767.385,   4924700.825],
    ])

ilrs_time = astropy.time.Time(ilrs_mjd + ilrs_secs / 86400, format='mjd', scale='utc')
x = ilrs_data[:,0]
y = ilrs_data[:,1]
z = ilrs_data[:,2]

# instrument location
topo_loc = astropy.coordinates.EarthLocation(
    lat=inst['latitude']*u.deg, lon=inst['longitude']*u.deg, height=inst['altitude']*u.m)

for meas in measurements:

    # measurements
    time_meas = meas['timestamp_measurement']
    ttime = astropy.time.Time(time_meas, format='datetime', scale='utc')
    trange = meas['range_value']

    # find closest point in time
    imin = abs(ilrs_time - ttime).argmin()

    # interpolate, 9 point
    tinterp = range(imin - 4, imin + 5)
    t0 = ilrs_time[tinterp][0]
    tl = (ilrs_time[tinterp] - t0).sec  # ilrs times
    teval = (ttime - t0).sec + np.array([-0.05, 0, 0.05])  # times to evaluate at
    px = scipy.interpolate.lagrange(tl, x[tinterp])
    py = scipy.interpolate.lagrange(tl, y[tinterp])
    pz = scipy.interpolate.lagrange(tl, z[tinterp])
    posi = np.array([px(teval), py(teval), pz(teval)])

    # position in ITRS
    dpos = astropy.coordinates.SkyCoord(
        x=posi[0, :], y=posi[1, :], z=posi[2, :],
        representation_type='cartesian',
        frame="itrs", unit="meter"
    )
    # get the Geocentric CIRS position
    cirs_geo = dpos.transform_to(CIRS()) 
    # get the Geocentric CIRS position of the observer
    obsrepr = topo_loc.get_itrs().transform_to(CIRS()).cartesian
    # take the difference
    cirs_repr = cirs_geo.cartesian - obsrepr
    # form the Topocentric CIRS position
    topocentric_cirs_frame = CIRS(location=topo_loc)
    cirs_topo = topocentric_cirs_frame.realize_frame(cirs_repr) 
    # position relative to site
    # We seem to need a time here
    t = astropy.time.Time('J2000')
    dpos_altaz = cirs_topo.transform_to(
        astropy.coordinates.AltAz(obstime=t, location=topo_loc)
    )
    rs = dpos_altaz.distance.m

    # residuals
    r_resid = trange - rs[1]

    print('\tCalculated measurement residual')
    print('\t\trange: %2.1fm' % r_resid)

  5. mkbrewer commented on Jun 17, 2022

    @mkbrewer
    mkbrewer
    on Jun 17, 2022
    Contributor

    A more straightforward (and faster) approach would be to just form the topocentric ITRS position and do the AltAz transform yourself:

    import numpy as np
import scipy.interpolate
import datetime
import dateutil.parser

import astropy.time, astropy.coordinates
from astropy.coordinates.matrix_utilities import rotation_matrix
import astropy.units as u
from astropy.utils import iers
iers.conf.iers_auto_url='https://maia.usno.navy.mil/ser7/finals2000A.all'


# Instrument:
inst = {
    'longitude': -85.479833,
    'latitude': 10.413333,
    'altitude': 23.0,
    }

measurements = [
    {
        'timestamp_measurement': dateutil.parser.parse('2022-06-08 08:18:45.103596+00'),
        'range_value': 1813723.24341158,
    },
    {
        'timestamp_measurement': dateutil.parser.parse('2022-06-08 08:18:44.983595+00'),
        'range_value':  1813538.3470108,
    },
    {
        'timestamp_measurement': dateutil.parser.parse('2022-06-08 08:18:44.926049+00'),
        'range_value':  1813449.3366306,
    },
]

ilrs_mjd = 59738
ilrs_secs = np.array([
    27900.00000, 28080.00000, 28260.00000, 28440.00000, 28620.00000,
    28800.00000, 28980.00000, 29160.00000, 29340.00000, 29520.00000,
    29700.00000, 29880.00000, 30060.00000, 30240.00000, 30420.00000,
    30600.00000, 30780.00000, 30960.00000, 31140.00000, 31320.00000,
    31500.00000, 31680.00000, 31860.00000, 32040.00000, 32220.00000,
    ])

ilrs_data = np.array([
    [-2741944.579,   -235689.755,  -7325387.829],
    [-2711657.122,  -1474364.598,  -7189654.207],
    [-2641154.833,  -2675083.522,  -6860757.255],
    [-2530855.876,  -3806852.586,  -6347386.337],
    [-2381441.653,  -4840435.015,  -5663181.834],
    [-2193903.565,  -5749124.122,  -4826398.133],
    [-1969601.044,  -6509460.381,  -3859440.304],
    [-1710326.586,  -7101871.185,  -2788282.056],
    [-1418372.400,  -7511212.894,  -1641776.504],
    [-1096593.106,  -7727194.734,   -450875.156],
    [ -748458.835,  -7744665.880,    752222.831],
    [ -378094.716,  -7563753.660,   1934963.517],
    [    9698.097,  -7189847.764,   3065344.786],
    [  409443.848,  -6633429.388,   4112807.145],
    [  815017.796,  -5909751.857,   5049077.091],
    [ 1219697.263,  -5038388.983,   5848936.440],
    [ 1616238.610,  -4042669.245,   6490897.504],
    [ 1996980.100,  -2949016.560,   6957766.909],
    [ 2353970.392,  -1786221.732,   7237086.985],
    [ 2679123.360,   -584665.510,   7321450.016],
    [ 2964396.995,    624488.873,   7208680.076],
    [ 3201992.320,   1810113.002,   6901878.965],
    [ 3384567.960,   2941919.905,   6409337.507],
    [ 3505465.893,   3991242.320,   5744316.504],
    [ 3558942.459,   4931767.385,   4924700.825],
    ])

ilrs_time = astropy.time.Time(ilrs_mjd + ilrs_secs / 86400, format='mjd', scale='utc')
x = ilrs_data[:,0]
y = ilrs_data[:,1]
z = ilrs_data[:,2]

# instrument location
topo_loc = astropy.coordinates.EarthLocation(
    lat=inst['latitude']*u.deg, lon=inst['longitude']*u.deg, height=inst['altitude']*u.m)

for meas in measurements:

    # measurements
    time_meas = meas['timestamp_measurement']
    ttime = astropy.time.Time(time_meas, format='datetime', scale='utc')
    trange = meas['range_value']

    # find closest point in time
    imin = abs(ilrs_time - ttime).argmin()

    # interpolate, 9 point
    tinterp = range(imin - 4, imin + 5)
    t0 = ilrs_time[tinterp][0]
    tl = (ilrs_time[tinterp] - t0).sec  # ilrs times
    teval = (ttime - t0).sec + np.array([-0.05, 0, 0.05])  # times to evaluate at
    px = scipy.interpolate.lagrange(tl, x[tinterp])
    py = scipy.interpolate.lagrange(tl, y[tinterp])
    pz = scipy.interpolate.lagrange(tl, z[tinterp])
    posi = np.array([px(teval), py(teval), pz(teval)])

    # position in ITRS
    dpos = astropy.coordinates.SkyCoord(
        x=posi[0, :], y=posi[1, :], z=posi[2, :],
        representation_type='cartesian',
        frame="itrs", unit="meter"
    )
    # form the Topocentric ITRS position
    topocentric_itrs_repr = dpos.cartesian - topo_loc.get_itrs().cartesian
    # form ITRS to AltAz matrix
    lon, lat, height = topo_loc.to_geodetic('WGS84')
    elong = lon.to_value(u.radian)
    elat = lat.to_value(u.radian)
    minus_x = np.eye(3)
    minus_x[0][0] *= -1.0
    r = (minus_x 
         @ rotation_matrix(np.pi / 2.0 - elat, 'y', unit=u.radian) 
         @ rotation_matrix(elong, 'z', unit=u.radian))
    # transform to AltAz
    altaz_repr = topocentric_itrs_repr.transform(r)
    altaz_frame = astropy.coordinates.AltAz(location=topo_loc)
    dpos_altaz = altaz_frame.realize_frame(altaz_repr) 
    rs = dpos_altaz.distance.m

    # residuals
    r_resid = trange - rs[1]

    print('\tCalculated measurement residual')
    print('\t\trange: %2.1fm' % r_resid)

  6. mkbrewer commented on Jun 17, 2022

    @mkbrewer
    mkbrewer
    on Jun 17, 2022
    Contributor

    Actually, it would be better to form the ITRS to AltAz matrix outside the loop.

  7. mkbrewer commented on Jun 18, 2022

    @mkbrewer
    mkbrewer
    on Jun 18, 2022
    Contributor

    Since this is a recurring issue, I will raise a separate issue to propose the above more straightforward approach. If I get the go ahead, I will submit a pull request that will solve this problem once and for all.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Metadata

Metadata

Assignees

No one assigned

    Labels

    Bugcoordinatestime

    Type

    No type

    Projects

    No projects

      Milestone

      No milestone

      Relationships

      None yet

      Development

      No branches or pull requests

      Issue actions